AMR for parabolic terms in 2D & 3D on TreeMeshes (#1629)

* Clean branch

* Un-Comment

* un-comment

* test coarsen

* remove redundancy

* Remove support for passive terms

* expand resize

* comments

* format

* Avoid code duplication

* Update src/callbacks_step/amr_dg1d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* comment

* comment & format

* Try to increase coverage

* Slightly more expressive names

* Apply suggestions from code review

* add specifier for 1d

* Structs for resizing parabolic helpers

* check if mortars are present

* reuse `reinitialize_containers!`

* resize calls for parabolic helpers

* update analysis callbacks

* Velocities for compr euler

* Init container

* correct copy-paste error

* resize each dim

* add dispatch

* Add AMR for shear layer

* USe only amr shear layer

* first steps towards p4est parabolic amr

* Add tests

* remove plots

* Format

* remove redundant line

* platform independent tests

* No need for different flux_viscous comps after adding container_viscous to p4est

* Laplace 3d

* Longer times to allow converage to hit coarsen!

* Increase testing of Laplace 3D

* Add tests for velocities

* remove comment

* Add comments

* Remove some specializations

* Add comments

* Use tuple for outer, fixed size datastruct for internal quantities

* Format

* Add comments

* Update examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/Trixi.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Move velocity into elixir

* remove tests

* Remove deprecated comments

* Add news

---------

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
Co-authored-by: Jesse Chan <[email protected]>

NEWS.md

@@ -14,6 +14,7 @@ for human readability.
 - Wetting and drying feature and examples for 1D and 2D shallow water equations
 - Implementation of the polytropic Euler equations in 2D
 - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh`
+- AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh`
 
 #### Changed
 

examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl → examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl

@@ -7,12 +7,20 @@ using Trixi
 
 # TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 1.0/3.0 * 10^(-3) # equivalent to Re = 3000
+mu() = 1.0/3.0 * 10^(-5) # equivalent to Re = 30,000
 
 equations = CompressibleEulerEquations2D(1.4)
 equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
                                                           Prandtl=prandtl_number())
 
+"""
+A compressible version of the double shear layer initial condition. Adapted from
+Brown and Minion (1995).
+
+- David L. Brown and Michael L. Minion (1995)
+  Performance of Under-resolved Two-Dimensional Incompressible Flow Simulations.
+  [DOI: 10.1006/jcph.1995.1205](https://doi.org/10.1006/jcph.1995.1205)
+"""
 function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D)
   # Shear layer parameters
   k = 80
@@ -22,8 +30,8 @@ function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquatio
   Ms = 0.1 # maximum Mach number
 
   rho = 1.0
-  v1  = x[2] <= 0.5 ? u0 * tanh(k*(x[2]*0.5 - 0.25)) : u0 * tanh(k*(0.75 -x[2]*0.5))
-  v2  = u0 * delta * sin(2*pi*(x[1]*0.5 + 0.25))
+  v1  = x[2] <= 0.5 ? u0 * tanh(k*(x[2] - 0.25)) : u0 * tanh(k*(0.75 -x[2]))
+  v2  = u0 * delta * sin(2*pi*(x[1]+ 0.25))
   p   = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms
 
   return prim2cons(SVector(rho, v1, v2, p), equations)
@@ -47,27 +55,43 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 2.0)
+tspan = (0.0, 1.0)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 50
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
-                                     extra_analysis_integrals=(energy_kinetic,
-                                                               energy_internal,
-                                                               enstrophy))
+analysis_interval = 2000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval=analysis_interval,)
 
+# This uses velocity-based AMR
+@inline function v1(u, equations::CompressibleEulerEquations2D)
+  rho, rho_v1, _, _ = u
+  return rho_v1 / rho
+end
+amr_indicator = IndicatorLöhner(semi, variable=v1)                                          
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 3,
+                                      med_level  = 5, med_threshold=0.2,
+                                      max_level  = 7, max_threshold=0.5)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=50,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.3)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
-                        alive_callback)
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
 
-time_int_tol = 1e-8
-sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
-            ode_default_options()..., callback=callbacks)
+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

examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl

@@ -0,0 +1,93 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+diffusivity() = 5.0e-4
+equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations)
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0, -1.0)
+coordinates_max = ( 1.0,  1.0,  1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=80_000)
+
+# Define initial condition
+function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation3D)
+  # Store translated coordinate for easy use of exact solution
+  x_trans = x - equation.advection_velocity * t
+
+  nu = diffusivity()
+  c = 1.0
+  A = 0.5
+  L = 2
+  f = 1/L
+  omega = 2 * pi * f
+  scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t)
+  return SVector(scalar)
+end
+initial_condition = initial_condition_diffusive_convergence_test
+
+# define periodic boundary conditions everywhere
+boundary_conditions = boundary_condition_periodic
+boundary_conditions_parabolic = boundary_condition_periodic
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions,
+                                                                  boundary_conditions_parabolic))
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(entropy,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=3,
+                                      med_level=4, med_threshold=1.2,
+                                      max_level=5, max_threshold=1.45)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+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);
+summary_callback() # print the timer summary

examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl

@@ -0,0 +1,91 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 5.0e-2
+advection_velocity = (1.0, 0.0, 0.0)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -0.5, -0.25) # minimum coordinates (min(x), min(y), min(z))
+coordinates_max = ( 0.0,  0.5,  0.25) # maximum coordinates (max(x), max(y), max(z))
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=3,
+                periodicity=false,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# Example setup taken from
+# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016).
+#   Robust DPG methods for transient convection-diffusion.
+#   In: Building bridges: connections and challenges in modern approaches
+#   to numerical partial differential equations.
+#   [DOI](https://doi.org/10.1007/978-3-319-41640-3_6).
+function initial_condition_eriksson_johnson(x, t, equations)
+  l = 4
+  epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt
+  lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
+  lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
+  r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
+  s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
+  u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) +
+      cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1))
+  return SVector{1}(u)
+end
+initial_condition = initial_condition_eriksson_johnson
+
+boundary_conditions = (; x_neg = BoundaryConditionDirichlet(initial_condition),
+                         y_neg = BoundaryConditionDirichlet(initial_condition),
+                         z_neg = boundary_condition_do_nothing,
+                         y_pos = BoundaryConditionDirichlet(initial_condition),
+                         x_pos = boundary_condition_do_nothing,
+                         z_pos = boundary_condition_do_nothing)
+
+boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions,
+                                                                  boundary_conditions_parabolic))
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+time_int_tol = 1.0e-11
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+
+# Print the timer summary
+summary_callback()

src/Trixi.jl

@@ -153,7 +153,7 @@ export AcousticPerturbationEquations2D,
        LinearizedEulerEquations2D,
        PolytropicEulerEquations2D
 
-export LaplaceDiffusion1D, LaplaceDiffusion2D,
+export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
        CompressibleNavierStokesDiffusion3D
 

src/callbacks_step/amr_dg1d.jl

@@ -83,30 +83,13 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
     # actually transferring the solution to the refined cells
     refine!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_refine)
 
-    # The remaining function only handles the necessary adaptation of the data structures
-    # for the parabolic part of the semidiscretization
-
-    # Get new list of leaf cells
-    leaf_cell_ids = local_leaf_cells(mesh.tree)
-
-    @unpack elements, viscous_container = cache_parabolic
-    resize!(elements, length(leaf_cell_ids))
-    init_elements!(elements, leaf_cell_ids, mesh, dg.basis)
-
     # Resize parabolic helper variables
+    @unpack viscous_container = cache_parabolic
     resize!(viscous_container, equations, dg, cache)
-
-    # re-initialize interfaces container
-    @unpack interfaces = cache_parabolic
-    resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids))
-    init_interfaces!(interfaces, elements, mesh)
-
-    # re-initialize boundaries container
-    @unpack boundaries = cache_parabolic
-    resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids))
-    init_boundaries!(boundaries, elements, mesh)
+    reinitialize_containers!(mesh, equations, dg, cache_parabolic)
 
     # Sanity check
+    @unpack interfaces = cache_parabolic
     if isperiodic(mesh.tree)
         @assert ninterfaces(interfaces)==1 * nelements(dg, cache_parabolic) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements")
     end
@@ -246,27 +229,13 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
     # actually transferring the solution to the coarsened cells
     coarsen!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_remove)
 
-    # Get new list of leaf cells
-    leaf_cell_ids = local_leaf_cells(mesh.tree)
-
-    @unpack elements, viscous_container = cache_parabolic
-    resize!(elements, length(leaf_cell_ids))
-    init_elements!(elements, leaf_cell_ids, mesh, dg.basis)
-
     # Resize parabolic helper variables
+    @unpack viscous_container = cache_parabolic
     resize!(viscous_container, equations, dg, cache)
-
-    # re-initialize interfaces container
-    @unpack interfaces = cache_parabolic
-    resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids))
-    init_interfaces!(interfaces, elements, mesh)
-
-    # re-initialize boundaries container
-    @unpack boundaries = cache_parabolic
-    resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids))
-    init_boundaries!(boundaries, elements, mesh)
+    reinitialize_containers!(mesh, equations, dg, cache_parabolic)
 
     # Sanity check
+    @unpack interfaces = cache_parabolic
     if isperiodic(mesh.tree)
         @assert ninterfaces(interfaces)==1 * nelements(dg, cache_parabolic) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements")
     end