AMR for 2D Parabolic P4est

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()

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()

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,

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

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,

src/callbacks_step/amr.jl

@@ -528,6 +528,103 @@ function copy_to_quad_iter_volume(info, user_data)
     return nothing
 end
 
+# specialized callback which includes the `cache_parabolic` argument
+function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh,
+                                     equations, dg::DG, cache, cache_parabolic,
+                                     semi,
+                                     t, iter;
+                                     only_refine = false, only_coarsen = false,
+                                     passive_args = ())
+    @unpack controller, adaptor = amr_callback
+
+    u = wrap_array(u_ode, mesh, equations, dg, cache)
+    lambda = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, cache,
+                                                          t = t, iter = iter)
+
+    @boundscheck begin
+        @assert axes(lambda)==(Base.OneTo(ncells(mesh)),) ("Indicator array (axes = $(axes(lambda))) and mesh cells (axes = $(Base.OneTo(ncells(mesh)))) have different axes")
+    end
+
+    # Copy controller value of each quad to the quad's user data storage
+    iter_volume_c = cfunction(copy_to_quad_iter_volume, Val(ndims(mesh)))
+
+    # The pointer to lambda will be interpreted as Ptr{Int} above
+    @assert lambda isa Vector{Int}
+    iterate_p4est(mesh.p4est, lambda; iter_volume_c = iter_volume_c)
+
+    @trixi_timeit timer() "refine" if !only_coarsen
+        # Refine mesh
+        refined_original_cells = @trixi_timeit timer() "mesh" refine!(mesh)
+
+        # Refine solver
+        @trixi_timeit timer() "solver" refine!(u_ode, adaptor, mesh, equations, dg,
+                                               cache, cache_parabolic,
+                                               refined_original_cells)
+        for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args
+            @trixi_timeit timer() "passive solver" refine!(p_u_ode, adaptor, p_mesh,
+                                                           p_equations,
+                                                           p_dg, p_cache,
+                                                           refined_original_cells)
+        end
+    else
+        # If there is nothing to refine, create empty array for later use
+        refined_original_cells = Int[]
+    end
+
+    @trixi_timeit timer() "coarsen" if !only_refine
+        # Coarsen mesh
+        coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh)
+
+        # coarsen solver
+        @trixi_timeit timer() "solver" coarsen!(u_ode, adaptor, mesh, equations, dg,
+                                                cache, cache_parabolic,
+                                                coarsened_original_cells)
+        for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args
+            @trixi_timeit timer() "passive solver" coarsen!(p_u_ode, adaptor, p_mesh,
+                                                            p_equations,
+                                                            p_dg, p_cache,
+                                                            coarsened_original_cells)
+        end
+    else
+        # If there is nothing to coarsen, create empty array for later use
+        coarsened_original_cells = Int[]
+    end
+
+    # Store whether there were any cells coarsened or refined and perform load balancing
+    has_changed = !isempty(refined_original_cells) || !isempty(coarsened_original_cells)
+    # Check if mesh changed on other processes
+    if mpi_isparallel()
+        has_changed = MPI.Allreduce!(Ref(has_changed), |, mpi_comm())[]
+    end
+
+    if has_changed # TODO: Taal decide, where shall we set this?
+        # don't set it to has_changed since there can be changes from earlier calls
+        mesh.unsaved_changes = true
+
+        if mpi_isparallel() && amr_callback.dynamic_load_balancing
+            @trixi_timeit timer() "dynamic load balancing" begin
+                global_first_quadrant = unsafe_wrap(Array,
+                                                    unsafe_load(mesh.p4est).global_first_quadrant,
+                                                    mpi_nranks() + 1)
+                old_global_first_quadrant = copy(global_first_quadrant)
+                partition!(mesh)
+                rebalance_solver!(u_ode, mesh, equations, dg, cache,
+                                  old_global_first_quadrant)
+            end
+        end
+
+        reinitialize_boundaries!(semi.boundary_conditions, cache)
+        # if the semidiscretization also stores parabolic boundary conditions, 
+        # reinitialize them after each refinement step as well.
+        if hasproperty(semi, :boundary_conditions_parabolic)
+            reinitialize_boundaries!(semi.boundary_conditions_parabolic, cache)
+        end
+    end
+
+    # Return true if there were any cells coarsened or refined, otherwise false
+    return has_changed
+end
+
 # 2D
 function cfunction(::typeof(copy_to_quad_iter_volume), ::Val{2})
     @cfunction(copy_to_quad_iter_volume, Cvoid,