Diffusive CFL P4est

examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl

@@ -94,6 +94,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, amr
 # run the simulation
 
 # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+ode_alg = RDPK3SpFSAL49()
 time_int_tol = 1.0e-11
-sol = solve(ode, dt = 1e-7, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+sol = solve(ode, dt = 1e-7, ode_alg; abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)

examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl

@@ -157,9 +157,11 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
+# Moderate number of threads (e.g. 4) advisable to speed things up
+ode_alg = RDPK3SpFSAL49(thread = Trixi.True())
 time_int_tol = 1e-7
-sol = solve(ode,
-            # Moderate number of threads (e.g. 4) advisable to speed things up
-            RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, ode_alg;
+            # not necessary, added for overwriting in tests
+            adaptive = true, dt = 1.4e-3,
             abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)

examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl

@@ -0,0 +1,89 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 1 / 17
+advection_velocity = (1.0, 0.0, 0.0)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations)
+
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -0.5, -0.5)
+coordinates_max = (0.0, 0.5, 0.5)
+
+# Affine type mapping to take the [-1,1]^3 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+function mapping(xi, eta, zeta)
+    y_ = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta))
+    x_ = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y_) * cos(0.5 * pi * zeta))
+    z_ = zeta + 1 / 6 * (cos(0.5 * pi * x_) * cos(pi * y_) * cos(0.5 * pi * zeta))
+
+    # Map from [-1, 1]^3 to [-1, -0.5, -0.5] x [0, 0.5, 0.5]
+    x = 0.5 * (x_ - 1)
+    y = 0.5 * y_
+    z = 0.5 * z_
+
+    return SVector(x, y, z)
+end
+
+trees_per_dimension = (5, 3, 3)
+mesh = P4estMesh(trees_per_dimension, polydeg = polydeg,
+                 mapping = mapping, 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() # NOTE: this requires epsilon <= 1/16 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(u)
+end
+initial_condition = initial_condition_eriksson_johnson
+
+boundary_conditions = boundary_condition_default(mesh,
+                                                 BoundaryConditionDirichlet(initial_condition))
+
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6,
+                                     cfl_diffusive = 0.25)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_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);

examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl

@@ -78,6 +78,10 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
+ode_alg = RDPK3SpFSAL49()
 time_int_tol = 1e-8
-sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+sol = solve(ode, ode_alg;
+            # not necessary, added for overwriting in tests
+            adaptive = true, dt = 3.913e-3,
+            abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)

src/callbacks_step/stepsize_dg1d.jl

@@ -8,8 +8,7 @@
 function max_dt(u, t, mesh::TreeMesh{1},
                 constant_speed::False, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
-    # e.g. for steady-state linear advection
+    # Avoid division by zero if the speed vanishes everywhere
     max_scaled_speed = nextfloat(zero(t))
 
     @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
@@ -33,28 +32,31 @@ function max_dt(u, t, mesh::TreeMesh{1},
                 constant_diffusivity::False, equations,
                 equations_parabolic::AbstractEquationsParabolic,
                 dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
-    max_scaled_speed = nextfloat(zero(t))
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
 
-    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
-        max_diffusivity_ = zero(max_scaled_speed)
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity_ = zero(max_scaled_diffusivity)
         for i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, element)
             diffusivity = max_diffusivity(u_node, equations_parabolic)
-            max_diffusivity_ = max(max_diffusivity_, diffusivity)
+            max_diffusivity_ = Base.max(max_diffusivity_, diffusivity)
         end
         inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx
-        max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_)
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          inv_jacobian^2 * max_diffusivity_)
     end
 
     # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2
