MPI Implementation for hyperbolic-parabolic system and the P4est Mesh

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 refinment 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 refinment bocks on the diagonal of the cube 
+# Used to explicitly include mortars in the parabolic MPI testscases
+
+function positional_rule(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 = IndicatorPositional(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)

src/Trixi.jl

@@ -321,7 +321,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, IndicatorPositional
 
 export PositivityPreservingLimiterZhangShu, EntropyBoundedLimiter
 

src/callbacks_step/amr.jl

@@ -633,6 +633,46 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh,
                 partition!(mesh)
                 rebalance_solver!(u_ode, mesh, equations, dg, cache,
                                   old_global_first_quadrant)
+
+                # Resync parabolic cache after repartitioning
+                @unpack viscous_container = cache_parabolic
+                resize!(viscous_container, equations, dg, cache)
+
+                if hasproperty(cache_parabolic, :elements)
+                    nelements = size(cache.elements.inverse_jacobian, ndims(mesh) + 1)
+                    resize!(cache_parabolic.elements, nelements)
+                    copyto!(cache_parabolic.elements.node_coordinates,
+                            cache.elements.node_coordinates)
+                    copyto!(cache_parabolic.elements.jacobian_matrix,
+                            cache.elements.jacobian_matrix)
+                    copyto!(cache_parabolic.elements.contravariant_vectors,
+                            cache.elements.contravariant_vectors)
+                    copyto!(cache_parabolic.elements.inverse_jacobian,
+                            cache.elements.inverse_jacobian)
+                    copyto!(cache_parabolic.elements.surface_flux_values,
+                            cache.elements.surface_flux_values)
+                end
+
+                if hasproperty(cache_parabolic, :interfaces)
+                    ninterfaces = size(cache.interfaces.neighbor_ids, 2)
+                    resize!(cache_parabolic.interfaces, ninterfaces)
+                    copyto!(cache_parabolic.interfaces.u, cache.interfaces.u)
+                    copyto!(cache_parabolic.interfaces.neighbor_ids,
+                            cache.interfaces.neighbor_ids)
+                    copyto!(cache_parabolic.interfaces.node_indices,
+                            cache.interfaces.node_indices)
+                end
+
+                if hasproperty(cache_parabolic, :boundaries)
+                    nboundaries = length(cache.boundaries.neighbor_ids)
+                    resize!(cache_parabolic.boundaries, nboundaries)
+                    copyto!(cache_parabolic.boundaries.u, cache.boundaries.u)
+                    copyto!(cache_parabolic.boundaries.neighbor_ids,
+                            cache.boundaries.neighbor_ids)
+                    copyto!(cache_parabolic.boundaries.node_indices,
+                            cache.boundaries.node_indices)
+                    copyto!(cache_parabolic.boundaries.name, cache.boundaries.name)
+                end
             end
         end
 

src/callbacks_step/amr_dg2d.jl

@@ -193,6 +193,37 @@ function refine!(u_ode::AbstractVector, adaptor,
     @unpack viscous_container = cache_parabolic
     resize!(viscous_container, equations, dg, cache)
 
+    # Keep parabolic topology containers in sync if they exist (currently only for P4estMesh thus the conditional
+    if hasproperty(cache_parabolic, :elements)
+        nelements = size(cache.elements.inverse_jacobian, ndims(mesh) + 1)
+        resize!(cache_parabolic.elements, nelements)
+
+        copyto!(cache_parabolic.elements.node_coordinates,
+                cache.elements.node_coordinates)
+        copyto!(cache_parabolic.elements.jacobian_matrix,
+                cache.elements.jacobian_matrix)
+        copyto!(cache_parabolic.elements.contravariant_vectors,
+                cache.elements.contravariant_vectors)
+        copyto!(cache_parabolic.elements.inverse_jacobian,
+                cache.elements.inverse_jacobian)
+        copyto!(cache_parabolic.elements.surface_flux_values,
+                cache.elements.surface_flux_values)
+    end
+    if hasproperty(cache_parabolic, :interfaces)
+        resize!(cache_parabolic.interfaces, size(cache.interfaces.neighbor_ids, 2))
+        copyto!(cache_parabolic.interfaces.u, cache.interfaces.u)
+        copyto!(cache_parabolic.interfaces.neighbor_ids, cache.interfaces.neighbor_ids)
+        copyto!(cache_parabolic.interfaces.node_indices, cache.interfaces.node_indices)
+    end
+
+    if hasproperty(cache_parabolic, :boundaries)
+        resize!(cache_parabolic.boundaries, length(cache.boundaries.neighbor_ids))
+        copyto!(cache_parabolic.boundaries.u, cache.boundaries.u)
+        copyto!(cache_parabolic.boundaries.neighbor_ids, cache.boundaries.neighbor_ids)
+        copyto!(cache_parabolic.boundaries.node_indices, cache.boundaries.node_indices)
+        copyto!(cache_parabolic.boundaries.name, cache.boundaries.name)
+    end
+
     return nothing
 end
 
@@ -387,6 +418,36 @@ function coarsen!(u_ode::AbstractVector, adaptor,
     # Resize parabolic helper variables
     @unpack viscous_container = cache_parabolic
     resize!(viscous_container, equations, dg, cache)
+    # Keep parabolic topology containers in sync if they exist (currently only for P4estMesh thus the conditional
+    if hasproperty(cache_parabolic, :elements)
+        nelements = size(cache.elements.inverse_jacobian, ndims(mesh) + 1)
+        resize!(cache_parabolic.elements, nelements)
+
+        copyto!(cache_parabolic.elements.node_coordinates,
+                cache.elements.node_coordinates)
+        copyto!(cache_parabolic.elements.jacobian_matrix,
+                cache.elements.jacobian_matrix)
+        copyto!(cache_parabolic.elements.contravariant_vectors,
+                cache.elements.contravariant_vectors)
+        copyto!(cache_parabolic.elements.inverse_jacobian,
+                cache.elements.inverse_jacobian)
+        copyto!(cache_parabolic.elements.surface_flux_values,
+                cache.elements.surface_flux_values)
+    end
+    if hasproperty(cache_parabolic, :interfaces)
+        resize!(cache_parabolic.interfaces, size(cache.interfaces.neighbor_ids, 2))
+        copyto!(cache_parabolic.interfaces.u, cache.interfaces.u)
+        copyto!(cache_parabolic.interfaces.neighbor_ids, cache.interfaces.neighbor_ids)
+        copyto!(cache_parabolic.interfaces.node_indices, cache.interfaces.node_indices)
+    end
+
+    if hasproperty(cache_parabolic, :boundaries)
+        resize!(cache_parabolic.boundaries, length(cache.boundaries.neighbor_ids))
+        copyto!(cache_parabolic.boundaries.u, cache.boundaries.u)
+        copyto!(cache_parabolic.boundaries.neighbor_ids, cache.boundaries.neighbor_ids)
+        copyto!(cache_parabolic.boundaries.node_indices, cache.boundaries.node_indices)
+        copyto!(cache_parabolic.boundaries.name, cache.boundaries.name)
+    end
 
     return nothing
 end

src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl

@@ -91,6 +91,7 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple,
 
     cache_parabolic = create_cache_parabolic(mesh, equations, solver,
                                              nelements(solver, cache), uEltype)
+    cache_parabolic = (; cache..., cache_parabolic...)
 
     _boundary_conditions_parabolic = digest_boundary_conditions(boundary_conditions_parabolic,
                                                                 mesh, solver, cache)
@@ -173,6 +174,12 @@ function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic)
         print(io, key)
     end
     print(io, "))")
+    print(io, ", cacheParabolic(")
+    for (idx, key) in enumerate(keys(semi.cache_parabolic))
+        idx > 1 && print(io, " ")
+        print(io, key)
+    end
+    print(io, "))")
     return nothing
 end
 

src/solvers/dgsem/indicators.jl

@@ -489,4 +489,44 @@
     end
     return nothing
 end
+
+"""
+    IndicatorPositional(f)
+
+Create an AMR indicator from a positional rule `f(x, t)`. The usecase for this rule is explict testing and or targeted "static" refinment 
+based on a functional rule. The function `f` is called with the element center `x::SVector{ndims,Float64}` and time `t`.
+"""
+struct IndicatorPositional{F, Cache} <: AbstractIndicator
+    rule::F
+    cache::Cache
+end
+
+# this method is used when the indicator is constructed as for AMR
+function IndicatorPositional(rule, semi::AbstractSemidiscretization)
+    cache = create_cache(IndicatorPositional, semi)
+    return IndicatorPositional{typeof(rule), typeof(cache)}(rule, cache)
+end
+
+function (positional::IndicatorPositional)(u, mesh, equations, dg, cache; t = 0.0,
+                                           iter = 0)
+    x = cache.elements.node_coordinates
+    @unpack alpha, center_threaded = positional.cache
+    resize!(alpha, nelements(dg, cache))
+    #### @threaded does not work with cfunction (positional.rule) thus Threads.@threads for now
+    Threads.@threads for element in Trixi.eachelement(dg, cache)
+        center = center_threaded[Threads.threadid()]
+        fill!(center, 0.0)
+        n = 0
+
+        for k in Trixi.eachnode(dg), j in Trixi.eachnode(dg), i in Trixi.eachnode(dg)
+            center[1] += x[1, i, j, k, element]
+            center[2] += x[2, i, j, k, element]
+            center[3] += x[3, i, j, k, element]
+            n += 1
+        end
+        center ./= n
+        alpha[element] = positional.rule(center, t)
+    end
+    return alpha
+end
 end # @muladd

src/solvers/dgsem_p4est/dg.jl

@@ -68,6 +68,52 @@ function create_cache(mesh::P4estMeshView, equations::AbstractEquations, dg::DG,
     return cache
 end
 
+function create_cache_parabolic(mesh::P4estMesh,
+                                equations_hyperbolic::AbstractEquations,
+                                dg::DG, n_elements, uEltype)
+    # Keep consistent with the baseline p4est cache
+    balance!(mesh)
+
+    elements = init_elements(mesh, equations_hyperbolic, dg.basis, uEltype)
+    interfaces = init_interfaces(mesh, equations_hyperbolic, dg.basis, elements)
+    boundaries = init_boundaries(mesh, equations_hyperbolic, dg.basis, elements)
+    viscous_container = init_viscous_container_3d(nvariables(equations_hyperbolic),
+                                                  nnodes(dg), n_elements,
+                                                  uEltype)
+
+    return (; elements, interfaces, boundaries, viscous_container)
+end
+
+function create_cache_parabolic(mesh::P4estMeshView,
+                                equations_hyperbolic::AbstractEquations,
+                                dg::DG, n_elements, uEltype)
+    balance!(mesh.parent)
+
+    elements_parent = init_elements(mesh.parent, equations_hyperbolic, dg.basis,
+                                    uEltype)
+    interfaces_parent = init_interfaces(mesh.parent, equations_hyperbolic, dg.basis,
+                                        elements_parent)
+    boundaries_parent = init_boundaries(mesh.parent, equations_hyperbolic, dg.basis,
+                                        elements_parent)
+    mortars_parent = init_mortars(mesh.parent, equations_hyperbolic, dg.basis,
+                                  elements_parent)
+
+    elements, interfaces, boundaries, mortars, neighbor_ids_parent = extract_p4est_mesh_view(elements_parent,
+                                                                                             interfaces_parent,
+                                                                                             boundaries_parent,
+                                                                                             mortars_parent,
+                                                                                             mesh,
+                                                                                             equations_hyperbolic,
+                                                                                             dg,
+                                                                                             uEltype)
+
+    viscous_container = init_viscous_container_3d(nvariables(equations_hyperbolic),
+                                                  nnodes(dg), n_elements,
+                                                  uEltype)
+
+    return (; elements, interfaces, boundaries, neighbor_ids_parent, viscous_container)
+end
+
 # Extract outward-pointing normal direction
 # (contravariant vector ±Ja^i, i = index)
 # Note that this vector is not normalized
@@ -91,6 +137,7 @@ include("dg_2d_parabolic.jl")
 
 include("dg_3d.jl")
 include("dg_3d_parabolic.jl")
+include("dg_3d_parabolic_parallel.jl")
 include("dg_parallel.jl")
 
 # Subcell limiters