Support for purely parabolic PDEs via SemidiscretizationParabolic

AUTHORS.md

@@ -38,6 +38,7 @@ are listed in alphabetical order:
 * Lucas Gemein
 * Sven Goldberg
 * Joshua Lampert
+* Tristan Montoya
 * Julia Odenthal
 * Sigrun Ortleb
 * Hendrik Ranocha

NEWS.md

@@ -8,6 +8,8 @@ for human readability.
 ## Changes in the v0.16 lifecycle
 
 #### Added
+- A new semidiscretization type `SemidiscretizationParabolic` has been added to support purely parabolic equations with no hyperbolic part.
+The new equation types `LinearDiffusionEquation1D` and `LinearDiffusionEquation2D` have been implemented to demonstrate this functionality ([#2874]).
 - GPU support extended to include AMD GPU with a buildkite workflow using `TRIXI_TEST=AMDGPU` ([#2834]).
 
 ## Changes when updating to v0.16 from v0.15.x

README.md

@@ -63,7 +63,7 @@ installation and postprocessing procedures. Its features include:
   * Hyperbolic diffusion equations for elliptic problems
   * Lattice-Boltzmann equations (D2Q9 and D3Q27 schemes)
   * Shallow water equations via [TrixiShallowWater.jl](https://github.com/trixi-framework/TrixiShallowWater.jl)
-  * Several scalar conservation laws (e.g., linear advection, Burgers' equation, LWR traffic flow)
+  * Several scalar conservation laws (e.g., linear advection, diffusion/heat equation, Burgers' equation, LWR traffic flow)
 * Multi-physics simulations
   * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics)
 * Shared-memory parallelization via multithreading

examples/tree_1d_dgsem/elixir_diffusion_ldg.jl

@@ -2,15 +2,13 @@ using OrdinaryDiffEqLowStorageRK
 using Trixi
 
 ###############################################################################
-# semidiscretization of the linear (advection) diffusion equation
+# semidiscretization of the (pure) linear diffusion equation
 
-advection_velocity = 0.0 # Note: This renders the equation mathematically purely parabolic
-equations = LinearScalarAdvectionEquation1D(advection_velocity)
 diffusivity() = 0.5
-equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)
+equations = LinearDiffusionEquation1D(diffusivity())
 
-# 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)
+# Create DG solver with polynomial degree = 3
+solver = DGSEM(polydeg = 3)
 
 coordinates_min = -convert(Float64, pi) # minimum coordinate
 coordinates_max = convert(Float64, pi) # maximum coordinate
@@ -36,17 +34,11 @@ function initial_condition_pure_diffusion_1d_convergence_test(x, t,
 end
 initial_condition = initial_condition_pure_diffusion_1d_convergence_test
 
-# define periodic boundary conditions everywhere
-boundary_conditions = boundary_condition_periodic
-boundary_conditions_parabolic = boundary_condition_periodic
-
 # A semidiscretization collects data structures and functions for the spatial discretization
 solver_parabolic = ParabolicFormulationLocalDG()
-semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
-                                             initial_condition,
-                                             solver; solver_parabolic,
-                                             boundary_conditions = (boundary_conditions,
-                                                                    boundary_conditions_parabolic))
+semi = SemidiscretizationParabolic(mesh, equations, initial_condition, solver;
+                                   solver_parabolic = solver_parabolic,
+                                   boundary_conditions = boundary_condition_periodic)
 
 ###############################################################################
 # ODE solvers, callbacks etc.

examples/tree_1d_dgsem/elixir_diffusion_ldg_amr_boundary_layer.jl

@@ -0,0 +1,89 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# Adaptive semidiscretization of the pure diffusion equation with mixed 
+# Dirichlet-Neumann BCs and an initial condition with a boundary layer.
+
+diffusivity = 0.25
+amplitude = 1.0
+boundary_layer_thickness = 0.01
+
+equations = LinearDiffusionEquation1D(diffusivity)
+
+solver = DGSEM(polydeg = 3)
+solver_parabolic = ParabolicFormulationLocalDG()
+
+mesh = TreeMesh((0.0,), (1.0,),
+                initial_refinement_level = 0,
+                periodicity = false,
+                n_cells_max = 30_000)
+
+function initial_condition_boundary_layer(x, t, equations)
+    return SVector(amplitude * (1 - exp(-x[1] / boundary_layer_thickness)) /
+                   (1 - exp(-1 / boundary_layer_thickness)))
+end
+
+initial_condition = initial_condition_boundary_layer
+boundary_condition_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))
+boundary_condition_neumann = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0))
+
+boundary_conditions = (; x_neg = boundary_condition_dirichlet,
+                       x_pos = boundary_condition_neumann)
+
+semi = SemidiscretizationParabolic(mesh, equations, initial_condition, solver;
+                                   solver_parabolic = solver_parabolic,
+                                   boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 200)
+
+alive_callback = AliveCallback(analysis_interval = 200)
+
+amr_indicator = IndicatorLöhner(semi, variable = first)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 3,
+                                      med_threshold = 0.01,
+                                      max_level = 6, max_threshold = 0.175)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 200,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = false)
+
+stepsize_callback = StepsizeCallback(cfl_parabolic = 0.05)
+
+# specify extra node variables to be saved in the `SaveSolutionCallback`
+extra_node_variables = (:dudx,)
+
+# note that using `get_node_variable` to access the gradient exposes Trixi.jl internals 
+# that are not part of the public API, so this usage is not guaranteed to be stable across 
+# releases
+function Trixi.get_node_variable(::Val{:dudx}, u, mesh, equations, dg, cache,
+                                 cache_parabolic)
+    return copy(@view cache_parabolic.parabolic_container.gradients[1, :, :])
+end
+
+# save the solution at the last time step, including the extra node variables
+save_solution = SaveSolutionCallback(dt = 2.0,
+                                     save_initial_solution = false,
+                                     save_final_solution = true,
+                                     solution_variables = cons2cons,
+                                     extra_node_variables = extra_node_variables)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        save_solution, amr_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = stepsize_callback(ode), adaptive = false,
+            ode_default_options()..., callback = callbacks, maxiters = 500_000)

