Conservative AMR (#2028)

* small typo fixes * Interpolate and project Ju instead of u to retain discrete conservation * add specialized routines for TreeMesh and P4estMesh. Hopefully some tests pass now * remove unused copy * remove false comment * simplify the method on the children * update 2d p4est tests * update 3d p4est tests * update 2d t8code tests * update 3d t8code tests * update mpi test values * fix formatting * remove unnecessary comment * add new 2d tests including conservation * add new 3d tests including conservation * update t8code 2d mpi test * remove unnecessary Union * add comments regarding how element ids are incremented during refine or coarsen * typo fixes * another typo fix * fix comment * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * combine GC.preserve blocks * fix formatting in elixir * simplify implementation and add if blocks for P4estMesh specialities * unify structure across P4estMesh and T8codeMesh * fix some of the bad formatting * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * add news * Update NEWS.md Co-authored-by: Andrew Winters <[email protected]> --------- Co-authored-by: Daniel Doehring <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Aug 19, 2024 · ba82d36 · ba82d36
1 parent 26ab5fa
commit ba82d36
Show file tree

Hide file tree

Showing 18 changed files with 980 additions and 129 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes in the v0.8 lifecycle
+
+#### Changed
+
+- The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product
+  of the Jacobian and the conserved variables instead of the conserved variables only
+  to make AMR fully conservative ([#2028]). This may change AMR results slightly.
+
 ## Changes when updating to v0.8 from v0.7.x
 
 #### Added

diff --git a/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,117 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+
+    r0 = 0.2
+    E = 1
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+# Get the DG approximation space
+
+# Activate the shock capturing + flux differencing
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+###############################################################################
+
+# Affine type mapping to take the [-1,1]^2 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+# Warping with the coefficient 0.2 is even more extreme.
+function mapping_twist(xi, eta)
+    y = eta + 0.125 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta)
+    x = xi + 0.125 * cos(0.5 * pi * xi) * cos(2.0 * pi * y)
+    return SVector(x, y)
+end
+
+# The mesh below can be made periodic
+# Create P4estMesh with 8 x 8 trees
+trees_per_dimension = (8, 8)
+mesh = P4estMesh(trees_per_dimension, polydeg = 4,
+                 mapping = mapping_twist,
+                 initial_refinement_level = 0,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 400
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.2,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.05,
+                                      max_level = 2, max_threshold = 0.1)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);#, maxiters=4);
+summary_callback() # print the timer summary
diff --git a/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,114 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+function initial_condition_weak_blast_wave(x, t,
+                                           equations::CompressibleEulerEquations3D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
+
+    r0 = 0.2
+    E = 1.0
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    v3 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 1.0,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# Setup a periodic mesh with 4 x 4 x 4 trees and 8 x 8 x 8 elements
+trees_per_dimension = (4, 4, 4)
+
+# Affine type mapping to take the [-1,1]^3 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+function mapping_twist(xi, eta, zeta)
+    y = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta))
+
+    x = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta))
+
+    z = zeta + 1 / 6 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta))
+
+    return SVector(x, y, z)
+end
+
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 2,
+                 mapping = mapping_twist,
+                 initial_refinement_level = 1,
+                 periodicity = true)
+
+# Create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 1,
+                                      med_level = 2, med_threshold = 0.05,
+                                      max_level = 3, max_threshold = 0.15)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 1,
+                           adapt_initial_condition = false,
+                           adapt_initial_condition_only_refine = false)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
@@ -61,7 +61,8 @@ amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first)
 amr_callback = AMRCallback(semi, amr_controller,
                            interval = 5,
                            adapt_initial_condition = true,
-                           adapt_initial_condition_only_refine = true)
+                           adapt_initial_condition_only_refine = true,
+                           dynamic_load_balancing = true)
 
 stepsize_callback = StepsizeCallback(cfl = 0.7)
 

diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,114 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+
+    r0 = 0.2
+    E = 1
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+# Get the DG approximation space
+
+# Activate the shock capturing + flux differencing
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# Affine type mapping to take the [-1,1]^2 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+# Warping with the coefficient 0.2 is even more extreme.
+function mapping_twist(xi, eta)
+    y = eta + 0.125 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta)
+    x = xi + 0.125 * cos(0.5 * pi * xi) * cos(2.0 * pi * y)
+    return SVector(x, y)
+end
+
+# The mesh below can be made periodic
+# Create T8codeMesh with 8 x 8 trees
+trees_per_dimension = (8, 8)
+mesh = T8codeMesh(trees_per_dimension, polydeg = 4,
+                  mapping = mapping_twist,
+                  initial_refinement_level = 0,
+                  periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 400
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.05,
+                                      max_level = 2, max_threshold = 0.1)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);#, maxiters=4);
+summary_callback() # print the timer summary
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)