From ec9091429cebcdeda6b22f2951c6a03b156a3e32 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 08:57:49 +0100 Subject: [PATCH 01/57] Add "GradientEquations2D" and start adding stuff to cache --- src/equations/equations.jl | 5 +++ src/equations/gradient_equations_2d.jl | 46 ++++++++++++++++++++++++++ src/solvers/dg/dg_2d.jl | 14 ++++++++ 3 files changed, 65 insertions(+) create mode 100644 src/equations/gradient_equations_2d.jl diff --git a/src/equations/equations.jl b/src/equations/equations.jl index c5165697de8..254a0f00442 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -48,6 +48,7 @@ function calcflux(u, orientation, equations) end # set sensible default values that may be overwritten by specific equations have_nonconservative_terms(::AbstractEquations) = Val(false) have_constant_speed(::AbstractEquations) = Val(false) +have_parabolic_terms(::AbstractEquations) = Val(false) default_analysis_errors(::AbstractEquations) = (:l2_error, :linf_error) default_analysis_integrals(::AbstractEquations) = (entropy_timederivative,) @@ -109,3 +110,7 @@ include("hyperbolic_diffusion_3d.jl") # Lattice-Boltzmann equation (advection part only) abstract type AbstractLatticeBoltzmannEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("lattice_boltzmann_2d.jl") + +# Gradient equations +abstract type AbstractGradientEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +include("lattice_boltzmann_2d.jl") diff --git a/src/equations/gradient_equations_2d.jl b/src/equations/gradient_equations_2d.jl new file mode 100644 index 00000000000..e4bf5e30b34 --- /dev/null +++ b/src/equations/gradient_equations_2d.jl @@ -0,0 +1,46 @@ + +@doc raw""" + GradientEquations2D + +The gradient equations +```math +q^d - \partial_d u = 0 +``` +in direction `d` in two space dimensions as required for, e.g., the Bassi & Rebay 1 (BR1) or the +local discontinuous Galerkin (LDG) schemes. +""" +struct GradientEquations2D{RealT<:Real, NVARS, Orientation} <: AbstractGradientEquations{2, NVARS} +end + +GradientEquations2D(::Type{RealT}, nvars, orientation) where RealT = GradientEquations2D{RealT, nvars, orientation}() + + +get_name(::GradientEquations2D) = "GradientEquations2D" +varnames(::typeof(cons2cons), equations::GradientEquations2D) = SVector(ntuple(v -> "gradient_"*string(v), nvariables(equations))) +varnames(::typeof(cons2prim), equations::GradientEquations2D) = varnames(cons2cons, equations) + +orient(::GradientEquations2D{RealT, NVARS, Orientation}) where {RealT, NVARS, Orientation} = Orientation + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_constant(x, t, equations::GradientEquations2D) + +A constant initial condition to test free-stream preservation. +""" +function initial_condition_constant(x, t, equation::GradientEquations2D) + return @SVector SVector(ntuple(v -> 0.0, nvariables(equations))) +end + + +# Pre-defined source terms should be implemented as +# function source_terms_WHATEVER(u, x, t, equations::GradientEquations2D) + + +# Calculate 1D flux in for a single point +@inline function calcflux(u, orientation, equation::GradientEquations2D) + if orientation == orient(equation) + return u + else + return @SVector SVector(ntuple(v -> 0.0, nvariables(equations))) + end +end diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index 8fcc2e6de45..aa655d85ada 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -23,6 +23,7 @@ function create_cache(mesh::TreeMesh{2}, equations::AbstractEquations{2}, # Add specialized parts of the cache required to compute the volume integral etc. cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg)...) cache = (;cache..., create_cache(mesh, equations, dg.mortar)...) + cache = (;cache..., create_cache(mesh, equations, dg, have_parabolic_terms(equations), cache)...) return cache end @@ -30,6 +31,19 @@ end # The methods below are specialized on the volume integral type # and called from the basic `create_cache` method at the top. +function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val{false}, cache) + NamedTuple() +end + +function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val{true}, cache) + gradients_x_ode = allocate_coefficients(mesh, equations, dg, cache) + gradients_y_ode = allocate_coefficients(mesh, equations, dg, cache) + gradients_x = wrap_array(gradients_x_ode, mesh, equations, dg, cache) + gradients_y = wrap_array(gradients_y_ode, mesh, equations, dg, cache) + + return (; gradients_x_ode, gradients_y_ode, gradients_x, gradients_y) +end + function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG) create_cache(mesh, have_nonconservative_terms(equations), equations, volume_integral, dg) end From 0422a0ac3419d5676b3b9fe6b8345c3ee107bfb9 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 14:38:57 +0100 Subject: [PATCH 02/57] Update gradient equations --- src/equations/gradient_equations_2d.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/equations/gradient_equations_2d.jl b/src/equations/gradient_equations_2d.jl index e4bf5e30b34..e3273441274 100644 --- a/src/equations/gradient_equations_2d.jl +++ b/src/equations/gradient_equations_2d.jl @@ -9,26 +9,26 @@ q^d - \partial_d u = 0 in direction `d` in two space dimensions as required for, e.g., the Bassi & Rebay 1 (BR1) or the local discontinuous Galerkin (LDG) schemes. """ -struct GradientEquations2D{RealT<:Real, NVARS, Orientation} <: AbstractGradientEquations{2, NVARS} +struct GradientEquations2D{RealT<:Real, NVARS} <: AbstractGradientEquations{2, NVARS} + orientation::Int end -GradientEquations2D(::Type{RealT}, nvars, orientation) where RealT = GradientEquations2D{RealT, nvars, orientation}() +GradientEquations2D(::Type{RealT}, nvars, orientation) where RealT = GradientEquations2D{RealT, nvars}(orientation) +GradientEquations2D(nvars, orientation) = GradientEquations2D(Float64, nvars, orientation) get_name(::GradientEquations2D) = "GradientEquations2D" varnames(::typeof(cons2cons), equations::GradientEquations2D) = SVector(ntuple(v -> "gradient_"*string(v), nvariables(equations))) varnames(::typeof(cons2prim), equations::GradientEquations2D) = varnames(cons2cons, equations) -orient(::GradientEquations2D{RealT, NVARS, Orientation}) where {RealT, NVARS, Orientation} = Orientation - # Set initial conditions at physical location `x` for time `t` """ initial_condition_constant(x, t, equations::GradientEquations2D) A constant initial condition to test free-stream preservation. """ -function initial_condition_constant(x, t, equation::GradientEquations2D) - return @SVector SVector(ntuple(v -> 0.0, nvariables(equations))) +function initial_condition_constant(x, t, equations::GradientEquations2D) + return @SVector SVector(ntuple(v -> zero(real(equations)), nvariables(equations))) end @@ -37,10 +37,10 @@ end # Calculate 1D flux in for a single point -@inline function calcflux(u, orientation, equation::GradientEquations2D) - if orientation == orient(equation) +@inline function calcflux(u, orientation, equations::GradientEquations2D) + if orientation == equations.orientation return u else - return @SVector SVector(ntuple(v -> 0.0, nvariables(equations))) + return @SVector SVector(ntuple(v -> zero(real(equations)), nvariables(equations))) end end From 59d52c601f3fe6a090492be8ecd66f31b2edee5e Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 14:58:18 +0100 Subject: [PATCH 03/57] Bug fixes for GradientEquations2D --- src/equations/equations.jl | 2 +- src/equations/gradient_equations_2d.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 254a0f00442..4f41769fcd0 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -113,4 +113,4 @@ include("lattice_boltzmann_2d.jl") # Gradient equations abstract type AbstractGradientEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end -include("lattice_boltzmann_2d.jl") +include("gradient_equations_2d.jl") diff --git a/src/equations/gradient_equations_2d.jl b/src/equations/gradient_equations_2d.jl index e3273441274..4bc7c4da184 100644 --- a/src/equations/gradient_equations_2d.jl +++ b/src/equations/gradient_equations_2d.jl @@ -28,7 +28,7 @@ varnames(::typeof(cons2prim), equations::GradientEquations2D) = varnames(cons2co A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::GradientEquations2D) - return @SVector SVector(ntuple(v -> zero(real(equations)), nvariables(equations))) + return SVector(ntuple(v -> zero(eltype(x)), nvariables(equations))) end @@ -41,6 +41,6 @@ end if orientation == equations.orientation return u else - return @SVector SVector(ntuple(v -> zero(real(equations)), nvariables(equations))) + return SVector(ntuple(v -> zero(eltype(u)), nvariables(equations))) end end From 01ad2ba614948a29c095184a539ff12726e26171 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 15:34:49 +0100 Subject: [PATCH 04/57] Add heat equation in 2D --- src/equations/equations.jl | 4 ++ src/equations/heat_equation_2d.jl | 86 +++++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+) create mode 100644 src/equations/heat_equation_2d.jl diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 4f41769fcd0..3c4bd0ef15b 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -114,3 +114,7 @@ include("lattice_boltzmann_2d.jl") # Gradient equations abstract type AbstractGradientEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("gradient_equations_2d.jl") + +# Heat equation +abstract type AbstractHeatEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +include("heat_equation_2d.jl") diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl new file mode 100644 index 00000000000..a80b431aa57 --- /dev/null +++ b/src/equations/heat_equation_2d.jl @@ -0,0 +1,86 @@ + +@doc raw""" + HeatEquation2D + +The linear scalar advection equation +```math +\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u = 0 +``` +in two space dimensions with constant velocity `a`. +""" +struct HeatEquation2D{RealT<:Real} <: AbstractHeatEquation{2, 1} + nu::RealT +end + +get_name(::HeatEquation2D) = "HeatEquation2D" +varnames(::typeof(cons2cons), ::HeatEquation2D) = SVector("scalar") +varnames(::typeof(cons2prim), ::HeatEquation2D) = SVector("scalar") + +have_parabolic_terms(::HeatEquations2D) = Val(true) + + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_constant(x, t, equations::HeatEquation2D) + +A constant initial condition to test free-stream preservation. +""" +function initial_condition_constant(x, t, equation::HeatEquation2D) + return @SVector [2.0] +end + + +""" + initial_condition_convergence_test(x, t, equations::HeatEquation2D) + +A smooth initial condition used for convergence tests. +""" +function initial_condition_convergence_test(x, t, equation::HeatEquation2D) + c = 1.0 + A = 0.5 + L = 2 + f = 1/L + omega = 2 * pi * f + scalar = c + A * sin(omega * sum(x)) + return @SVector [scalar] +end + + +# Pre-defined source terms should be implemented as +# function source_terms_WHATEVER(u, x, t, equations::HeatEquation2D) + + +# Calculate 1D flux in for a single point +@inline function calcflux(u, orientation, equation::HeatEquation2D) + return zero(u) +end + + +function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::HeatEquation2D) + return zero(u) +end + + + +@inline have_constant_speed(::HeatEquation2D) = Val(true) + +@inline function max_abs_speeds(equation::HeatEquation2D) + return NaN +end + + +# Convert conservative variables to primitive +@inline cons2prim(u, equation::HeatEquation2D) = u + +# Convert conservative variables to entropy variables +@inline cons2entropy(u, equation::HeatEquation2D) = u + + +# Calculate entropy for a conservative state `cons` +@inline entropy(u::Real, ::HeatEquation2D) = 0.5 * u^2 +@inline entropy(u, equation::HeatEquation2D) = entropy(u[1], equation) + + +# Calculate total energy for a conservative state `cons` +@inline energy_total(u::Real, ::HeatEquation2D) = 0.5 * u^2 +@inline energy_total(u, equation::HeatEquation2D) = energy_total(u[1], equation) From 9a5396a1443ef980c760c63f0f3d0518fb944079 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 15:35:51 +0100 Subject: [PATCH 05/57] Fix sign in gradient equations to calculate gradient and not negative gradient --- src/equations/gradient_equations_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/gradient_equations_2d.jl b/src/equations/gradient_equations_2d.jl index 4bc7c4da184..c227985d336 100644 --- a/src/equations/gradient_equations_2d.jl +++ b/src/equations/gradient_equations_2d.jl @@ -39,7 +39,7 @@ end # Calculate 1D flux in for a single point @inline function calcflux(u, orientation, equations::GradientEquations2D) if orientation == equations.orientation - return u + return -u else return SVector(ntuple(v -> zero(eltype(u)), nvariables(equations))) end From d3bcf3c107ab48ab988ee512a272890664e2a49d Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 15:36:11 +0100 Subject: [PATCH 06/57] Hack in the calculation of the Laplacian in DG 2D --- src/solvers/dg/dg_2d.jl | 41 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index aa655d85ada..12407ec5479 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -36,12 +36,22 @@ function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val end function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val{true}, cache) - gradients_x_ode = allocate_coefficients(mesh, equations, dg, cache) - gradients_y_ode = allocate_coefficients(mesh, equations, dg, cache) - gradients_x = wrap_array(gradients_x_ode, mesh, equations, dg, cache) - gradients_y = wrap_array(gradients_y_ode, mesh, equations, dg, cache) + solver_grad_n = DGSEM(polydeg(dg), flux_central) - return (; gradients_x_ode, gradients_y_ode, gradients_x, gradients_y) + equations_grad_x = GradientEquations2D(nvariables(equations), 1) + semi_grad_x = SemidiscretizationHyperbolic(mesh, equations_grad_x, + initial_condition_constant, solver_grad_x) + + equations_grad_y = GradientEquations2D(nvariables(equations), 2) + semi_grad_y = SemidiscretizationHyperbolic(mesh, equations_grad_y, + initial_condition_constant, solver_grad_y) + + u_ode_grad_x = compute_coefficients(nothing, semi_grad_x) + u_ode_grad_y = compute_coefficients(nothing, semi_grad_y) + u_ode_grad_xx = compute_coefficients(nothing, semi_grad_x) + u_ode_grad_yy = compute_coefficients(nothing, semi_grad_y) + + return (; semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy) end function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG) @@ -178,6 +188,27 @@ function rhs!(du::AbstractArray{<:Any,4}, u, t, # Apply Jacobian from mapping to reference element @timeit_debug timer() "Jacobian" apply_jacobian!(du, equations, dg, cache) + # Add Laplacian + if have_parabolic_terms(equations) == Val(true) + # Compute gradients + @unpack semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy = cache + rhs!(u_ode_grad_x, vec(u), semi_grad_x, t) + rhs!(u_ode_grad_y, vec(u), semi_grad_y, t) + + rhs!(u_ode_grad_xx, u_ode_grad_x, semi_grad_x, t) + rhs!(u_ode_grad_yy, u_ode_grad_y, semi_grad_y, t) + + grad_xx = wrap_array(u_ode_grad_xx, semi_grad_x) + grad_yy = wrap_array(u_ode_grad_yy, semi_grad_y) + + laplacian = grad_xx + laplacian .+= grad_yy + + # Add Laplacian to rhs + nu = 1.2 + du .+= nu .* laplacian + end + # Calculate source terms @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache) From cf6cc6a4438ef6df788d6c254dc3ba54fc3a4305 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 16:15:20 +0100 Subject: [PATCH 07/57] Add elixir for heat equation --- examples/2d/elixir_heat_basic.jl | 62 ++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 examples/2d/elixir_heat_basic.jl diff --git a/examples/2d/elixir_heat_basic.jl b/examples/2d/elixir_heat_basic.jl new file mode 100644 index 00000000000..af5c0096960 --- /dev/null +++ b/examples/2d/elixir_heat_basic.jl @@ -0,0 +1,62 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +nu = 1.2e-2 +equations = HeatEquation2D(nu) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(3, flux_lax_friedrichs) + +coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1, 1) # maximum coordinates (max(x), max(y)) + +# Create a uniformely refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=2, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +initial_condition = initial_condition_convergence_test +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + + +############################################################################### +# 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 SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step +# stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +# Print the timer summary +summary_callback() From d21e1c445a7f51304ef17d9b7bac641e332791bf Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 11 Dec 2020 16:15:55 +0100 Subject: [PATCH 08/57] Bug fixes --- src/Trixi.jl | 3 ++- src/equations/heat_equation_2d.jl | 7 ++++--- src/solvers/dg/dg_2d.jl | 6 +++--- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index bcb0da45b4a..a2f2bbe350b 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -67,7 +67,8 @@ export CompressibleEulerEquations1D, CompressibleEulerEquations2D, CompressibleE IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D, LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D, - LatticeBoltzmannEquations2D + LatticeBoltzmannEquations2D, + HeatEquation2D export flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_upwind, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_kennedy_gruber, flux_shima_etal diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index a80b431aa57..a21717c15c2 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -16,7 +16,7 @@ get_name(::HeatEquation2D) = "HeatEquation2D" varnames(::typeof(cons2cons), ::HeatEquation2D) = SVector("scalar") varnames(::typeof(cons2prim), ::HeatEquation2D) = SVector("scalar") -have_parabolic_terms(::HeatEquations2D) = Val(true) +have_parabolic_terms(::HeatEquation2D) = Val(true) # Set initial conditions at physical location `x` for time `t` @@ -36,12 +36,13 @@ end A smooth initial condition used for convergence tests. """ function initial_condition_convergence_test(x, t, equation::HeatEquation2D) + @unpack nu = equation c = 1.0 A = 0.5 L = 2 f = 1/L omega = 2 * pi * f - scalar = c + A * sin(omega * sum(x)) + scalar = c + A * sin(omega * sum(x)) * exp(-2 * nu * omega^2 * t) return @SVector [scalar] end @@ -57,7 +58,7 @@ end function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::HeatEquation2D) - return zero(u) + return zero(u_ll) end diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index 12407ec5479..008ee3b92a2 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -40,11 +40,11 @@ function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val equations_grad_x = GradientEquations2D(nvariables(equations), 1) semi_grad_x = SemidiscretizationHyperbolic(mesh, equations_grad_x, - initial_condition_constant, solver_grad_x) + initial_condition_constant, solver_grad_n) equations_grad_y = GradientEquations2D(nvariables(equations), 2) semi_grad_y = SemidiscretizationHyperbolic(mesh, equations_grad_y, - initial_condition_constant, solver_grad_y) + initial_condition_constant, solver_grad_n) u_ode_grad_x = compute_coefficients(nothing, semi_grad_x) u_ode_grad_y = compute_coefficients(nothing, semi_grad_y) @@ -205,7 +205,7 @@ function rhs!(du::AbstractArray{<:Any,4}, u, t, laplacian .+= grad_yy # Add Laplacian to rhs - nu = 1.2 + @unpack nu = equations du .+= nu .* laplacian end From 9eefbab1926179927c60fbdc81236bc926235772 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Sun, 13 Dec 2020 07:30:25 +0100 Subject: [PATCH 09/57] Add advection-diffusion equation --- examples/2d/elixir_advdiff_basic.jl | 67 +++++++++++ src/Trixi.jl | 3 +- src/equations/equations.jl | 4 + .../linear_advection_diffusion_2d.jl | 110 ++++++++++++++++++ 4 files changed, 183 insertions(+), 1 deletion(-) create mode 100644 examples/2d/elixir_advdiff_basic.jl create mode 100644 src/equations/linear_advection_diffusion_2d.jl diff --git a/examples/2d/elixir_advdiff_basic.jl b/examples/2d/elixir_advdiff_basic.jl new file mode 100644 index 00000000000..5c8bf9e8d23 --- /dev/null +++ b/examples/2d/elixir_advdiff_basic.jl @@ -0,0 +1,67 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advectionvelocity = (0.2, 0.6) +nu = 1.2e-2 +equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(4, flux_lax_friedrichs) + +coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1, 1) # maximum coordinates (max(x), max(y)) +refinement_patches = ( + (type="box", coordinates_min=(0.0, -1.0), coordinates_max=(1.0, 1.0)), +) + +# Create a uniformely refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=2, + refinement_patches=refinement_patches, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +initial_condition = initial_condition_convergence_test +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + + +############################################################################### +# 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 SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step +# stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index a2f2bbe350b..71b2048f6dd 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -68,7 +68,8 @@ export CompressibleEulerEquations1D, CompressibleEulerEquations2D, CompressibleE HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D, LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D, LatticeBoltzmannEquations2D, - HeatEquation2D + HeatEquation2D, + LinearAdvectionDiffusionEquation2D export flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_upwind, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_kennedy_gruber, flux_shima_etal diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 3c4bd0ef15b..df2a80b9550 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -118,3 +118,7 @@ include("gradient_equations_2d.jl") # Heat equation abstract type AbstractHeatEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("heat_equation_2d.jl") + +# Linear scalar advection-diffusion equation +abstract type AbstractLinearAdvectionDiffusionEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +include("linear_advection_diffusion_2d.jl") diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl new file mode 100644 index 00000000000..53c6d35279d --- /dev/null +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -0,0 +1,110 @@ + +@doc raw""" + LinearAdvectionDiffusionEquation2D + +The linear scalar advection equation +```math +\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u = 0 +``` +in two space dimensions with constant velocity `a`. +""" +struct LinearAdvectionDiffusionEquation2D{RealT<:Real} <: AbstractLinearAdvectionDiffusionEquation{2, 1} + advectionvelocity::SVector{2, RealT} + nu::RealT +end + +function LinearAdvectionDiffusionEquation2D(a::NTuple{2,<:Real}, nu::Real) + LinearAdvectionDiffusionEquation2D(SVector(a), nu) +end + +function LinearAdvectionDiffusionEquation2D(a1::Real, a2::Real, nu::Real) + a1, a2, nu = promote(a1, a2, nu) + LinearAdvectionDiffusionEquation2D(SVector(a1, a2), nu) +end + +get_name(::LinearAdvectionDiffusionEquation2D) = "LinearAdvectionDiffusionEquation2D" +varnames(::typeof(cons2cons), ::LinearAdvectionDiffusionEquation2D) = SVector("scalar") +varnames(::typeof(cons2prim), ::LinearAdvectionDiffusionEquation2D) = SVector("scalar") + +have_parabolic_terms(::LinearAdvectionDiffusionEquation2D) = Val(true) + +# Calculates translated coordinates `x` for a periodic domain +function x_trans_periodic_2d(x, domain_length = SVector(2, 2), center = SVector(0, 0)) + x_normalized = x .- center + x_shifted = x_normalized .% domain_length + x_offset = ((x_shifted .< -0.5*domain_length) - (x_shifted .> 0.5*domain_length)) .* domain_length + return center + x_shifted + x_offset +end + + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_constant(x, t, equations::LinearAdvectionDiffusionEquation2D) + +A constant initial condition to test free-stream preservation. +""" +function initial_condition_constant(x, t, equation::LinearAdvectionDiffusionEquation2D) + return @SVector [2.0] +end + + +""" + initial_condition_convergence_test(x, t, equations::LinearAdvectionDiffusionEquation2D) + +A smooth initial condition used for convergence tests. +""" +function initial_condition_convergence_test(x, t, equation::LinearAdvectionDiffusionEquation2D) + # Store translated coordinate for easy use of exact solution + x_trans = x - equation.advectionvelocity * t + + @unpack nu = equation + c = 1.0 + A = 0.5 + L = 2 + f = 1/L + omega = 2 * pi * f + scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) + return @SVector [scalar] +end + + +# Pre-defined source terms should be implemented as +# function source_terms_WHATEVER(u, x, t, equations::LinearAdvectionDiffusionEquation2D) + + +# Calculate 1D flux in for a single point +@inline function calcflux(u, orientation, equation::LinearAdvectionDiffusionEquation2D) + a = equation.advectionvelocity[orientation] + return a * u +end + + +function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::LinearAdvectionDiffusionEquation2D) + a = equation.advectionvelocity[orientation] + return 0.5 * ( a * (u_ll + u_rr) - abs(a) * (u_rr - u_ll) ) +end + + + +@inline have_constant_speed(::LinearAdvectionDiffusionEquation2D) = Val(true) + +@inline function max_abs_speeds(equation::LinearAdvectionDiffusionEquation2D) + return abs.(equation.advectionvelocity) +end + + +# Convert conservative variables to primitive +@inline cons2prim(u, equation::LinearAdvectionDiffusionEquation2D) = u + +# Convert conservative variables to entropy variables +@inline cons2entropy(u, equation::LinearAdvectionDiffusionEquation2D) = u + + +# Calculate entropy for a conservative state `cons` +@inline entropy(u::Real, ::LinearAdvectionDiffusionEquation2D) = 0.5 * u^2 +@inline entropy(u, equation::LinearAdvectionDiffusionEquation2D) = entropy(u[1], equation) + + +# Calculate total energy for a conservative state `cons` +@inline energy_total(u::Real, ::LinearAdvectionDiffusionEquation2D) = 0.5 * u^2 +@inline energy_total(u, equation::LinearAdvectionDiffusionEquation2D) = energy_total(u[1], equation) From de8618b19bb43f5295e36321c8451c0bac792194 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Sun, 17 Jan 2021 07:53:54 +0100 Subject: [PATCH 10/57] Add test for linear advection-diffusion --- test/test_examples_2d.jl | 3 +++ test/test_examples_2d_advdiff.jl | 19 +++++++++++++++++++ 2 files changed, 22 insertions(+) create mode 100644 test/test_examples_2d_advdiff.jl diff --git a/test/test_examples_2d.jl b/test/test_examples_2d.jl index a90977a568d..494e1ebea87 100644 --- a/test/test_examples_2d.jl +++ b/test/test_examples_2d.jl @@ -19,6 +19,9 @@ isdir(outdir) && rm(outdir, recursive=true) # Linear advection include("test_examples_2d_advection.jl") + # Linear advection-diffusion + include("test_examples_2d_advdiff.jl") + # Hyperbolic diffusion include("test_examples_2d_hypdiff.jl") diff --git a/test/test_examples_2d_advdiff.jl b/test/test_examples_2d_advdiff.jl new file mode 100644 index 00000000000..df6bd8de548 --- /dev/null +++ b/test/test_examples_2d_advdiff.jl @@ -0,0 +1,19 @@ +module TestExamples2DAdvDiff + +using Test +using Trixi + +include("test_trixi.jl") + +# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory +EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") + +@testset "Scalar advection-diffusion" begin + @testset "elixir_advdiff_basic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_basic.jl"), + l2 = [6.707801888154751e-5], + linf = [0.0005912932758251888]) + end +end + +end # module From 4265f1d38c15b57b453007e398721a1b30052ec7 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Sun, 17 Jan 2021 07:54:12 +0100 Subject: [PATCH 11/57] Remove duplicate function --- src/equations/linear_advection_diffusion_2d.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index 53c6d35279d..77eef096adb 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -28,14 +28,6 @@ varnames(::typeof(cons2prim), ::LinearAdvectionDiffusionEquation2D) = SVector("s have_parabolic_terms(::LinearAdvectionDiffusionEquation2D) = Val(true) -# Calculates translated coordinates `x` for a periodic domain -function x_trans_periodic_2d(x, domain_length = SVector(2, 2), center = SVector(0, 0)) - x_normalized = x .- center - x_shifted = x_normalized .% domain_length - x_offset = ((x_shifted .< -0.5*domain_length) - (x_shifted .> 0.5*domain_length)) .* domain_length - return center + x_shifted + x_offset -end - # Set initial conditions at physical location `x` for time `t` """ From 926b17300a52d0552e518fb7378b5c461d18d457 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 07:23:41 +0100 Subject: [PATCH 12/57] Add test for heat equation --- test/test_examples_2d.jl | 3 +++ test/test_examples_2d_heat.jl | 19 +++++++++++++++++++ 2 files changed, 22 insertions(+) create mode 100644 test/test_examples_2d_heat.jl diff --git a/test/test_examples_2d.jl b/test/test_examples_2d.jl index 494e1ebea87..b11ebb6aa46 100644 --- a/test/test_examples_2d.jl +++ b/test/test_examples_2d.jl @@ -22,6 +22,9 @@ isdir(outdir) && rm(outdir, recursive=true) # Linear advection-diffusion include("test_examples_2d_advdiff.jl") + # Heat equation + include("test_examples_2d_heat.jl") + # Hyperbolic diffusion include("test_examples_2d_hypdiff.jl") diff --git a/test/test_examples_2d_heat.jl b/test/test_examples_2d_heat.jl new file mode 100644 index 00000000000..76de3ed8342 --- /dev/null +++ b/test/test_examples_2d_heat.jl @@ -0,0 +1,19 @@ +module TestExamples2DHeat + +using Test +using Trixi + +include("test_trixi.jl") + +# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory +EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") + +@testset "Heat equation" begin + @testset "elixir_heat_basic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_basic.jl"), + l2 = [0.0008071867035363602], + linf = [0.004961971226386641]) + end +end + +end # module From 7eaeb20ce7e1814195aaeb6590a46489932b8ec9 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 07:54:32 +0100 Subject: [PATCH 13/57] Add semi parabolic aux vars --- src/Trixi.jl | 4 +- .../semidiscretization_parabolic_auxvars.jl | 184 ++++++++++++++++++ 2 files changed, 187 insertions(+), 1 deletion(-) create mode 100644 src/semidiscretization/semidiscretization_parabolic_auxvars.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index 71b2048f6dd..6301ae665dc 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -53,6 +53,7 @@ include("mesh/mesh.jl") include("solvers/dg/dg.jl") include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") +include("semidiscretization/semidiscretization_parabolic_auxvars.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") @@ -118,7 +119,8 @@ export DG, export nelements, nnodes, nvariables, eachelement, eachnode, eachvariable -export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate +export SemidiscretizationHyperbolic, SemidiscretizationParabolicAuxVars, + semidiscretize, compute_coefficients, integrate export SemidiscretizationEulerGravity, ParametersEulerGravity, timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N! diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl new file mode 100644 index 00000000000..c2b6630deb8 --- /dev/null +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -0,0 +1,184 @@ + +# TODO: Taal refactor, Mesh<:AbstractMesh{NDIMS}, Equations<:AbstractEquations{NDIMS} ? +""" + SemidiscretizationParabolicAuxVars + +A struct containing everything needed to describe a spatial semidiscretization +of a parabolic conservation law. +""" +struct SemidiscretizationParabolicAuxVars{Mesh, Equations, InitialCondition, BoundaryConditions, + SourceTerms, Solver, Cache} <: AbstractSemidiscretization + + mesh::Mesh + equations::Equations + + # This guy is a bit messy since we abuse it as some kind of "exact solution" + # although this doesn't really exist... + initial_condition::InitialCondition + + boundary_conditions::BoundaryConditions + source_terms::SourceTerms + solver::Solver + cache::Cache + performance_counter::PerformanceCounter + + function SemidiscretizationParabolicAuxVars{Mesh, Equations, InitialCondition, BoundaryConditions, + SourceTerms, Solver, Cache}( + mesh::Mesh, equations::Equations, + initial_condition::InitialCondition, boundary_conditions::BoundaryConditions, + source_terms::SourceTerms, + solver::Solver, cache::Cache) where {Mesh, Equations, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache} + @assert ndims(mesh) == ndims(equations) + @assert have_parabolic_terms(equations) == Val(true) "Cannot create parabolic semidiscretization for purely hyperbolic system" + + performance_counter = PerformanceCounter() + + new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache, performance_counter) + end +end + +""" + SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic) + +Construct a semidiscretization of a parabolic PDE. +""" +function SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + RealT=real(solver)) + + cache = create_cache(mesh, equations, solver, RealT) + _boundary_conditions = digest_boundary_conditions(boundary_conditions) + + solver_auxvars = create_solver_auxvars(solver) + + equations_grad_x = GradientEquations2D(nvariables(equations), 1) + semi_grad_x = SemidiscretizationHyperbolic(mesh, equations_grad_x, + initial_condition_constant, solver_auxvars) + + equations_grad_y = GradientEquations2D(nvariables(equations), 2) + semi_grad_y = SemidiscretizationHyperbolic(mesh, equations_grad_y, + initial_condition_constant, solver_auxvars) + + u_ode_grad_x = compute_coefficients(nothing, semi_grad_x) + u_ode_grad_y = compute_coefficients(nothing, semi_grad_y) + u_ode_grad_xx = compute_coefficients(nothing, semi_grad_x) + u_ode_grad_yy = compute_coefficients(nothing, semi_grad_y) + + cache = (; cache..., semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy) + + SemidiscretizationParabolicAuxVars{typeof(mesh), typeof(equations), typeof(initial_condition), typeof(_boundary_conditions), typeof(source_terms), typeof(solver), typeof(cache)}( + mesh, equations, initial_condition, _boundary_conditions, source_terms, solver, cache) +end + + +function Base.show(io::IO, semi::SemidiscretizationParabolicAuxVars) + print(io, "SemidiscretizationParabolicAuxVars(") + 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, ", cache(") + for (idx,key) in enumerate(keys(semi.cache)) + idx > 1 && print(io, " ") + print(io, key) + end + print(io, "))") +end + +function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationParabolicAuxVars) + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationParabolicAuxVars") + summary_line(io, "#spatial dimensions", ndims(semi.equations)) + summary_line(io, "mesh", semi.mesh) + summary_line(io, "equations", typeof(semi.equations).name) + summary_line(io, "initial condition", semi.initial_condition) + summary_line(io, "boundary conditions", 2*ndims(semi)) + if (semi.boundary_conditions isa Tuple || + semi.boundary_conditions isa NamedTuple || + semi.boundary_conditions isa AbstractArray) + bcs = semi.boundary_conditions + else + bcs = collect(semi.boundary_conditions for _ in 1:(2*ndims(semi))) + end + summary_line(increment_indent(io), "negative x", bcs[1]) + summary_line(increment_indent(io), "positive x", bcs[2]) + if ndims(semi) > 1 + summary_line(increment_indent(io), "negative y", bcs[3]) + summary_line(increment_indent(io), "positive y", bcs[4]) + end + if ndims(semi) > 2 + summary_line(increment_indent(io), "negative z", bcs[5]) + summary_line(increment_indent(io), "positive z", bcs[6]) + end + summary_line(io, "source terms", semi.source_terms) + summary_line(io, "solver", typeof(semi.solver).name) + summary_line(io, "total #DOFs", ndofs(semi)) + summary_footer(io) + end +end + + +@inline Base.ndims(semi::SemidiscretizationParabolicAuxVars) = ndims(semi.mesh) + +@inline nvariables(semi::SemidiscretizationParabolicAuxVars) = nvariables(semi.equations) + +@inline Base.real(semi::SemidiscretizationParabolicAuxVars) = real(semi.solver) + + +@inline function mesh_equations_solver_cache(semi::SemidiscretizationParabolicAuxVars) + @unpack mesh, equations, solver, cache = semi + return mesh, equations, solver, cache +end + + +function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationParabolicAuxVars) + @unpack mesh, equations, initial_condition, solver, cache = semi + u = wrap_array(u_ode, mesh, equations, solver, cache) + + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) +end + +function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationParabolicAuxVars) + @unpack mesh, equations, initial_condition, solver, cache = semi + + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) +end + + +function compute_coefficients(t, semi::SemidiscretizationParabolicAuxVars) + compute_coefficients(semi.initial_condition, t, semi) +end + +function compute_coefficients!(u_ode::AbstractVector, t, semi::SemidiscretizationParabolicAuxVars) + compute_coefficients!(u_ode, semi.initial_condition, t, semi) +end + + +function rhs!(du_ode, u_ode, semi::SemidiscretizationParabolicAuxVars, t) + @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi + + u = wrap_array(u_ode, mesh, equations, solver, cache) + du = wrap_array(du_ode, mesh, equations, solver, cache) + + # Compute gradients + @unpack semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y = cache + rhs!(u_ode_grad_x, u_ode, semi_grad_x, t) + rhs!(u_ode_grad_y, u_ode, semi_grad_y, t) + gradients = (wrap_array(u_ode_grad_x, mesh, equations, solver, cache), + wrap_array(u_ode_grad_y, mesh, equations, solver, cache)) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @timeit_debug timer() "rhs!" rhs!(du, u, gradients, t, mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache) + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end From 8d1941f05acb1fd978f038fae0429ad97de71aa6 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 07:56:11 +0100 Subject: [PATCH 14/57] Add parabolic flux function to heat equation --- src/equations/heat_equation_2d.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index a21717c15c2..7aaad8fa06c 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -51,12 +51,18 @@ end # function source_terms_WHATEVER(u, x, t, equations::HeatEquation2D) -# Calculate 1D flux in for a single point +# Calculate hyperbolic 1D flux in axis `orientation` for a single point @inline function calcflux(u, orientation, equation::HeatEquation2D) return zero(u) end +# Calculate parabolic 1D flux in axis `orientation` for a single point +@inline function calcflux(u, gradients, orientation, equation::HeatEquation2D) + return gradients[orientation] +end + + function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::HeatEquation2D) return zero(u_ll) end From 5d768397efdfdf8e0047b2492d6768a44c97eacf Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 07:57:16 +0100 Subject: [PATCH 15/57] Allow creation of solvers for auxiliary variables --- src/solvers/dg/dg.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index 79a77280d8c..39dadea6670 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -274,6 +274,13 @@ function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, end +function create_solver_auxvars(dg::DGSEM) + @unpack basis = dg + + return DGSEM(basis, flux_central, VolumeIntegralWeakForm(), MortarL2(basis)) +end + + # indicators used for shock-capturing and AMR include("indicators.jl") From 4b4b2cde87325c8273b8a15eb95a3fe6aa4502ac Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 07:57:44 +0100 Subject: [PATCH 16/57] Switch heat equation to parabolic semi --- examples/2d/elixir_heat_basic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/2d/elixir_heat_basic.jl b/examples/2d/elixir_heat_basic.jl index af5c0096960..1baf319e46e 100644 --- a/examples/2d/elixir_heat_basic.jl +++ b/examples/2d/elixir_heat_basic.jl @@ -21,7 +21,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # A semidiscretization collects data structures and functions for the spatial discretization initial_condition = initial_condition_convergence_test -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver) ############################################################################### From a79f8978d4a77000ce55b007b0b0c4e90bb68c9f Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 07:58:19 +0100 Subject: [PATCH 17/57] Heat equation test still passes --- src/solvers/dg/dg_2d.jl | 113 +++++++++++++++++++++++++++++++++------- 1 file changed, 94 insertions(+), 19 deletions(-) diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index 008ee3b92a2..8f57407a663 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -36,22 +36,24 @@ function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val end function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val{true}, cache) - solver_grad_n = DGSEM(polydeg(dg), flux_central) - - equations_grad_x = GradientEquations2D(nvariables(equations), 1) - semi_grad_x = SemidiscretizationHyperbolic(mesh, equations_grad_x, - initial_condition_constant, solver_grad_n) - - equations_grad_y = GradientEquations2D(nvariables(equations), 2) - semi_grad_y = SemidiscretizationHyperbolic(mesh, equations_grad_y, - initial_condition_constant, solver_grad_n) - - u_ode_grad_x = compute_coefficients(nothing, semi_grad_x) - u_ode_grad_y = compute_coefficients(nothing, semi_grad_y) - u_ode_grad_xx = compute_coefficients(nothing, semi_grad_x) - u_ode_grad_yy = compute_coefficients(nothing, semi_grad_y) - - return (; semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy) + return NamedTuple() + +# solver_grad_n = DGSEM(polydeg(dg), flux_central) +# +# equations_grad_x = GradientEquations2D(nvariables(equations), 1) +# semi_grad_x = SemidiscretizationHyperbolic(mesh, equations_grad_x, +# initial_condition_constant, solver_grad_n) +# +# equations_grad_y = GradientEquations2D(nvariables(equations), 2) +# semi_grad_y = SemidiscretizationHyperbolic(mesh, equations_grad_y, +# initial_condition_constant, solver_grad_n) +# +# u_ode_grad_x = compute_coefficients(nothing, semi_grad_x) +# u_ode_grad_y = compute_coefficients(nothing, semi_grad_y) +# u_ode_grad_xx = compute_coefficients(nothing, semi_grad_x) +# u_ode_grad_yy = compute_coefficients(nothing, semi_grad_y) +# +# return (; semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy) end function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG) @@ -188,13 +190,55 @@ function rhs!(du::AbstractArray{<:Any,4}, u, t, # Apply Jacobian from mapping to reference element @timeit_debug timer() "Jacobian" apply_jacobian!(du, equations, dg, cache) + # Calculate source terms + @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache) + + return nothing +end + +function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, + mesh::TreeMesh{2}, equations, + initial_condition, boundary_conditions, source_terms, + dg::DG, cache) + # Reset du + @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du)) + + # Calculate volume integral + @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) + + # Prolong solution to interfaces + @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(cache, u, equations, dg) + + # Calculate interface fluxes + @timeit_debug timer() "interface flux" calc_interface_flux!(cache.elements.surface_flux_values, + have_nonconservative_terms(equations), equations, + dg, cache) + + # Prolong solution to boundaries + @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, u, equations, dg) + + # Calculate boundary fluxes + @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, equations, dg) + + # Prolong solution to mortars + @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, u, equations, dg.mortar, dg) + + # Calculate mortar fluxes + @timeit_debug timer() "mortar flux" calc_mortar_flux!(cache.elements.surface_flux_values, + have_nonconservative_terms(equations), equations, + dg.mortar, dg, cache) + + # Calculate surface integrals + @timeit_debug timer() "surface integral" calc_surface_integral!(du, equations, dg, cache) + + # Apply Jacobian from mapping to reference element + @timeit_debug timer() "Jacobian" apply_jacobian!(du, equations, dg, cache) + # Add Laplacian if have_parabolic_terms(equations) == Val(true) # Compute gradients @unpack semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy = cache - rhs!(u_ode_grad_x, vec(u), semi_grad_x, t) - rhs!(u_ode_grad_y, vec(u), semi_grad_y, t) - rhs!(u_ode_grad_xx, u_ode_grad_x, semi_grad_x, t) rhs!(u_ode_grad_yy, u_ode_grad_y, semi_grad_y, t) @@ -244,6 +288,37 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, end +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, gradients, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralWeakForm, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + + Threads.@threads for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + grad_x_node = get_node_vars(gradients[1], equations, dg, i, j, element) + grad_y_node = get_node_vars(gradients[2], equations, dg, i, j, element) + grad_node = (grad_x_node, grad_y_node) + + flux1 = calcflux(u_node, grad_node, 1, equations) + for ii in eachnode(dg) + integral_contribution = derivative_dhat[ii, i] * flux1 + add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element) + end + + flux2 = calcflux(u_node, grad_node, 2, equations) + for jj in eachnode(dg) + integral_contribution = derivative_dhat[jj, j] * flux2 + add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element) + end + end + end + + return nothing +end + + # Calculate 2D twopoint flux (element version) @inline function calcflux_twopoint!(f1, f2, u::AbstractArray{<:Any,4}, element, volume_flux, equations, dg::DG, cache) From 7dfa562fb40e9ac0c8997e9bc5bc8ef5d46f9461 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 11:45:32 +0100 Subject: [PATCH 18/57] Add parabolic central flux --- src/equations/equations.jl | 12 +++++++++++- src/equations/heat_equation_2d.jl | 18 +++++------------- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index df2a80b9550..014153eef71 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -56,10 +56,12 @@ default_analysis_integrals(::AbstractEquations) = (entropy_timederivative,) """ flux_central(u_ll, u_rr, orientation, equations::AbstractEquations) + flux_central(u_ll, u_rr, gradients_ll, gradients_rr, orientation, equations::AbstractEquations) The classical central numerical flux `f((u_ll) + f(u_rr)) / 2`. When this flux is used as volume flux, the discretization is equivalent to the classical weak form -DG method (except floating point errors). +DG method (except floating point errors). The second version call the parabolic flux, which requires +the gradients as additional arguments. """ @inline function flux_central(u_ll, u_rr, orientation, equations::AbstractEquations) # Calculate regular 1D fluxes @@ -69,6 +71,14 @@ DG method (except floating point errors). # Average regular fluxes return 0.5 * (f_ll + f_rr) end +@inline function flux_central(u_ll, u_rr, gradients_ll, gradients_rr, orientation, equations::AbstractEquations) + # Calculate regular 1D fluxes + f_ll = calcflux(u_ll, gradients_ll, orientation, equations) + f_rr = calcflux(u_rr, gradients_rr, orientation, equations) + + # Average regular fluxes + return 0.5 * (f_ll + f_rr) +end @inline cons2cons(u, ::AbstractEquations) = u diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index 7aaad8fa06c..19c418027dd 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -29,6 +29,10 @@ function initial_condition_constant(x, t, equation::HeatEquation2D) return @SVector [2.0] end +function initial_condition_linear_xy(x, t, equation::HeatEquation2D) + return @SVector [2*x[1] + 3*x[2]] +end + """ initial_condition_convergence_test(x, t, equations::HeatEquation2D) @@ -51,24 +55,12 @@ end # function source_terms_WHATEVER(u, x, t, equations::HeatEquation2D) -# Calculate hyperbolic 1D flux in axis `orientation` for a single point -@inline function calcflux(u, orientation, equation::HeatEquation2D) - return zero(u) -end - - # Calculate parabolic 1D flux in axis `orientation` for a single point @inline function calcflux(u, gradients, orientation, equation::HeatEquation2D) - return gradients[orientation] + return equation.nu*gradients[orientation] end -function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::HeatEquation2D) - return zero(u_ll) -end - - - @inline have_constant_speed(::HeatEquation2D) = Val(true) @inline function max_abs_speeds(equation::HeatEquation2D) From 1e7a811b3c08d7aef0e0fb0fcb26c9c4543f6f36 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 11:46:21 +0100 Subject: [PATCH 19/57] Add functionality for DG to handle parabolic system --- src/solvers/dg/dg.jl | 4 +- src/solvers/dg/dg_2d.jl | 91 +++++++++++++++++++++++++++++++---------- 2 files changed, 71 insertions(+), 24 deletions(-) diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index 39dadea6670..08aef51c8d1 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -275,9 +275,9 @@ end function create_solver_auxvars(dg::DGSEM) - @unpack basis = dg + @unpack basis, mortar = dg - return DGSEM(basis, flux_central, VolumeIntegralWeakForm(), MortarL2(basis)) + return DGSEM(basis, flux_central, VolumeIntegralWeakForm(), mortar) end diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index 8f57407a663..7bc713cd6ff 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -199,21 +199,21 @@ end function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, mesh::TreeMesh{2}, equations, initial_condition, boundary_conditions, source_terms, - dg::DG, cache) + dg::DG, cache, cache_gradients) # Reset du @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du)) # Calculate volume integral - @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, have_nonconservative_terms(equations), equations, + @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, gradients, have_nonconservative_terms(equations), equations, dg.volume_integral, dg, cache) # Prolong solution to interfaces - @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(cache, u, equations, dg) + @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(cache, cache_gradients, u, gradients, equations, dg) # Calculate interface fluxes @timeit_debug timer() "interface flux" calc_interface_flux!(cache.elements.surface_flux_values, have_nonconservative_terms(equations), equations, - dg, cache) + dg, cache, cache_gradients) # Prolong solution to boundaries @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, u, equations, dg) @@ -236,22 +236,22 @@ function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, @timeit_debug timer() "Jacobian" apply_jacobian!(du, equations, dg, cache) # Add Laplacian - if have_parabolic_terms(equations) == Val(true) - # Compute gradients - @unpack semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy = cache - rhs!(u_ode_grad_xx, u_ode_grad_x, semi_grad_x, t) - rhs!(u_ode_grad_yy, u_ode_grad_y, semi_grad_y, t) + # if have_parabolic_terms(equations) == Val(true) + # # Compute gradients + # @unpack semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy = cache + # rhs!(u_ode_grad_xx, u_ode_grad_x, semi_grad_x, t) + # rhs!(u_ode_grad_yy, u_ode_grad_y, semi_grad_y, t) - grad_xx = wrap_array(u_ode_grad_xx, semi_grad_x) - grad_yy = wrap_array(u_ode_grad_yy, semi_grad_y) + # grad_xx = wrap_array(u_ode_grad_xx, semi_grad_x) + # grad_yy = wrap_array(u_ode_grad_yy, semi_grad_y) - laplacian = grad_xx - laplacian .+= grad_yy + # laplacian = grad_xx + # laplacian .+= grad_yy - # Add Laplacian to rhs - @unpack nu = equations - du .+= nu .* laplacian - end + # # Add Laplacian to rhs + # @unpack nu = equations + # du .+= nu .* laplacian + # end # Calculate source terms @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache) @@ -297,17 +297,17 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, gradients, Threads.@threads for element in eachelement(dg, cache) for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, element) - grad_x_node = get_node_vars(gradients[1], equations, dg, i, j, element) - grad_y_node = get_node_vars(gradients[2], equations, dg, i, j, element) - grad_node = (grad_x_node, grad_y_node) + gradients_x_node = get_node_vars(gradients[1], equations, dg, i, j, element) + gradients_y_node = get_node_vars(gradients[2], equations, dg, i, j, element) + gradients_node = (gradients_x_node, gradients_y_node) - flux1 = calcflux(u_node, grad_node, 1, equations) + flux1 = calcflux(u_node, gradients_node, 1, equations) for ii in eachnode(dg) integral_contribution = derivative_dhat[ii, i] * flux1 add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element) end - flux2 = calcflux(u_node, grad_node, 2, equations) + flux2 = calcflux(u_node, gradients_node, 2, equations) for jj in eachnode(dg) integral_contribution = derivative_dhat[jj, j] * flux2 add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element) @@ -689,6 +689,14 @@ function prolong2interfaces!(cache, u::AbstractArray{<:Any,4}, equations, dg::DG return nothing end +function prolong2interfaces!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, equations, dg::DG) + prolong2interfaces!(cache, u, equations, dg) + prolong2interfaces!(cache_gradients[1], gradients[1], equations, dg) + prolong2interfaces!(cache_gradients[2], gradients[2], equations, dg) + + return nothing +end + function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, nonconservative_terms::Val{false}, equations, dg::DG, cache) @@ -722,6 +730,45 @@ function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, return nothing end +function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, + nonconservative_terms::Val{false}, equations, + dg::DG, cache, cache_gradients) + @unpack surface_flux = dg + @unpack u, neighbor_ids, orientations = cache.interfaces + gradients_x = cache_gradients[1].interfaces.u + gradients_y = cache_gradients[2].interfaces.u + + Threads.@threads for interface in eachinterface(dg, cache) + # Get neighboring elements + left_id = neighbor_ids[1, interface] + right_id = neighbor_ids[2, interface] + + # Determine interface direction with respect to elements: + # orientation = 1: left -> 2, right -> 1 + # orientation = 2: left -> 4, right -> 3 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface) + gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_x, equations, dg, i, interface) + gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_y, equations, dg, i, interface) + gradients_ll = (gradients_x_ll, gradients_y_ll) + gradients_rr = (gradients_x_rr, gradients_y_rr) + flux = surface_flux(u_ll, u_rr, gradients_ll, gradients_rr, orientations[interface], equations) + + # Copy flux to left and right element storage + for v in eachvariable(equations) + surface_flux_values[v, i, left_direction, left_id] = flux[v] + surface_flux_values[v, i, right_direction, right_id] = flux[v] + end + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, nonconservative_terms::Val{true}, equations, dg::DG, cache) From 98f619fdc30a9b135ba286c2ffa02ebe6209b809 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 11:47:05 +0100 Subject: [PATCH 20/57] Pass gradients cache to rhs! --- src/semidiscretization/semidiscretization_parabolic_auxvars.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index c2b6630deb8..69f9e898dc7 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -173,10 +173,11 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationParabolicAuxVars, t) rhs!(u_ode_grad_y, u_ode, semi_grad_y, t) gradients = (wrap_array(u_ode_grad_x, mesh, equations, solver, cache), wrap_array(u_ode_grad_y, mesh, equations, solver, cache)) + cache_gradients = (semi_grad_x.cache, semi_grad_y.cache) # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() - @timeit_debug timer() "rhs!" rhs!(du, u, gradients, t, mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache) + @timeit_debug timer() "rhs!" rhs!(du, u, gradients, t, mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache, cache_gradients) runtime = time_ns() - time_start put!(semi.performance_counter, runtime) From d748522dfb23012ae5edbbd1be14460916e014a3 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 11:47:12 +0100 Subject: [PATCH 21/57] Use central flux --- examples/2d/elixir_heat_basic.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/2d/elixir_heat_basic.jl b/examples/2d/elixir_heat_basic.jl index 1baf319e46e..f878032c295 100644 --- a/examples/2d/elixir_heat_basic.jl +++ b/examples/2d/elixir_heat_basic.jl @@ -8,8 +8,8 @@ using Trixi nu = 1.2e-2 equations = HeatEquation2D(nu) -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(3, flux_lax_friedrichs) +# Create DG solver with polynomial degree = 3 and the central flux as surface flux +solver = DGSEM(3, flux_central) coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y)) coordinates_max = ( 1, 1) # maximum coordinates (max(x), max(y)) @@ -56,7 +56,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); + save_everystep=false, callback=callbacks, maxiters=1e5); # Print the timer summary summary_callback() From c9b813b6f70c05ab42a89158bfb7016d3f3e5c56 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 13:03:20 +0100 Subject: [PATCH 22/57] Fix sign of the flux --- src/equations/heat_equation_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index 19c418027dd..79c8df342c7 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -57,7 +57,7 @@ end # Calculate parabolic 1D flux in axis `orientation` for a single point @inline function calcflux(u, gradients, orientation, equation::HeatEquation2D) - return equation.nu*gradients[orientation] + return -equation.nu*gradients[orientation] end From 7c7e295125daa690333946cc196c27d2dca6eb91 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 17:44:48 +0100 Subject: [PATCH 23/57] Add parabolic flux to advection-diffusion equation --- src/equations/linear_advection_diffusion_2d.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index 77eef096adb..214871540c9 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -64,12 +64,17 @@ end # function source_terms_WHATEVER(u, x, t, equations::LinearAdvectionDiffusionEquation2D) -# Calculate 1D flux in for a single point +# Calculate hyperbolic 1D flux in axis `orientation` for a single point @inline function calcflux(u, orientation, equation::LinearAdvectionDiffusionEquation2D) a = equation.advectionvelocity[orientation] return a * u end +# Calculate parabolic 1D flux in axis `orientation` for a single point +@inline function calcflux(u, gradients, orientation, equation::LinearAdvectionDiffusionEquation2D) + return -equation.nu*gradients[orientation] +end + function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::LinearAdvectionDiffusionEquation2D) a = equation.advectionvelocity[orientation] From a9d11c7cba10d4184934521e44ed67ae7ff30319 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 17:45:05 +0100 Subject: [PATCH 24/57] Add hyperbolic-parabolic semidiscretization --- src/Trixi.jl | 2 + ...semidiscretization_hyperbolic_parabolic.jl | 123 ++++++++++++++++++ 2 files changed, 125 insertions(+) create mode 100644 src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index 6301ae665dc..dcc0bb8163a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -54,6 +54,7 @@ include("solvers/dg/dg.jl") include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") include("semidiscretization/semidiscretization_parabolic_auxvars.jl") +include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") @@ -120,6 +121,7 @@ export nelements, nnodes, nvariables, eachelement, eachnode, eachvariable export SemidiscretizationHyperbolic, SemidiscretizationParabolicAuxVars, + SemidiscretizationHyperbolicParabolic, semidiscretize, compute_coefficients, integrate export SemidiscretizationEulerGravity, ParametersEulerGravity, diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl new file mode 100644 index 00000000000..a8ff2ca1e04 --- /dev/null +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -0,0 +1,123 @@ + +# TODO: Taal refactor, Mesh<:AbstractMesh{NDIMS}, Equations<:AbstractEquations{NDIMS} ? +""" + SemidiscretizationHyperbolicParabolic + +A struct containing everything needed to describe a spatial semidiscretization +of a hyperbolic-parabolic conservation law. +""" +struct SemidiscretizationHyperbolicParabolic{SemiHyperbolic, SemiParabolic, Cache} <: AbstractSemidiscretization + semi_hyperbolic::SemiHyperbolic + semi_parabolic::SemiParabolic + cache::Cache + performance_counter::PerformanceCounter + + function SemidiscretizationHyperbolicParabolic{SemiHyperbolic, SemiParabolic, Cache}( + semi_hyperbolic::SemiHyperbolic, semi_parabolic::SemiParabolic, + cache::Cache) where {SemiHyperbolic, SemiParabolic, Cache} + + performance_counter = PerformanceCounter() + + new(semi_hyperbolic, semi_parabolic, cache, performance_counter) + end +end + +""" + SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic) + +Construct a semidiscretization of a hyperbolic-parabolic PDE by combining the purely hyperbolic and +purely parabolic semidiscretizations. +""" +function SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic) + cache = (; du_ode_parabolic=Vector{real(semi_parabolic)}()) + + SemidiscretizationHyperbolicParabolic{typeof(semi_hyperbolic), typeof(semi_parabolic), typeof(cache)}( + semi_hyperbolic, semi_parabolic, cache) +end + + +function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) + print(io, "SemidiscretizationHyperbolicParabolic(") + print(io, semi.semi_hyperbolic) + print(io, ", ", semi.semi_parabolic) + print(io, ", cache(") + for (idx,key) in enumerate(keys(semi.cache)) + idx > 1 && print(io, " ") + print(io, key) + end + print(io, "))") +end + +function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperbolicParabolic) + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationHyperbolicParabolic") + summary_line(io, "hyperbolic system", semi.semi_hyperbolic) + summary_line(io, "parabolic system", semi.semi_parabolic) + summary_footer(io) + end +end + + +@inline Base.ndims(semi::SemidiscretizationHyperbolicParabolic) = ndims(semi.semi_hyperbolic.mesh) + +@inline nvariables(semi::SemidiscretizationHyperbolicParabolic) = nvariables(semi.semi_hyperbolic.equations) + +@inline Base.real(semi::SemidiscretizationHyperbolicParabolic) = real(semi.semi_hyperbolic.solver) + + +@inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolicParabolic) + @unpack mesh, equations, solver, cache = semi.semi_hyperbolic + return mesh, equations, solver, cache +end + + +function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationHyperbolicParabolic) + @unpack mesh, equations, initial_condition, solver, cache = semi.semi_hyperbolic + u = wrap_array(u_ode, mesh, equations, solver, cache) + + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) +end + +function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationHyperbolicParabolic) + @unpack mesh, equations, initial_condition, solver, cache = semi.semi_hyperbolic + + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) +end + + +function compute_coefficients(t, semi::SemidiscretizationHyperbolicParabolic) + compute_coefficients(semi.semi_hyperbolic.initial_condition, t, semi.semi_hyperbolic) +end + +function compute_coefficients!(u_ode::AbstractVector, t, semi::SemidiscretizationHyperbolicParabolic) + compute_coefficients!(u_ode, semi.semi_hyperbolic.initial_condition, t, semi.semi_hyperbolic) +end + + +function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) + @unpack du_ode_parabolic = semi.cache + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + + # Resize the storage for the parabolic time derivative in case it was changed + resize!(du_ode_parabolic, length(du_ode)) + + @timeit_debug timer() "rhs!" begin + # Compute hyperbolic time derivative + rhs!(du_ode, u_ode, semi.semi_hyperbolic, t) + + # Compute parabolic time derivative + rhs!(du_ode_parabolic, u_ode, semi.semi_parabolic, t) + + # Add parabolic update to hyperbolic time derivative + du_ode .+= du_ode_parabolic + end + + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end From 317a76bc2d4c3cda796c66032128afb55aff034a Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 17:48:25 +0100 Subject: [PATCH 25/57] Change adv-diff elixir such that it works again (without mortars yet) --- examples/2d/elixir_advdiff_basic.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/examples/2d/elixir_advdiff_basic.jl b/examples/2d/elixir_advdiff_basic.jl index 5c8bf9e8d23..dffb579739e 100644 --- a/examples/2d/elixir_advdiff_basic.jl +++ b/examples/2d/elixir_advdiff_basic.jl @@ -10,7 +10,8 @@ nu = 1.2e-2 equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(4, flux_lax_friedrichs) +polydeg = 3 # 4 +solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs) coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y)) coordinates_max = ( 1, 1) # maximum coordinates (max(x), max(y)) @@ -21,12 +22,17 @@ refinement_patches = ( # Create a uniformely refined mesh with periodic boundaries mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=2, - refinement_patches=refinement_patches, + # refinement_patches=refinement_patches, n_cells_max=30_000) # set maximum capacity of tree data structure # A semidiscretization collects data structures and functions for the spatial discretization initial_condition = initial_condition_convergence_test -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +semi_hyperbolic = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver_hyperbolic) + +solver_parabolic = DGSEM(polydeg, flux_central) +semi_parabolic = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver_parabolic) + +semi = SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic) ############################################################################### From 3aca3e20014e44307e490f7b190f42534250b5f5 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 17:52:33 +0100 Subject: [PATCH 26/57] Move parabolic implementaiton to different file --- src/solvers/dg/dg.jl | 1 + src/solvers/dg/dg_2d.jl | 167 ------------------------------ src/solvers/dg/dg_2d_parabolic.jl | 126 ++++++++++++++++++++++ 3 files changed, 127 insertions(+), 167 deletions(-) create mode 100644 src/solvers/dg/dg_2d_parabolic.jl diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index 08aef51c8d1..739fb8c9a2b 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -296,6 +296,7 @@ include("dg_1d.jl") include("containers_2d.jl") include("dg_2d.jl") include("dg_2d_parallel.jl") +include("dg_2d_parabolic.jl") # 3D DG implementation include("containers_3d.jl") diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index 7bc713cd6ff..8fcc2e6de45 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -23,7 +23,6 @@ function create_cache(mesh::TreeMesh{2}, equations::AbstractEquations{2}, # Add specialized parts of the cache required to compute the volume integral etc. cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg)...) cache = (;cache..., create_cache(mesh, equations, dg.mortar)...) - cache = (;cache..., create_cache(mesh, equations, dg, have_parabolic_terms(equations), cache)...) return cache end @@ -31,31 +30,6 @@ end # The methods below are specialized on the volume integral type # and called from the basic `create_cache` method at the top. -function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val{false}, cache) - NamedTuple() -end - -function create_cache(mesh::TreeMesh{2}, equations, dg::DG, parabolic_terms::Val{true}, cache) - return NamedTuple() - -# solver_grad_n = DGSEM(polydeg(dg), flux_central) -# -# equations_grad_x = GradientEquations2D(nvariables(equations), 1) -# semi_grad_x = SemidiscretizationHyperbolic(mesh, equations_grad_x, -# initial_condition_constant, solver_grad_n) -# -# equations_grad_y = GradientEquations2D(nvariables(equations), 2) -# semi_grad_y = SemidiscretizationHyperbolic(mesh, equations_grad_y, -# initial_condition_constant, solver_grad_n) -# -# u_ode_grad_x = compute_coefficients(nothing, semi_grad_x) -# u_ode_grad_y = compute_coefficients(nothing, semi_grad_y) -# u_ode_grad_xx = compute_coefficients(nothing, semi_grad_x) -# u_ode_grad_yy = compute_coefficients(nothing, semi_grad_y) -# -# return (; semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy) -end - function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG) create_cache(mesh, have_nonconservative_terms(equations), equations, volume_integral, dg) end @@ -196,69 +170,6 @@ function rhs!(du::AbstractArray{<:Any,4}, u, t, return nothing end -function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, - mesh::TreeMesh{2}, equations, - initial_condition, boundary_conditions, source_terms, - dg::DG, cache, cache_gradients) - # Reset du - @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du)) - - # Calculate volume integral - @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, gradients, have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) - - # Prolong solution to interfaces - @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(cache, cache_gradients, u, gradients, equations, dg) - - # Calculate interface fluxes - @timeit_debug timer() "interface flux" calc_interface_flux!(cache.elements.surface_flux_values, - have_nonconservative_terms(equations), equations, - dg, cache, cache_gradients) - - # Prolong solution to boundaries - @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, u, equations, dg) - - # Calculate boundary fluxes - @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, equations, dg) - - # Prolong solution to mortars - @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, u, equations, dg.mortar, dg) - - # Calculate mortar fluxes - @timeit_debug timer() "mortar flux" calc_mortar_flux!(cache.elements.surface_flux_values, - have_nonconservative_terms(equations), equations, - dg.mortar, dg, cache) - - # Calculate surface integrals - @timeit_debug timer() "surface integral" calc_surface_integral!(du, equations, dg, cache) - - # Apply Jacobian from mapping to reference element - @timeit_debug timer() "Jacobian" apply_jacobian!(du, equations, dg, cache) - - # Add Laplacian - # if have_parabolic_terms(equations) == Val(true) - # # Compute gradients - # @unpack semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy = cache - # rhs!(u_ode_grad_xx, u_ode_grad_x, semi_grad_x, t) - # rhs!(u_ode_grad_yy, u_ode_grad_y, semi_grad_y, t) - - # grad_xx = wrap_array(u_ode_grad_xx, semi_grad_x) - # grad_yy = wrap_array(u_ode_grad_yy, semi_grad_y) - - # laplacian = grad_xx - # laplacian .+= grad_yy - - # # Add Laplacian to rhs - # @unpack nu = equations - # du .+= nu .* laplacian - # end - - # Calculate source terms - @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache) - - return nothing -end - function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms::Val{false}, equations, @@ -288,37 +199,6 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, end -function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, gradients, - nonconservative_terms::Val{false}, equations, - volume_integral::VolumeIntegralWeakForm, - dg::DGSEM, cache) - @unpack derivative_dhat = dg.basis - - Threads.@threads for element in eachelement(dg, cache) - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, element) - gradients_x_node = get_node_vars(gradients[1], equations, dg, i, j, element) - gradients_y_node = get_node_vars(gradients[2], equations, dg, i, j, element) - gradients_node = (gradients_x_node, gradients_y_node) - - flux1 = calcflux(u_node, gradients_node, 1, equations) - for ii in eachnode(dg) - integral_contribution = derivative_dhat[ii, i] * flux1 - add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element) - end - - flux2 = calcflux(u_node, gradients_node, 2, equations) - for jj in eachnode(dg) - integral_contribution = derivative_dhat[jj, j] * flux2 - add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element) - end - end - end - - return nothing -end - - # Calculate 2D twopoint flux (element version) @inline function calcflux_twopoint!(f1, f2, u::AbstractArray{<:Any,4}, element, volume_flux, equations, dg::DG, cache) @@ -689,14 +569,6 @@ function prolong2interfaces!(cache, u::AbstractArray{<:Any,4}, equations, dg::DG return nothing end -function prolong2interfaces!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, equations, dg::DG) - prolong2interfaces!(cache, u, equations, dg) - prolong2interfaces!(cache_gradients[1], gradients[1], equations, dg) - prolong2interfaces!(cache_gradients[2], gradients[2], equations, dg) - - return nothing -end - function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, nonconservative_terms::Val{false}, equations, dg::DG, cache) @@ -730,45 +602,6 @@ function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, return nothing end -function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, - nonconservative_terms::Val{false}, equations, - dg::DG, cache, cache_gradients) - @unpack surface_flux = dg - @unpack u, neighbor_ids, orientations = cache.interfaces - gradients_x = cache_gradients[1].interfaces.u - gradients_y = cache_gradients[2].interfaces.u - - Threads.@threads for interface in eachinterface(dg, cache) - # Get neighboring elements - left_id = neighbor_ids[1, interface] - right_id = neighbor_ids[2, interface] - - # Determine interface direction with respect to elements: - # orientation = 1: left -> 2, right -> 1 - # orientation = 2: left -> 4, right -> 3 - left_direction = 2 * orientations[interface] - right_direction = 2 * orientations[interface] - 1 - - for i in eachnode(dg) - # Call pointwise Riemann solver - u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface) - gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_x, equations, dg, i, interface) - gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_y, equations, dg, i, interface) - gradients_ll = (gradients_x_ll, gradients_y_ll) - gradients_rr = (gradients_x_rr, gradients_y_rr) - flux = surface_flux(u_ll, u_rr, gradients_ll, gradients_rr, orientations[interface], equations) - - # Copy flux to left and right element storage - for v in eachvariable(equations) - surface_flux_values[v, i, left_direction, left_id] = flux[v] - surface_flux_values[v, i, right_direction, right_id] = flux[v] - end - end - end - - return nothing -end - function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, nonconservative_terms::Val{true}, equations, dg::DG, cache) diff --git a/src/solvers/dg/dg_2d_parabolic.jl b/src/solvers/dg/dg_2d_parabolic.jl new file mode 100644 index 00000000000..44b0ccb83b7 --- /dev/null +++ b/src/solvers/dg/dg_2d_parabolic.jl @@ -0,0 +1,126 @@ +# This file collects all methods that have been updated to work with parabolic systems of equations + +function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, + mesh::TreeMesh{2}, equations, + initial_condition, boundary_conditions, source_terms, + dg::DG, cache, cache_gradients) + # Reset du + @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du)) + + # Calculate volume integral + @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, gradients, have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) + + # Prolong solution to interfaces + @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(cache, cache_gradients, u, gradients, equations, dg) + + # Calculate interface fluxes + @timeit_debug timer() "interface flux" calc_interface_flux!(cache.elements.surface_flux_values, + have_nonconservative_terms(equations), equations, + dg, cache, cache_gradients) + + # Prolong solution to boundaries + @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, u, equations, dg) + + # Calculate boundary fluxes + @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, equations, dg) + + # Prolong solution to mortars + @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, u, equations, dg.mortar, dg) + + # Calculate mortar fluxes + @timeit_debug timer() "mortar flux" calc_mortar_flux!(cache.elements.surface_flux_values, + have_nonconservative_terms(equations), equations, + dg.mortar, dg, cache) + + # Calculate surface integrals + @timeit_debug timer() "surface integral" calc_surface_integral!(du, equations, dg, cache) + + # Apply Jacobian from mapping to reference element + @timeit_debug timer() "Jacobian" apply_jacobian!(du, equations, dg, cache) + + # Calculate source terms + @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache) + + return nothing +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, gradients, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralWeakForm, + dg::DGSEM, cache) + @unpack derivative_dhat = dg.basis + + Threads.@threads for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + gradients_x_node = get_node_vars(gradients[1], equations, dg, i, j, element) + gradients_y_node = get_node_vars(gradients[2], equations, dg, i, j, element) + gradients_node = (gradients_x_node, gradients_y_node) + + flux1 = calcflux(u_node, gradients_node, 1, equations) + for ii in eachnode(dg) + integral_contribution = derivative_dhat[ii, i] * flux1 + add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element) + end + + flux2 = calcflux(u_node, gradients_node, 2, equations) + for jj in eachnode(dg) + integral_contribution = derivative_dhat[jj, j] * flux2 + add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element) + end + end + end + + return nothing +end + + +function prolong2interfaces!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, equations, dg::DG) + prolong2interfaces!(cache, u, equations, dg) + prolong2interfaces!(cache_gradients[1], gradients[1], equations, dg) + prolong2interfaces!(cache_gradients[2], gradients[2], equations, dg) + + return nothing +end + + +function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, + nonconservative_terms::Val{false}, equations, + dg::DG, cache, cache_gradients) + @unpack surface_flux = dg + @unpack u, neighbor_ids, orientations = cache.interfaces + gradients_x = cache_gradients[1].interfaces.u + gradients_y = cache_gradients[2].interfaces.u + + Threads.@threads for interface in eachinterface(dg, cache) + # Get neighboring elements + left_id = neighbor_ids[1, interface] + right_id = neighbor_ids[2, interface] + + # Determine interface direction with respect to elements: + # orientation = 1: left -> 2, right -> 1 + # orientation = 2: left -> 4, right -> 3 + left_direction = 2 * orientations[interface] + right_direction = 2 * orientations[interface] - 1 + + for i in eachnode(dg) + # Call pointwise Riemann solver + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface) + gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_x, equations, dg, i, interface) + gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_y, equations, dg, i, interface) + gradients_ll = (gradients_x_ll, gradients_y_ll) + gradients_rr = (gradients_x_rr, gradients_y_rr) + flux = surface_flux(u_ll, u_rr, gradients_ll, gradients_rr, orientations[interface], equations) + + # Copy flux to left and right element storage + for v in eachvariable(equations) + surface_flux_values[v, i, left_direction, left_id] = flux[v] + surface_flux_values[v, i, right_direction, right_id] = flux[v] + end + end + end + + return nothing +end From 9e996770ef0da66dfb6bbd76d8a97b668e1501d4 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 18:23:50 +0100 Subject: [PATCH 27/57] Implement mortar fluxes --- src/solvers/dg/dg_2d_parabolic.jl | 69 +++++++++++++++++++++++++++++-- 1 file changed, 66 insertions(+), 3 deletions(-) diff --git a/src/solvers/dg/dg_2d_parabolic.jl b/src/solvers/dg/dg_2d_parabolic.jl index 44b0ccb83b7..ebe285bdb60 100644 --- a/src/solvers/dg/dg_2d_parabolic.jl +++ b/src/solvers/dg/dg_2d_parabolic.jl @@ -26,12 +26,13 @@ function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, equations, dg) # Prolong solution to mortars - @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, u, equations, dg.mortar, dg) + @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, cache_gradients, u, gradients, + equations, dg.mortar, dg) # Calculate mortar fluxes @timeit_debug timer() "mortar flux" calc_mortar_flux!(cache.elements.surface_flux_values, have_nonconservative_terms(equations), equations, - dg.mortar, dg, cache) + dg.mortar, dg, cache, cache_gradients) # Calculate surface integrals @timeit_debug timer() "surface integral" calc_surface_integral!(du, equations, dg, cache) @@ -77,7 +78,8 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, gradients, end -function prolong2interfaces!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, equations, dg::DG) +function prolong2interfaces!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, + equations, dg::DG) prolong2interfaces!(cache, u, equations, dg) prolong2interfaces!(cache_gradients[1], gradients[1], equations, dg) prolong2interfaces!(cache_gradients[2], gradients[2], equations, dg) @@ -124,3 +126,64 @@ function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, return nothing end + + +function prolong2mortars!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, + equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) + prolong2mortars!(cache, u, equations, mortar_l2, dg) + prolong2mortars!(cache_gradients[1], gradients[1], equations, mortar_l2, dg) + prolong2mortars!(cache_gradients[2], gradients[2], equations, mortar_l2, dg) + + return nothing +end + + +function calc_mortar_flux!(surface_flux_values::AbstractArray{<:Any,4}, + nonconservative_terms::Val{false}, equations, + mortar_l2::LobattoLegendreMortarL2, dg::DG, cache, cache_gradients) + @unpack neighbor_ids, u_lower, u_upper, orientations = cache.mortars + @unpack fstar_upper_threaded, fstar_lower_threaded = cache + gradients_x_upper = cache_gradients[1].mortars.u_upper + gradients_x_lower = cache_gradients[1].mortars.u_lower + gradients_y_upper = cache_gradients[2].mortars.u_upper + gradients_y_lower = cache_gradients[2].mortars.u_lower + gradients_upper = (gradients_x_upper, gradients_y_upper) + gradients_lower = (gradients_x_lower, gradients_y_lower) + + Threads.@threads for mortar in eachmortar(dg, cache) + # Choose thread-specific pre-allocated container + fstar_upper = fstar_upper_threaded[Threads.threadid()] + fstar_lower = fstar_lower_threaded[Threads.threadid()] + + # Calculate fluxes + orientation = orientations[mortar] + calc_fstar!(fstar_upper, equations, dg, u_upper, gradients_upper, mortar, orientation) + calc_fstar!(fstar_lower, equations, dg, u_lower, gradients_lower, mortar, orientation) + + mortar_fluxes_to_elements!(surface_flux_values, equations, mortar_l2, dg, cache, + mortar, fstar_upper, fstar_lower) + end + + return nothing +end + + +@inline function calc_fstar!(destination::AbstractArray{<:Any,2}, equations, dg::DGSEM, + u_interfaces, gradients_interfaces, interface, orientation) + @unpack surface_flux = dg + + for i in eachnode(dg) + # Call pointwise two-point numerical flux function + u_ll, u_rr = get_surface_node_vars(u_interfaces, equations, dg, i, interface) + gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_interfaces[1], equations, dg, i, interface) + gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_interfaces[2], equations, dg, i, interface) + gradients_ll = (gradients_x_ll, gradients_y_ll) + gradients_rr = (gradients_x_rr, gradients_y_rr) + flux = surface_flux(u_ll, u_rr, gradients_x_ll, gradients_y_ll, orientation, equations) + + # Copy flux to left and right element storage + set_node_vars!(destination, flux, equations, dg, i) + end + + return nothing +end From 3efdb0d935a21023f30e22fa5ade08dd1b31fbf3 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 18:24:33 +0100 Subject: [PATCH 28/57] re-enable mortars for adv-div; convergence tests works but different results in comparison to "naive" implementation --- examples/2d/elixir_advdiff_basic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/2d/elixir_advdiff_basic.jl b/examples/2d/elixir_advdiff_basic.jl index dffb579739e..7af928a63b7 100644 --- a/examples/2d/elixir_advdiff_basic.jl +++ b/examples/2d/elixir_advdiff_basic.jl @@ -10,7 +10,7 @@ nu = 1.2e-2 equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -polydeg = 3 # 4 +polydeg = 4 solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs) coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y)) @@ -22,7 +22,7 @@ refinement_patches = ( # Create a uniformely refined mesh with periodic boundaries mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=2, - # refinement_patches=refinement_patches, + refinement_patches=refinement_patches, n_cells_max=30_000) # set maximum capacity of tree data structure # A semidiscretization collects data structures and functions for the spatial discretization From 44d1bd99a8bbbe94e897880af8c789ae80b65e19 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 18:29:46 +0100 Subject: [PATCH 29/57] Remove `have_parabolic_terms` at it is currently not needed --- src/equations/equations.jl | 1 - src/equations/heat_equation_2d.jl | 2 -- src/equations/linear_advection_diffusion_2d.jl | 2 -- src/semidiscretization/semidiscretization_parabolic_auxvars.jl | 1 - 4 files changed, 6 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 014153eef71..7cedc65ec42 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -48,7 +48,6 @@ function calcflux(u, orientation, equations) end # set sensible default values that may be overwritten by specific equations have_nonconservative_terms(::AbstractEquations) = Val(false) have_constant_speed(::AbstractEquations) = Val(false) -have_parabolic_terms(::AbstractEquations) = Val(false) default_analysis_errors(::AbstractEquations) = (:l2_error, :linf_error) default_analysis_integrals(::AbstractEquations) = (entropy_timederivative,) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index 79c8df342c7..e173cb8f9c0 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -16,8 +16,6 @@ get_name(::HeatEquation2D) = "HeatEquation2D" varnames(::typeof(cons2cons), ::HeatEquation2D) = SVector("scalar") varnames(::typeof(cons2prim), ::HeatEquation2D) = SVector("scalar") -have_parabolic_terms(::HeatEquation2D) = Val(true) - # Set initial conditions at physical location `x` for time `t` """ diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index 214871540c9..33cfe6a79fc 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -26,8 +26,6 @@ get_name(::LinearAdvectionDiffusionEquation2D) = "LinearAdvectionDiffusionEquati varnames(::typeof(cons2cons), ::LinearAdvectionDiffusionEquation2D) = SVector("scalar") varnames(::typeof(cons2prim), ::LinearAdvectionDiffusionEquation2D) = SVector("scalar") -have_parabolic_terms(::LinearAdvectionDiffusionEquation2D) = Val(true) - # Set initial conditions at physical location `x` for time `t` """ diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index 69f9e898dc7..8ea9c16d78d 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -29,7 +29,6 @@ struct SemidiscretizationParabolicAuxVars{Mesh, Equations, InitialCondition, Bou source_terms::SourceTerms, solver::Solver, cache::Cache) where {Mesh, Equations, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache} @assert ndims(mesh) == ndims(equations) - @assert have_parabolic_terms(equations) == Val(true) "Cannot create parabolic semidiscretization for purely hyperbolic system" performance_counter = PerformanceCounter() From d477c137e44b3ddf2b0933f623003b50778ba787 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 19:51:03 +0100 Subject: [PATCH 30/57] Add convenience method `SemidiscretizationHyperbolicParabolicBR1` --- src/Trixi.jl | 2 +- ...semidiscretization_hyperbolic_parabolic.jl | 31 +++++++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index dcc0bb8163a..0696932b6be 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -121,7 +121,7 @@ export nelements, nnodes, nvariables, eachelement, eachnode, eachvariable export SemidiscretizationHyperbolic, SemidiscretizationParabolicAuxVars, - SemidiscretizationHyperbolicParabolic, + SemidiscretizationHyperbolicParabolic, SemidiscretizationHyperbolicParabolicBR1, semidiscretize, compute_coefficients, integrate export SemidiscretizationEulerGravity, ParametersEulerGravity, diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index a8ff2ca1e04..95b83552013 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -121,3 +121,34 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) return nothing end + + +""" + SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver::DGSEM; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic) + +Convenience function to construct a semidiscretization of a hyperbolic-parabolic PDE with the BR1 +scheme for a DGSEM solver. The passed arguments must correspond to the hyperbolic setup, including +the `solver`. Internally, a matching solver for the parabolic system is created and both a +hyperbolic and a parabolic semidiscretization are created and passed to a +`SemidiscretizationHyperbolicParabolic`. +""" +function SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver::DGSEM; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + RealT=real(solver)) + semi_hyperbolic = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + source_terms=source_terms, + boundary_conditions=boundary_conditions, + RealT=RealT) + + solver_parabolic = DGSEM(solver.basis, flux_central, solver.volume_integral, solver.mortar) + semi_parabolic = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, + solver_parabolic, + source_terms=source_terms, + boundary_conditions=boundary_conditions, + RealT=RealT) + + return SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic) +end From 3125aa53a983f28ab636dc9d45a8c7bfcc92d591 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 19:51:38 +0100 Subject: [PATCH 31/57] Add convenience method `SemidiscretizationHyperbolicParabolicBR1` --- examples/2d/elixir_advdiff_basic.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/examples/2d/elixir_advdiff_basic.jl b/examples/2d/elixir_advdiff_basic.jl index 7af928a63b7..ec1e6cbdcef 100644 --- a/examples/2d/elixir_advdiff_basic.jl +++ b/examples/2d/elixir_advdiff_basic.jl @@ -27,12 +27,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # A semidiscretization collects data structures and functions for the spatial discretization initial_condition = initial_condition_convergence_test -semi_hyperbolic = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver_hyperbolic) - -solver_parabolic = DGSEM(polydeg, flux_central) -semi_parabolic = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver_parabolic) - -semi = SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic) +semi = SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver_hyperbolic) ############################################################################### From 7e27661eb5fde4ed25e42a6c752838f4ca1f2033 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 19:52:28 +0100 Subject: [PATCH 32/57] Change test reference values although they should not change; however, we get EOCs of ~5.37 and 5.15 for L2/Linf, so it cannot be all bad --- test/test_examples_2d_advdiff.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/test_examples_2d_advdiff.jl b/test/test_examples_2d_advdiff.jl index df6bd8de548..6710b7606d1 100644 --- a/test/test_examples_2d_advdiff.jl +++ b/test/test_examples_2d_advdiff.jl @@ -11,8 +11,12 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") @testset "Scalar advection-diffusion" begin @testset "elixir_advdiff_basic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_basic.jl"), - l2 = [6.707801888154751e-5], - linf = [0.0005912932758251888]) + # FIXME: These values changed from the naive implementation for the parabolic terms to the + # "proper" implementaton by O(1e-6): Find out why and if it is OK + # l2 = [6.707801888154751e-5], + # linf = [0.0005912932758251888]) + l2 = [6.447564451432685e-5], + linf = [0.0005850692252580281]) end end From a343eca6d10a52749715c829c3ee499cf98f6b87 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 20:21:13 +0100 Subject: [PATCH 33/57] Fix call syntax in semidiscretizations for `calc_error_norms` --- .../semidiscretization_hyperbolic_parabolic.jl | 8 ++++---- .../semidiscretization_parabolic_auxvars.jl | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 95b83552013..e1930fb74d8 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -73,17 +73,17 @@ end end -function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationHyperbolicParabolic) +function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis) @unpack mesh, equations, initial_condition, solver, cache = semi.semi_hyperbolic u = wrap_array(u_ode, mesh, equations, solver, cache) - calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end -function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationHyperbolicParabolic) +function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis) @unpack mesh, equations, initial_condition, solver, cache = semi.semi_hyperbolic - calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index 8ea9c16d78d..c87db31be31 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -137,17 +137,17 @@ end end -function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationParabolicAuxVars) +function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationParabolicAuxVars, cache_analysis) @unpack mesh, equations, initial_condition, solver, cache = semi u = wrap_array(u_ode, mesh, equations, solver, cache) - calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end -function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationParabolicAuxVars) +function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationParabolicAuxVars, cache_analysis) @unpack mesh, equations, initial_condition, solver, cache = semi - calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache) + calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end From fab783c47b8bcf4593dce0736edbab17afce31cf Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 18 Jan 2021 20:21:44 +0100 Subject: [PATCH 34/57] Add Gauss pulse IC to heat equation/advdiff equation --- src/equations/heat_equation_2d.jl | 15 +++++++++++---- src/equations/linear_advection_diffusion_2d.jl | 11 +++++++++++ 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index e173cb8f9c0..b3605fb69b5 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -27,10 +27,6 @@ function initial_condition_constant(x, t, equation::HeatEquation2D) return @SVector [2.0] end -function initial_condition_linear_xy(x, t, equation::HeatEquation2D) - return @SVector [2*x[1] + 3*x[2]] -end - """ initial_condition_convergence_test(x, t, equations::HeatEquation2D) @@ -49,6 +45,17 @@ function initial_condition_convergence_test(x, t, equation::HeatEquation2D) end +""" + initial_condition_gauss(x, t, equation::HeatEquation2D) + +A Gaussian pulse used together with +[`boundary_condition_gauss`](@ref). +""" +function initial_condition_gauss(x, t, equation::HeatEquation2D) + return @SVector [exp(-(x[1]^2 + x[2]^2))] +end + + # Pre-defined source terms should be implemented as # function source_terms_WHATEVER(u, x, t, equations::HeatEquation2D) diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index 33cfe6a79fc..417611adc52 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -58,6 +58,17 @@ function initial_condition_convergence_test(x, t, equation::LinearAdvectionDiffu end +""" + initial_condition_gauss(x, t, equation::LinearAdvectionDiffusionEquation2D) + +A Gaussian pulse used together with +[`boundary_condition_gauss`](@ref). +""" +function initial_condition_gauss(x, t, equation::LinearAdvectionDiffusionEquation2D) + return @SVector [exp(-(x[1]^2 + x[2]^2))] +end + + # Pre-defined source terms should be implemented as # function source_terms_WHATEVER(u, x, t, equations::LinearAdvectionDiffusionEquation2D) From fa9b8f0900d947ddd8854a975257d855e7eaad1e Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 06:08:53 +0100 Subject: [PATCH 35/57] Remove unnecessary variables --- .../semidiscretization_parabolic_auxvars.jl | 32 +++++++++---------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index c87db31be31..61f7dcba1f1 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -53,20 +53,18 @@ function SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver_auxvars = create_solver_auxvars(solver) - equations_grad_x = GradientEquations2D(nvariables(equations), 1) - semi_grad_x = SemidiscretizationHyperbolic(mesh, equations_grad_x, - initial_condition_constant, solver_auxvars) + equations_gradients_x = GradientEquations2D(nvariables(equations), 1) + semi_gradients_x = SemidiscretizationHyperbolic(mesh, equations_gradients_x, + initial_condition_constant, solver_auxvars) - equations_grad_y = GradientEquations2D(nvariables(equations), 2) - semi_grad_y = SemidiscretizationHyperbolic(mesh, equations_grad_y, - initial_condition_constant, solver_auxvars) + equations_gradients_y = GradientEquations2D(nvariables(equations), 2) + semi_gradients_y = SemidiscretizationHyperbolic(mesh, equations_gradients_y, + initial_condition_constant, solver_auxvars) - u_ode_grad_x = compute_coefficients(nothing, semi_grad_x) - u_ode_grad_y = compute_coefficients(nothing, semi_grad_y) - u_ode_grad_xx = compute_coefficients(nothing, semi_grad_x) - u_ode_grad_yy = compute_coefficients(nothing, semi_grad_y) + u_ode_gradients_x = compute_coefficients(nothing, semi_gradients_x) + u_ode_gradients_y = compute_coefficients(nothing, semi_gradients_y) - cache = (; cache..., semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y, u_ode_grad_xx, u_ode_grad_yy) + cache = (; cache..., semi_gradients_x, u_ode_gradients_x, semi_gradients_y, u_ode_gradients_y) SemidiscretizationParabolicAuxVars{typeof(mesh), typeof(equations), typeof(initial_condition), typeof(_boundary_conditions), typeof(source_terms), typeof(solver), typeof(cache)}( mesh, equations, initial_condition, _boundary_conditions, source_terms, solver, cache) @@ -167,12 +165,12 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationParabolicAuxVars, t) du = wrap_array(du_ode, mesh, equations, solver, cache) # Compute gradients - @unpack semi_grad_x, u_ode_grad_x, semi_grad_y, u_ode_grad_y = cache - rhs!(u_ode_grad_x, u_ode, semi_grad_x, t) - rhs!(u_ode_grad_y, u_ode, semi_grad_y, t) - gradients = (wrap_array(u_ode_grad_x, mesh, equations, solver, cache), - wrap_array(u_ode_grad_y, mesh, equations, solver, cache)) - cache_gradients = (semi_grad_x.cache, semi_grad_y.cache) + @unpack semi_gradients_x, u_ode_gradients_x, semi_gradients_y, u_ode_gradients_y = cache + rhs!(u_ode_gradients_x, u_ode, semi_gradients_x, t) + rhs!(u_ode_gradients_y, u_ode, semi_gradients_y, t) + gradients = (wrap_array(u_ode_gradients_x, mesh, equations, solver, cache), + wrap_array(u_ode_gradients_y, mesh, equations, solver, cache)) + cache_gradients = (semi_gradients_x.cache, semi_gradients_y.cache) # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() From 291884a457f06a7d29213ff538c2aad0a1cdafb1 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 06:40:11 +0100 Subject: [PATCH 36/57] Fix bug in mortar code \o/ --- src/solvers/dg/dg_2d_parabolic.jl | 2 +- test/test_examples_2d_advdiff.jl | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/solvers/dg/dg_2d_parabolic.jl b/src/solvers/dg/dg_2d_parabolic.jl index ebe285bdb60..f6a093e47e4 100644 --- a/src/solvers/dg/dg_2d_parabolic.jl +++ b/src/solvers/dg/dg_2d_parabolic.jl @@ -179,7 +179,7 @@ end gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_interfaces[2], equations, dg, i, interface) gradients_ll = (gradients_x_ll, gradients_y_ll) gradients_rr = (gradients_x_rr, gradients_y_rr) - flux = surface_flux(u_ll, u_rr, gradients_x_ll, gradients_y_ll, orientation, equations) + flux = surface_flux(u_ll, u_rr, gradients_ll, gradients_rr, orientation, equations) # Copy flux to left and right element storage set_node_vars!(destination, flux, equations, dg, i) diff --git a/test/test_examples_2d_advdiff.jl b/test/test_examples_2d_advdiff.jl index 6710b7606d1..df6bd8de548 100644 --- a/test/test_examples_2d_advdiff.jl +++ b/test/test_examples_2d_advdiff.jl @@ -11,12 +11,8 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") @testset "Scalar advection-diffusion" begin @testset "elixir_advdiff_basic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_basic.jl"), - # FIXME: These values changed from the naive implementation for the parabolic terms to the - # "proper" implementaton by O(1e-6): Find out why and if it is OK - # l2 = [6.707801888154751e-5], - # linf = [0.0005912932758251888]) - l2 = [6.447564451432685e-5], - linf = [0.0005850692252580281]) + l2 = [6.707801888154751e-5], + linf = [0.0005912932758251888]) end end From 0b5984c178ae1f66649ac5b6d11f36aa7efad8d2 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 06:54:35 +0100 Subject: [PATCH 37/57] Make AMR work --- examples/2d/elixir_advdiff_amr.jl | 72 +++++++++++++++++++ ...semidiscretization_hyperbolic_parabolic.jl | 27 +++++++ .../semidiscretization_parabolic_auxvars.jl | 23 ++++++ test/test_examples_2d_advdiff.jl | 6 ++ 4 files changed, 128 insertions(+) create mode 100644 examples/2d/elixir_advdiff_amr.jl diff --git a/examples/2d/elixir_advdiff_amr.jl b/examples/2d/elixir_advdiff_amr.jl new file mode 100644 index 00000000000..ee56b13e641 --- /dev/null +++ b/examples/2d/elixir_advdiff_amr.jl @@ -0,0 +1,72 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advectionvelocity = (1.1, 0.6) +nu = 1.2e-2 +equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +polydeg = 3 +solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs) + +coordinates_min = (-5, -5) +coordinates_max = ( 5, 5) + +# Create a uniformely refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +initial_condition = initial_condition_gauss +semi = SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver_hyperbolic) + + +############################################################################### +# 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 SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), + base_level=4, + med_level=5, med_threshold=0.1, + max_level=6, max_threshold=0.6) +amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + +# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step +# stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index e1930fb74d8..1aaf9c131fb 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -152,3 +152,30 @@ function SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condi return SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic) end + + +@inline function (amr_callback::AMRCallback)(u_ode::AbstractVector, + semi::SemidiscretizationHyperbolicParabolic, + t, iter; kwargs...) + # Get semidiscretizations + @unpack semi_hyperbolic, semi_parabolic = semi + @unpack semi_gradients_x, semi_gradients_y = semi_parabolic.cache + + # Perform AMR based on hyperbolic solver + has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi_hyperbolic)..., t, iter; + kwargs...) + + if has_changed + # Reinitialize data structures for all secondary solvers (no solution data conversion needed) + reinitialize_containers!(mesh_equations_solver_cache(semi_parabolic)...) + reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_x)...) + reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_y)...) + + # Resize storage for auxiliary variables (= gradients) + @unpack u_ode_gradients_x, u_ode_gradients_y = semi_parabolic.cache + resize!(u_ode_gradients_x, length(u_ode)) + resize!(u_ode_gradients_y, length(u_ode)) + end + + return has_changed +end diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index 61f7dcba1f1..be18b2ecbd1 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -180,3 +180,26 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationParabolicAuxVars, t) return nothing end + + +@inline function (amr_callback::AMRCallback)(u_ode::AbstractVector, + semi::SemidiscretizationParabolicAuxVars, + t, iter; + kwargs...) + # Perform AMR based on main solver + has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi)..., t, iter; kwargs...) + + if has_changed + # Reinitialize data structures for all secondary solvers (no solution data conversion needed) + @unpack semi_gradients_x, semi_gradients_y = semi.cache + reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_x)...) + reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_y)...) + + # Resize storage for auxiliary variables (= gradients) + @unpack u_ode_gradients_x, u_ode_gradients_y = semi.cache + resize!(u_ode_gradients_x, length(u_ode)) + resize!(u_ode_gradients_y, length(u_ode)) + end + + return has_changed +end diff --git a/test/test_examples_2d_advdiff.jl b/test/test_examples_2d_advdiff.jl index df6bd8de548..18460b0faf6 100644 --- a/test/test_examples_2d_advdiff.jl +++ b/test/test_examples_2d_advdiff.jl @@ -14,6 +14,12 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") l2 = [6.707801888154751e-5], linf = [0.0005912932758251888]) end + + @testset "elixir_advdiff_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_amr.jl"), + l2 = [0.1282304080765362], + linf = [0.8335673252435933]) + end end end # module From 3e545da547d47b857a22d82fbcc6fd61ff59c06d Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 07:22:37 +0100 Subject: [PATCH 38/57] Allow changing plot properties for mesh plots --- src/visualization/plot_recipes.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/visualization/plot_recipes.jl b/src/visualization/plot_recipes.jl index debcc56d174..8fbf2c0802f 100644 --- a/src/visualization/plot_recipes.jl +++ b/src/visualization/plot_recipes.jl @@ -239,9 +239,9 @@ end legend --> :none # Set series properties - seriestype := :path - linecolor := :black - linewidth := 1 + seriestype --> :path + linecolor --> :black + linewidth --> 1 # Return data for plotting mesh_vertices_x, mesh_vertices_y From 8486e60db63c1aa9a67e3836566dc9cde0620ea3 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 07:24:12 +0100 Subject: [PATCH 39/57] Make AMR test more fun... --- examples/2d/elixir_advdiff_amr.jl | 16 ++++++++-------- test/test_examples_2d_advdiff.jl | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/2d/elixir_advdiff_amr.jl b/examples/2d/elixir_advdiff_amr.jl index ee56b13e641..d33049d203f 100644 --- a/examples/2d/elixir_advdiff_amr.jl +++ b/examples/2d/elixir_advdiff_amr.jl @@ -5,20 +5,20 @@ using Trixi ############################################################################### # semidiscretization of the linear advection equation -advectionvelocity = (1.1, 0.6) -nu = 1.2e-2 +advectionvelocity = (1.8, 0.6) +nu = 2.4e-1 equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux polydeg = 3 solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs) -coordinates_min = (-5, -5) -coordinates_max = ( 5, 5) +coordinates_min = (-4, -4) +coordinates_max = ( 4, 4) # Create a uniformely refined mesh with periodic boundaries mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=4, + initial_refinement_level=3, n_cells_max=30_000) # set maximum capacity of tree data structure # A semidiscretization collects data structures and functions for the spatial discretization @@ -45,9 +45,9 @@ save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim) amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), - base_level=4, - med_level=5, med_threshold=0.1, - max_level=6, max_threshold=0.6) + base_level=3, + med_level=4, med_threshold=0.1, + max_level=5, max_threshold=0.5) amr_callback = AMRCallback(semi, amr_controller, interval=5, adapt_initial_condition=true, diff --git a/test/test_examples_2d_advdiff.jl b/test/test_examples_2d_advdiff.jl index 18460b0faf6..7b9f89a09bf 100644 --- a/test/test_examples_2d_advdiff.jl +++ b/test/test_examples_2d_advdiff.jl @@ -17,8 +17,8 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") @testset "elixir_advdiff_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_amr.jl"), - l2 = [0.1282304080765362], - linf = [0.8335673252435933]) + l2 = [0.16503675954115254], + linf = [0.9236968722827653]) end end From 6af4cd42aa87e1b7a3a7575f347ca38b41b1307d Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 07:30:29 +0100 Subject: [PATCH 40/57] Add test for heat equation AMR --- examples/2d/elixir_heat_amr.jl | 71 ++++++++++++++++++++++++++++++++++ test/test_examples_2d_heat.jl | 6 +++ 2 files changed, 77 insertions(+) create mode 100644 examples/2d/elixir_heat_amr.jl diff --git a/examples/2d/elixir_heat_amr.jl b/examples/2d/elixir_heat_amr.jl new file mode 100644 index 00000000000..1ea24343cfe --- /dev/null +++ b/examples/2d/elixir_heat_amr.jl @@ -0,0 +1,71 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +nu = 4.8e-1 +equations = HeatEquation2D(nu) + +# Create DG solver with polynomial degree = 3 and the central flux as surface flux +solver = DGSEM(3, flux_central) + +coordinates_min = (-4, -4) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 4, 4) # maximum coordinates (max(x), max(y)) + +# Create a uniformely refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +initial_condition = initial_condition_gauss +semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver) + + +############################################################################### +# 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 SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), + base_level=3, + med_level=4, med_threshold=0.05, + max_level=5, max_threshold=0.3) +amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + +# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step +# stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks, maxiters=1e5); + +# Print the timer summary +summary_callback() diff --git a/test/test_examples_2d_heat.jl b/test/test_examples_2d_heat.jl index 76de3ed8342..4300a5d03af 100644 --- a/test/test_examples_2d_heat.jl +++ b/test/test_examples_2d_heat.jl @@ -14,6 +14,12 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") l2 = [0.0008071867035363602], linf = [0.004961971226386641]) end + + @testset "elixir_heat_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_amr.jl"), + l2 = [0.08891001607226263], + linf = [0.657532205699553]) + end end end # module From 2d737b734bb02f8f25a035bc62ae1ee9dcef8fd0 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 07:42:54 +0100 Subject: [PATCH 41/57] Update include order --- src/Trixi.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 7a088b6016e..5468b3388f2 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -56,11 +56,11 @@ include("mesh/mesh.jl") include("solvers/dg/dg.jl") include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") -include("semidiscretization/semidiscretization_parabolic_auxvars.jl") -include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") +include("semidiscretization/semidiscretization_parabolic_auxvars.jl") +include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("time_integration/time_integration.jl") # `trixi_include` and special elixirs such as `convergence_test` From a0a8f21e0476442b6b1addd5c3a174f5e50861a6 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 09:15:32 +0100 Subject: [PATCH 42/57] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- .../semidiscretization_hyperbolic_parabolic.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 1aaf9c131fb..55785f224a3 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -68,20 +68,19 @@ end @inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolicParabolic) - @unpack mesh, equations, solver, cache = semi.semi_hyperbolic - return mesh, equations, solver, cache + return mesh_equations_solver_cache(semi.semi_hyperbolic) end function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis) - @unpack mesh, equations, initial_condition, solver, cache = semi.semi_hyperbolic + mesh, equations, initial_condition, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic) u = wrap_array(u_ode, mesh, equations, solver, cache) calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis) - @unpack mesh, equations, initial_condition, solver, cache = semi.semi_hyperbolic + mesh, equations, initial_condition, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic) calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end From 4cc23f435849d0d12d0c9d3360fca00eeb416a1a Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 09:16:18 +0100 Subject: [PATCH 43/57] Fix docstring for heat equation --- src/equations/heat_equation_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index b3605fb69b5..1dfca6f0fa3 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -2,11 +2,11 @@ @doc raw""" HeatEquation2D -The linear scalar advection equation +The heat equation ```math -\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u = 0 +\partial_t u - \nu (\partial^2_1 u + \partial^2_2 u) = 0 ``` -in two space dimensions with constant velocity `a`. +in two space dimensions with constant viscosity ``\nu``. """ struct HeatEquation2D{RealT<:Real} <: AbstractHeatEquation{2, 1} nu::RealT From 7716534b9869ab59e1b8f46e59379f560f4f6106 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 09:17:21 +0100 Subject: [PATCH 44/57] Remove unused function --- src/equations/heat_equation_2d.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index 1dfca6f0fa3..b0b86a8e681 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -68,10 +68,6 @@ end @inline have_constant_speed(::HeatEquation2D) = Val(true) -@inline function max_abs_speeds(equation::HeatEquation2D) - return NaN -end - # Convert conservative variables to primitive @inline cons2prim(u, equation::HeatEquation2D) = u From 379c0599fe44b96b44672d6fbc4cadf9b6673a57 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 09:19:15 +0100 Subject: [PATCH 45/57] Fix docstring for advection-diffusion equation --- src/equations/linear_advection_diffusion_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index 417611adc52..f3bc4ee4595 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -2,11 +2,11 @@ @doc raw""" LinearAdvectionDiffusionEquation2D -The linear scalar advection equation +The linear advection-diffusion equation ```math -\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u = 0 +\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u - \nu (\partial^2_1 u + \partial^2_2 u)= 0 ``` -in two space dimensions with constant velocity `a`. +in two space dimensions with constant velocity `a` and constant viscosity ``\nu``. """ struct LinearAdvectionDiffusionEquation2D{RealT<:Real} <: AbstractLinearAdvectionDiffusionEquation{2, 1} advectionvelocity::SVector{2, RealT} From 79129395e477476162c77c2dce0684e43bde4614 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 10:11:22 +0100 Subject: [PATCH 46/57] Add comment as suggested by @andrewwinters5000 --- src/solvers/dg/dg_2d_parabolic.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solvers/dg/dg_2d_parabolic.jl b/src/solvers/dg/dg_2d_parabolic.jl index f6a093e47e4..12f4cf0974d 100644 --- a/src/solvers/dg/dg_2d_parabolic.jl +++ b/src/solvers/dg/dg_2d_parabolic.jl @@ -11,7 +11,7 @@ function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, gradients, have_nonconservative_terms(equations), equations, dg.volume_integral, dg, cache) - # Prolong solution to interfaces + # Prolong solution and gradients to interfaces @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(cache, cache_gradients, u, gradients, equations, dg) # Calculate interface fluxes @@ -19,13 +19,13 @@ function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, have_nonconservative_terms(equations), equations, dg, cache, cache_gradients) - # Prolong solution to boundaries + # Prolong solution and gradients to boundaries @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, u, equations, dg) # Calculate boundary fluxes @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, equations, dg) - # Prolong solution to mortars + # Prolong solution and gradients to mortars @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, cache_gradients, u, gradients, equations, dg.mortar, dg) From 50242e0b7aac004bebbad691b9222efad7eb3d77 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 10:21:45 +0100 Subject: [PATCH 47/57] Fix bug with mesh_equations_solver_cache use --- .../semidiscretization_hyperbolic_parabolic.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 55785f224a3..948122f75e6 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -73,14 +73,17 @@ end function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis) - mesh, equations, initial_condition, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic) + @unpack initial_condition = semi.semi_hyperbolic + u = wrap_array(u_ode, mesh, equations, solver, cache) calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis) - mesh, equations, initial_condition, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic) + @unpack initial_condition = semi.semi_hyperbolic calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis) end From 665e99f0885d222e70a137e2c90c4545c1073c95 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 19 Jan 2021 10:22:05 +0100 Subject: [PATCH 48/57] Make weak form the default for the BR1 scheme --- .../semidiscretization_hyperbolic_parabolic.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 948122f75e6..726b308e6c1 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -139,13 +139,14 @@ hyperbolic and a parabolic semidiscretization are created and passed to a function SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver::DGSEM; source_terms=nothing, boundary_conditions=boundary_condition_periodic, - RealT=real(solver)) + RealT=real(solver), + volume_integral_parabolic=VolumeIntegralWeakForm()) semi_hyperbolic = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; source_terms=source_terms, boundary_conditions=boundary_conditions, RealT=RealT) - solver_parabolic = DGSEM(solver.basis, flux_central, solver.volume_integral, solver.mortar) + solver_parabolic = DGSEM(solver.basis, flux_central, volume_integral_parabolic, solver.mortar) semi_parabolic = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver_parabolic, source_terms=source_terms, From 5769eae95d62e6fd97cee305dd5785b8b877d27c Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 2 Feb 2021 16:33:08 +0100 Subject: [PATCH 49/57] Add very early support for boundary conditions for parabolic problems (convergence test seems OK) --- examples/2d/elixir_heat_nonperiodic.jl | 70 ++++++++++++ src/Trixi.jl | 1 + src/equations/gradient_equations_2d.jl | 19 ++++ src/equations/heat_equation_2d.jl | 18 ++++ .../semidiscretization_parabolic_auxvars.jl | 6 +- src/solvers/dg/dg_2d_parabolic.jl | 102 +++++++++++++++++- test/test_examples_2d_heat.jl | 6 ++ 7 files changed, 218 insertions(+), 4 deletions(-) create mode 100644 examples/2d/elixir_heat_nonperiodic.jl diff --git a/examples/2d/elixir_heat_nonperiodic.jl b/examples/2d/elixir_heat_nonperiodic.jl new file mode 100644 index 00000000000..c66d9c684de --- /dev/null +++ b/examples/2d/elixir_heat_nonperiodic.jl @@ -0,0 +1,70 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +nu = 1.2e-2 +equations = HeatEquation2D(nu) + +# Create DG solver with polynomial degree = 3 and the central flux as surface flux +solver = DGSEM(3, flux_central) + +coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1, 1) # maximum coordinates (max(x), max(y)) + +# Create a uniformely refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=2, + n_cells_max=30_000, # set maximum capacity of tree data structure + periodicity=(false, true)) + +# A semidiscretization collects data structures and functions for the spatial discretization +initial_condition = initial_condition_sin_x + +# 1 => -x, 2 => +x, 3 => -y, 4 => +y as usual for orientations +boundary_conditions = (x_neg=boundary_condition_sin_x, + x_pos=boundary_condition_sin_x, + y_neg=boundary_condition_periodic, + y_pos=boundary_condition_periodic) +semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver, + 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 SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step +# stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks, maxiters=1e5); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index 5468b3388f2..b64e4241c27 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -113,6 +113,7 @@ export initial_condition_sedov_self_gravity, boundary_condition_sedov_self_gravi export initial_condition_eoc_test_coupled_euler_gravity, source_terms_eoc_test_coupled_euler_gravity, source_terms_eoc_test_euler export initial_condition_lid_driven_cavity, boundary_condition_lid_driven_cavity export initial_condition_couette_steady, initial_condition_couette_unsteady, boundary_condition_couette +export initial_condition_sin_x, boundary_condition_sin_x export cons2cons, cons2prim, cons2macroscopic export density, pressure, density_pressure, velocity diff --git a/src/equations/gradient_equations_2d.jl b/src/equations/gradient_equations_2d.jl index c227985d336..33bace780a6 100644 --- a/src/equations/gradient_equations_2d.jl +++ b/src/equations/gradient_equations_2d.jl @@ -32,6 +32,25 @@ function initial_condition_constant(x, t, equations::GradientEquations2D) end +function boundary_condition_sin_x(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::GradientEquations2D) + if orientation == equations.orientation + # u_boundary = initial_condition_sin_x(x, t, equation) + nu = 1.2e-2 + omega = pi + scalar = sin(omega * x[1]) * exp(-nu * omega^2 * t) + u_boundary = SVector(scalar) + + return -u_boundary + else + return SVector(ntuple(v -> zero(eltype(u_inner)), nvariables(equations))) + end + + return flux +end + + # Pre-defined source terms should be implemented as # function source_terms_WHATEVER(u, x, t, equations::GradientEquations2D) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index b0b86a8e681..bed41f9751c 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -56,6 +56,24 @@ function initial_condition_gauss(x, t, equation::HeatEquation2D) end +function initial_condition_sin_x(x, t, equation::HeatEquation2D) + @unpack nu = equation + + omega = pi + scalar = sin(omega * x[1]) * exp(-nu * omega^2 * t) + + return @SVector [scalar] +end + + +function boundary_condition_sin_x(u_inner, gradients_inner, orientation, direction, x, t, + surface_flux_function, + equation::HeatEquation2D) + # OBS! This is specifically implemented for BR1 + return calcflux(u_inner, gradients_inner, orientation, equation) +end + + # Pre-defined source terms should be implemented as # function source_terms_WHATEVER(u, x, t, equations::HeatEquation2D) diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index be18b2ecbd1..7499b7204c3 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -55,11 +55,13 @@ function SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, equations_gradients_x = GradientEquations2D(nvariables(equations), 1) semi_gradients_x = SemidiscretizationHyperbolic(mesh, equations_gradients_x, - initial_condition_constant, solver_auxvars) + initial_condition_constant, solver_auxvars, + boundary_conditions=boundary_conditions) equations_gradients_y = GradientEquations2D(nvariables(equations), 2) semi_gradients_y = SemidiscretizationHyperbolic(mesh, equations_gradients_y, - initial_condition_constant, solver_auxvars) + initial_condition_constant, solver_auxvars, + boundary_conditions=boundary_conditions) u_ode_gradients_x = compute_coefficients(nothing, semi_gradients_x) u_ode_gradients_y = compute_coefficients(nothing, semi_gradients_y) diff --git a/src/solvers/dg/dg_2d_parabolic.jl b/src/solvers/dg/dg_2d_parabolic.jl index 12f4cf0974d..9dbf0b6f5d4 100644 --- a/src/solvers/dg/dg_2d_parabolic.jl +++ b/src/solvers/dg/dg_2d_parabolic.jl @@ -20,10 +20,10 @@ function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t, dg, cache, cache_gradients) # Prolong solution and gradients to boundaries - @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, u, equations, dg) + @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, cache_gradients, u, gradients, equations, dg) # Calculate boundary fluxes - @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, equations, dg) + @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, cache_gradients, t, boundary_conditions, equations, dg) # Prolong solution and gradients to mortars @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, cache_gradients, u, gradients, @@ -128,6 +128,104 @@ function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4}, end +function prolong2boundaries!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, + equations, dg::DG) + prolong2boundaries!(cache, u, equations, dg) + prolong2boundaries!(cache_gradients[1], gradients[1], equations, dg) + prolong2boundaries!(cache_gradients[2], gradients[2], equations, dg) + + return nothing +end + + +# TODO: Taal dimension agnostic +function calc_boundary_flux!(cache, cache_gradients, t, + boundary_condition::BoundaryConditionPeriodic, + equations::AbstractEquations{2}, dg::DG) + @assert isempty(eachboundary(dg, cache)) +end + +# TODO: Taal dimension agnostic +function calc_boundary_flux!(cache, cache_gradients, t, boundary_condition, + equations::AbstractEquations{2}, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + for direction in eachindex(firsts) + calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_condition, + equations, dg, cache, cache_gradients, + direction, firsts[direction], lasts[direction]) + end +end + +function calc_boundary_flux!(cache, cache_gradients, t, + boundary_conditions::Union{NamedTuple,Tuple}, + equations::AbstractEquations{2}, dg::DG) + @unpack surface_flux_values = cache.elements + @unpack n_boundaries_per_direction = cache.boundaries + + # Calculate indices + lasts = accumulate(+, n_boundaries_per_direction) + firsts = lasts - n_boundaries_per_direction .+ 1 + + # Calc boundary fluxes in each direction + calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1], + equations, dg, cache, cache_gradients, 1, firsts[1], lasts[1]) + calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2], + equations, dg, cache, cache_gradients, 2, firsts[2], lasts[2]) + calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[3], + equations, dg, cache, cache_gradients, 3, firsts[3], lasts[3]) + calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[4], + equations, dg, cache, cache_gradients, 4, firsts[4], lasts[4]) +end + +function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t, + boundary_condition, equations, dg::DG, + cache, cache_gradients, + direction, first_boundary, last_boundary) + @unpack surface_flux = dg + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + gradients_x = cache_gradients[1].boundaries.u + gradients_y = cache_gradients[2].boundaries.u + + Threads.@threads for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for i in eachnode(dg) + # Get boundary flux + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, boundary) + gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_x, equations, dg, i, boundary) + gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_y, equations, dg, i, boundary) + gradients_ll = (gradients_x_ll, gradients_y_ll) + gradients_rr = (gradients_x_rr, gradients_y_rr) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + u_inner = u_ll + gradients_inner = gradients_ll + else # Element is on the right, boundary on the left + u_inner = u_rr + gradients_inner = gradients_rr + end + x = get_node_coords(node_coordinates, equations, dg, i, boundary) + flux = boundary_condition(u_inner, gradients_inner, orientations[boundary], direction, x, t, + surface_flux, equations) + + # Copy flux to left and right element storage + for v in eachvariable(equations) + surface_flux_values[v, i, direction, neighbor] = flux[v] + end + end + end + + return nothing +end + + function prolong2mortars!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients, equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) prolong2mortars!(cache, u, equations, mortar_l2, dg) diff --git a/test/test_examples_2d_heat.jl b/test/test_examples_2d_heat.jl index 4300a5d03af..295adcec1d3 100644 --- a/test/test_examples_2d_heat.jl +++ b/test/test_examples_2d_heat.jl @@ -20,6 +20,12 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") l2 = [0.08891001607226263], linf = [0.657532205699553]) end + + @testset "elixir_heat_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_nonperiodic.jl"), + l2 = [0.0012846393772784545], + linf = [0.005583936479563221]) + end end end # module From 8991ea4fd32bc03864568b3ff9ed75e79181d121 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 2 Feb 2021 17:01:48 +0100 Subject: [PATCH 50/57] Add source terms for parabolic system (convtest seems ok) --- examples/2d/elixir_heat_source_terms.jl | 63 +++++++++++++++++++++++++ src/equations/heat_equation_2d.jl | 24 +++++++++- test/test_examples_2d_heat.jl | 6 +++ 3 files changed, 91 insertions(+), 2 deletions(-) create mode 100644 examples/2d/elixir_heat_source_terms.jl diff --git a/examples/2d/elixir_heat_source_terms.jl b/examples/2d/elixir_heat_source_terms.jl new file mode 100644 index 00000000000..3aedfcee641 --- /dev/null +++ b/examples/2d/elixir_heat_source_terms.jl @@ -0,0 +1,63 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +nu = 1.2e-2 +equations = HeatEquation2D(nu) + +# Create DG solver with polynomial degree = 3 and the central flux as surface flux +solver = DGSEM(3, flux_central) + +coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1, 1) # maximum coordinates (max(x), max(y)) + +# Create a uniformely refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=2, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# A semidiscretization collects data structures and functions for the spatial discretization +initial_condition = initial_condition_poisson_periodic +semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver, + source_terms=source_terms_poisson_periodic) + + +############################################################################### +# 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 SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step +# stepsize_callback = StepsizeCallback(cfl=1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks, maxiters=1e5); + +# Print the timer summary +summary_callback() diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index bed41f9751c..746011f7067 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -74,8 +74,28 @@ function boundary_condition_sin_x(u_inner, gradients_inner, orientation, directi end -# Pre-defined source terms should be implemented as -# function source_terms_WHATEVER(u, x, t, equations::HeatEquation2D) +function initial_condition_poisson_periodic(x, t, equation::HeatEquation2D) + # elliptic equation: -νΔϕ = f + # depending on initial constant state, c, for phi this converges to the solution ϕ + c + @unpack nu = equation + + phi = sin(2.0*pi*x[1])*sin(2.0*pi*x[2])*(1 - exp(-8*nu*pi^2*t)) + return @SVector [phi] +end + + +@inline function source_terms_poisson_periodic(u, x, t, equation::HeatEquation2D) + # elliptic equation: -νΔϕ = f + # analytical solution: phi = sin(2πx)*sin(2πy) and f = -8νπ^2 sin(2πx)*sin(2πy) + C = -8 * equation.nu * pi^2 + + x1, x2 = x + tmp1 = sinpi(2 * x1) + tmp2 = sinpi(2 * x2) + du1 = -C*tmp1*tmp2 + + return SVector(du1) +end # Calculate parabolic 1D flux in axis `orientation` for a single point diff --git a/test/test_examples_2d_heat.jl b/test/test_examples_2d_heat.jl index 295adcec1d3..f75159d1c73 100644 --- a/test/test_examples_2d_heat.jl +++ b/test/test_examples_2d_heat.jl @@ -26,6 +26,12 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") l2 = [0.0012846393772784545], linf = [0.005583936479563221]) end + + @testset "elixir_heat_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_source_terms.jl"), + l2 = [0.014470380266085592], + linf = [0.04433923071172761]) + end end end # module From 8645e7d26769c86ce989d8f961a97cf1bdfbe06c Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 2 Feb 2021 18:04:42 +0100 Subject: [PATCH 51/57] Implement stepsize calculation for parabolic systems in the StepsizeCallback --- examples/2d/elixir_heat_amr.jl | 8 ++--- examples/2d/elixir_heat_basic.jl | 11 +++--- examples/2d/elixir_heat_nonperiodic.jl | 11 +++--- examples/2d/elixir_heat_source_terms.jl | 11 +++--- src/callbacks_step/stepsize.jl | 19 +++++----- src/callbacks_step/stepsize_dg1d.jl | 8 ++--- src/callbacks_step/stepsize_dg2d.jl | 36 +++++++++++++------ src/callbacks_step/stepsize_dg3d.jl | 8 ++--- src/equations/heat_equation_2d.jl | 7 +++- .../semidiscretization_euler_gravity.jl | 20 +++++++++-- .../semidiscretization_hyperbolic.jl | 12 +++++++ .../semidiscretization_parabolic_auxvars.jl | 11 ++++++ test/test_examples_2d_heat.jl | 16 ++++----- 13 files changed, 115 insertions(+), 63 deletions(-) diff --git a/examples/2d/elixir_heat_amr.jl b/examples/2d/elixir_heat_amr.jl index 1ea24343cfe..00981bafbb3 100644 --- a/examples/2d/elixir_heat_amr.jl +++ b/examples/2d/elixir_heat_amr.jl @@ -52,11 +52,11 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition_only_refine=true) # The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl=1.6) +stepsize_callback = StepsizeCallback(cfl=0.5) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback, + stepsize_callback) ############################################################################### @@ -64,7 +64,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_ # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks, maxiters=1e5); # Print the timer summary diff --git a/examples/2d/elixir_heat_basic.jl b/examples/2d/elixir_heat_basic.jl index f878032c295..18d0be4cf01 100644 --- a/examples/2d/elixir_heat_basic.jl +++ b/examples/2d/elixir_heat_basic.jl @@ -27,8 +27,8 @@ semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, so ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.0 -tspan = (0.0, 1.0) +# Create ODE problem with time span from 0.0 to 50.0 +tspan = (0.0, 50.0) ode = semidiscretize(semi, tspan); # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup @@ -43,11 +43,10 @@ save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim) # The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl=1.6) +stepsize_callback = StepsizeCallback(cfl=0.5) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) ############################################################################### @@ -55,7 +54,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks, maxiters=1e5); # Print the timer summary diff --git a/examples/2d/elixir_heat_nonperiodic.jl b/examples/2d/elixir_heat_nonperiodic.jl index c66d9c684de..d6c63b90bca 100644 --- a/examples/2d/elixir_heat_nonperiodic.jl +++ b/examples/2d/elixir_heat_nonperiodic.jl @@ -35,8 +35,8 @@ semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, so ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.0 -tspan = (0.0, 1.0) +# Create ODE problem with time span from 0.0 to 100.0 +tspan = (0.0, 100.0) ode = semidiscretize(semi, tspan); # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup @@ -51,11 +51,10 @@ save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim) # The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl=1.6) +stepsize_callback = StepsizeCallback(cfl=0.5) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) ############################################################################### @@ -63,7 +62,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks, maxiters=1e5); # Print the timer summary diff --git a/examples/2d/elixir_heat_source_terms.jl b/examples/2d/elixir_heat_source_terms.jl index 3aedfcee641..e3b3a87944b 100644 --- a/examples/2d/elixir_heat_source_terms.jl +++ b/examples/2d/elixir_heat_source_terms.jl @@ -28,8 +28,8 @@ semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, so ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to 1.0 -tspan = (0.0, 1.0) +# Create ODE problem with time span from 0.0 to 100.0 +tspan = (0.0, 100.0) ode = semidiscretize(semi, tspan); # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup @@ -44,11 +44,10 @@ save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim) # The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl=1.6) +stepsize_callback = StepsizeCallback(cfl=0.5) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) ############################################################################### @@ -56,7 +55,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks, maxiters=1e5); # Print the timer summary diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 5a836bbb9c7..9bcbe1baa65 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -1,11 +1,13 @@ """ - StepsizeCallback(; cfl=1.0) + StepsizeCallback(; cfl) Set the time step size according to a CFL condition with CFL number `cfl` -if the time integration method isn't adaptive itself. +if the time integration method isn't adaptive itself. When using a semidiscretization that wraps +multiple solvers, you might need to provide a tuple of CFL numbers. """ mutable struct StepsizeCallback{RealT} + # FIXME: RealT -> CFLNumber cfl_number::RealT end @@ -30,7 +32,7 @@ function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{Condition,Af end -function StepsizeCallback(; cfl::Real=1.0) +function StepsizeCallback(; cfl) stepsize_callback = StepsizeCallback(cfl) @@ -58,13 +60,10 @@ end t = integrator.t u_ode = integrator.u semi = integrator.p - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @unpack cfl_number = stepsize_callback - u = wrap_array(u_ode, mesh, equations, solver, cache) - dt = @timeit_debug timer() "calculate dt" cfl_number * max_dt(u, t, mesh, - have_constant_speed(equations), equations, - solver, cache) + dt = @timeit_debug timer() "calculate dt" max_dt(u_ode, t, cfl_number, semi) + set_proposed_dt!(integrator, dt) integrator.opts.dtmax = dt integrator.dtcache = dt @@ -86,10 +85,8 @@ function (cb::DiscreteCallback{Condition,Affect!})(ode::ODEProblem) where {Condi u_ode = ode.u0 t = first(ode.tspan) semi = ode.p - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - u = wrap_array(u_ode, mesh, equations, solver, cache) - return cfl_number * max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) + return max_dt(u_ode, t, cfl_number, semi) end diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl index 65973311bad..1940451b8d7 100644 --- a/src/callbacks_step/stepsize_dg1d.jl +++ b/src/callbacks_step/stepsize_dg1d.jl @@ -1,6 +1,6 @@ -function max_dt(u::AbstractArray{<:Any,3}, t, mesh::TreeMesh{1}, - constant_speed::Val{false}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,3}, t, mesh::TreeMesh{1}, + constant_speed::Val{false}, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -20,8 +20,8 @@ function max_dt(u::AbstractArray{<:Any,3}, t, mesh::TreeMesh{1}, end -function max_dt(u::AbstractArray{<:Any,3}, t, mesh::TreeMesh{1}, - constant_speed::Val{true}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,3}, t, mesh::TreeMesh{1}, + constant_speed::Val{true}, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 56c453a88c0..2d275c425a1 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -1,6 +1,6 @@ -function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, - constant_speed::Val{false}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, + constant_speed::Val{false}, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -21,8 +21,8 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, end -function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, - constant_speed::Val{true}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, + constant_speed::Val{true}, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -37,13 +37,13 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, end -function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2}, - constant_speed::Val{false}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2}, + constant_speed::Val{false}, equations, dg::DG, cache) # call the method accepting a general `mesh::TreeMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. - dt = invoke(max_dt, + dt = invoke(max_dt_hyperbolic, Tuple{typeof(u), typeof(t), TreeMesh{2}, typeof(constant_speed), typeof(equations), typeof(dg), typeof(cache)}, u, t, mesh, constant_speed, equations, dg, cache) @@ -53,13 +53,13 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2}, end -function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2}, - constant_speed::Val{true}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2}, + constant_speed::Val{true}, equations, dg::DG, cache) # call the method accepting a general `mesh::TreeMesh{2}` # TODO: MPI, we should improve this; maybe we should dispatch on `u` # and create some MPI array type, overloading broadcasting and mapreduce etc. # Then, this specific array type should also work well with DiffEq etc. - dt = invoke(max_dt, + dt = invoke(max_dt_hyperbolic, Tuple{typeof(u), typeof(t), TreeMesh{2}, typeof(constant_speed), typeof(equations), typeof(dg), typeof(cache)}, u, t, mesh, constant_speed, equations, dg, cache) @@ -67,3 +67,19 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2}, return dt end + + +function max_dt_parabolic(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2}, + constant_speed::Val{true}, equations, dg::DG, cache) + # to avoid a division by zero if the diffusion vanishes everywhere, + # e.g. for steady-state linear advection + max_scaled_diffusion = nextfloat(zero(t)) + + for element in eachelement(dg, cache) + max_λ1, max_λ2 = max_abs_diffusions(equations) + inv_jacobian = cache.elements.inverse_jacobian[element] + max_scaled_diffusion = max(max_scaled_diffusion, inv_jacobian^2 * (max_λ1 + max_λ2)) + end + + return 2^2 / (nnodes(dg)^2 * max_scaled_diffusion) +end diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl index a038fabbfa2..60750eabab6 100644 --- a/src/callbacks_step/stepsize_dg3d.jl +++ b/src/callbacks_step/stepsize_dg3d.jl @@ -1,6 +1,6 @@ -function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3}, - constant_speed::Val{false}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3}, + constant_speed::Val{false}, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -22,8 +22,8 @@ function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3}, end -function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3}, - constant_speed::Val{true}, equations, dg::DG, cache) +function max_dt_hyperbolic(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3}, + constant_speed::Val{true}, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index 746011f7067..740ade89324 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -104,8 +104,13 @@ end end -@inline have_constant_speed(::HeatEquation2D) = Val(true) +@inline have_constant_diffusion(::HeatEquation2D) = Val(true) +# FIXME: Find a better name than `max_abs_diffusions` or `max_abs_speeds_viscous` +@inline function max_abs_diffusions(equation::HeatEquation2D) + @unpack nu = equation + return (nu, nu) +end # Convert conservative variables to primitive @inline cons2prim(u, equation::HeatEquation2D) = u diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index 5a1f6bff8bb..fda51e711cf 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -165,6 +165,17 @@ end end +function max_dt(u_ode::AbstractVector, t, cfl_number::Real, semi::SemidiscretizationEulerGravity) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + u = wrap_array(u_ode, mesh, equations, solver, cache) + + dt = cfl_number * max_dt_hyperbolic(u, t, mesh, have_constant_speed(equations), equations, solver, + cache) + + return dt +end + + function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t) @unpack semi_euler, semi_gravity, cache = semi @@ -229,9 +240,12 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode::AbstractVe # iterate gravity solver until convergence or maximum number of iterations are reached @unpack equations = semi_gravity while !finalstep - dt = @timeit_debug timer() "calculate dt" cfl * max_dt(u_gravity, t, semi_gravity.mesh, - have_constant_speed(equations), equations, - semi_gravity.solver, semi_gravity.cache) + dt = @timeit_debug timer() "calculate dt" cfl * max_dt_hyperbolic(u_gravity, t, + semi_gravity.mesh, + have_constant_speed(equations), + equations, + semi_gravity.solver, + semi_gravity.cache) # evolve solution by one pseudo-time step time_start = time_ns() diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 737f278e440..0a430cdd1c3 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -164,6 +164,18 @@ function compute_coefficients!(u_ode::AbstractVector, t, semi::Semidiscretizatio end +function max_dt(u_ode::AbstractVector, t, cfl_number::Real, semi::SemidiscretizationHyperbolic) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + u = wrap_array(u_ode, mesh, equations, solver, cache) + + # FIXME: add deprecation for max_dt -> max_dt_hyperbolic + dt = cfl_number * max_dt_hyperbolic(u, t, mesh, have_constant_speed(equations), equations, solver, + cache) + + return dt +end + + function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t) @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index 7499b7204c3..703e9b35bec 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -160,6 +160,17 @@ function compute_coefficients!(u_ode::AbstractVector, t, semi::Semidiscretizatio end +function max_dt(u_ode::AbstractVector, t, cfl_number::Real, semi::SemidiscretizationParabolicAuxVars) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + u = wrap_array(u_ode, mesh, equations, solver, cache) + + dt = cfl_number * max_dt_parabolic(u, t, mesh, have_constant_diffusion(equations), equations, + solver, cache) + + return dt +end + + function rhs!(du_ode, u_ode, semi::SemidiscretizationParabolicAuxVars, t) @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi diff --git a/test/test_examples_2d_heat.jl b/test/test_examples_2d_heat.jl index f75159d1c73..0d5d0388377 100644 --- a/test/test_examples_2d_heat.jl +++ b/test/test_examples_2d_heat.jl @@ -11,26 +11,26 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") @testset "Heat equation" begin @testset "elixir_heat_basic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_basic.jl"), - l2 = [0.0008071867035363602], - linf = [0.004961971226386641]) + l2 = [1.2700756869568872e-8], + linf = [6.629198634477973e-8]) end @testset "elixir_heat_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_amr.jl"), - l2 = [0.08891001607226263], - linf = [0.657532205699553]) + l2 = [0.08891001605654104], + linf = [0.6575322065283009]) end @testset "elixir_heat_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_nonperiodic.jl"), - l2 = [0.0012846393772784545], - linf = [0.005583936479563221]) + l2 = [2.050593060856897e-8], + linf = [6.631714299464696e-8]) end @testset "elixir_heat_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_source_terms.jl"), - l2 = [0.014470380266085592], - linf = [0.04433923071172761]) + l2 = [0.025350269628084208], + linf = [0.09030986815695541]) end end From 6231077ba9d5cf01f9eeb33801f81b3805a70342 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 2 Feb 2021 18:15:45 +0100 Subject: [PATCH 52/57] Implement stepsize calculation for hyperbolic-parabolic systems --- examples/2d/elixir_advdiff_amr.jl | 8 +++++--- examples/2d/elixir_advdiff_basic.jl | 8 ++++---- src/equations/linear_advection_diffusion_2d.jl | 7 +++++++ .../semidiscretization_hyperbolic_parabolic.jl | 10 ++++++++++ test/test_examples_2d_advdiff.jl | 8 ++++---- 5 files changed, 30 insertions(+), 11 deletions(-) diff --git a/examples/2d/elixir_advdiff_amr.jl b/examples/2d/elixir_advdiff_amr.jl index d33049d203f..ce0a4907e38 100644 --- a/examples/2d/elixir_advdiff_amr.jl +++ b/examples/2d/elixir_advdiff_amr.jl @@ -54,10 +54,12 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition_only_refine=true) # The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl=1.6) +# The first CFL number is for the hyperbolic system, the second for the parabolic system +stepsize_callback = StepsizeCallback(cfl=(1.6, 0.5)) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback, + stepsize_callback) ############################################################################### @@ -65,7 +67,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_ # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks); # Print the timer summary diff --git a/examples/2d/elixir_advdiff_basic.jl b/examples/2d/elixir_advdiff_basic.jl index ec1e6cbdcef..ac5b44f648c 100644 --- a/examples/2d/elixir_advdiff_basic.jl +++ b/examples/2d/elixir_advdiff_basic.jl @@ -49,11 +49,11 @@ save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim) # The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl=1.6) +# The first CFL number is for the hyperbolic system, the second for the parabolic system +stepsize_callback = StepsizeCallback(cfl=(1.0, 0.1)) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) ############################################################################### @@ -61,7 +61,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks); # Print the timer summary diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index f3bc4ee4595..2f2d05d2712 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -98,6 +98,13 @@ end return abs.(equation.advectionvelocity) end +@inline have_constant_diffusion(::LinearAdvectionDiffusionEquation2D) = Val(true) + +@inline function max_abs_diffusions(equation::LinearAdvectionDiffusionEquation2D) + @unpack nu = equation + return (nu, nu) +end + # Convert conservative variables to primitive @inline cons2prim(u, equation::LinearAdvectionDiffusionEquation2D) = u diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 726b308e6c1..19a3d98aa5d 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -98,6 +98,16 @@ function compute_coefficients!(u_ode::AbstractVector, t, semi::Semidiscretizatio end +function max_dt(u_ode::AbstractVector, t, cfl_number, + semi::SemidiscretizationHyperbolicParabolic) + @unpack semi_hyperbolic, semi_parabolic = semi + dt = min(max_dt(u_ode, t, cfl_number[1], semi_hyperbolic), + max_dt(u_ode, t, cfl_number[2], semi_parabolic)) + + return dt +end + + function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) @unpack du_ode_parabolic = semi.cache diff --git a/test/test_examples_2d_advdiff.jl b/test/test_examples_2d_advdiff.jl index 7b9f89a09bf..e29c6806718 100644 --- a/test/test_examples_2d_advdiff.jl +++ b/test/test_examples_2d_advdiff.jl @@ -11,14 +11,14 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d") @testset "Scalar advection-diffusion" begin @testset "elixir_advdiff_basic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_basic.jl"), - l2 = [6.707801888154751e-5], - linf = [0.0005912932758251888]) + l2 = [6.707802391539727e-5], + linf = [0.0005912935421852339]) end @testset "elixir_advdiff_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_amr.jl"), - l2 = [0.16503675954115254], - linf = [0.9236968722827653]) + l2 = [0.16503674544814984], + linf = [0.9236893624634693]) end end From a0e46948835ba907d60c25e7e9f3f2d14f586dc1 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Mon, 8 Feb 2021 14:40:52 +0100 Subject: [PATCH 53/57] Add "experimental code" warning to be able to change the implementation in a breaking way without bumping the minor version --- src/equations/heat_equation_2d.jl | 3 +++ src/equations/linear_advection_diffusion_2d.jl | 3 +++ .../semidiscretization_hyperbolic_parabolic.jl | 9 +++++++++ .../semidiscretization_parabolic_auxvars.jl | 10 ++++++++-- 4 files changed, 23 insertions(+), 2 deletions(-) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index 740ade89324..f1875bab71f 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -7,6 +7,9 @@ The heat equation \partial_t u - \nu (\partial^2_1 u + \partial^2_2 u) = 0 ``` in two space dimensions with constant viscosity ``\nu``. + +!!! warning "Experimental code" + This system of equations is experimental and can change any time. """ struct HeatEquation2D{RealT<:Real} <: AbstractHeatEquation{2, 1} nu::RealT diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index 2f2d05d2712..528cea96f4e 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -7,6 +7,9 @@ The linear advection-diffusion equation \partial_t u + a_1 \partial_1 u + a_2 \partial_2 u - \nu (\partial^2_1 u + \partial^2_2 u)= 0 ``` in two space dimensions with constant velocity `a` and constant viscosity ``\nu``. + +!!! warning "Experimental code" + This system of equations is experimental and can change any time. """ struct LinearAdvectionDiffusionEquation2D{RealT<:Real} <: AbstractLinearAdvectionDiffusionEquation{2, 1} advectionvelocity::SVector{2, RealT} diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 19a3d98aa5d..a92876a829e 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -5,6 +5,9 @@ A struct containing everything needed to describe a spatial semidiscretization of a hyperbolic-parabolic conservation law. + +!!! warning "Experimental code" + This semidiscretization is experimental and can change any time. """ struct SemidiscretizationHyperbolicParabolic{SemiHyperbolic, SemiParabolic, Cache} <: AbstractSemidiscretization semi_hyperbolic::SemiHyperbolic @@ -27,6 +30,9 @@ end Construct a semidiscretization of a hyperbolic-parabolic PDE by combining the purely hyperbolic and purely parabolic semidiscretizations. + +!!! warning "Experimental code" + This semidiscretization is experimental and can change any time. """ function SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic) cache = (; du_ode_parabolic=Vector{real(semi_parabolic)}()) @@ -145,6 +151,9 @@ scheme for a DGSEM solver. The passed arguments must correspond to the hyperboli the `solver`. Internally, a matching solver for the parabolic system is created and both a hyperbolic and a parabolic semidiscretization are created and passed to a `SemidiscretizationHyperbolicParabolic`. + +!!! warning "Experimental code" + This semidiscretization is experimental and can change any time. """ function SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver::DGSEM; source_terms=nothing, diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl index 703e9b35bec..0d7f8a53543 100644 --- a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl +++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl @@ -5,6 +5,9 @@ A struct containing everything needed to describe a spatial semidiscretization of a parabolic conservation law. + +!!! warning "Experimental code" + This semidiscretization is experimental and can change any time. """ struct SemidiscretizationParabolicAuxVars{Mesh, Equations, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache} <: AbstractSemidiscretization @@ -38,10 +41,13 @@ end """ SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver; - source_terms=nothing, - boundary_conditions=boundary_condition_periodic) + source_terms=nothing, + boundary_conditions=boundary_condition_periodic) Construct a semidiscretization of a parabolic PDE. + +!!! warning "Experimental code" + This semidiscretization is experimental and can change any time. """ function SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver; source_terms=nothing, From 8f9657afb3c8e9448e651b9a41f4b9dcfef98089 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 17 Feb 2021 10:12:36 +0100 Subject: [PATCH 54/57] Small comment fixes --- src/equations/gradient_equations_2d.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/equations/gradient_equations_2d.jl b/src/equations/gradient_equations_2d.jl index 33bace780a6..83639a7769a 100644 --- a/src/equations/gradient_equations_2d.jl +++ b/src/equations/gradient_equations_2d.jl @@ -8,6 +8,9 @@ q^d - \partial_d u = 0 ``` in direction `d` in two space dimensions as required for, e.g., the Bassi & Rebay 1 (BR1) or the local discontinuous Galerkin (LDG) schemes. + +!!! warning "Experimental code" + This system of equations is experimental and can change any time. """ struct GradientEquations2D{RealT<:Real, NVARS} <: AbstractGradientEquations{2, NVARS} orientation::Int @@ -25,7 +28,7 @@ varnames(::typeof(cons2prim), equations::GradientEquations2D) = varnames(cons2co """ initial_condition_constant(x, t, equations::GradientEquations2D) -A constant initial condition to test free-stream preservation. +A constant initial condition (zero). """ function initial_condition_constant(x, t, equations::GradientEquations2D) return SVector(ntuple(v -> zero(eltype(x)), nvariables(equations))) @@ -50,16 +53,13 @@ function boundary_condition_sin_x(u_inner, orientation, direction, x, t, return flux end - -# Pre-defined source terms should be implemented as -# function source_terms_WHATEVER(u, x, t, equations::GradientEquations2D) - - # Calculate 1D flux in for a single point @inline function calcflux(u, orientation, equations::GradientEquations2D) if orientation == equations.orientation + # In the direction specified in the constructor, the flux is equal to the current state. return -u else + # Otherwise, the flux is zero. return SVector(ntuple(v -> zero(eltype(u)), nvariables(equations))) end end From c10d7966f96d4725bda6828db8bf22593b17a9a3 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Sat, 3 Apr 2021 08:58:34 +0200 Subject: [PATCH 55/57] Replace @SVector by SVector --- src/equations/heat_equation_2d.jl | 10 +++++----- src/equations/linear_advection_diffusion_2d.jl | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl index f1875bab71f..1c15db699f4 100644 --- a/src/equations/heat_equation_2d.jl +++ b/src/equations/heat_equation_2d.jl @@ -27,7 +27,7 @@ varnames(::typeof(cons2prim), ::HeatEquation2D) = SVector("scalar") A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::HeatEquation2D) - return @SVector [2.0] + return SVector(2.0) end @@ -44,7 +44,7 @@ function initial_condition_convergence_test(x, t, equation::HeatEquation2D) f = 1/L omega = 2 * pi * f scalar = c + A * sin(omega * sum(x)) * exp(-2 * nu * omega^2 * t) - return @SVector [scalar] + return SVector(scalar) end @@ -55,7 +55,7 @@ A Gaussian pulse used together with [`boundary_condition_gauss`](@ref). """ function initial_condition_gauss(x, t, equation::HeatEquation2D) - return @SVector [exp(-(x[1]^2 + x[2]^2))] + return SVector(exp(-(x[1]^2 + x[2]^2))) end @@ -65,7 +65,7 @@ function initial_condition_sin_x(x, t, equation::HeatEquation2D) omega = pi scalar = sin(omega * x[1]) * exp(-nu * omega^2 * t) - return @SVector [scalar] + return SVector(scalar) end @@ -83,7 +83,7 @@ function initial_condition_poisson_periodic(x, t, equation::HeatEquation2D) @unpack nu = equation phi = sin(2.0*pi*x[1])*sin(2.0*pi*x[2])*(1 - exp(-8*nu*pi^2*t)) - return @SVector [phi] + return SVector(phi) end diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl index 528cea96f4e..d2474fab836 100644 --- a/src/equations/linear_advection_diffusion_2d.jl +++ b/src/equations/linear_advection_diffusion_2d.jl @@ -37,7 +37,7 @@ varnames(::typeof(cons2prim), ::LinearAdvectionDiffusionEquation2D) = SVector("s A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::LinearAdvectionDiffusionEquation2D) - return @SVector [2.0] + return SVector(2.0) end @@ -57,7 +57,7 @@ function initial_condition_convergence_test(x, t, equation::LinearAdvectionDiffu f = 1/L omega = 2 * pi * f scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) - return @SVector [scalar] + return SVector(scalar) end @@ -68,7 +68,7 @@ A Gaussian pulse used together with [`boundary_condition_gauss`](@ref). """ function initial_condition_gauss(x, t, equation::LinearAdvectionDiffusionEquation2D) - return @SVector [exp(-(x[1]^2 + x[2]^2))] + return SVector(exp(-(x[1]^2 + x[2]^2))) end From 2d8736bbc626d3ecebbec1d56aa7631b47fa9dd0 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Sat, 3 Apr 2021 09:08:34 +0200 Subject: [PATCH 56/57] Deprecate `max_dt` for `max_dt_hyperbolic` --- NEWS.md | 1 + src/callbacks_step/stepsize.jl | 4 ++++ src/semidiscretization/semidiscretization_hyperbolic.jl | 1 - 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 873b8c1c162..6417045e45f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -37,5 +37,6 @@ for human readability. `cons2prim` - `varnames_cons(equations)` → `varnames(cons2cons, equations)` - `varnames_prim(equations)` → `varnames(cons2prim, equations)` +- `max_dt(u, t, mesh, ...)` → `max_dt_hyperbolic(u, t, mesh, ...)` ([#378](https://github.com/trixi-framework/Trixi.jl/pull/378)) #### Removed diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 3f72bd30ec7..659fac3ca2e 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -94,6 +94,10 @@ function (cb::DiscreteCallback{Condition,Affect!})(ode::ODEProblem) where {Condi end +# FIXME: Deprecations introduced in v0.3 +@deprecate max_dt(u::AbstractArray, t, mesh::AbstractMesh, constant_speed, equations, dg::DG, cache) max_dt_hyperbolic(u, t, mesh, constant_speed, equations, dg, cache) + + include("stepsize_dg1d.jl") include("stepsize_dg2d.jl") include("stepsize_dg3d.jl") diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 230beab1cee..76f4750f8a6 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -195,7 +195,6 @@ function max_dt(u_ode::AbstractVector, t, cfl_number::Real, semi::Semidiscretiza mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) - # FIXME: add deprecation for max_dt -> max_dt_hyperbolic dt = cfl_number * max_dt_hyperbolic(u, t, mesh, have_constant_speed(equations), equations, solver, cache) From 5113a59ca770503bec9a6ecb8f0763bed80802a2 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Sat, 3 Apr 2021 09:20:36 +0200 Subject: [PATCH 57/57] Make CFL number an arbitray type to support multiple CFL numbers --- src/callbacks_step/stepsize.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index 659fac3ca2e..88f697ea82a 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -6,9 +6,8 @@ Set the time step size according to a CFL condition with CFL number `cfl` if the time integration method isn't adaptive itself. When using a semidiscretization that wraps multiple solvers, you might need to provide a tuple of CFL numbers. """ -mutable struct StepsizeCallback{RealT} - # FIXME: RealT -> CFLNumber - cfl_number::RealT +mutable struct StepsizeCallback{CFLNumber} + cfl_number::CFLNumber end