-    return 4 / (nnodes(dg) * max_scaled_speed)
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
 end
 
 function max_dt(u, t, mesh::TreeMesh{1},
                 constant_speed::True, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -71,32 +73,32 @@ function max_dt(u, t, mesh::TreeMesh{1},
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
-function max_dt(u, t, mesh::TreeMesh,
+function max_dt(u, t, mesh::TreeMesh, # for all dimensions
                 constant_diffusivity::True, equations,
                 equations_parabolic::AbstractEquationsParabolic,
                 dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
-    max_scaled_speed = nextfloat(zero(t))
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
 
     diffusivity = max_diffusivity(equations_parabolic)
 
-    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
         inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx
         # Note: For the currently supported parabolic equations
-        # Diffusion & Navier-Stokes, we only have one diffusivity,
+        # Diffusion & Navier-Stokes, we only have one isotropic diffusivity,
         # so this is valid for 1D, 2D and 3D.
-        max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * diffusivity)
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          inv_jacobian^2 * diffusivity)
     end
 
     # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2
-    return 4 / (nnodes(dg) * max_scaled_speed)
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
 end
 
 function max_dt(u, t, mesh::StructuredMesh{1},
                 constant_speed::False, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
-    # e.g. for steady-state linear advection
+    # Avoid division by zero if the speed vanishes everywhere
     max_scaled_speed = nextfloat(zero(t))
 
     @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
@@ -123,7 +125,7 @@ end
 function max_dt(u, t, mesh::StructuredMesh{1},
                 constant_speed::True, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 

src/callbacks_step/stepsize_dg2d.jl

@@ -7,7 +7,7 @@
 
 function max_dt(u, t, mesh::TreeMesh{2},
                 constant_speed::False, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -33,7 +33,7 @@ function max_dt(u, t, mesh::TreeMesh{2},
                 constant_diffusivity::False, equations,
                 equations_parabolic::AbstractEquationsParabolic,
                 dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
+    # Avoid division by zero if the diffusivity vanishes everywhere
     max_scaled_speed = nextfloat(zero(t))
 
     @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
@@ -55,7 +55,7 @@ end
 
 function max_dt(u, t, mesh::TreeMesh{2},
                 constant_speed::True, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -114,7 +114,7 @@ function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
                             T8codeMesh{2}, StructuredMeshView{2}},
                 constant_speed::False, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -148,13 +148,59 @@ function max_dt(u, t,
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
+function max_dt(u, t,
+                mesh::P4estMesh{2}, # Parabolic terms currently only for `TreeMesh` and `P4estMesh`
+                constant_diffusivity::False, equations,
+                equations_parabolic::AbstractEquationsParabolic,
+                dg::DG, cache)
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
+
+    @unpack contravariant_vectors, inverse_jacobian = cache.elements
+
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity1 = max_diffusivity2 = zero(max_scaled_diffusivity)
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            diffusivity = max_diffusivity(u_node, equations_parabolic)
+
+            # Local diffusivity transformed to the reference element
+            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors,
+                                                  i, j, element)
+            # The metric terms are squared due to the second derivatives in the parabolic terms,
+            # which lead to application of the chain rule (coordinate transformation) twice.
+            # See for instance also the implementation in FLEXI
+            # https://github.com/flexi-framework/flexi/blob/e980e8635e6605daf906fc176c9897314a9cd9b1/src/equations/navierstokes/calctimestep.f90#L75-L77
+            # or FLUXO:
+            # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L130-L132
+            diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2)
+            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors,
+                                                  i, j, element)
+            diffusivity2_transformed = diffusivity * (Ja21^2 + Ja22^2)
+
+            inv_jacobian = abs(inverse_jacobian[i, j, element])
+
+            max_diffusivity1 = Base.max(max_diffusivity1,
+                                        diffusivity1_transformed * inv_jacobian^2)
+            max_diffusivity2 = Base.max(max_diffusivity2,
+                                        diffusivity2_transformed * inv_jacobian^2)
+        end
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          max_diffusivity1 + max_diffusivity2)
+    end
+
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
+end
+
 function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
                             P4estMeshView{2}, T8codeMesh{2}, StructuredMeshView{2}},
                 constant_speed::True, equations, dg::DG, cache)
     @unpack contravariant_vectors, inverse_jacobian = cache.elements
 
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -182,6 +228,51 @@ function max_dt(u, t,
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
+function max_dt(u, t,
+                mesh::P4estMesh{2}, # Parabolic terms currently only for `TreeMesh` and `P4estMesh`
+                constant_diffusivity::True, equations,
+                equations_parabolic::AbstractEquationsParabolic,
+                dg::DG, cache)
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
+
+    diffusivity = max_diffusivity(equations_parabolic)
+
+    @unpack contravariant_vectors, inverse_jacobian = cache.elements
+
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity1 = max_diffusivity2 = zero(max_scaled_diffusivity)
+        for j in eachnode(dg), i in eachnode(dg)
+            # Local diffusivity transformed to the reference element
+            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors,
+                                                  i, j, element)
+            # The metric terms are squared due to the second derivatives in the parabolic terms,
+            # which lead to application of the chain rule (coordinate transformation) twice.
+            # See for instance also the implementation in FLEXI
+            # https://github.com/flexi-framework/flexi/blob/e980e8635e6605daf906fc176c9897314a9cd9b1/src/equations/navierstokes/calctimestep.f90#L75-L77
+            # or FLUXO:
+            # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L130-L132
+            diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2)
+            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors,
+                                                  i, j, element)
+            diffusivity2_transformed = diffusivity * (Ja21^2 + Ja22^2)
+
+            inv_jacobian = abs(inverse_jacobian[i, j, element])
+
+            max_diffusivity1 = Base.max(max_diffusivity1,
+                                        diffusivity1_transformed * inv_jacobian^2)
+            max_diffusivity2 = Base.max(max_diffusivity2,
+                                        diffusivity2_transformed * inv_jacobian^2)
+        end
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          max_diffusivity1 + max_diffusivity2)
+    end
+
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
+end
+
 function max_dt(u, t, mesh::ParallelP4estMesh{2},
                 constant_speed::False, equations, dg::DG, cache)
     # call the method accepting a general `mesh::P4estMesh{2}`