diff --git a/AUTHORS.md b/AUTHORS.md index eb56a6aa0ba..b54b8c45f41 100644 --- a/AUTHORS.md +++ b/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 diff --git a/NEWS.md b/NEWS.md index 1e560b8b287..de2de6efba6 100644 --- a/NEWS.md +++ b/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]). - A new AMR indicator `IndicatorNodalFunction` is introduced, which allows AMR depending on the solution, space, and time. This can be useful, for example, for testing AMR implementations, but also when the solution behavior is known a priori ([#2881]). - GPU support extended to include AMD GPU with a buildkite workflow using `TRIXI_TEST=AMDGPU` ([#2834]). diff --git a/README.md b/README.md index 096efbfd016..9d929b4969c 100644 --- a/README.md +++ b/README.md @@ -64,6 +64,7 @@ installation and postprocessing procedures. Its features include: * 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) + * Linear diffusion/heat equation * Multi-physics simulations * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics) * Shared-memory parallelization via multithreading diff --git a/docs/src/index.md b/docs/src/index.md index 4792e97dc57..1a047c23fb4 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -59,6 +59,7 @@ installation and postprocessing procedures. Its features include: * 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) + * Linear diffusion/heat equation * Multi-physics simulations * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics) * Shared-memory parallelization via multithreading diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl index 8f79de2d3c4..3f3b060838e 100644 --- a/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl +++ b/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. diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg_amr_boundary_layer.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg_amr_boundary_layer.jl new file mode 100644 index 00000000000..3fa9683d7c8 --- /dev/null +++ b/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) diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg_dirichlet.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg_dirichlet.jl new file mode 100644 index 00000000000..bb273463c12 --- /dev/null +++ b/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) diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl index f52b62738cd..5270c9e4dbe 100644 --- a/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl +++ b/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. diff --git a/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl b/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl index 51530f3286f..af3c12aeab4 100644 --- a/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl +++ b/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) diff --git a/src/Trixi.jl b/src/Trixi.jl index 1d8098ac283..b2031942e18 100644 --- a/src/Trixi.jl +++ b/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 diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 6931ec2a645..1fc3a591df1 100644 --- a/src/callbacks_step/amr.jl +++ b/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) diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl index 68110a3b257..62c1973d169 100644 --- a/src/callbacks_step/save_solution.jl +++ b/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) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 9c2b9020cff..a3258a147ef 100644 --- a/src/callbacks_step/stepsize.jl +++ b/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) diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 11d56ba2525..f57ed103247 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -11,4 +11,6 @@ include("laplace_diffusion.jl") include("laplace_diffusion_entropy_variables.jl") +include("linear_diffusion_equation.jl") + include("compressible_navier_stokes.jl") diff --git a/src/equations/laplace_diffusion_1d.jl b/src/equations/laplace_diffusion_1d.jl index c340c673a2e..007d080097d 100644 --- a/src/equations/laplace_diffusion_1d.jl +++ b/src/equations/laplace_diffusion_1d.jl @@ -3,6 +3,9 @@ `LaplaceDiffusion1D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` with diffusivity ``\kappa`` applied to each solution component defined by `equations`. +This is intended for use as the parabolic part of a hyperbolic-parabolic system, where the +hyperbolic part is defined by `equations`. For a purely parabolic diffusion equation +without any hyperbolic part, see [`LinearDiffusionEquation1D`](@ref). """ struct LaplaceDiffusion1D{E, N, T} <: AbstractLaplaceDiffusion{1, N} diffusivity::T diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index c243f3de364..a740cbda31b 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -3,6 +3,9 @@ `LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` with diffusivity ``\kappa`` applied to each solution component defined by `equations`. +This is intended for use as the parabolic part of a hyperbolic-parabolic system, where the +hyperbolic part is defined by `equations`. For a purely parabolic diffusion equation +without any hyperbolic part, see [`LinearDiffusionEquation2D`](@ref). """ struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusion{2, N} diffusivity::T diff --git a/src/equations/linear_diffusion_equation.jl b/src/equations/linear_diffusion_equation.jl new file mode 100644 index 00000000000..41550c3628c --- /dev/null +++ b/src/equations/linear_diffusion_equation.jl @@ -0,0 +1,2 @@ +include("linear_diffusion_equation_1d.jl") +include("linear_diffusion_equation_2d.jl") diff --git a/src/equations/linear_diffusion_equation_1d.jl b/src/equations/linear_diffusion_equation_1d.jl new file mode 100644 index 00000000000..288709decc4 --- /dev/null +++ b/src/equations/linear_diffusion_equation_1d.jl @@ -0,0 +1,39 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + LinearDiffusionEquation1D(diffusivity) + +The linear diffusion equation (or heat equation) in one space dimension with constant `diffusivity` ``\kappa``: +```math +\partial_t u = \partial_1 \left( \kappa \partial_1 u \right). +``` +Unlike [`LaplaceDiffusion1D`](@ref), which represents the parabolic part of a +hyperbolic-parabolic equation, `LinearDiffusionEquation1D` represents a purely parabolic +equation without any hyperbolic part. +""" +struct LinearDiffusionEquation1D{RealT <: Real} <: AbstractLaplaceDiffusion{1, 1} + diffusivity::RealT +end + +varnames(::typeof(cons2cons), ::LinearDiffusionEquation1D) = ("scalar",) +varnames(::typeof(cons2prim), ::LinearDiffusionEquation1D) = ("scalar",) +varnames(::typeof(cons2entropy), ::LinearDiffusionEquation1D) = ("scalar",) + +@inline cons2prim(u, equations::LinearDiffusionEquation1D) = u +@inline cons2entropy(u, equations::LinearDiffusionEquation1D) = u + +@inline entropy(u::Real, ::LinearDiffusionEquation1D) = 0.5f0 * u^2 +@inline entropy(u, equations::LinearDiffusionEquation1D) = entropy(u[1], equations) + +@inline function flux(u, gradients, orientation::Integer, + equations::LinearDiffusionEquation1D) + dudx, = gradients + # orientation == 1 + return equations.diffusivity * dudx +end +end # @muladd diff --git a/src/equations/linear_diffusion_equation_2d.jl b/src/equations/linear_diffusion_equation_2d.jl new file mode 100644 index 00000000000..31bec47bc0c --- /dev/null +++ b/src/equations/linear_diffusion_equation_2d.jl @@ -0,0 +1,43 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + LinearDiffusionEquation2D(diffusivity) + +The linear diffusion equation (or heat equation) in two space dimensions with constant `diffusivity` ``\kappa``: +```math +\partial_t u = \partial_1 \left( \kappa \partial_1 u \right) + + \partial_2 \left( \kappa \partial_2 u \right). +``` +Unlike [`LaplaceDiffusion2D`](@ref), which represents the parabolic part of a +hyperbolic-parabolic equation, `LinearDiffusionEquation2D` represents a purely parabolic +equation without any hyperbolic part. +""" +struct LinearDiffusionEquation2D{RealT <: Real} <: AbstractLaplaceDiffusion{2, 1} + diffusivity::RealT +end + +varnames(::typeof(cons2cons), ::LinearDiffusionEquation2D) = ("scalar",) +varnames(::typeof(cons2prim), ::LinearDiffusionEquation2D) = ("scalar",) +varnames(::typeof(cons2entropy), ::LinearDiffusionEquation2D) = ("scalar",) + +@inline cons2prim(u, equations::LinearDiffusionEquation2D) = u +@inline cons2entropy(u, equations::LinearDiffusionEquation2D) = u + +@inline entropy(u::Real, ::LinearDiffusionEquation2D) = 0.5f0 * u^2 +@inline entropy(u, equations::LinearDiffusionEquation2D) = entropy(u[1], equations) + +@inline function flux(u, gradients, orientation::Integer, + equations::LinearDiffusionEquation2D) + dudx, dudy = gradients + if orientation == 1 + return SVector(equations.diffusivity * dudx) + else # if orientation == 2 + return SVector(equations.diffusivity * dudy) + end +end +end # @muladd diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index 65425213327..97dbe524fb1 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -65,6 +65,9 @@ function integrate(u_ode, semi::AbstractSemidiscretization; normalize = true) return integrate(cons2cons, u_ode, semi; normalize = normalize) end +# Select the right-hand side function corresponding to the semidiscretization `semi`. +@inline default_rhs(::AbstractSemidiscretization) = rhs! + """ calc_error_norms([func=(u_node,equations)->u_node,] u_ode, t, analyzer, semi::AbstractSemidiscretization, cache_analysis) @@ -123,18 +126,19 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan; end u0_ode = compute_coefficients(first(tspan), semi) # Invoke initial condition + rhs_semi! = default_rhs(semi) # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using # mpi_isparallel() && MPI.Barrier(mpi_comm()) # See https://github.com/trixi-framework/Trixi.jl/issues/328 - iip = true # is-inplace, i.e., we modify a vector when calling rhs! - specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi) + iip = true # is-inplace, i.e., we modify a vector when calling `rhs_semi!` + specialize = SciMLBase.FullSpecialize # specialize on `rhs_semi!` and parameters (semi) # Check if Jacobian prototype is provided for sparse Jacobian if jac_prototype !== nothing # Convert `jac_prototype` to real type, as seen here: # https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Declaring-a-Sparse-Jacobian-with-Automatic-Sparsity-Detection - ode = SciMLBase.ODEFunction(rhs!, + ode = SciMLBase.ODEFunction(rhs_semi!, jac_prototype = convert.(eltype(u0_ode), jac_prototype), colorvec = colorvec) # coloring vector is optional @@ -142,9 +146,9 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan; return ODEProblem{iip, specialize}(ode, u0_ode, tspan, semi) else # We could also construct an `ODEFunction` explicitly without the Jacobian here, - # but we stick to the lean direct in-place function `rhs!` and + # but we stick to the lean direct in-place function `rhs_semi!` and # let OrdinaryDiffEq.jl handle the rest - return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) + return ODEProblem{iip, specialize}(rhs_semi!, u0_ode, tspan, semi) end end @@ -177,18 +181,19 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan, end u0_ode = load_restart_file(semi, restart_file) # Load initial condition from restart file + rhs_semi! = default_rhs(semi) # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using # mpi_isparallel() && MPI.Barrier(mpi_comm()) # See https://github.com/trixi-framework/Trixi.jl/issues/328 - iip = true # is-inplace, i.e., we modify a vector when calling rhs! - specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi) + iip = true # is-inplace, i.e., we modify a vector when calling `rhs_semi!` + specialize = SciMLBase.FullSpecialize # specialize on `rhs_semi!` and parameters (semi) # Check if Jacobian prototype is provided for sparse Jacobian if jac_prototype !== nothing # Convert `jac_prototype` to real type, as seen here: # https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Declaring-a-Sparse-Jacobian-with-Automatic-Sparsity-Detection - ode = SciMLBase.ODEFunction(rhs!, + ode = SciMLBase.ODEFunction(rhs_semi!, jac_prototype = convert.(eltype(u0_ode), jac_prototype), colorvec = colorvec) # coloring vector is optional @@ -196,9 +201,9 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan, return ODEProblem{iip, specialize}(ode, u0_ode, tspan, semi) else # We could also construct an `ODEFunction` explicitly without the Jacobian here, - # but we stick to the lean direct in-place function `rhs!` and + # but we stick to the lean direct in-place function `rhs_semi!` and # let OrdinaryDiffEq.jl handle the rest - return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) + return ODEProblem{iip, specialize}(rhs_semi!, u0_ode, tspan, semi) end end @@ -236,6 +241,33 @@ function compute_coefficients!(u_ode, func, t, semi::AbstractSemidiscretization) mesh_equations_solver_cache(semi)...) end +# Helper function to compute linear structure for a given semidiscretization +# `semi` and applied RHS `apply_rhs!` +function _linear_structure_from_rhs(semi::AbstractSemidiscretization, apply_rhs!) + # allocate memory + u_ode = allocate_coefficients(mesh_equations_solver_cache(semi)...) + du_ode = similar(u_ode) + + # get the right hand side from boundary conditions and optional source terms + u_ode .= zero(eltype(u_ode)) + apply_rhs!(du_ode, u_ode) + b = -du_ode + + # Create a copy of `b` used internally to extract the linear part of `semi`. + # This is necessary to get everything correct when the user updates the + # returned vector `b`. + b_tmp = copy(b) + + # wrap the linear operator + A = LinearMap(length(u_ode), ismutating = true) do dest, src + apply_rhs!(dest, src) + @. dest += b_tmp + return dest + end + + return A, b +end + """ linear_structure(semi::AbstractSemidiscretization; t0 = zero(real(semi))) @@ -270,28 +302,11 @@ function linear_structure(semi::AbstractSemidiscretization; throw(ArgumentError("`linear_structure` expects linear equations.")) end - # allocate memory - u_ode = allocate_coefficients(mesh_equations_solver_cache(semi)...) - du_ode = similar(u_ode) - - # get the right hand side from boundary conditions and optional source terms - u_ode .= zero(eltype(u_ode)) - rhs!(du_ode, u_ode, semi, t0) - b = -du_ode - - # Create a copy of `b` used internally to extract the linear part of `semi`. - # This is necessary to get everything correct when the user updates the - # returned vector `b`. - b_tmp = copy(b) - - # wrap the linear operator - A = LinearMap(length(u_ode), ismutating = true) do dest, src - rhs!(dest, src, semi, t0) - @. dest += b_tmp - return dest + apply_rhs! = function (dest, src) + return rhs!(dest, src, semi, t0) end - return A, b + return _linear_structure_from_rhs(semi, apply_rhs!) end """ @@ -311,9 +326,10 @@ function jacobian_fd(semi::AbstractSemidiscretization; du0_ode = similar(u_ode) dup_ode = similar(u_ode) dum_ode = similar(u_ode) + rhs_semi! = default_rhs(semi) # compute residual of linearization state - rhs!(du0_ode, u_ode, semi, t0) + rhs_semi!(du0_ode, u_ode, semi, t0) # initialize Jacobian matrix J = zeros(eltype(u_ode), length(u_ode), length(u_ode)) @@ -331,11 +347,11 @@ function jacobian_fd(semi::AbstractSemidiscretization; # plus fluctuation u_ode[idx] = u0_ode[idx] + epsilon - rhs!(dup_ode, u_ode, semi, t0) + rhs_semi!(dup_ode, u_ode, semi, t0) # minus fluctuation u_ode[idx] = u0_ode[idx] - epsilon - rhs!(dum_ode, u_ode, semi, t0) + rhs_semi!(dum_ode, u_ode, semi, t0) # restore linearization state u_ode[idx] = u0_ode[idx] @@ -374,11 +390,12 @@ end function _jacobian_ad_forward(semi, t0, u0_ode, du_ode, config) new_semi = remake(semi, uEltype = eltype(config)) + rhs_semi! = default_rhs(new_semi) # Create anonymous function passed as first argument to `ForwardDiff.jacobian` to match # `ForwardDiff.jacobian(f!, y::AbstractArray, x::AbstractArray, # cfg::JacobianConfig = JacobianConfig(f!, y, x), check=Val{true}())` J = ForwardDiff.jacobian(du_ode, u0_ode, config) do du_ode, u_ode - return Trixi.rhs!(du_ode, u_ode, new_semi, t0) + return rhs_semi!(du_ode, u_ode, new_semi, t0) end return J @@ -409,6 +426,7 @@ end function _jacobian_ad_forward_structarrays(semi, t0, u0_ode_plain, du_ode_plain, config) new_semi = remake(semi, uEltype = eltype(config)) + rhs_semi! = default_rhs(new_semi) # Create anonymous function passed as first argument to `ForwardDiff.jacobian` to match # `ForwardDiff.jacobian(f!, y::AbstractArray, x::AbstractArray, # cfg::JacobianConfig = JacobianConfig(f!, y, x), check=Val{true}())` @@ -424,7 +442,7 @@ function _jacobian_ad_forward_structarrays(semi, t0, u0_ode_plain, du_ode_plain, :, v), nvariables(semi))) - return Trixi.rhs!(du_ode, u_ode, new_semi, t0) + return rhs_semi!(du_ode, u_ode, new_semi, t0) end return J @@ -447,11 +465,12 @@ end function _jacobian_ad_forward_staticarrays(semi, t0, u0_ode_plain, du_ode_plain, config) new_semi = remake(semi, uEltype = eltype(config)) + rhs_semi! = default_rhs(new_semi) J = ForwardDiff.jacobian(du_ode_plain, u0_ode_plain, config) do du_ode_plain, u_ode_plain u_ode = reinterpret(SVector{nvariables(semi), eltype(config)}, u_ode_plain) du_ode = reinterpret(SVector{nvariables(semi), eltype(config)}, du_ode_plain) - return Trixi.rhs!(du_ode, u_ode, new_semi, t0) + return rhs_semi!(du_ode, u_ode, new_semi, t0) end return J diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 69c5e8fbc8d..460c1a196f1 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -432,37 +432,17 @@ function linear_structure(semi::SemidiscretizationHyperbolicParabolic; throw(ArgumentError("`linear_structure` expects linear equations.")) end - # allocate memory - u_ode = allocate_coefficients(mesh_equations_solver_cache(semi)...) - du_ode = similar(u_ode) - - # get the right hand side from boundary conditions and optional source terms - u_ode .= zero(eltype(u_ode)) - rhs!(du_ode, u_ode, semi, t0) - b = -du_ode - - # Repeat for parabolic part - rhs_parabolic!(du_ode, u_ode, semi, t0) - @. b -= du_ode - - # Create a copy of `b` used internally to extract the linear part of `semi`. - # This is necessary to get everything correct when the user updates the - # returned vector `b`. - b_tmp = copy(b) - # additional storage for parabolic part - dest_para = similar(du_ode) + dest_para = allocate_coefficients(mesh_equations_solver_cache(semi)...) - # wrap the linear operator - A = LinearMap(length(u_ode), ismutating = true) do dest, src + apply_rhs! = function (dest, src) rhs!(dest, src, semi, t0) rhs_parabolic!(dest_para, src, semi, t0) - - @. dest += dest_para + b_tmp + @. dest += dest_para return dest end - return A, b + return _linear_structure_from_rhs(semi, apply_rhs!) end function _jacobian_ad_forward(semi::SemidiscretizationHyperbolicParabolic, t0, u0_ode, @@ -519,29 +499,11 @@ function linear_structure_parabolic(semi::SemidiscretizationHyperbolicParabolic; throw(ArgumentError("`linear_structure_parabolic` expects equations with constant diffusive terms.")) end - # allocate memory - u_ode = allocate_coefficients(mesh_equations_solver_cache(semi)...) - du_ode = similar(u_ode) - - # get the parabolic right hand side from boundary conditions and optional source terms - u_ode .= zero(eltype(u_ode)) - rhs_parabolic!(du_ode, u_ode, semi, t0) - b = -du_ode - - # Create a copy of `b` used internally to extract the linear part of `semi`. - # This is necessary to get everything correct when the user updates the - # returned vector `b`. - b_tmp = copy(b) - - # wrap the linear operator - A = LinearMap(length(u_ode), ismutating = true) do dest, src - rhs_parabolic!(dest, src, semi, t0) - - @. dest += b_tmp - return dest + apply_rhs! = function (dest, src) + return rhs_parabolic!(dest, src, semi, t0) end - return A, b + return _linear_structure_from_rhs(semi, apply_rhs!) end """ diff --git a/src/semidiscretization/semidiscretization_parabolic.jl b/src/semidiscretization/semidiscretization_parabolic.jl new file mode 100644 index 00000000000..4c663c70381 --- /dev/null +++ b/src/semidiscretization/semidiscretization_parabolic.jl @@ -0,0 +1,236 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +""" + SemidiscretizationParabolic + +A struct containing everything needed to describe a spatial semidiscretization of a purely +parabolic PDE. +""" +mutable struct SemidiscretizationParabolic{Mesh, Equations, InitialCondition, + BoundaryConditions, SourceTerms, + Solver, SolverParabolic, + Cache, CacheParabolic} <: + AbstractSemidiscretization + mesh::Mesh + equations::Equations + + const initial_condition::InitialCondition + + const boundary_conditions::BoundaryConditions + const source_terms::SourceTerms + + const solver::Solver + const solver_parabolic::SolverParabolic + + cache::Cache + cache_parabolic::CacheParabolic + + performance_counter::PerformanceCounter +end + +""" + SemidiscretizationParabolic(mesh, equations, initial_condition, solver; + solver_parabolic=default_parabolic_solver(), + source_terms=nothing, + boundary_conditions, + RealT=real(solver), + uEltype=RealT) + +Construct a semidiscretization of a purely parabolic PDE. + +Boundary conditions must be provided explicitly either as a `NamedTuple` or as a +single boundary condition that gets applied to all boundaries. +""" +function SemidiscretizationParabolic(mesh, equations::AbstractEquationsParabolic, + initial_condition, solver; + solver_parabolic = default_parabolic_solver(), + source_terms = nothing, + boundary_conditions, + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT = real(solver), uEltype = RealT) + @assert ndims(mesh) == ndims(equations) + + cache = create_cache(mesh, equations, solver, RealT, uEltype) + _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, + cache) + check_periodicity_mesh_boundary_conditions(mesh, _boundary_conditions) + + cache_parabolic = create_cache_parabolic(mesh, equations, solver, + nelements(solver, cache), uEltype) + + performance_counter = PerformanceCounter() + + return SemidiscretizationParabolic{typeof(mesh), typeof(equations), + typeof(initial_condition), + typeof(_boundary_conditions), + typeof(source_terms), + typeof(solver), + typeof(solver_parabolic), + typeof(cache), + typeof(cache_parabolic)}(mesh, equations, + initial_condition, + _boundary_conditions, + source_terms, + solver, + solver_parabolic, + cache, + cache_parabolic, + performance_counter) +end + +# @eval due to @muladd +@eval Adapt.@adapt_structure(SemidiscretizationParabolic) + +# Create a new semidiscretization but change some parameters compared to the input. +function remake(semi::SemidiscretizationParabolic; uEltype = real(semi.solver), + mesh = semi.mesh, + equations = semi.equations, + initial_condition = semi.initial_condition, + solver = semi.solver, + solver_parabolic = semi.solver_parabolic, + source_terms = semi.source_terms, + boundary_conditions = semi.boundary_conditions) + return SemidiscretizationParabolic(mesh, equations, initial_condition, solver; + solver_parabolic = solver_parabolic, + source_terms = source_terms, + boundary_conditions = boundary_conditions, + uEltype = uEltype) +end + +function Base.show(io::IO, semi::SemidiscretizationParabolic) + @nospecialize semi # reduce precompilation time + + print(io, "SemidiscretizationParabolic(") + print(io, semi.mesh) + print(io, ", ", semi.equations) + print(io, ", ", semi.initial_condition) + print(io, ", ", semi.boundary_conditions) + print(io, ", ", semi.source_terms) + print(io, ", ", semi.solver) + print(io, ", ", semi.solver_parabolic) + print(io, ", cache(") + for (idx, key) in enumerate(keys(semi.cache)) + idx > 1 && print(io, " ") + print(io, key) + end + print(io, "))") + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationParabolic) + @nospecialize semi # reduce precompilation time + + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationParabolic") + summary_line(io, "#spatial dimensions", ndims(semi.equations)) + summary_line(io, "mesh", semi.mesh) + summary_line(io, "equations", semi.equations |> typeof |> nameof) + summary_line(io, "initial condition", semi.initial_condition) + summary_line(io, "source terms", semi.source_terms) + summary_line(io, "solver", semi.solver |> typeof |> nameof) + summary_line(io, "parabolic solver", semi.solver_parabolic |> typeof |> nameof) + summary_line(io, "total #DOFs per field", ndofsglobal(semi)) + summary_footer(io) + end +end + +@inline Base.ndims(semi::SemidiscretizationParabolic) = ndims(semi.mesh) + +@inline function nvariables(semi::SemidiscretizationParabolic) + return nvariables(semi.equations) +end + +@inline Base.real(semi::SemidiscretizationParabolic) = real(semi.solver) + +@inline function mesh_equations_solver_cache(semi::SemidiscretizationParabolic) + @unpack mesh, equations, solver, cache = semi + return mesh, equations, solver, cache +end + +function calc_error_norms(func, u_ode, t, analyzer, + semi::SemidiscretizationParabolic, cache_analysis) + @unpack mesh, equations, initial_condition, solver, cache = semi + u = wrap_array(u_ode, mesh, equations, solver, cache) + + return calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, + solver, cache, cache_analysis) +end + +function compute_coefficients(t, semi::SemidiscretizationParabolic) + return compute_coefficients(semi.initial_condition, t, semi) +end + +function compute_coefficients!(u_ode, t, semi::SemidiscretizationParabolic) + return compute_coefficients!(u_ode, semi.initial_condition, t, semi) +end + +# Required for storing `extra_node_variables` in the `SaveSolutionCallback`. +# Not to be confused with `get_node_vars` which returns the variables of the simulated equation. +function get_node_variables!(node_variables, u_ode, + semi::SemidiscretizationParabolic) + return get_node_variables!(node_variables, u_ode, + mesh_equations_solver_cache(semi)..., + semi.cache_parabolic) +end + +""" + linear_structure(semi::SemidiscretizationParabolic; + t0 = zero(real(semi))) + +Wraps the right-hand side operator of the parabolic semidiscretization `semi` +at time `t0` as an affine-linear operator given by a linear operator `A` +and a vector `b`: +```math +\\partial_t u(t) = A u(t) - b. +``` +Works only for equations for which `have_constant_diffusivity(equations) == True()`. +As in [`linear_structure(semi::AbstractSemidiscretization)`](@ref), returns `A` and `b`, +with `A` represented as a matrix-free `LinearMap`. +""" +function linear_structure(semi::SemidiscretizationParabolic; + t0 = zero(real(semi))) + if have_constant_diffusivity(semi.equations) == False() + throw(ArgumentError("`linear_structure` expects equations with constant diffusive terms.")) + end + + apply_rhs! = function (dest, src) + return rhs_parabolic!(dest, src, semi, t0) + end + + return _linear_structure_from_rhs(semi, apply_rhs!) +end + +# For a purely parabolic semidiscretization, the right-hand side is `rhs_parabolic!` +# instead of the default `rhs!`. +@inline default_rhs(::SemidiscretizationParabolic) = rhs_parabolic! + +function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationParabolic, t) + @unpack mesh, equations, boundary_conditions, source_terms, solver, solver_parabolic, cache, cache_parabolic = semi + + u = wrap_array(u_ode, mesh, equations, solver, cache) + du = wrap_array(du_ode, mesh, equations, solver, cache) + backend = trixi_backend(u) + + time_start = time_ns() + @trixi_timeit_ext backend timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, + equations, + boundary_conditions, + source_terms, + solver, + solver_parabolic, + cache, + cache_parabolic) + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end +end # @muladd diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 911213d1c84..6d2c7207719 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -38,7 +38,22 @@ function get_node_variables!(node_variables, u_ode, mesh, equations, return nothing end -# Version for parabolic-extended equations + +# Version for purely parabolic equations (adds cache_parabolic). +function get_node_variables!(node_variables, u_ode, mesh, equations, + dg, cache, cache_parabolic) + if !isempty(node_variables) + u = wrap_array(u_ode, mesh, equations, dg, cache) + for var in keys(node_variables) + node_variables[var] = get_node_variable(Val(var), u, mesh, equations, + dg, cache, cache_parabolic) + end + end + + return nothing +end + +# Version for hyperbolic-parabolic equations (adds equations_parabolic and cache_parabolic). function get_node_variables!(node_variables, u_ode, mesh, equations, dg, cache, equations_parabolic, cache_parabolic) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 9fef07fed7a..f16cd664326 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -61,28 +61,6 @@ end @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) end -@trixi_testset "TreeMesh1D: elixir_diffusion_ldg.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", - "elixir_diffusion_ldg.jl"), - initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, - l2=[9.235894939144276e-6], linf=[5.402550135213957e-5]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) -end - -@trixi_testset "TreeMesh1D: elixir_diffusion_ldg_newton_krylov.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", - "elixir_diffusion_ldg_newton_krylov.jl"), - atol_lin_solve=1e-11, rtol_lin_solve=1e-10, - atol_ode_solve=1e-10, rtol_ode_solve=1e-9, - l2=[4.14999791227157e-6], linf=[2.424658410971059e-5]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) - @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) -end - @trixi_testset "TreeMesh1D: elixir_advection_diffusion_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", "elixir_advection_diffusion_restart.jl"), @@ -526,6 +504,96 @@ end end end +@testset "SemidiscretizationParabolic (1D)" begin +#! format: noindent + +@trixi_testset "TreeMesh1D: elixir_diffusion_ldg.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_diffusion_ldg.jl"), + initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, + l2=[9.235894939144276e-6], linf=[5.402550135213957e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) +end + +@trixi_testset "TreeMesh1D: elixir_diffusion_ldg_newton_krylov.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_diffusion_ldg_newton_krylov.jl"), + atol_lin_solve=1e-11, rtol_lin_solve=1e-10, + atol_ode_solve=1e-10, rtol_ode_solve=1e-9, + l2=[4.14999791227157e-6], linf=[2.424658410971059e-5], + # Relax error tols to avoid stochastic CI failures + atol=1e-12) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) +end + +@trixi_testset "TreeMesh1D: elixir_diffusion_ldg_amr_boundary_layer.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_diffusion_ldg_amr_boundary_layer.jl"), + l2=[0.5881457102264551], linf=[0.9302621795999283]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + + # Test `show` method not exercised in elixirs + @trixi_test_nowarn show(IOContext(stdout, :compact => true), MIME"text/plain"(), + semi) + + # Test basic semidiscretization dispatches + @test ndims(semi) == ndims(semi.mesh) + @test nvariables(semi) == nvariables(semi.equations) + @test real(semi) == real(semi.solver) + + # Test that `remake` works for `SemidiscretizationParabolic` + semi_remade = remake(semi) + @test semi_remade isa SemidiscretizationParabolic + @test semi_remade !== semi + @test semi_remade.mesh === semi.mesh + @test Trixi.ndofsglobal(semi_remade) == Trixi.ndofsglobal(semi) +end + +@trixi_testset "TreeMesh1D consistency check: elixir_diffusion_ldg_dirichlet.jl" begin + # Run the Dirichlet-Dirichlet elixir (uses `SemidiscretizationParabolic`) + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_diffusion_ldg_dirichlet.jl"), + tspan=(0.0, 0.1), + analysis_callback=AnalysisCallback(semi, + interval = 100, + extra_analysis_errors = (:l2_error_primitive, + :linf_error_primitive), + extra_analysis_integrals = (entropy,)), + l2=[2.3481439150004898e-6], linf=[2.4576876189230656e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) + + # Store reference solution for comparison + reference_solution = copy(sol.u[end]) + + # Run again using an advection-diffusion equation with advection velocity zero + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_diffusion_ldg_dirichlet.jl"), + tspan=(0.0, 0.1), + equations=LinearScalarAdvectionEquation1D(0.0), + semi=SemidiscretizationHyperbolicParabolic(mesh, + (equations, + LaplaceDiffusion1D(diffusivity(), + equations)), + initial_condition, + solver; + solver_parabolic = solver_parabolic, + boundary_conditions = (boundary_conditions, + boundary_conditions))) + # Check if the solutions for `SemidiscretizationParabolic` match those from + # `SemidiscretizationHyperbolicParabolic` using the same Float64 tolerance defaults as + # `@test_trixi_include` in TrixiTest.jl. + @test sol.u[end]≈reference_solution atol=500 * eps(Float64) rtol=sqrt(eps(Float64)) +end +end + # Clean up afterwards: delete Trixi output directory @test_nowarn isdir(outdir) && rm(outdir, recursive = true) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index b10876dc263..ac99445476e 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -393,18 +393,6 @@ end @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) end -@trixi_testset "TreeMesh2D: elixir_diffusion_steady_state_linear_map.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", - "elixir_diffusion_steady_state_linear_map.jl"), - l2=[2.9029827892716424e-5], linf=[0.0003022506331279151], - # Relax error tols to avoid stochastic CI failures - atol=1e-10) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) - @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) -end - @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), @@ -1367,6 +1355,27 @@ end end end +@testset "SemidiscretizationParabolic (2D)" begin +#! format: noindent + +@trixi_testset "TreeMesh2D: elixir_diffusion_steady_state_linear_map.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", + "elixir_diffusion_steady_state_linear_map.jl"), + tspan=(0.0, 1.0e-4), + analysis_callback=AnalysisCallback(semi, + interval = 1, + extra_analysis_errors = (:l2_error_primitive, + :linf_error_primitive), + extra_analysis_integrals = (entropy,)), + l2=[2.9029827892716424e-5], linf=[0.0003022506331279151], + # Relax error tols to avoid stochastic CI failures + atol=1e-10) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) +end +end + # Clean up afterwards: delete Trixi.jl output directory @test_nowarn isdir(outdir) && rm(outdir, recursive = true) diff --git a/test/test_type.jl b/test/test_type.jl index 06d9c4b1192..ef4d7453612 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1863,6 +1863,20 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Linear Diffusion Equation" begin + for RealT in (Float32, Float64) + u = SVector(one(RealT)) + + equations_1d = LinearDiffusionEquation1D(RealT(0.1)) + @test eltype(@inferred cons2prim(u, equations_1d)) == RealT + @test eltype(@inferred cons2entropy(u, equations_1d)) == RealT + + equations_2d = LinearDiffusionEquation2D(RealT(0.1)) + @test eltype(@inferred cons2prim(u, equations_2d)) == RealT + @test eltype(@inferred cons2entropy(u, equations_2d)) == RealT + end + end + @timed_testset "Laplace Diffusion 2D" begin for RealT in (Float32, Float64) equations = LinearScalarAdvectionEquation2D(RealT(1), RealT(1))