AMR for 2D Parabolic P4est (#1691)

* 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 * look into p4est amr parabolic * platform independent tests * Working on mortars for parabolic p4est * print analysis every step * No need for different flux_viscous comps after adding container_viscous to p4est * Laplace 3d * coarsen also for p4est mesh * Eliminate unused stuff * Longer times to allow converage to hit coarsen! * Increase testing of Laplace 3D * Add tests for velocities * Continue working on AMR p4est 2d * Note possible remedy for paraboic p4est mortars * remove comment * remove some comments * clean up attempts * add elixir for amr testing * adding commented out mortar routines in 2D * Adding Mortar to 2d dg parabolic term * remove testing snippet * fix comments * add more arguments for dispatch * add some temporary todo notes * some updates for AP and KS * specialize mortar_fluxes_to_elements * BUGFIX: apply_jacobian_parabolic! was incorrect for P4estMesh * fixed rhs_parabolic! for mortars * more changes to elixir * indexing bug * comments * Adding the example for nonperiodic BCs with amr * hopefully this fixes AMR boundaries for parabolic terms * add elixir * Example with non periodic bopundary conditions * remove cruft * clean up merge nonsense * restore examples * restore project * restore * more similar * more similar to standard case * correct amr file * delete doubled tests * remove unused * remove duplicate test * format * BCs are not periodic * straighten examples, provide tests * format * fmt * reduce test time intervals * shorten test * test vals updated ode * cache_parabolic does not require mortars * add couple comments * Update src/callbacks_step/amr.jl Co-authored-by: Jesse Chan <[email protected]> * Update src/callbacks_step/amr.jl Co-authored-by: Jesse Chan <[email protected]> * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl Co-authored-by: Jesse Chan <[email protected]> * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl Co-authored-by: Jesse Chan <[email protected]> * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl Co-authored-by: Jesse Chan <[email protected]> * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl Co-authored-by: Jesse Chan <[email protected]> * fmt * add news entry * Update src/callbacks_step/amr.jl Co-authored-by: Erik Faulhaber <[email protected]> * typo --------- Co-authored-by: Michael Schlottke-Lakemper <[email protected]> Co-authored-by: Jesse Chan <[email protected]> Co-authored-by: Jesse Chan <[email protected]> Co-authored-by: Ahmad Peyvan <[email protected]> Co-authored-by: Erik Faulhaber <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Nov 21, 2023 · 3c1762b · 3c1762b
1 parent 9317f7c
commit 3c1762b
Show file tree

Hide file tree

Showing 11 changed files with 690 additions and 21 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -7,6 +7,7 @@ for human readability.
 ## Changes when updating to v0.6 from v0.5.x
 
 #### Added
+- AMR for hyperbolic-parabolic equations on 2D `P4estMesh`
 
 #### Changed
 

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl
@@ -0,0 +1,98 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 5.0e-2
+advection_velocity = (1.0, 0.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+equations_parabolic = LaplaceDiffusion2D(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) # minimum coordinates (min(x), min(y))
+coordinates_max = (0.0, 0.5) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 2,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = false)
+
+# 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 = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
+                           :y_neg => BoundaryConditionDirichlet(initial_condition),
+                           :y_pos => BoundaryConditionDirichlet(initial_condition),
+                           :x_pos => boundary_condition_do_nothing)
+
+boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
+                                     :x_pos => BoundaryConditionDirichlet(initial_condition),
+                                     :y_neg => BoundaryConditionDirichlet(initial_condition),
+                                     :y_pos => 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 = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 1,
+                                      med_level = 2, med_threshold = 0.9,
+                                      max_level = 3, max_threshold = 1.0)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 50)
+
+# 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, amr_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, dt = 1e-7, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+            ode_default_options()..., callback = callbacks)
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_amr.jl
@@ -0,0 +1,83 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+advection_velocity = (1.5, 1.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+diffusivity() = 5.0e-2
+equations_parabolic = LaplaceDiffusion2D(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, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 1,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max)
+
+# Define initial condition
+function initial_condition_diffusive_convergence_test(x, t,
+                                                      equation::LinearScalarAdvectionEquation2D)
+    # 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
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# 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)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 1,
+                                      med_level = 2, med_threshold = 1.25,
+                                      max_level = 3, max_threshold = 1.45)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 20)
+
+# 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, amr_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()
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
@@ -41,7 +41,6 @@ heat_bc = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
 boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)
 
-# define periodic boundary conditions everywhere
 boundary_conditions = Dict(:x_neg => boundary_condition_slip_wall,
                            :y_neg => boundary_condition_slip_wall,
                            :y_pos => boundary_condition_slip_wall,

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl
@@ -0,0 +1,94 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 0.001
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+                                                          Prandtl = prandtl_number())
+
+# 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, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh
+trees_per_dimension = (6, 6)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 2,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, false))
+
+function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D)
+    Ma = 0.1
+    rho = 1.0
+    u, v = 0.0, 0.0
+    p = 1.0 / (Ma^2 * equations.gamma)
+    return prim2cons(SVector(rho, u, v, p), equations)
+end
+initial_condition = initial_condition_cavity
+
+# BC types
+velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0))
+velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0))
+heat_bc = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
+boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)
+
+boundary_conditions = Dict(:x_neg => boundary_condition_slip_wall,
+                           :y_neg => boundary_condition_slip_wall,
+                           :y_pos => boundary_condition_slip_wall,
+                           :x_pos => boundary_condition_slip_wall)
+
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_cavity,
+                                     :y_neg => boundary_condition_cavity,
+                                     :y_pos => boundary_condition_lid,
+                                     :x_pos => boundary_condition_cavity)
+
+# 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, 25.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval = 2000)
+analysis_interval = 2000
+analysis_callback = AnalysisCallback(semi, 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.005,
+                                      max_level = 2, max_threshold = 0.01)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 50,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, amr_callback)
+# callbacks = CallbackSet(summary_callback, alive_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)
+
+summary_callback() # print the timer summary
diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
@@ -40,7 +40,6 @@ heat_bc = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
 boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)
 
-# define periodic boundary conditions everywhere
 boundary_conditions = boundary_condition_slip_wall
 
 boundary_conditions_parabolic = (; x_neg = boundary_condition_cavity,