examples/tree_1d_dgsem/elixir_diffusion_ldg_dirichlet.jl

@@ -0,0 +1,60 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the pure diffusion equation
+
+diffusivity() = 0.5
+equations = LinearDiffusionEquation1D(diffusivity())
+
+# Create DG solver with polynomial degree = 3
+solver = DGSEM(polydeg = 3)
+solver_parabolic = ParabolicFormulationLocalDG()
+
+# Create a uniformly refined mesh with nonperiodic boundaries
+mesh = TreeMesh(0.0, 1.0,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000, # set maximum capacity of tree data structure
+                periodicity = false)
+
+function analytical_solution(x, t, equations)
+    scalar = sinpi(x[1]) * exp(-diffusivity() * pi^2 * t)
+    return SVector(scalar)
+end
+initial_condition = analytical_solution
+
+boundary_conditions = (; x_neg = BoundaryConditionDirichlet(initial_condition),
+                       x_pos = BoundaryConditionDirichlet(initial_condition))
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationParabolic(mesh, equations, initial_condition, solver;
+                                   solver_parabolic = solver_parabolic,
+                                   boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+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_callback = AnalysisCallback(semi, interval = 100)
+
+# The AliveCallback prints short status information in regular intervals
+alive_callback = AliveCallback(analysis_interval = 100)
+
+# 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)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+# For CI purposes, we use fixed time-stepping for this elixir.
+sol = solve(ode, RDPK3SpFSAL35(); dt = 1.0e-4, adaptive = false,
+            ode_default_options()..., callback = callbacks)

examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl

@@ -5,15 +5,12 @@ using LinearSolve # For Jacobian-free Newton-Krylov (GMRES) solver
 using ADTypes # For automatic differentiation via finite differences
 
 ###############################################################################
-# semidiscretization of the linear (advection) diffusion equation
+# semidiscretization of the pure diffusion equation
 
-advection_velocity = 0.0 # Note: This renders the equation mathematically purely parabolic
-equations = LinearScalarAdvectionEquation1D(advection_velocity)
 diffusivity() = 0.5
-equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)
+equations = LinearDiffusionEquation1D(diffusivity())
 
