MPI P4est parabolic 3D

examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr_mortarTestcase.jl

@@ -0,0 +1,106 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+#This testcases uses a artificial refinement rule to explicitly test refinement and unrifinement. 
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+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_parabolic) -> SVector(1.0, 0.0))
+velocity_bc_cavity = NoSlip((x, t, equations_parabolic) -> SVector(0.0, 0.0))
+heat_bc = Adiabatic((x, t, equations_parabolic) -> 0.0)
+boundary_condition_lid = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
+boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)
+
+boundary_conditions = (; 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 = (; 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)
+
+function positional_rule(u, x, t)
+    xtar = sin(π * t)
+    if ((x[1] < xtar) && (x[2] < xtar))
+        return 1.0   # refine
+    else
+        return 0.0   # keep coarse
+    end
+end
+
+amr_indicator = IndicatorNodalFunction(positional_rule, semi)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator;
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.5,
+                                      max_level = 2, max_threshold = 0.9)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 50,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+callbacks = CallbackSet(save_solution, summary_callback, alive_callback, analysis_callback,
+                        amr_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)

examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr_mortar.jl

@@ -0,0 +1,129 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+# Why: This testcases adds uneven refinement in the taylor green vortex (TGV) to specifically test (MPI) mortar treatment in the parabolic part.
+
+###############################################################################
+# Physics
+
+prandtl_number() = 0.72
+mu = 6.25e-4  # Re ≈ 1600
+
+equations = CompressibleEulerEquations3D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
+                                                          Prandtl = prandtl_number())
+
+###############################################################################
+# Initial condition (Taylor-Green vortex)
+
+function initial_condition_taylor_green_vortex(x, t,
+                                               equations::CompressibleEulerEquations3D)
+    A = 1.0
+    Ms = 0.1
+
+    rho = 1.0
+    v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3])
+    v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
+    v3 = 0.0
+
+    p = (A / Ms)^2 * rho / equations.gamma
+    p += 1.0 / 16.0 * A^2 * rho *
+         (cos(2x[1]) * cos(2x[3]) +
+          2cos(2x[2]) + 2cos(2x[1]) +
+          cos(2x[2]) * cos(2x[3]))
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+initial_condition = initial_condition_taylor_green_vortex
+
+###############################################################################
+# DG setup
+
+volume_flux = flux_ranocha
+
+solver = DGSEM(polydeg = 3,
+               surface_flux = flux_lax_friedrichs,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Mesh
+
+coordinates_min = (-1.0, -1.0, -1.0) .* pi
+coordinates_max = (1.0, 1.0, 1.0) .* pi
+
+trees_per_dimension = (2, 2, 2)
+
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3,
+                 coordinates_min = coordinates_min,
+                 coordinates_max = coordinates_max,
+                 periodicity = (true, true, true),
+                 initial_refinement_level = 0)
+
+###############################################################################
+# Semidiscretization
+
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition,
+                                             solver;
+                                             boundary_conditions = (boundary_condition_periodic,
+                                                                    boundary_condition_periodic))
+
+###############################################################################
+# AMR: positional rule → creates 2 refinement bocks on the diagonal of the cube 
+# Used to explicitly include mortars in the parabolic MPI testscases
+
+function positional_rule(u, x, t)
+    if ((x[1] < 0) && (x[2] < 0) && (x[3] < 0)) ||
+       ((x[1] > 0) && (x[2] > 0) && (x[3] > 0))
+        return 1.0   # refine
+    else
+        return 0.0   # keep coarse
+    end
+end
+
+amr_indicator = IndicatorNodalFunction(positional_rule, semi)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator;
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.5,
+                                      max_level = 2, max_threshold = 0.9)
+
+# refine only once to get a static mortar mesh
+amr_callback = AMRCallback(semi,
+                           amr_controller;
+                           interval = typemax(Int),                 # no further AMR
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+###############################################################################
+# Callbacks 
+analysis_interval = 50
+analysis_callback = AnalysisCallback(semi,
+                                     interval = analysis_interval,  # effectively disabled for short runs
+                                     save_analysis = true,
+                                     extra_analysis_integrals = (energy_kinetic,
+                                                                 energy_internal,
+                                                                 enstrophy))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(analysis_callback,
+                        alive_callback,
+                        amr_callback)
+
+###############################################################################
+# Time integration
+
+tspan = (0.0, 5)
+
+ode = semidiscretize(semi, tspan)
+
+sol = solve(ode,
+            RDPK3SpFSAL49();
+            abstol = 1e-8,
+            reltol = 1e-8,
+            ode_default_options()...,
+            callback = callbacks)

examples/tree_1d_dgsem/elixir_advection_amr.jl

@@ -39,7 +39,9 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+rule(u, x, t) = u[1]
+
+amr_controller = ControllerThreeLevel(semi, IndicatorNodalFunction(rule, semi),
                                       base_level = 4,
                                       med_level = 5, med_threshold = 0.1,
                                       max_level = 6, max_threshold = 0.6)

src/Trixi.jl

@@ -322,7 +322,7 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!
 
 export ControllerThreeLevel, ControllerThreeLevelCombined,
-       IndicatorLöhner, IndicatorLoehner, IndicatorMax
+       IndicatorLöhner, IndicatorLoehner, IndicatorMax, IndicatorNodalFunction
 
 export PositivityPreservingLimiterZhangShu, EntropyBoundedLimiter
 

src/callbacks_step/amr.jl

@@ -633,6 +633,8 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh,
                 partition!(mesh)
                 rebalance_solver!(u_ode, mesh, equations, dg, cache,
                                   old_global_first_quadrant)
+                @unpack parabolic_container = cache_parabolic
+                resize!(parabolic_container, equations, dg, cache)
             end
         end
 

src/callbacks_step/analysis_surface_integral_2d.jl

@@ -250,7 +250,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t,
             j_node += j_node_step
         end
     end
-    return surface_integral
+    return distribute_surface_integral(surface_integral, mesh)
 end
 
 # 2D version of the `analyze` function for `AnalysisSurfaceIntegral` of viscous, i.e.,
@@ -320,6 +320,19 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t,
             j_node += j_node_step
         end
     end
-    return surface_integral
+    return distribute_surface_integral(surface_integral, mesh)
+end
+
+# Serial/default: do nothing
+distribute_surface_integral(val, mesh) = val
+
+# Parallel: sum over all ranks
+function distribute_surface_integral(val,
+                                     mesh::Union{P4estMeshParallel{2},
+                                                 T8codeMeshParallel{2}})
+    comm = MPI.COMM_WORLD
+    buf = [val]
+    MPI.Allreduce!(buf, MPI.SUM, comm)
+    return buf[1]
 end
 end # muladd

src/solvers/dgsem/indicators.jl

@@ -489,4 +489,42 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
     return nothing
 end
+
+"""
+    IndicatorNodalFunction(f)
+
+Create an AMR indicator from a solution, space and time dependent indicator function `f(u, x, t)`.
+The function `f` is evaluated at the nodal points. The maximum of `f`` over all nodes in each element is used as indicator for the element.
+The function can be solution independent allowing for user-defined mesh refinement/coarsening varying in space and time, similar to the `refinement_patches` keyword for the [`Treemesh`](@ref).
+"""
+struct IndicatorNodalFunction{F, Cache} <: AbstractIndicator
+    indicator_function::F
+    cache::Cache
+end
+
+# this method is used when the indicator is constructed as for AMR
+function IndicatorNodalFunction(indicator_function, semi::AbstractSemidiscretization)
+    cache = create_cache(IndicatorNodalFunction, semi)
+    return IndicatorNodalFunction{typeof(indicator_function), typeof(cache)}(indicator_function,
+                                                                             cache)
+end
+
+function Base.show(io::IO, indicator::IndicatorNodalFunction)
+    @nospecialize indicator # reduce precompilation time
+    print(io, "IndicatorNodalFunction")
+    print(io, indicator.indicator_function)
+    return nothing
+end
+
+function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorNodalFunction)
+    @nospecialize indicator # reduce precompilation time
+    if get(io, :compact, false)
+        show(io, indicator)
+    else
+        setup = [
+            "Indicator Function" => indicator.indicator_function
+        ]
+        summary_box(io, "IndicatorNodalFunction", setup)
+    end
+end
 end # @muladd