-# surface flux does not matter for pure diffusion problem
-solver = DGSEM(polydeg = 3, surface_flux = flux_central)
+solver = DGSEM(polydeg = 3)
 
 coordinates_min = -convert(Float64, pi)
 coordinates_max = convert(Float64, pi)
@@ -34,11 +31,9 @@ end
 initial_condition = initial_condition_pure_diffusion_1d_convergence_test
 
 solver_parabolic = ParabolicFormulationLocalDG()
-semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
-                                             initial_condition,
-                                             solver; solver_parabolic,
-                                             boundary_conditions = (boundary_condition_periodic,
-                                                                    boundary_condition_periodic))
+semi = SemidiscretizationParabolic(mesh, equations, initial_condition, solver;
+                                   solver_parabolic = solver_parabolic,
+                                   boundary_conditions = boundary_condition_periodic)
 
 ###############################################################################
 # ODE solvers, callbacks etc.

examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl

@@ -3,14 +3,9 @@ using Trixi
 ###############################################################################
 
 # Build pure diffusion (Laplace) operator
-advection_velocity = (0, 0)
-equations = LinearScalarAdvectionEquation2D(advection_velocity)
 diffusivity() = 1
-equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
-
-# The hyperbolic flux does not matter for this example since
-# the hyperbolic part is zero.
-solver = DGSEM(polydeg = 5, surface_flux = flux_central)
+equations = LinearDiffusionEquation2D(diffusivity())
+solver = DGSEM(polydeg = 5)
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (1.0, 1.0)
@@ -38,19 +33,15 @@ function bc_sin(x, t, equations)
 end
 bc_sin_dirichlet = BoundaryConditionDirichlet(bc_sin)
 
-# Same boundary conditions for hyperbolic and parabolic part
 boundary_conditions = (; x_neg = bc_homogeneous_dirichlet,
                        y_neg = bc_sin_dirichlet,
                        y_pos = bc_sin_dirichlet,
                        x_pos = bc_homogeneous_dirichlet)
 
-# `solver_parabolic = ParabolicFormulationLocalDG()` strictly required for elliptic/diffusion-dominated problem
-semi = SemidiscretizationHyperbolicParabolic(mesh,
-                                             (equations, equations_parabolic),
-                                             initial_condition, solver;
-                                             solver_parabolic = ParabolicFormulationLocalDG(),
-                                             boundary_conditions = (boundary_conditions,
-                                                                    boundary_conditions))
+# Build the parabolic semidiscretization with a local DG formulation for diffusion
+semi = SemidiscretizationParabolic(mesh, equations, initial_condition, solver;
+                                   solver_parabolic = ParabolicFormulationLocalDG(),
+                                   boundary_conditions = boundary_conditions)
 
 # Note that `linear_structure` does not access the `initial_condition`/steady-state solution
 A_map, b = linear_structure(semi)

src/Trixi.jl

@@ -145,6 +145,7 @@ include("solvers/solvers.jl")
 include("equations/equations_parabolic.jl") # these depend on parabolic solver types
 include("semidiscretization/semidiscretization.jl")
 include("semidiscretization/semidiscretization_hyperbolic.jl")
+include("semidiscretization/semidiscretization_parabolic.jl")
 include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
 include("semidiscretization/semidiscretization_coupled.jl")
@@ -186,7 +187,8 @@ export AcousticPerturbationEquations2D,
 export NonIdealCompressibleEulerEquations1D, NonIdealCompressibleEulerEquations2D
 export IdealGas, VanDerWaals, PengRobinson
 
-export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
+export LinearDiffusionEquation1D, LinearDiffusionEquation2D,
+       LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        LaplaceDiffusionEntropyVariables1D, LaplaceDiffusionEntropyVariables2D,
        LaplaceDiffusionEntropyVariables3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
@@ -295,6 +297,8 @@ export nelements, nnodes, nvariables,
 
 export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate
 
+export SemidiscretizationParabolic
+
 export SemidiscretizationHyperbolicParabolic
 export have_constant_diffusivity, max_diffusivity
 

src/callbacks_step/amr.jl

@@ -213,7 +213,8 @@ end
 end
 
 @inline function (amr_callback::AMRCallback)(u_ode::AbstractVector,
-                                             semi::SemidiscretizationHyperbolicParabolic,
+                                             semi::Union{SemidiscretizationHyperbolicParabolic,
+                                                         SemidiscretizationParabolic},
                                              t, iter;
                                              kwargs...)
     # Note that we don't `wrap_array` the vector `u_ode` to be able to `resize!`
@@ -381,13 +382,13 @@ end
 function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh,
                                      equations, dg::DG,
                                      cache, cache_parabolic,
-                                     semi::SemidiscretizationHyperbolicParabolic,
+                                     semi::Union{SemidiscretizationHyperbolicParabolic,
+                                                 SemidiscretizationParabolic},
                                      t, iter;
                                      only_refine = false, only_coarsen = false)
     @unpack controller, adaptor = amr_callback
 
     u = wrap_array(u_ode, mesh, equations, dg, cache)
-    # Indicator kept based on hyperbolic variables
     lambda = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, cache,
                                                           t = t, iter = iter)
 

src/callbacks_step/save_solution.jl

@@ -34,7 +34,16 @@ end
 and must return an array of dimension
 `(ntuple(_ -> n_nodes, ndims(mesh))..., n_elements)`.
 
-For parabolic-hyperbolic equations `equations_parabolic` and `cache_parabolic` must be added:
+For purely parabolic equations, `cache_parabolic` must be added:
+```julia
+function get_node_variable(::Val{symbol}, u, mesh, equations, dg, cache,
+                           cache_parabolic)
+    # Implementation goes here
+end
+```
+
+For hyperbolic-parabolic equations, `equations_parabolic` and `cache_parabolic` must be 
+added:
 ```julia
 function get_node_variable(::Val{symbol}, u, mesh, equations, dg, cache,
                            equations_parabolic, cache_parabolic)

src/callbacks_step/stepsize.jl

@@ -18,7 +18,8 @@ a ramp-up of the time step.
 
 One can additionally supply a parabolic CFL number `cfl_parabolic` to
 limit the admissible timestep also respecting parabolic restrictions.
-This is only applicable for semidiscretizations of type [`SemidiscretizationHyperbolicParabolic`](@ref).
+This is only applicable for semidiscretizations of type
+[`SemidiscretizationHyperbolicParabolic`](@ref) and [`SemidiscretizationParabolic`](@ref).
 To enable checking for parabolic timestep restrictions, provide a value greater than zero for `cfl_parabolic`.
 By default, `cfl_parabolic` is set to zero which means that only the hyperbolic CFL number `cfl` is considered.
 The keyword argument `cfl_parabolic` must be either a `Real` number, corresponding to a constant
@@ -152,6 +153,17 @@ function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic,
                   solver, cache)
 end
 
+# Case for a purely parabolic semidiscretization
+function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic,
+                      semi::SemidiscretizationParabolic)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+    return cfl_parabolic(t) * max_dt(u, t, mesh,
+                  have_constant_diffusivity(equations), equations,
+                  equations, solver, cache)
+end
+
 # For Euler-Acoustic simulations with `EulerAcousticsCouplingCallback`
 function calculate_dt(u_ode, t, cfl_hyperbolic::Real, cfl_parabolic::Real,
                       semi::AbstractSemidiscretization)