From 7d3925bf636e79935c231aa639377ddb31bc83e6 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 17 May 2022 15:38:19 -0500 Subject: [PATCH 01/86] minor formatting --- src/solvers/dgmulti/dg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index e4d25920eac..e23c5e8b84a 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -56,8 +56,8 @@ end # interface with semidiscretization_hyperbolic wrap_array(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode wrap_array_native(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode -function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}, mesh::DGMultiMesh, - dg::DGMulti, cache) where {Keys,ValueTypes<:NTuple{N,Any}} where {N} +function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys, ValueTypes}, mesh::DGMultiMesh, + dg::DGMulti, cache) where {Keys, ValueTypes<:NTuple{N, Any}} where {N} return boundary_conditions end From 96e16ed2813e0cce5203534cc6169f5bc0bec17c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 19 May 2022 16:00:31 -0500 Subject: [PATCH 02/86] adding multi-term semidiscretization --- src/Trixi.jl | 3 + .../semidiscretization_multi_term.jl | 182 ++++++++++++++++++ 2 files changed, 185 insertions(+) create mode 100644 src/semidiscretization/semidiscretization_multi_term.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index 3460b1761cd..f7e06e7d8ca 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -106,6 +106,7 @@ include("meshes/meshes.jl") include("solvers/solvers.jl") include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") +include("semidiscretization/semidiscretization_multi_term.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") @@ -183,6 +184,8 @@ export nelements, nnodes, nvariables, export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate +export SemidiscretizationMultiTerm + export SemidiscretizationEulerAcoustics export SemidiscretizationEulerGravity, ParametersEulerGravity, diff --git a/src/semidiscretization/semidiscretization_multi_term.jl b/src/semidiscretization/semidiscretization_multi_term.jl new file mode 100644 index 00000000000..d2aeea7ac37 --- /dev/null +++ b/src/semidiscretization/semidiscretization_multi_term.jl @@ -0,0 +1,182 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +# @muladd begin + + +""" + SemidiscretizationMultiTerm + +A struct containing everything needed to describe a spatial semidiscretization +of a conservation law with multiple equation terms (e.g., a hyperbolic and parabolic equation). +""" +struct SemidiscretizationMultiTerm{NTerms, Mesh, Equations <: NTuple{NTerms}, 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 SemidiscretizationMultiTerm{NTerms, 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 {NTerms, Mesh, Equations<:NTuple{NTerms}, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache} + @assert ndims(mesh) == ndims(first(equations)) + + # TODO: check that all equations have the same number of variables and same dimension + + performance_counter = PerformanceCounter() + + new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache, performance_counter) + end +end + +""" + SemidiscretizationMultiTerm(mesh, equations, initial_condition, solver; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + RealT=real(solver), + uEltype=RealT, + initial_cache=NamedTuple()) + +Construct a semidiscretization of a hyperbolic PDE. +""" +function SemidiscretizationMultiTerm(mesh, equations, initial_condition, solver; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT=real(solver), uEltype=RealT, + initial_cache=NamedTuple()) + + cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...) + _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) + + SemidiscretizationMultiTerm{length(equations), 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 + + +# Create a new semidiscretization but change some parameters compared to the input. +# `Base.similar` follows a related concept but would require us to `copy` the `mesh`, +# which would impact the performance. Instead, `SciMLBase.remake` has exactly the +# semantics we want to use here. In particular, it allows us to re-use mutable parts, +# e.g. `remake(semi).mesh === semi.mesh`. +function remake(semi::SemidiscretizationMultiTerm; uEltype=real(semi.solver), + mesh=semi.mesh, + equations=semi.equations, + initial_condition=semi.initial_condition, + solver=semi.solver, + source_terms=semi.source_terms, + boundary_conditions=semi.boundary_conditions + ) + # TODO: Which parts do we want to `remake`? At least the solver needs some + # special care if shock-capturing volume integrals are used (because of + # the indicators and their own caches...). + SemidiscretizationMultiTerm( + mesh, equations, initial_condition, solver; source_terms, boundary_conditions, uEltype) +end + +function Base.show(io::IO, semi::SemidiscretizationMultiTerm) + @nospecialize semi # reduce precompilation time + + print(io, "SemidiscretizationMultiTerm(") + print(io, semi.mesh) + for equation in semi.equations + print(io, ", ", equation) + end + 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::SemidiscretizationMultiTerm) + @nospecialize semi # reduce precompilation time + + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationMultiTerm") + summary_line(io, "#spatial dimensions", ndims(first(semi.equations))) + summary_line(io, "mesh", semi.mesh) + for (i, equation) in enumerate(semi.equations) + summary_line(io, "equation $i", equation |> typeof |> nameof) + end + summary_line(io, "initial condition", semi.initial_condition) + + # print_boundary_conditions(io, semi) + + summary_line(io, "source terms", semi.source_terms) + summary_line(io, "solver", semi.solver |> typeof |> nameof) + summary_line(io, "total #DOFs", ndofs(semi)) + summary_footer(io) + end +end + +@inline Base.ndims(semi::SemidiscretizationMultiTerm) = ndims(semi.mesh) + +# we assume all equations have the same number of variables +@inline nvariables(semi::SemidiscretizationMultiTerm) = nvariables(first(semi.equations)) +@inline nvariables(equations::NTuple) = nvariables(first(equations)) + +@inline Base.real(semi::SemidiscretizationMultiTerm) = real(semi.solver) + + +# to reuse existing functionality, we dispatch based on the first set of equations +@inline function mesh_equations_solver_cache(semi::SemidiscretizationMultiTerm) + @unpack mesh, equations, solver, cache = semi + return mesh, first(equations), solver, cache +end + +function calc_error_norms(func, u_ode, t, analyzer, semi::SemidiscretizationMultiTerm, cache_analysis) + @unpack mesh, equations, initial_condition, solver, cache = semi + u = wrap_array(u_ode, mesh, first(equations), solver, cache) + + calc_error_norms(func, u, t, analyzer, mesh, first(equations), initial_condition, solver, cache, cache_analysis) +end + +function compute_coefficients(t, semi::SemidiscretizationMultiTerm) + # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl` + compute_coefficients(semi.initial_condition, t, semi) +end + +function compute_coefficients!(u_ode, t, semi::SemidiscretizationMultiTerm) + compute_coefficients!(u_ode, semi.initial_condition, t, semi) +end + +function rhs!(du_ode, u_ode, semi::SemidiscretizationMultiTerm, t) + @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi + + u = wrap_array(u_ode, mesh, first(equations), solver, cache) + du = wrap_array(du_ode, mesh, first(equations), solver, cache) + + time_start = time_ns() + for (i, equation) in enumerate(equations) + @trixi_timeit timer() "$(i)th rhs!" rhs!(du, u, t, mesh, equation, initial_condition, boundary_conditions, source_terms, solver, cache) + end + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end + + +# end # @muladd From d320529c58add2509a3c1586ce60028c0840f599 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 19 May 2022 17:21:49 -0500 Subject: [PATCH 03/86] fixing bug and adding test --- .../semidiscretization_multi_term.jl | 21 ++++++++----- test/test_dgmulti_2d.jl | 30 +++++++++++++++++++ 2 files changed, 44 insertions(+), 7 deletions(-) diff --git a/src/semidiscretization/semidiscretization_multi_term.jl b/src/semidiscretization/semidiscretization_multi_term.jl index d2aeea7ac37..ff6f347c486 100644 --- a/src/semidiscretization/semidiscretization_multi_term.jl +++ b/src/semidiscretization/semidiscretization_multi_term.jl @@ -2,8 +2,7 @@ # Since these FMAs can increase the performance of many numerical algorithms, # we need to opt-in explicitly. # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -# @muladd begin - +@muladd begin """ SemidiscretizationMultiTerm @@ -34,7 +33,9 @@ struct SemidiscretizationMultiTerm{NTerms, Mesh, Equations <: NTuple{NTerms}, In solver::Solver, cache::Cache) where {NTerms, Mesh, Equations<:NTuple{NTerms}, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache} @assert ndims(mesh) == ndims(first(equations)) - # TODO: check that all equations have the same number of variables and same dimension + # check that all equations have the same number of variables and same dimension + @assert all(ndims.(equations) .== ndims(first(equations))) + @assert all(nvariables.(equations) .== nvariables(first(equations))) performance_counter = PerformanceCounter() @@ -133,7 +134,6 @@ end @inline Base.ndims(semi::SemidiscretizationMultiTerm) = ndims(semi.mesh) -# we assume all equations have the same number of variables @inline nvariables(semi::SemidiscretizationMultiTerm) = nvariables(first(semi.equations)) @inline nvariables(equations::NTuple) = nvariables(first(equations)) @@ -167,10 +167,18 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationMultiTerm, t) u = wrap_array(u_ode, mesh, first(equations), solver, cache) du = wrap_array(du_ode, mesh, first(equations), solver, cache) + du_single_term = similar(du) + # compute RHS for each term/equation time_start = time_ns() + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, solver, cache) for (i, equation) in enumerate(equations) - @trixi_timeit timer() "$(i)th rhs!" rhs!(du, u, t, mesh, equation, initial_condition, boundary_conditions, source_terms, solver, cache) + @trixi_timeit timer() "`rhs!` number $(i)" rhs!(du_single_term, u, t, mesh, equation, + initial_condition, boundary_conditions, source_terms, + solver, cache) + @threaded for i in eachindex(du) + du[i] = du[i] + du_single_term[i] + end end runtime = time_ns() - time_start put!(semi.performance_counter, runtime) @@ -178,5 +186,4 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationMultiTerm, t) return nothing end - -# end # @muladd +end # @muladd diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index d406be33c3e..760d5e17545 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -12,6 +12,36 @@ outdir = "out" isdir(outdir) && rm(outdir, recursive=true) @testset "DGMulti 2D" begin + + @trixi_testset "SemidiscretizationMultiTerm" begin + dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_central), + volume_integral = VolumeIntegralWeakForm()) + mesh = DGMultiMesh(dg, cells_per_dimension=(4, 4), + periodicity=(true, true)) + + initial_condition = initial_condition_gauss + + # split linear advection with vector = (1, 1) into two terms with + # advection vectors (1, 0) and (0, 1). + equations = LinearScalarAdvectionEquation2D(1.0, 1.0) + two_equations = (LinearScalarAdvectionEquation2D(1.0, 0.0), + LinearScalarAdvectionEquation2D(0.0, 1.0)) + + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) + semi_multi = SemidiscretizationMultiTerm(mesh, two_equations, initial_condition, dg) + + ode = semidiscretize(semi, (0.0, 0.01)) + ode_multi = semidiscretize(semi_multi, (0.0, 0.01)) + + du = similar(ode.u0) + du_multi = similar(ode_multi.u0) + Trixi.rhs!(du, ode.u0, semi, 0.0) + Trixi.rhs!(du_multi, ode_multi.u0, semi_multi, 0.0) + + @test norm(norm.(du .- du_multi)) < 100 * eps() + end + @trixi_testset "elixir_euler_weakform.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), cells_per_dimension = (4, 4), From 84006ba6d6e2567b9c379edfed02dee6e8ebaa8b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 20 May 2022 11:19:50 -0500 Subject: [PATCH 04/86] changing name --- ...semidiscretization_hyperbolic_parabolic.jl | 236 ++++++++++++++++++ .../semidiscretization_multi_term.jl | 189 -------------- test/test_dgmulti_2d.jl | 20 +- 3 files changed, 248 insertions(+), 197 deletions(-) create mode 100644 src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl delete mode 100644 src/semidiscretization/semidiscretization_multi_term.jl diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl new file mode 100644 index 00000000000..5231697c60c --- /dev/null +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -0,0 +1,236 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + + +""" + SemidiscretizationHyperbolicParabolic + +A struct containing everything needed to describe a spatial semidiscretization +of a mixed hyperbolic-parabolic conservation law. +""" + +struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations, InitialCondition, + BoundaryConditions, ParabolicBoundaryConditions, + SourceTerms, Solver, Cache, ParabolicCache} <: AbstractSemidiscretization + + mesh::Mesh + + equations::Equations + parabolic_equations::ParabolicEquations + + # 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 + parabolic_boundary_conditions::ParabolicBoundaryConditions + + source_terms::SourceTerms + + solver::Solver + # TODO: do we want to introduce `parabolic_solver` for future specialization? + + cache::Cache + parabolic_cache::ParabolicCache + + performance_counter::PerformanceCounter + + function SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations, InitialCondition, BoundaryConditions, ParabolicBoundaryConditions, SourceTerms, Solver, Cache, ParabolicCache}( + mesh::Mesh, equations::Equations, parabolic_equations::ParabolicEquations, initial_condition::InitialCondition, + boundary_conditions::BoundaryConditions, parabolic_boundary_conditions::ParabolicBoundaryConditions, + source_terms::SourceTerms, solver::Solver, cache::Cache, parabolic_cache::ParabolicCache) where {Mesh, Equations, ParabolicEquations, InitialCondition, BoundaryConditions, ParabolicBoundaryConditions, SourceTerms, Solver, Cache, ParabolicCache} + @assert ndims(mesh) == ndims(equations) + + # Todo: assert nvariables(equations)==nvariables(parabolic_equations) + + performance_counter = PerformanceCounter() + + new(mesh, equations, parabolic_equations, initial_condition, + boundary_conditions, parabolic_boundary_conditions, + source_terms, solver, cache, parabolic_cache, performance_counter) + end +end + +""" + SemidiscretizationHyperbolicParabolic(mesh, equations, initial_condition, solver; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + RealT=real(solver), + uEltype=RealT, + initial_cache=NamedTuple()) + +Construct a semidiscretization of a hyperbolic PDE. +""" +function SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, + initial_condition, solver; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + parabolic_boundary_conditions=boundary_condition_periodic, + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT=real(solver), uEltype=RealT, + initial_cache=NamedTuple(), + initial_parabolic_cache=NamedTuple()) + + cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...) + _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) + _parabolic_boundary_conditions = digest_boundary_conditions(parabolic_boundary_conditions, mesh, solver, cache) + + parabolic_cache = (; create_cache(mesh, parabolic_equations, solver, RealT, uEltype)..., + initial_parabolic_cache...) + + SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(parabolic_equations), + typeof(initial_condition), typeof(_boundary_conditions), typeof(_parabolic_boundary_conditions), + typeof(source_terms), typeof(solver), typeof(cache), typeof(parabolic_cache)}( + mesh, equations, parabolic_equations, initial_condition, + _boundary_conditions, _parabolic_boundary_conditions, source_terms, + solver, cache, parabolic_cache) +end + + +# Create a new semidiscretization but change some parameters compared to the input. +# `Base.similar` follows a related concept but would require us to `copy` the `mesh`, +# which would impact the performance. Instead, `SciMLBase.remake` has exactly the +# semantics we want to use here. In particular, it allows us to re-use mutable parts, +# e.g. `remake(semi).mesh === semi.mesh`. +function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.solver), + mesh=semi.mesh, + equations=semi.equations, + parabolic_equations=semi.parabolic_equations, + initial_condition=semi.initial_condition, + solver=semi.solver, + source_terms=semi.source_terms, + boundary_conditions=semi.boundary_conditions, + parabolic_boundary_conditions=semi.parabolic_boundary_conditions + ) + # TODO: Which parts do we want to `remake`? At least the solver needs some + # special care if shock-capturing volume integrals are used (because of + # the indicators and their own caches...). + SemidiscretizationHyperbolicParabolic( + mesh, equations, parabolic_equations, initial_condition, solver; source_terms, boundary_conditions, parabolic_boundary_conditions, uEltype) +end + +function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) + @nospecialize semi # reduce precompilation time + + print(io, "SemidiscretizationHyperbolicParabolic(") + print(io, semi.mesh) + print(io, ", ", semi.equations) + print(io, ", ", semi.parabolic_equations) + print(io, ", ", semi.initial_condition) + print(io, ", ", semi.boundary_conditions) + print(io, ", ", semi.parabolic_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::SemidiscretizationHyperbolicParabolic) + @nospecialize semi # reduce precompilation time + + if get(io, :compact, false) + show(io, semi) + else + summary_header(io, "SemidiscretizationHyperbolicParabolic") + summary_line(io, "#spatial dimensions", ndims(semi.equations)) + summary_line(io, "mesh", semi.mesh) + summary_line(io, "equations", semi.equations |> typeof |> nameof) + summary_line(io, "parabolic_equations", semi.parabolic_equations |> typeof |> nameof) + summary_line(io, "initial condition", semi.initial_condition) + + # print_boundary_conditions(io, semi) + + summary_line(io, "source terms", semi.source_terms) + summary_line(io, "solver", semi.solver |> typeof |> nameof) + summary_line(io, "total #DOFs", ndofs(semi)) + summary_footer(io) + end +end + +@inline Base.ndims(semi::SemidiscretizationHyperbolicParabolic) = ndims(semi.mesh) + +@inline nvariables(semi::SemidiscretizationHyperbolicParabolic) = nvariables(semi.equations) + +@inline Base.real(semi::SemidiscretizationHyperbolicParabolic) = real(semi.solver) + +# retain dispatch on hyperbolic equations only +@inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolicParabolic) + @unpack mesh, equations, solver, cache = semi + return mesh, equations, solver, cache +end + + +function calc_error_norms(func, u_ode, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, 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, cache_analysis) +end + + +function compute_coefficients(t, semi::SemidiscretizationHyperbolicParabolic) + # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl` + compute_coefficients(semi.initial_condition, t, semi) +end + +function compute_coefficients!(u_ode, t, semi::SemidiscretizationHyperbolicParabolic) + compute_coefficients!(u_ode, semi.initial_condition, t, semi) +end + +""" + semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan) + +Wrap the semidiscretization `semi` as a Split ODE problem in the time interval `tspan` +that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). +""" +function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan) + u0_ode = compute_coefficients(first(tspan), semi) + # TODO: MPI, do we want to synchonize loading and print debug statements, e.g. using + # mpi_isparallel() && MPI.Barrier(mpi_comm()) + # See https://github.com/trixi-framework/Trixi.jl/issues/328 + return SplitODEProblem(rhs!, rhs_parabolic!, u0_ode, tspan, semi) +end + +function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, 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) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @trixi_timeit timer() "rhs!" rhs!(du, u, t, mesh, equations, initial_condition, + boundary_conditions, source_terms, solver, cache) + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end + +function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) + @unpack mesh, parabolic_equations, initial_condition, parabolic_boundary_conditions, source_terms, solver, parabolic_cache = semi + + u = wrap_array(u_ode, mesh, parabolic_equations, solver, parabolic_cache) + du = wrap_array(du_ode, mesh, parabolic_equations, solver, parabolic_cache) + + # TODO: Taal decide, do we need to pass the mesh? + time_start = time_ns() + @trixi_timeit timer() "parabolic rhs!" rhs!(du, u, t, mesh, parabolic_equations, initial_condition, + parabolic_boundary_conditions, source_terms, + solver, parabolic_cache) + runtime = time_ns() - time_start + put!(semi.performance_counter, runtime) + + return nothing +end + + +end # @muladd diff --git a/src/semidiscretization/semidiscretization_multi_term.jl b/src/semidiscretization/semidiscretization_multi_term.jl deleted file mode 100644 index ff6f347c486..00000000000 --- a/src/semidiscretization/semidiscretization_multi_term.jl +++ /dev/null @@ -1,189 +0,0 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin - -""" - SemidiscretizationMultiTerm - -A struct containing everything needed to describe a spatial semidiscretization -of a conservation law with multiple equation terms (e.g., a hyperbolic and parabolic equation). -""" -struct SemidiscretizationMultiTerm{NTerms, Mesh, Equations <: NTuple{NTerms}, 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 SemidiscretizationMultiTerm{NTerms, 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 {NTerms, Mesh, Equations<:NTuple{NTerms}, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache} - @assert ndims(mesh) == ndims(first(equations)) - - # check that all equations have the same number of variables and same dimension - @assert all(ndims.(equations) .== ndims(first(equations))) - @assert all(nvariables.(equations) .== nvariables(first(equations))) - - performance_counter = PerformanceCounter() - - new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache, performance_counter) - end -end - -""" - SemidiscretizationMultiTerm(mesh, equations, initial_condition, solver; - source_terms=nothing, - boundary_conditions=boundary_condition_periodic, - RealT=real(solver), - uEltype=RealT, - initial_cache=NamedTuple()) - -Construct a semidiscretization of a hyperbolic PDE. -""" -function SemidiscretizationMultiTerm(mesh, equations, initial_condition, solver; - source_terms=nothing, - boundary_conditions=boundary_condition_periodic, - # `RealT` is used as real type for node locations etc. - # while `uEltype` is used as element type of solutions etc. - RealT=real(solver), uEltype=RealT, - initial_cache=NamedTuple()) - - cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...) - _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) - - SemidiscretizationMultiTerm{length(equations), 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 - - -# Create a new semidiscretization but change some parameters compared to the input. -# `Base.similar` follows a related concept but would require us to `copy` the `mesh`, -# which would impact the performance. Instead, `SciMLBase.remake` has exactly the -# semantics we want to use here. In particular, it allows us to re-use mutable parts, -# e.g. `remake(semi).mesh === semi.mesh`. -function remake(semi::SemidiscretizationMultiTerm; uEltype=real(semi.solver), - mesh=semi.mesh, - equations=semi.equations, - initial_condition=semi.initial_condition, - solver=semi.solver, - source_terms=semi.source_terms, - boundary_conditions=semi.boundary_conditions - ) - # TODO: Which parts do we want to `remake`? At least the solver needs some - # special care if shock-capturing volume integrals are used (because of - # the indicators and their own caches...). - SemidiscretizationMultiTerm( - mesh, equations, initial_condition, solver; source_terms, boundary_conditions, uEltype) -end - -function Base.show(io::IO, semi::SemidiscretizationMultiTerm) - @nospecialize semi # reduce precompilation time - - print(io, "SemidiscretizationMultiTerm(") - print(io, semi.mesh) - for equation in semi.equations - print(io, ", ", equation) - end - 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::SemidiscretizationMultiTerm) - @nospecialize semi # reduce precompilation time - - if get(io, :compact, false) - show(io, semi) - else - summary_header(io, "SemidiscretizationMultiTerm") - summary_line(io, "#spatial dimensions", ndims(first(semi.equations))) - summary_line(io, "mesh", semi.mesh) - for (i, equation) in enumerate(semi.equations) - summary_line(io, "equation $i", equation |> typeof |> nameof) - end - summary_line(io, "initial condition", semi.initial_condition) - - # print_boundary_conditions(io, semi) - - summary_line(io, "source terms", semi.source_terms) - summary_line(io, "solver", semi.solver |> typeof |> nameof) - summary_line(io, "total #DOFs", ndofs(semi)) - summary_footer(io) - end -end - -@inline Base.ndims(semi::SemidiscretizationMultiTerm) = ndims(semi.mesh) - -@inline nvariables(semi::SemidiscretizationMultiTerm) = nvariables(first(semi.equations)) -@inline nvariables(equations::NTuple) = nvariables(first(equations)) - -@inline Base.real(semi::SemidiscretizationMultiTerm) = real(semi.solver) - - -# to reuse existing functionality, we dispatch based on the first set of equations -@inline function mesh_equations_solver_cache(semi::SemidiscretizationMultiTerm) - @unpack mesh, equations, solver, cache = semi - return mesh, first(equations), solver, cache -end - -function calc_error_norms(func, u_ode, t, analyzer, semi::SemidiscretizationMultiTerm, cache_analysis) - @unpack mesh, equations, initial_condition, solver, cache = semi - u = wrap_array(u_ode, mesh, first(equations), solver, cache) - - calc_error_norms(func, u, t, analyzer, mesh, first(equations), initial_condition, solver, cache, cache_analysis) -end - -function compute_coefficients(t, semi::SemidiscretizationMultiTerm) - # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl` - compute_coefficients(semi.initial_condition, t, semi) -end - -function compute_coefficients!(u_ode, t, semi::SemidiscretizationMultiTerm) - compute_coefficients!(u_ode, semi.initial_condition, t, semi) -end - -function rhs!(du_ode, u_ode, semi::SemidiscretizationMultiTerm, t) - @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi - - u = wrap_array(u_ode, mesh, first(equations), solver, cache) - du = wrap_array(du_ode, mesh, first(equations), solver, cache) - du_single_term = similar(du) - - # compute RHS for each term/equation - time_start = time_ns() - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, solver, cache) - for (i, equation) in enumerate(equations) - @trixi_timeit timer() "`rhs!` number $(i)" rhs!(du_single_term, u, t, mesh, equation, - initial_condition, boundary_conditions, source_terms, - solver, cache) - @threaded for i in eachindex(du) - du[i] = du[i] + du_single_term[i] - end - end - runtime = time_ns() - time_start - put!(semi.performance_counter, runtime) - - return nothing -end - -end # @muladd diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 760d5e17545..fa08861f39e 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -13,7 +13,7 @@ isdir(outdir) && rm(outdir, recursive=true) @testset "DGMulti 2D" begin - @trixi_testset "SemidiscretizationMultiTerm" begin + @trixi_testset "SemidiscretizationHyperbolicParabolic" begin dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), surface_integral = SurfaceIntegralWeakForm(flux_central), volume_integral = VolumeIntegralWeakForm()) @@ -25,21 +25,25 @@ isdir(outdir) && rm(outdir, recursive=true) # split linear advection with vector = (1, 1) into two terms with # advection vectors (1, 0) and (0, 1). equations = LinearScalarAdvectionEquation2D(1.0, 1.0) - two_equations = (LinearScalarAdvectionEquation2D(1.0, 0.0), - LinearScalarAdvectionEquation2D(0.0, 1.0)) + equation1 = LinearScalarAdvectionEquation2D(1.0, 0.0) + equation2 = LinearScalarAdvectionEquation2D(0.0, 1.0) # our "parabolic" equation semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) - semi_multi = SemidiscretizationMultiTerm(mesh, two_equations, initial_condition, dg) + semi_hyp_par = SemidiscretizationHyperbolicParabolic(mesh, equation1, equation2, initial_condition, dg) ode = semidiscretize(semi, (0.0, 0.01)) - ode_multi = semidiscretize(semi_multi, (0.0, 0.01)) + ode_hyp_par = semidiscretize(semi_hyp_par, (0.0, 0.01)) du = similar(ode.u0) - du_multi = similar(ode_multi.u0) + du_hyp = similar(ode_hyp_par.u0) + du_par = similar(ode_hyp_par.u0) Trixi.rhs!(du, ode.u0, semi, 0.0) - Trixi.rhs!(du_multi, ode_multi.u0, semi_multi, 0.0) + Trixi.rhs!(du_hyp, ode_hyp_par.u0, semi_hyp_par, 0.0) + Trixi.rhs_parabolic!(du_par, ode_hyp_par.u0, semi_hyp_par, 0.0) - @test norm(norm.(du .- du_multi)) < 100 * eps() + # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs + # for `equation` with SemidiscretizationHyperbolic. + @test norm(norm.(du .- (du_hyp .+ du_par))) < 100 * eps() end @trixi_testset "elixir_euler_weakform.jl" begin From 93c718eb00bcc86a1835b39da47bfea313987ae1 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 20 May 2022 11:20:05 -0500 Subject: [PATCH 05/86] changing name, importing SplitODEFunction --- src/Trixi.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index f7e06e7d8ca..752f0612990 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -26,7 +26,8 @@ using SparseArrays: AbstractSparseMatrix, AbstractSparseMatrixCSC, sparse, dropt using Reexport: @reexport using SciMLBase: CallbackSet, DiscreteCallback, - ODEProblem, ODESolution, ODEFunction + ODEProblem, ODESolution, ODEFunction, + SplitODEProblem import SciMLBase: get_du, get_tmp_cache, u_modified!, AbstractODEIntegrator, init, step!, check_error, get_proposed_dt, set_proposed_dt!, @@ -106,7 +107,7 @@ include("meshes/meshes.jl") include("solvers/solvers.jl") include("semidiscretization/semidiscretization.jl") include("semidiscretization/semidiscretization_hyperbolic.jl") -include("semidiscretization/semidiscretization_multi_term.jl") +include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") @@ -184,7 +185,7 @@ export nelements, nnodes, nvariables, export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate -export SemidiscretizationMultiTerm +export SemidiscretizationHyperbolicParabolic export SemidiscretizationEulerAcoustics From 7dc7b2d6f8d2ee6af362286fb41252b203aa33ff Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 20 May 2022 12:31:15 -0500 Subject: [PATCH 06/86] fixing test by adding norm --- test/test_dgmulti_2d.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index fa08861f39e..1cf5193be74 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -2,6 +2,7 @@ module TestExamplesDGMulti2D using Test using Trixi +using LinearAlgebra: norm include("test_trixi.jl") From 1772269092d6189457982b709f8a85f2e0994c93 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 07:22:42 -0500 Subject: [PATCH 07/86] trying another test fix --- test/test_dgmulti_2d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 1cf5193be74..ec5af9ac81d 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -2,7 +2,6 @@ module TestExamplesDGMulti2D using Test using Trixi -using LinearAlgebra: norm include("test_trixi.jl") @@ -44,7 +43,7 @@ isdir(outdir) && rm(outdir, recursive=true) # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs # for `equation` with SemidiscretizationHyperbolic. - @test norm(norm.(du .- (du_hyp .+ du_par))) < 100 * eps() + @test all(du .≈ (du_hyp .+ du_par)) end @trixi_testset "elixir_euler_weakform.jl" begin From fae867b728e50334543d152600339a18f266808b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 08:19:43 -0500 Subject: [PATCH 08/86] adding norm back --- test/test_dgmulti_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index ec5af9ac81d..1cf5193be74 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -2,6 +2,7 @@ module TestExamplesDGMulti2D using Test using Trixi +using LinearAlgebra: norm include("test_trixi.jl") @@ -43,7 +44,7 @@ isdir(outdir) && rm(outdir, recursive=true) # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs # for `equation` with SemidiscretizationHyperbolic. - @test all(du .≈ (du_hyp .+ du_par)) + @test norm(norm.(du .- (du_hyp .+ du_par))) < 100 * eps() end @trixi_testset "elixir_euler_weakform.jl" begin From cd16b365f3f9b66034022549b1b3390ef48d8fc7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 08:21:21 -0500 Subject: [PATCH 09/86] adjusting tol --- test/test_dgmulti_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 1cf5193be74..c9af7d14ebc 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -44,7 +44,7 @@ isdir(outdir) && rm(outdir, recursive=true) # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs # for `equation` with SemidiscretizationHyperbolic. - @test norm(norm.(du .- (du_hyp .+ du_par))) < 100 * eps() + @test norm(norm.(du .- (du_hyp .+ du_par))) < 10 * length(du) * eps() end @trixi_testset "elixir_euler_weakform.jl" begin From d6690ba6c203cbd955f6b1360e3f7093b4d79fa0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 12:40:01 -0500 Subject: [PATCH 10/86] try to fix test again --- test/test_dgmulti_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index c9af7d14ebc..97d94bdcb62 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -2,7 +2,7 @@ module TestExamplesDGMulti2D using Test using Trixi -using LinearAlgebra: norm +using Trixi.LinearAlgebra.norm include("test_trixi.jl") From fe0203af4608cb8012b74ab894d5822271060524 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 13:01:41 -0500 Subject: [PATCH 11/86] syntax error --- test/test_dgmulti_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 97d94bdcb62..9fb40635839 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -2,7 +2,7 @@ module TestExamplesDGMulti2D using Test using Trixi -using Trixi.LinearAlgebra.norm +using Trixi.LinearAlgebra: norm include("test_trixi.jl") From 64e6ebf1aa504a090e01a0631ab771bddb4c98d1 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 14:27:41 -0500 Subject: [PATCH 12/86] removing norm frm test --- test/test_dgmulti_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 9fb40635839..511d5fdda5c 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -2,7 +2,6 @@ module TestExamplesDGMulti2D using Test using Trixi -using Trixi.LinearAlgebra: norm include("test_trixi.jl") @@ -44,7 +43,8 @@ isdir(outdir) && rm(outdir, recursive=true) # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs # for `equation` with SemidiscretizationHyperbolic. - @test norm(norm.(du .- (du_hyp .+ du_par))) < 10 * length(du) * eps() + du_diff = du .- (du_hyp .+ du_par) + @test sum(abs.(getindex.(du_diff, 1))) < 10 * length(du) * eps() end @trixi_testset "elixir_euler_weakform.jl" begin From dd60ebf5770c712beaa4a790832770d855e17b5d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 22:49:07 -0500 Subject: [PATCH 13/86] adding cache to parabolic terms (enable reuse of inviscid cache) --- .../semidiscretization_hyperbolic_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 5231697c60c..d218d306f12 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -225,7 +225,7 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol time_start = time_ns() @trixi_timeit timer() "parabolic rhs!" rhs!(du, u, t, mesh, parabolic_equations, initial_condition, parabolic_boundary_conditions, source_terms, - solver, parabolic_cache) + solver, cache, parabolic_cache) runtime = time_ns() - time_start put!(semi.performance_counter, runtime) From c141fc17fc78765ea82e3d8ff2667aec9e2a1c39 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 22:49:50 -0500 Subject: [PATCH 14/86] generalize dgmulti iterators --- src/solvers/dgmulti/dg.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index e23c5e8b84a..371b8c96da5 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -41,7 +41,7 @@ end # iteration over all elements in a mesh @inline ndofs(mesh::DGMultiMesh, dg::DGMulti, cache) = dg.basis.Np * mesh.md.num_elements -@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(mesh.md.num_elements) +@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, other_args...) = eachelement(mesh, dg) # iteration over quantities in a single element @inline nnodes(basis::RefElemData) = basis.Np @@ -50,8 +50,8 @@ end # iteration over quantities over the entire mesh (dofs, quad nodes, face nodes). @inline each_dof_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(ndofs(mesh, dg, cache)) -@inline each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nq * mesh.md.num_elements) -@inline each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nfq * mesh.md.num_elements) +@inline each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nq * mesh.md.num_elements) +@inline each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nfq * mesh.md.num_elements) # interface with semidiscretization_hyperbolic wrap_array(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode From 943450f7313c30834be63c0705969bbc3e8a2e09 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 22:53:14 -0500 Subject: [PATCH 15/86] fix typo --- src/solvers/dgmulti/dg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 371b8c96da5..9a213ce51d5 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -41,7 +41,7 @@ end # iteration over all elements in a mesh @inline ndofs(mesh::DGMultiMesh, dg::DGMulti, cache) = dg.basis.Np * mesh.md.num_elements -@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, other_args...) = eachelement(mesh, dg) +@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(mesh.md.num_elements) # iteration over quantities in a single element @inline nnodes(basis::RefElemData) = basis.Np From 093f27f7d079d62d5deef437975b7493b55acc83 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 23:06:55 -0500 Subject: [PATCH 16/86] generalizing reset_du! for dgmulti --- src/solvers/dgmulti/dg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 9a213ce51d5..27cb61b698a 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -67,7 +67,7 @@ function allocate_nested_array(uEltype, nvars, array_dimensions, dg) return StructArray{SVector{nvars, uEltype}}(ntuple(_->zeros(uEltype, array_dimensions...), nvars)) end -function reset_du!(du, dg::DGMulti, cache) +function reset_du!(du, dg::DGMulti, other_args...) @threaded for i in eachindex(du) du[i] = zero(eltype(du)) end From ace2ce1d48eb36ed95a18f4e5357a50d7d92ad19 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 23:36:23 -0500 Subject: [PATCH 17/86] WIP gradient and divergence computations --- src/Trixi.jl | 3 + src/solvers/dgmulti/diffusion.jl | 147 +++++++++++++++++++++++++++++++ 2 files changed, 150 insertions(+) create mode 100644 src/solvers/dgmulti/diffusion.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index c7a3de79d7d..b7aa4b6e2f6 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -114,6 +114,9 @@ include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") include("time_integration/time_integration.jl") +# TODO: move this around +include("solvers/dgmulti/diffusion.jl") + # `trixi_include` and special elixirs such as `convergence_test` include("auxiliary/special_elixirs.jl") diff --git a/src/solvers/dgmulti/diffusion.jl b/src/solvers/dgmulti/diffusion.jl new file mode 100644 index 00000000000..e7aabfa5e7a --- /dev/null +++ b/src/solvers/dgmulti/diffusion.jl @@ -0,0 +1,147 @@ + +abstract type AbstractParabolicEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end + +struct ScalarDiffusion2D{T} <: AbstractParabolicEquations{2, 1} + diffusivity::T +end + +# no orientation here since the flux is vector-valued +function flux(u, grad_u, equations::ScalarDiffusion) + dudx, dudy = grad_u + return equations.diffusivity * dudx, equations.diffusivity * dudy +end + +function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, + dg::DGMultiWeakForm, RealT, uEltype) + nvars = nvariables(equations) + + # u_parabolic stores "transformed" variables for + @unpack md = mesh + u_transformed = allocate_nested_array(uEltype, nvars, size(md.xq), dg) + u_grad = ntuple(_ -> similar(u_transformed), ndims(mesh)) + u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) + flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) + return (; u_transformed, u_grad, viscous_flux=similar.(u_grad), u_face_values, + flux_face_values, local_flux_values = similar(flux_face_values[:, 1])) +end + +# Transform variables prior to taking the gradient. Defaults to doing nothing. +# TODO: can we avoid copying data? +function transform_variables(u_transformed, u, equations) + @threaded for i in eachindex(u) + u_transformed[i] = u[i] + end +end + +# interpolates from solution coefficients to face quadrature points +function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, + surface_integral, dg::DGMulti, cache) + rd = dg.basis + apply_to_each_field(mul_by!(rd.Vf), u_face_values, u) +end + +function calc_gradient_surface_integral(u_grad, u, flux_face_values, + mesh, equations::AbstractParabolicEquations, + dg::DGMulti, cache, parabolic_cache) + @unpack local_flux_values = parabolic_cache + @threaded for e in eachelement(mesh, dg) + for dim in eachdim(mesh) + for i in eachindex(local_flux_values) + # compute [u] * (nx, ny, nz) + local_flux_values[i] = flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] + end + apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(u_grad[dim], :, e), local_flux_values) + end + end +end + +function calc_gradient!(u_grad, u::StructArray, mesh::DGMultiMesh, + equations::AbstractParabolicEquations, + boundary_conditions, dg::DGMulti, cache, parabolic_cache) + + for dim in 1:length(u_grad) + reset_du!(u_grad[dim], dg) + end + + # compute volume contributions to gradients + @threaded for e in eachelement(mesh, dg) + for i in eachdim(mesh), j in eachdim(mesh) + dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine + StructArrays.foreachfield(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), + view(u_grad[i], :, e), view(u, :, e)) + end + end + + prolong2interfaces!(cache.u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) + + # compute fluxes at interfaces + @unpack u_face_values, flux_face_values = cache + @unpack mapM, mapP, Jf = mesh.md + @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) + idM, idP = mapM[face_node_index], mapP[face_node_index] + uM = u_face_values[idM] + # compute flux if node is not a boundary node + if idM != idP + uP = u_face_values[idP] + flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? + end + end + + # compute surface contributions + calc_gradient_surface_integral(u_grad, u, flux_face_values, + mesh, equations, dg, cache, parabolic_cache) + + calc_gradient_boundary_integral!(u_grad, u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) + + for dim in eachdim(mesh) + invert_jacobian!(u_grad[dim], mesh, equations, dg, cache) + end +end + +# This routine differs from the one in dgmulti/dg.jl in that we do not negate the result. +function invert_jacobian!(du, mesh::DGMultiMesh, equations::AbstractParabolicEquations, + dg::DGMulti, cache) + @threaded for i in Base.OneTo(length(du)) + du[i] = du[i] * cache.invJ[i] + end +end + +# do nothing for periodic domains +function calc_gradient_boundary_integral!(du, u, mesh, equations, ::BoundaryConditionPeriodic, + dg, cache, parabolic_cache) + return nothing +end + +function calc_viscous_fluxes!(u_flux, u, grad_u, mesh::DGMultiMesh, + equations::AbstractParabolicEquations, + dg::DGMulti, cache, parabolic_cache) + @threaded for i in eachindex(u) + u_flux[i] = flux(u, grad_u, equations) + end +end + +# function calc_divergence!(du, u::StructArray, mesh::DGMultiMesh, +# equations::AbstractParabolicEquations, +# boundary_conditions, dg::DGMulti, cache, parabolic_cache) +# calc_divergence_volume_integral!(du, u, mesh, equations, dg, cache, parabolic_cache) +# calc_divergence_surface_integral!(du, u, mesh, equations, dg, cache, parabolic_cache) +# calc_divergence_boundary_integral!(du, u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) +# end + +# # assumptions: parabolic terms are of the form div(f(u, grad(u))) and +# # will be discretized in first order form +# # - compute grad(u) +# # - compute f(u, grad(u)) +# # - compute div(u) +# # boundary conditions will be applied to both grad(u) and div(u). +# function rhs!(du, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, +# initial_condition, boundary_conditions, source_terms, +# dg::DGMulti, cache, parabolic_cache) +# @unpack u_transformed, grad_u, viscous_flux = parabolic_cache +# transform_variables!(u_transformed, u, equations) +# calc_gradient!(grad_u, u_transformed, mesh, equations, +# boundary_conditions, dg, cache, parabolic_cache) +# calc_viscous_fluxes!(viscous_flux, u_transformed, grad_u, +# mesh, equations, dg, cache, parabolic_cache) +# calc_divergence!(du, grad_u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) +# end From 9133d8456e570a4009a5e5f6073d43ad5c824dee Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sat, 21 May 2022 23:49:02 -0500 Subject: [PATCH 18/86] fix name --- src/solvers/dgmulti/diffusion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/diffusion.jl b/src/solvers/dgmulti/diffusion.jl index e7aabfa5e7a..cf7d1990aa3 100644 --- a/src/solvers/dgmulti/diffusion.jl +++ b/src/solvers/dgmulti/diffusion.jl @@ -6,7 +6,7 @@ struct ScalarDiffusion2D{T} <: AbstractParabolicEquations{2, 1} end # no orientation here since the flux is vector-valued -function flux(u, grad_u, equations::ScalarDiffusion) +function flux(u, grad_u, equations::ScalarDiffusion2D) dudx, dudy = grad_u return equations.diffusivity * dudx, equations.diffusivity * dudy end From ade76e3bce77ac110de47a29ca25371890bb0707 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 00:20:27 -0500 Subject: [PATCH 19/86] more generalization --- src/solvers/dgmulti/dg.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 27cb61b698a..9952c84b989 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -40,16 +40,16 @@ end @inline eachdim(mesh) = Base.OneTo(ndims(mesh)) # iteration over all elements in a mesh -@inline ndofs(mesh::DGMultiMesh, dg::DGMulti, cache) = dg.basis.Np * mesh.md.num_elements +@inline ndofs(mesh::DGMultiMesh, dg::DGMulti, other_args...) = dg.basis.Np * mesh.md.num_elements @inline eachelement(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(mesh.md.num_elements) # iteration over quantities in a single element @inline nnodes(basis::RefElemData) = basis.Np -@inline each_face_node(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nfq) -@inline each_quad_node(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nq) +@inline each_face_node(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nfq) +@inline each_quad_node(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nq) # iteration over quantities over the entire mesh (dofs, quad nodes, face nodes). -@inline each_dof_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(ndofs(mesh, dg, cache)) +@inline each_dof_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(ndofs(mesh, dg, other_args...)) @inline each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nq * mesh.md.num_elements) @inline each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nfq * mesh.md.num_elements) @@ -353,9 +353,10 @@ end # Todo: DGMulti. Specialize for modal DG on curved meshes using WADG -function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache) +# inverts Jacobian and scales by -1.0 +function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1.0) @threaded for i in each_dof_global(mesh, dg, cache) - du[i] *= -cache.invJ[i] + du[i] *= scaling * cache.invJ[i] end end From 54a4b7ec99e55e5f0526b05b413206216261c210 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 00:20:41 -0500 Subject: [PATCH 20/86] divergence + remove repeated code --- src/solvers/dgmulti/diffusion.jl | 129 +++++++++++++++++++++---------- 1 file changed, 90 insertions(+), 39 deletions(-) diff --git a/src/solvers/dgmulti/diffusion.jl b/src/solvers/dgmulti/diffusion.jl index cf7d1990aa3..45abff9810c 100644 --- a/src/solvers/dgmulti/diffusion.jl +++ b/src/solvers/dgmulti/diffusion.jl @@ -20,9 +20,16 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, u_transformed = allocate_nested_array(uEltype, nvars, size(md.xq), dg) u_grad = ntuple(_ -> similar(u_transformed), ndims(mesh)) u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) + grad_u_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) + + viscous_flux = similar.(u_grad) + local_viscous_flux = ntuple(_ -> similar(u_transformed, size(md.xq, 1)), ndims(mesh)) + flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) - return (; u_transformed, u_grad, viscous_flux=similar.(u_grad), u_face_values, - flux_face_values, local_flux_values = similar(flux_face_values[:, 1])) + return (; u_transformed, u_grad, + viscous_flux, local_viscous_flux, + u_face_values, grad_u_face_values, + flux_face_values, local_flux_values=similar(flux_face_values[:, 1])) end # Transform variables prior to taking the gradient. Defaults to doing nothing. @@ -94,15 +101,7 @@ function calc_gradient!(u_grad, u::StructArray, mesh::DGMultiMesh, calc_gradient_boundary_integral!(u_grad, u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) for dim in eachdim(mesh) - invert_jacobian!(u_grad[dim], mesh, equations, dg, cache) - end -end - -# This routine differs from the one in dgmulti/dg.jl in that we do not negate the result. -function invert_jacobian!(du, mesh::DGMultiMesh, equations::AbstractParabolicEquations, - dg::DGMulti, cache) - @threaded for i in Base.OneTo(length(du)) - du[i] = du[i] * cache.invJ[i] + invert_jacobian!(u_grad[dim], mesh, equations, dg, cache; scaling=1.0) end end @@ -112,36 +111,88 @@ function calc_gradient_boundary_integral!(du, u, mesh, equations, ::BoundaryCond return nothing end -function calc_viscous_fluxes!(u_flux, u, grad_u, mesh::DGMultiMesh, +function calc_viscous_fluxes!(u_flux, u, u_grad, mesh::DGMultiMesh, equations::AbstractParabolicEquations, dg::DGMulti, cache, parabolic_cache) - @threaded for i in eachindex(u) - u_flux[i] = flux(u, grad_u, equations) + + @unpack local_viscous_flux = parabolic_cache + for dim in eachdim(mesh) + reset_du!(local_viscous_flux[dim], dg) + end + + @threaded for e in eachelement(mesh, dg) + # interpolate to quadrature points + for dim in eachdim(mesh) + StructArrays.foreachfield(mul_by!(dg.basis.Vq), local_viscous_flux[dim], view(u_grad[dim], :, e)) + end + + # compute viscous flux at quad points + for i in eachindex(local_viscous_flux) + u_grad_i = getindex.(local_viscous_flux, i) # TODO: check if this allocates. Shouldn't for tuples or SVector... + setindex!.(local_viscous_flux, flux(u[i,e], u_grad_i, equations), i) + end + + # project back to the DG approximation space + for dim in eachdim(mesh) + StructArrays.foreachfield(mul_by!(dg.basis.Pq), view(u_flux[dim], :, e), local_viscous_flux[dim]) + end + end +end + +function calc_divergence!(du, u::StructArray, u_flux, mesh::DGMultiMesh, + equations::AbstractParabolicEquations, + boundary_conditions, dg::DGMulti, cache, parabolic_cache) + + reset_du!(du, dg) + + # compute volume contributions to divergence + @threaded for e in eachelement(mesh, dg) + for i in eachdim(mesh), j in eachdim(mesh) + dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine + StructArrays.foreachfield(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), + view(du, :, e), view(u_flux[i], :, e)) + end end + + # prolong2interfaces!(cache.u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) + + # # compute fluxes at interfaces + # @unpack u_face_values, flux_face_values = cache + # @unpack mapM, mapP, Jf = mesh.md + # @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) + # idM, idP = mapM[face_node_index], mapP[face_node_index] + # uM = u_face_values[idM] + # # compute flux if node is not a boundary node + # if idM != idP + # uP = u_face_values[idP] + # flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? + # end + # end + + # calc_divergence_boundary_integral!(du, u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) + + invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0) end -# function calc_divergence!(du, u::StructArray, mesh::DGMultiMesh, -# equations::AbstractParabolicEquations, -# boundary_conditions, dg::DGMulti, cache, parabolic_cache) -# calc_divergence_volume_integral!(du, u, mesh, equations, dg, cache, parabolic_cache) -# calc_divergence_surface_integral!(du, u, mesh, equations, dg, cache, parabolic_cache) -# calc_divergence_boundary_integral!(du, u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) -# end - -# # assumptions: parabolic terms are of the form div(f(u, grad(u))) and -# # will be discretized in first order form -# # - compute grad(u) -# # - compute f(u, grad(u)) -# # - compute div(u) -# # boundary conditions will be applied to both grad(u) and div(u). -# function rhs!(du, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, -# initial_condition, boundary_conditions, source_terms, -# dg::DGMulti, cache, parabolic_cache) -# @unpack u_transformed, grad_u, viscous_flux = parabolic_cache -# transform_variables!(u_transformed, u, equations) -# calc_gradient!(grad_u, u_transformed, mesh, equations, -# boundary_conditions, dg, cache, parabolic_cache) -# calc_viscous_fluxes!(viscous_flux, u_transformed, grad_u, -# mesh, equations, dg, cache, parabolic_cache) -# calc_divergence!(du, grad_u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) -# end +# assumptions: parabolic terms are of the form div(f(u, grad(u))) and +# will be discretized in first order form +# - compute grad(u) +# - compute f(u, grad(u)) +# - compute div(u) +# boundary conditions will be applied to both grad(u) and div(u). +function rhs!(du, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, + initial_condition, boundary_conditions, source_terms, + dg::DGMulti, cache, parabolic_cache) + + @unpack u_transformed, grad_u, viscous_flux = parabolic_cache + transform_variables!(u_transformed, u, equations) + + calc_gradient!(grad_u, u_transformed, mesh, equations, + boundary_conditions, dg, cache, parabolic_cache) + + calc_viscous_fluxes!(viscous_flux, u_transformed, grad_u, + mesh, equations, dg, cache, parabolic_cache) + + calc_divergence!(du, grad_u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) + +end From 681c81805dbf1466f894a15821a1f87b51d69b95 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 08:02:03 -0500 Subject: [PATCH 21/86] adding cache --- .../semidiscretization_hyperbolic_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index d218d306f12..8cfe1ea688c 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -216,7 +216,7 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) end function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) - @unpack mesh, parabolic_equations, initial_condition, parabolic_boundary_conditions, source_terms, solver, parabolic_cache = semi + @unpack mesh, parabolic_equations, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache = semi u = wrap_array(u_ode, mesh, parabolic_equations, solver, parabolic_cache) du = wrap_array(du_ode, mesh, parabolic_equations, solver, parabolic_cache) From 421a3177aa92c7a3bdafbda82fbb17add3f54689 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 13:26:17 -0500 Subject: [PATCH 22/86] remove test --- test/test_dgmulti_2d.jl | 66 ++++++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 511d5fdda5c..1c95a20e6a2 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -13,39 +13,39 @@ isdir(outdir) && rm(outdir, recursive=true) @testset "DGMulti 2D" begin - @trixi_testset "SemidiscretizationHyperbolicParabolic" begin - dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), - surface_integral = SurfaceIntegralWeakForm(flux_central), - volume_integral = VolumeIntegralWeakForm()) - mesh = DGMultiMesh(dg, cells_per_dimension=(4, 4), - periodicity=(true, true)) - - initial_condition = initial_condition_gauss - - # split linear advection with vector = (1, 1) into two terms with - # advection vectors (1, 0) and (0, 1). - equations = LinearScalarAdvectionEquation2D(1.0, 1.0) - equation1 = LinearScalarAdvectionEquation2D(1.0, 0.0) - equation2 = LinearScalarAdvectionEquation2D(0.0, 1.0) # our "parabolic" equation - - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) - semi_hyp_par = SemidiscretizationHyperbolicParabolic(mesh, equation1, equation2, initial_condition, dg) - - ode = semidiscretize(semi, (0.0, 0.01)) - ode_hyp_par = semidiscretize(semi_hyp_par, (0.0, 0.01)) - - du = similar(ode.u0) - du_hyp = similar(ode_hyp_par.u0) - du_par = similar(ode_hyp_par.u0) - Trixi.rhs!(du, ode.u0, semi, 0.0) - Trixi.rhs!(du_hyp, ode_hyp_par.u0, semi_hyp_par, 0.0) - Trixi.rhs_parabolic!(du_par, ode_hyp_par.u0, semi_hyp_par, 0.0) - - # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs - # for `equation` with SemidiscretizationHyperbolic. - du_diff = du .- (du_hyp .+ du_par) - @test sum(abs.(getindex.(du_diff, 1))) < 10 * length(du) * eps() - end + # @trixi_testset "SemidiscretizationHyperbolicParabolic" begin + # dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), + # surface_integral = SurfaceIntegralWeakForm(flux_central), + # volume_integral = VolumeIntegralWeakForm()) + # mesh = DGMultiMesh(dg, cells_per_dimension=(4, 4), + # periodicity=(true, true)) + + # initial_condition = initial_condition_gauss + + # # split linear advection with vector = (1, 1) into two terms with + # # advection vectors (1, 0) and (0, 1). + # equations = LinearScalarAdvectionEquation2D(1.0, 1.0) + # equation1 = LinearScalarAdvectionEquation2D(1.0, 0.0) + # equation2 = LinearScalarAdvectionEquation2D(0.0, 1.0) # our "parabolic" equation + + # semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) + # semi_hyp_par = SemidiscretizationHyperbolicParabolic(mesh, equation1, equation2, initial_condition, dg) + + # ode = semidiscretize(semi, (0.0, 0.01)) + # ode_hyp_par = semidiscretize(semi_hyp_par, (0.0, 0.01)) + + # du = similar(ode.u0) + # du_hyp = similar(ode_hyp_par.u0) + # du_par = similar(ode_hyp_par.u0) + # Trixi.rhs!(du, ode.u0, semi, 0.0) + # Trixi.rhs!(du_hyp, ode_hyp_par.u0, semi_hyp_par, 0.0) + # Trixi.rhs_parabolic!(du_par, ode_hyp_par.u0, semi_hyp_par, 0.0) + + # # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs + # # for `equation` with SemidiscretizationHyperbolic. + # du_diff = du .- (du_hyp .+ du_par) + # @test sum(abs.(getindex.(du_diff, 1))) < 10 * length(du) * eps() + # end @trixi_testset "elixir_euler_weakform.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), From d0d708336d8c8d8a89c1b1161aff15507ad6e795 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 14:06:20 -0500 Subject: [PATCH 23/86] update diffusion test --- test/test_dgmulti_2d.jl | 72 ++++++++++++++++++++++------------------- 1 file changed, 39 insertions(+), 33 deletions(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 1c95a20e6a2..5d50fb6749f 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -13,39 +13,45 @@ isdir(outdir) && rm(outdir, recursive=true) @testset "DGMulti 2D" begin - # @trixi_testset "SemidiscretizationHyperbolicParabolic" begin - # dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), - # surface_integral = SurfaceIntegralWeakForm(flux_central), - # volume_integral = VolumeIntegralWeakForm()) - # mesh = DGMultiMesh(dg, cells_per_dimension=(4, 4), - # periodicity=(true, true)) - - # initial_condition = initial_condition_gauss - - # # split linear advection with vector = (1, 1) into two terms with - # # advection vectors (1, 0) and (0, 1). - # equations = LinearScalarAdvectionEquation2D(1.0, 1.0) - # equation1 = LinearScalarAdvectionEquation2D(1.0, 0.0) - # equation2 = LinearScalarAdvectionEquation2D(0.0, 1.0) # our "parabolic" equation - - # semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) - # semi_hyp_par = SemidiscretizationHyperbolicParabolic(mesh, equation1, equation2, initial_condition, dg) - - # ode = semidiscretize(semi, (0.0, 0.01)) - # ode_hyp_par = semidiscretize(semi_hyp_par, (0.0, 0.01)) - - # du = similar(ode.u0) - # du_hyp = similar(ode_hyp_par.u0) - # du_par = similar(ode_hyp_par.u0) - # Trixi.rhs!(du, ode.u0, semi, 0.0) - # Trixi.rhs!(du_hyp, ode_hyp_par.u0, semi_hyp_par, 0.0) - # Trixi.rhs_parabolic!(du_par, ode_hyp_par.u0, semi_hyp_par, 0.0) - - # # Test that the sum of the rhs for `equations1` and `equations2` is the same as the rhs - # # for `equation` with SemidiscretizationHyperbolic. - # du_diff = du .- (du_hyp .+ du_par) - # @test sum(abs.(getindex.(du_diff, 1))) < 10 * length(du) * eps() - # end + @trixi_testset "SemidiscretizationHyperbolicParabolic" begin + using Test + + dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_central), + volume_integral = VolumeIntegralWeakForm()) + mesh = DGMultiMesh(dg, cells_per_dimension=(2, 2)) + + # test with polynomial initial condition x^2 * y + initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2]) + + equations = LinearScalarAdvectionEquation2D(1.0, 1.0) + equations_parabolic = Trixi.ScalarDiffusion2D(1.0) + + semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) + ode = semidiscretize(semi, (0.0, 0.01)) + + @unpack cache, parabolic_cache, parabolic_equations = semi + @unpack u_grad = parabolic_cache + for dim in 1:length(u_grad) + fill!(u_grad[dim], zero(eltype(u_grad[dim]))) + end + + # pass in BoundaryConditionPeriodic() to fake a "do nothing" BC + Trixi.calc_gradient!(u_grad, ode.u0, mesh, parabolic_equations, Trixi.BoundaryConditionPeriodic(), dg, cache, parabolic_cache) + @unpack x, y = mesh.md + @test norm(getindex.(u_grad[1], 1) - 2 * x .* y) < 1e3 * eps() + @test norm(getindex.(u_grad[2], 1) - x.^2) < 1e3 * eps() + + u_flux = similar.(u_grad) + Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, parabolic_equations, dg, cache, parabolic_cache) + @test u_flux[1] ≈ u_grad[1] + @test u_flux[2] ≈ u_grad[2] + + du = similar(ode.u0) + # pass in BoundaryConditionPeriodic() to fake a "do nothing" BC + Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, parabolic_equations, Trixi.BoundaryConditionPeriodic(), dg, cache, parabolic_cache) + @test norm(-getindex.(du, 1) - 2 * y) < 1e3 * eps() + end @trixi_testset "elixir_euler_weakform.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), From a9fde8a86e397699512b2f72b32faf1cad2c2e8f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 14:06:33 -0500 Subject: [PATCH 24/86] update diffusion rhs --- src/solvers/dgmulti/diffusion.jl | 57 ++++++++++++++++++++++---------- 1 file changed, 39 insertions(+), 18 deletions(-) diff --git a/src/solvers/dgmulti/diffusion.jl b/src/solvers/dgmulti/diffusion.jl index 45abff9810c..f3d1fdea2a1 100644 --- a/src/solvers/dgmulti/diffusion.jl +++ b/src/solvers/dgmulti/diffusion.jl @@ -43,8 +43,7 @@ end # interpolates from solution coefficients to face quadrature points function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, surface_integral, dg::DGMulti, cache) - rd = dg.basis - apply_to_each_field(mul_by!(rd.Vf), u_face_values, u) + apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) end function calc_gradient_surface_integral(u_grad, u, flux_face_values, @@ -111,6 +110,12 @@ function calc_gradient_boundary_integral!(du, u, mesh, equations, ::BoundaryCond return nothing end +# do nothing for periodic domains +function calc_divergence_boundary_integral!(du, u, u_flux, mesh, equations, ::BoundaryConditionPeriodic, + dg, cache, parabolic_cache) + return nothing +end + function calc_viscous_fluxes!(u_flux, u, u_grad, mesh::DGMultiMesh, equations::AbstractParabolicEquations, dg::DGMulti, cache, parabolic_cache) @@ -154,24 +159,40 @@ function calc_divergence!(du, u::StructArray, u_flux, mesh::DGMultiMesh, end end - # prolong2interfaces!(cache.u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) + # interpolates from solution coefficients to face quadrature points + @unpack grad_u_face_values = parabolic_cache + for dim in eachdim(mesh) + prolong2interfaces!(grad_u_face_values[dim], u_flux[dim], mesh, equations, + dg.surface_integral, dg, cache) + end + + # compute fluxes at interfaces + @unpack grad_u_face_values, flux_face_values = parabolic_cache + @unpack mapM, mapP, nxyzJ, Jf = mesh.md + @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) + idM, idP = mapM[face_node_index], mapP[face_node_index] + + # compute f(u, ∇u) ⋅ n + flux_face_value = zero(eltype(flux_face_values)) + for dim in eachdim(mesh) + uM = grad_u_face_values[dim][idM] + # compute flux if node is not a boundary node + if idM != idP + uP = grad_u_face_values[dim][idP] + # TODO: use strong/weak formulation? + flux_face_value = flux_face_value + 0.5 * (uP - uM) * nxyzJ[dim][face_node_index] + end + end + flux_face_values[idM] = flux_face_value + end - # # compute fluxes at interfaces - # @unpack u_face_values, flux_face_values = cache - # @unpack mapM, mapP, Jf = mesh.md - # @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) - # idM, idP = mapM[face_node_index], mapP[face_node_index] - # uM = u_face_values[idM] - # # compute flux if node is not a boundary node - # if idM != idP - # uP = u_face_values[idP] - # flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? - # end - # end + # surface contributions + apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, flux_face_values) - # calc_divergence_boundary_integral!(du, u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) + calc_divergence_boundary_integral!(du, u, u_flux, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) - invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0) + # multiply by -1 when adding to the RHS + invert_jacobian!(du, mesh, equations, dg, cache; scaling=-1.0) end # assumptions: parabolic terms are of the form div(f(u, grad(u))) and @@ -193,6 +214,6 @@ function rhs!(du, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, calc_viscous_fluxes!(viscous_flux, u_transformed, grad_u, mesh, equations, dg, cache, parabolic_cache) - calc_divergence!(du, grad_u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) + calc_divergence!(du, u_transformed, grad_u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) end From 827defbb015addf94c9177b5e4afa97d4b64252b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 15:15:25 -0500 Subject: [PATCH 25/86] moving routines around --- src/Trixi.jl | 3 +-- src/equations/equations.jl | 2 ++ src/equations/scalar_diffusion_2d.jl | 11 ++++++++++ src/solvers/dgmulti.jl | 3 +++ .../dgmulti/{diffusion.jl => dg_parabolic.jl} | 22 ++++--------------- 5 files changed, 21 insertions(+), 20 deletions(-) create mode 100644 src/equations/scalar_diffusion_2d.jl rename src/solvers/dgmulti/{diffusion.jl => dg_parabolic.jl} (92%) diff --git a/src/Trixi.jl b/src/Trixi.jl index b7aa4b6e2f6..b8a3d75bb29 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -114,8 +114,7 @@ include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") include("time_integration/time_integration.jl") -# TODO: move this around -include("solvers/dgmulti/diffusion.jl") +include("solvers/dgmulti/dg_parabolic.jl") # `trixi_include` and special elixirs such as `convergence_test` include("auxiliary/special_elixirs.jl") diff --git a/src/equations/equations.jl b/src/equations/equations.jl index dc3bfffba30..c43bb88fa9d 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -329,5 +329,7 @@ include("lattice_boltzmann_3d.jl") abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("acoustic_perturbation_2d.jl") +abstract type AbstractParabolicEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +include("scalar_diffusion_2d.jl") end # @muladd diff --git a/src/equations/scalar_diffusion_2d.jl b/src/equations/scalar_diffusion_2d.jl new file mode 100644 index 00000000000..1b980a48de0 --- /dev/null +++ b/src/equations/scalar_diffusion_2d.jl @@ -0,0 +1,11 @@ + +struct ScalarDiffusion2D{T} <: AbstractParabolicEquations{2, 1} + diffusivity::T +end + +# no orientation here since the flux is vector-valued +function flux(u, grad_u, equations::ScalarDiffusion2D) + dudx, dudy = grad_u + return equations.diffusivity * dudx, equations.diffusivity * dudy +end + diff --git a/src/solvers/dgmulti.jl b/src/solvers/dgmulti.jl index d5c8f12d395..318a11b678e 100644 --- a/src/solvers/dgmulti.jl +++ b/src/solvers/dgmulti.jl @@ -9,3 +9,6 @@ include("dgmulti/sbp.jl") # specialization of DGMulti to specific equations include("dgmulti/flux_differencing_compressible_euler.jl") + +# parabolic terms for DGMulti solvers +include("dgmulti/dg_parabolic.jl") \ No newline at end of file diff --git a/src/solvers/dgmulti/diffusion.jl b/src/solvers/dgmulti/dg_parabolic.jl similarity index 92% rename from src/solvers/dgmulti/diffusion.jl rename to src/solvers/dgmulti/dg_parabolic.jl index f3d1fdea2a1..073d6079590 100644 --- a/src/solvers/dgmulti/diffusion.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -1,21 +1,8 @@ - -abstract type AbstractParabolicEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end - -struct ScalarDiffusion2D{T} <: AbstractParabolicEquations{2, 1} - diffusivity::T -end - -# no orientation here since the flux is vector-valued -function flux(u, grad_u, equations::ScalarDiffusion2D) - dudx, dudy = grad_u - return equations.diffusivity * dudx, equations.diffusivity * dudy -end - function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, dg::DGMultiWeakForm, RealT, uEltype) nvars = nvariables(equations) - # u_parabolic stores "transformed" variables for + # u_parabolic stores "transformed" variables for computing the gradient @unpack md = mesh u_transformed = allocate_nested_array(uEltype, nvars, size(md.xq), dg) u_grad = ntuple(_ -> similar(u_transformed), ndims(mesh)) @@ -26,10 +13,9 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, local_viscous_flux = ntuple(_ -> similar(u_transformed, size(md.xq, 1)), ndims(mesh)) flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) - return (; u_transformed, u_grad, - viscous_flux, local_viscous_flux, - u_face_values, grad_u_face_values, - flux_face_values, local_flux_values=similar(flux_face_values[:, 1])) + return (; u_transformed, u_grad, viscous_flux, local_viscous_flux, + u_face_values, grad_u_face_values, flux_face_values, + local_flux_values=similar(flux_face_values[:, 1])) end # Transform variables prior to taking the gradient. Defaults to doing nothing. From a0a4b8bbfe117b13d8c93823561184825c80138c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 15:39:45 -0500 Subject: [PATCH 26/86] fix multiple includes --- src/Trixi.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index b8a3d75bb29..c7a3de79d7d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -114,8 +114,6 @@ include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") include("time_integration/time_integration.jl") -include("solvers/dgmulti/dg_parabolic.jl") - # `trixi_include` and special elixirs such as `convergence_test` include("auxiliary/special_elixirs.jl") From f4f181e8f3ecf953447657a03da08acf2aa669c7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 19:49:09 -0500 Subject: [PATCH 27/86] removing time dependence for now from viscous rhs --- .../semidiscretization_hyperbolic_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 8cfe1ea688c..1fe667b5284 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -223,7 +223,7 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() - @trixi_timeit timer() "parabolic rhs!" rhs!(du, u, t, mesh, parabolic_equations, initial_condition, + @trixi_timeit timer() "parabolic rhs!" rhs!(du, u, mesh, parabolic_equations, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache) runtime = time_ns() - time_start From 333bed3f64c17909d94e484b5ce625b6e97a87be Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 19:49:18 -0500 Subject: [PATCH 28/86] export ScalarDiffusion2D --- src/Trixi.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Trixi.jl b/src/Trixi.jl index c7a3de79d7d..c435dc20dfc 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -133,6 +133,8 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D +export ScalarDiffusion2D + export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, From 2a1a6b402c641a8c9e3d12104df1806be61d20b8 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 19:49:24 -0500 Subject: [PATCH 29/86] fix test --- test/test_dgmulti_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 5d50fb6749f..32f24fef774 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -50,7 +50,7 @@ isdir(outdir) && rm(outdir, recursive=true) du = similar(ode.u0) # pass in BoundaryConditionPeriodic() to fake a "do nothing" BC Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, parabolic_equations, Trixi.BoundaryConditionPeriodic(), dg, cache, parabolic_cache) - @test norm(-getindex.(du, 1) - 2 * y) < 1e3 * eps() + @test norm(getindex.(du, 1) - 2 * y) < 1e3 * eps() end @trixi_testset "elixir_euler_weakform.jl" begin From a909b5416afcdd86385f3895b4d93f94b7e5fcf2 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 19:49:44 -0500 Subject: [PATCH 30/86] working periodic advection RHS --- .../dgmulti_2d/elixir_advection_diffusion.jl | 30 +++++ src/solvers/dgmulti/dg_parabolic.jl | 105 +++++++++++------- 2 files changed, 96 insertions(+), 39 deletions(-) create mode 100644 examples/dgmulti_2d/elixir_advection_diffusion.jl diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl new file mode 100644 index 00000000000..c1b5f559054 --- /dev/null +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -0,0 +1,30 @@ + +using Trixi, OrdinaryDiffEq + +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +equations = LinearScalarAdvectionEquation2D(1.0, 1.0) +parabolic_equations = ScalarDiffusion2D(1e-2) +initial_condition = (x, t, equations) -> SVector(exp(-100 * (x[1]^2 + x[2]^2))) + +mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) +semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg) +# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +tol = 1e-6 +sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 073d6079590..09347567080 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -4,23 +4,24 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, # u_parabolic stores "transformed" variables for computing the gradient @unpack md = mesh - u_transformed = allocate_nested_array(uEltype, nvars, size(md.xq), dg) + u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) u_grad = ntuple(_ -> similar(u_transformed), ndims(mesh)) u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) grad_u_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) viscous_flux = similar.(u_grad) - local_viscous_flux = ntuple(_ -> similar(u_transformed, size(md.xq, 1)), ndims(mesh)) - flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) - return (; u_transformed, u_grad, viscous_flux, local_viscous_flux, - u_face_values, grad_u_face_values, flux_face_values, - local_flux_values=similar(flux_face_values[:, 1])) + + local_viscous_flux_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] + local_flux_values_threaded = [similar(flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] + + return (; u_transformed, u_grad, viscous_flux, u_face_values, grad_u_face_values, flux_face_values, + local_viscous_flux_threaded, local_flux_values_threaded) end # Transform variables prior to taking the gradient. Defaults to doing nothing. # TODO: can we avoid copying data? -function transform_variables(u_transformed, u, equations) +function transform_variables!(u_transformed, u, equations) @threaded for i in eachindex(u) u_transformed[i] = u[i] end @@ -35,8 +36,9 @@ end function calc_gradient_surface_integral(u_grad, u, flux_face_values, mesh, equations::AbstractParabolicEquations, dg::DGMulti, cache, parabolic_cache) - @unpack local_flux_values = parabolic_cache + @unpack local_flux_values_threaded = parabolic_cache @threaded for e in eachelement(mesh, dg) + local_flux_values = local_flux_values_threaded[Threads.threadid()] for dim in eachdim(mesh) for i in eachindex(local_flux_values) # compute [u] * (nx, ny, nz) @@ -59,17 +61,18 @@ function calc_gradient!(u_grad, u::StructArray, mesh::DGMultiMesh, @threaded for e in eachelement(mesh, dg) for i in eachdim(mesh), j in eachdim(mesh) dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine - StructArrays.foreachfield(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), - view(u_grad[i], :, e), view(u, :, e)) + apply_to_each_field(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), + view(u_grad[i], :, e), view(u, :, e)) end end - prolong2interfaces!(cache.u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) + prolong2interfaces!(cache.u_face_values, u, + mesh, equations, dg.surface_integral, dg, cache) # compute fluxes at interfaces @unpack u_face_values, flux_face_values = cache @unpack mapM, mapP, Jf = mesh.md - @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) + @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] # compute flux if node is not a boundary node @@ -91,46 +94,66 @@ function calc_gradient!(u_grad, u::StructArray, mesh::DGMultiMesh, end # do nothing for periodic domains -function calc_gradient_boundary_integral!(du, u, mesh, equations, ::BoundaryConditionPeriodic, - dg, cache, parabolic_cache) +function calc_gradient_boundary_integral!(du, u, mesh, equations::AbstractParabolicEquations, + ::BoundaryConditionPeriodic, dg, cache, parabolic_cache) return nothing end # do nothing for periodic domains -function calc_divergence_boundary_integral!(du, u, u_flux, mesh, equations, ::BoundaryConditionPeriodic, - dg, cache, parabolic_cache) +function calc_divergence_boundary_integral!(du, u, viscous_flux, mesh, equations::AbstractParabolicEquations, + ::BoundaryConditionPeriodic, dg, cache, parabolic_cache) return nothing end -function calc_viscous_fluxes!(u_flux, u, u_grad, mesh::DGMultiMesh, +function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, equations::AbstractParabolicEquations, dg::DGMulti, cache, parabolic_cache) - @unpack local_viscous_flux = parabolic_cache for dim in eachdim(mesh) - reset_du!(local_viscous_flux[dim], dg) + reset_du!(viscous_flux[dim], dg) end + # @threaded for i in eachindex(u) + # for dim in eachdim(mesh) + # viscous_flux[dim][i] = equations.diffusivity * u_grad[dim][i] + # end + # end + + @unpack local_viscous_flux_threaded = parabolic_cache + local_u_values_threaded = parabolic_cache.local_flux_values_threaded + @threaded for e in eachelement(mesh, dg) - # interpolate to quadrature points + + # reset local storage for each element + local_viscous_flux = local_viscous_flux_threaded[Threads.threadid()] + local_u_values = local_u_values_threaded[Threads.threadid()] + fill!(local_u_values, zero(eltype(local_u_values))) for dim in eachdim(mesh) - StructArrays.foreachfield(mul_by!(dg.basis.Vq), local_viscous_flux[dim], view(u_grad[dim], :, e)) + fill!(local_viscous_flux[dim], zero(eltype(local_viscous_flux[dim]))) + end + + # interpolate u and gradient to quadrature points, store in `local_viscous_flux` + apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) # TODO: can we avoid this when we don't need it? + for dim in eachdim(mesh) + apply_to_each_field(mul_by!(dg.basis.Vq), local_viscous_flux[dim], view(u_grad[dim], :, e)) end # compute viscous flux at quad points - for i in eachindex(local_viscous_flux) + for i in eachindex(local_u_values) + u_i = local_u_values[i] u_grad_i = getindex.(local_viscous_flux, i) # TODO: check if this allocates. Shouldn't for tuples or SVector... - setindex!.(local_viscous_flux, flux(u[i,e], u_grad_i, equations), i) + viscous_flux_i = flux(u_i, u_grad_i, equations) + setindex!.(local_viscous_flux, viscous_flux_i, i) end # project back to the DG approximation space for dim in eachdim(mesh) - StructArrays.foreachfield(mul_by!(dg.basis.Pq), view(u_flux[dim], :, e), local_viscous_flux[dim]) + apply_to_each_field(mul_by!(dg.basis.Pq), view(viscous_flux[dim], :, e), local_viscous_flux[dim]) end end end -function calc_divergence!(du, u::StructArray, u_flux, mesh::DGMultiMesh, +function calc_divergence!(du, u::StructArray, viscous_flux, mesh::DGMultiMesh, equations::AbstractParabolicEquations, boundary_conditions, dg::DGMulti, cache, parabolic_cache) @@ -140,21 +163,21 @@ function calc_divergence!(du, u::StructArray, u_flux, mesh::DGMultiMesh, @threaded for e in eachelement(mesh, dg) for i in eachdim(mesh), j in eachdim(mesh) dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine - StructArrays.foreachfield(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), - view(du, :, e), view(u_flux[i], :, e)) + apply_to_each_field(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), + view(du, :, e), view(viscous_flux[i], :, e)) end end # interpolates from solution coefficients to face quadrature points @unpack grad_u_face_values = parabolic_cache for dim in eachdim(mesh) - prolong2interfaces!(grad_u_face_values[dim], u_flux[dim], mesh, equations, + prolong2interfaces!(grad_u_face_values[dim], viscous_flux[dim], mesh, equations, dg.surface_integral, dg, cache) end # compute fluxes at interfaces @unpack grad_u_face_values, flux_face_values = parabolic_cache - @unpack mapM, mapP, nxyzJ, Jf = mesh.md + @unpack mapM, mapP, nxyzJ = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) idM, idP = mapM[face_node_index], mapP[face_node_index] @@ -175,10 +198,9 @@ function calc_divergence!(du, u::StructArray, u_flux, mesh::DGMultiMesh, # surface contributions apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, flux_face_values) - calc_divergence_boundary_integral!(du, u, u_flux, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) + calc_divergence_boundary_integral!(du, u, viscous_flux, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) - # multiply by -1 when adding to the RHS - invert_jacobian!(du, mesh, equations, dg, cache; scaling=-1.0) + invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0) end # assumptions: parabolic terms are of the form div(f(u, grad(u))) and @@ -187,19 +209,24 @@ end # - compute f(u, grad(u)) # - compute div(u) # boundary conditions will be applied to both grad(u) and div(u). -function rhs!(du, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, +function rhs!(du, u, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEquations, initial_condition, boundary_conditions, source_terms, dg::DGMulti, cache, parabolic_cache) - @unpack u_transformed, grad_u, viscous_flux = parabolic_cache - transform_variables!(u_transformed, u, equations) + reset_du!(du, dg) - calc_gradient!(grad_u, u_transformed, mesh, equations, + @unpack u_transformed, u_grad, viscous_flux = parabolic_cache + transform_variables!(u_transformed, u, parabolic_equations) + + calc_gradient!(u_grad, u_transformed, mesh, parabolic_equations, boundary_conditions, dg, cache, parabolic_cache) - calc_viscous_fluxes!(viscous_flux, u_transformed, grad_u, - mesh, equations, dg, cache, parabolic_cache) + calc_viscous_fluxes!(viscous_flux, u_transformed, u_grad, + mesh, parabolic_equations, dg, cache, parabolic_cache) + + calc_divergence!(du, u_transformed, viscous_flux, mesh, parabolic_equations, + boundary_conditions, dg, cache, parabolic_cache) - calc_divergence!(du, u_transformed, grad_u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) + return nothing end From 3316faba77d46a074abb5eec7866cc138ea547cc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 19:58:09 -0500 Subject: [PATCH 31/86] cleanup and comments --- examples/dgmulti_2d/elixir_advection_diffusion.jl | 9 ++++++--- src/solvers/dgmulti/dg_parabolic.jl | 8 ++++---- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index c1b5f559054..565f6c9d4cb 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -7,13 +7,16 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( equations = LinearScalarAdvectionEquation2D(1.0, 1.0) parabolic_equations = ScalarDiffusion2D(1e-2) -initial_condition = (x, t, equations) -> SVector(exp(-100 * (x[1]^2 + x[2]^2))) + +function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) + return SVector(exp(-100 * (x[1]^2 + x[2]^2))) +end +initial_condition = initial_condition_sharp_gaussian mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg) -# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) -tspan = (0.0, 0.5) +tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 09347567080..5967fff01e6 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -204,10 +204,10 @@ function calc_divergence!(du, u::StructArray, viscous_flux, mesh::DGMultiMesh, end # assumptions: parabolic terms are of the form div(f(u, grad(u))) and -# will be discretized in first order form -# - compute grad(u) -# - compute f(u, grad(u)) -# - compute div(u) +# will be discretized first order form as follows: +# 1. compute grad(u) +# 2. compute f(u, grad(u)) +# 3. compute div(u) # boundary conditions will be applied to both grad(u) and div(u). function rhs!(du, u, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEquations, initial_condition, boundary_conditions, source_terms, From 3fe88d7c3c5549600b1ab94eb462c129efd89b5b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 21:20:02 -0500 Subject: [PATCH 32/86] making periodic advection-diffusion example --- .../elixir_advection_diffusion_periodic.jl | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl new file mode 100644 index 00000000000..565f6c9d4cb --- /dev/null +++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl @@ -0,0 +1,33 @@ + +using Trixi, OrdinaryDiffEq + +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +equations = LinearScalarAdvectionEquation2D(1.0, 1.0) +parabolic_equations = ScalarDiffusion2D(1e-2) + +function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) + return SVector(exp(-100 * (x[1]^2 + x[2]^2))) +end +initial_condition = initial_condition_sharp_gaussian + +mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) +semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg) + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +tol = 1e-6 +sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, callback=callbacks) +summary_callback() # print the timer summary From 341de95dd2155ef24b19ef7e2dffdc64031db371 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 21:20:12 -0500 Subject: [PATCH 33/86] trim whitespace --- src/equations/scalar_diffusion_2d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/equations/scalar_diffusion_2d.jl b/src/equations/scalar_diffusion_2d.jl index 1b980a48de0..26cef6692b2 100644 --- a/src/equations/scalar_diffusion_2d.jl +++ b/src/equations/scalar_diffusion_2d.jl @@ -1,4 +1,3 @@ - struct ScalarDiffusion2D{T} <: AbstractParabolicEquations{2, 1} diffusivity::T end From fbf27c35ab91d97dd82642016901eac67f8c7b8a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 22:54:41 -0500 Subject: [PATCH 34/86] adding t back in to parabolic rhs add t back in --- .../semidiscretization_hyperbolic_parabolic.jl | 2 +- src/solvers/dgmulti/dg_parabolic.jl | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 1fe667b5284..8cfe1ea688c 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -223,7 +223,7 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() - @trixi_timeit timer() "parabolic rhs!" rhs!(du, u, mesh, parabolic_equations, initial_condition, + @trixi_timeit timer() "parabolic rhs!" rhs!(du, u, t, mesh, parabolic_equations, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache) runtime = time_ns() - time_start diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 5967fff01e6..a6c5201a725 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -49,7 +49,7 @@ function calc_gradient_surface_integral(u_grad, u, flux_face_values, end end -function calc_gradient!(u_grad, u::StructArray, mesh::DGMultiMesh, +function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, equations::AbstractParabolicEquations, boundary_conditions, dg::DGMulti, cache, parabolic_cache) @@ -153,7 +153,7 @@ function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, end end -function calc_divergence!(du, u::StructArray, viscous_flux, mesh::DGMultiMesh, +function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh, equations::AbstractParabolicEquations, boundary_conditions, dg::DGMulti, cache, parabolic_cache) @@ -209,7 +209,7 @@ end # 2. compute f(u, grad(u)) # 3. compute div(u) # boundary conditions will be applied to both grad(u) and div(u). -function rhs!(du, u, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEquations, +function rhs!(du, u, t, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEquations, initial_condition, boundary_conditions, source_terms, dg::DGMulti, cache, parabolic_cache) @@ -218,13 +218,13 @@ function rhs!(du, u, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEq @unpack u_transformed, u_grad, viscous_flux = parabolic_cache transform_variables!(u_transformed, u, parabolic_equations) - calc_gradient!(u_grad, u_transformed, mesh, parabolic_equations, + calc_gradient!(u_grad, u_transformed, t, mesh, parabolic_equations, boundary_conditions, dg, cache, parabolic_cache) calc_viscous_fluxes!(viscous_flux, u_transformed, u_grad, mesh, parabolic_equations, dg, cache, parabolic_cache) - calc_divergence!(du, u_transformed, viscous_flux, mesh, parabolic_equations, + calc_divergence!(du, u_transformed, t, viscous_flux, mesh, parabolic_equations, boundary_conditions, dg, cache, parabolic_cache) return nothing From b4bf34af684c6bcb965ac0bfd5cd39c9127cb5ff Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 22:57:15 -0500 Subject: [PATCH 35/86] refactoring boundary condition code --- src/solvers/dgmulti/dg_parabolic.jl | 43 +++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 9 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index a6c5201a725..388b5594cbc 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -78,30 +78,54 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, # compute flux if node is not a boundary node if idM != idP uP = u_face_values[idP] + # {u} - u = 0.5 * (uP - uM) flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? end end + calc_boundary_flux!(flux_face_values, u, t, Gradient(), boundary_conditions, + mesh, equations, dg, cache, parabolic_cache) + # compute surface contributions calc_gradient_surface_integral(u_grad, u, flux_face_values, mesh, equations, dg, cache, parabolic_cache) - calc_gradient_boundary_integral!(u_grad, u, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) - for dim in eachdim(mesh) invert_jacobian!(u_grad[dim], mesh, equations, dg, cache; scaling=1.0) end end +# operator types used for dispatch on boundary fluxes +struct Gradient end +struct Divergence end + # do nothing for periodic domains -function calc_gradient_boundary_integral!(du, u, mesh, equations::AbstractParabolicEquations, - ::BoundaryConditionPeriodic, dg, cache, parabolic_cache) +function calc_boundary_flux!(flux, u, t, operator_type, ::BoundaryConditionPeriodic, + mesh, equations::AbstractParabolicEquations, dg::DGMulti, + cache, parabolic_cache) return nothing end -# do nothing for periodic domains -function calc_divergence_boundary_integral!(du, u, viscous_flux, mesh, equations::AbstractParabolicEquations, - ::BoundaryConditionPeriodic, dg, cache, parabolic_cache) +# "lispy tuple programming" instead of for loop for type stability +function calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions, + mesh, equations, dg::DGMulti, cache, parabolic_cache) + # peel off first boundary condition + calc_single_boundary_flux!(flux, u, t, operator_type, first(boundary_conditions), first(keys(boundary_conditions)), + mesh, equations, dg, cache, parabolic_cache) + + # recurse on the remainder of the boundary conditions + calc_boundary_flux!(flux, u, t, operator_type, Base.tail(boundary_conditions), + mesh, equations, dg, cache, parabolic_cache) +end + +# terminate recursion +calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{(),Tuple{}}, + mesh, equations, dg::DGMulti, cache, parabolic_cache) = nothing + +# TODO: finish +function calc_single_boundary_flux!(flux, u, t, operator_type, boundary_condition, + boundary_key, mesh, equations, dg::DGMulti, + cache, parabolic_cache) return nothing end @@ -195,11 +219,12 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh flux_face_values[idM] = flux_face_value end + calc_boundary_flux!(flux_face_values, u, t, Divergence(), boundary_conditions, + mesh, equations, dg, cache, parabolic_cache) + # surface contributions apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, flux_face_values) - calc_divergence_boundary_integral!(du, u, viscous_flux, mesh, equations, boundary_conditions, dg, cache, parabolic_cache) - invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0) end From 5f35752dbb75ba577f638c4a7a57a095bd5bc441 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 22:57:45 -0500 Subject: [PATCH 36/86] update elixir with BCs --- .../dgmulti_2d/elixir_advection_diffusion.jl | 37 +++++++++++++------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index 565f6c9d4cb..557a4f0d671 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -6,28 +6,41 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) -parabolic_equations = ScalarDiffusion2D(1e-2) +parabolic_equations = ScalarDiffusion2D(1e-3) -function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) - return SVector(exp(-100 * (x[1]^2 + x[2]^2))) -end -initial_condition = initial_condition_sharp_gaussian +initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) +initial_condition = initial_condition_zero -mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) -semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg) +# example where we tag two separate boundary segments of the mesh +inflow(x, tol=50*eps()) = abs(x[1]+1) inflow, :outflow => rest_of_boundary) -tspan = (0.0, 2.0) +mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16), is_on_boundary=is_on_boundary) + +# define inviscid boundary conditions +boundary_condition_inflow = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)) +boundary_conditions = (; :inflow => boundary_condition_inflow) + +# define viscous boundary conditions +boundary_condition_wall = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) +parabolic_boundary_conditions = (; boundary_condition_wall) + +semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg; + boundary_conditions)#, parabolic_boundary_conditions) + +tspan = (0.0, .25) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() alive_callback = AliveCallback(alive_interval=10) -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) -callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) +callbacks = CallbackSet(summary_callback, alive_callback) ############################################################################### # run the simulation tol = 1e-6 -sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, callback=callbacks) +tsave = LinRange(tspan..., 10) +sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, + saveat=tsave, callback=callbacks) summary_callback() # print the timer summary From 00cdd2b10a52603893c22605510d53ca21a2a627 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 22:59:39 -0500 Subject: [PATCH 37/86] change default DGMulti interface flux behavior same as in https://github.com/trixi-framework/Trixi.jl/pull/1139 --- src/solvers/dgmulti/dg.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 9952c84b989..f4bbcc4d6e1 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -208,13 +208,9 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm, # inner (idM -> minus) and outer (idP -> plus) indices idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] - - # compute flux if node is not a boundary node - if idM != idP - uP = u_face_values[idP] - normal = SVector{NDIMS}(getindex.(nxyzJ, idM)) / Jf[idM] - flux_face_values[idM] = surface_flux(uM, uP, normal, equations) * Jf[idM] - end + uP = u_face_values[idP] + normal = SVector{NDIMS}(getindex.(nxyzJ, idM)) / Jf[idM] + flux_face_values[idM] = surface_flux(uM, uP, normal, equations) * Jf[idM] end end From 2083bb282595743b969bf2ede913174bf31a7ab4 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 23:54:29 -0500 Subject: [PATCH 38/86] cleanup --- src/solvers/dgmulti/dg_parabolic.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 388b5594cbc..4bab291aef5 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -137,12 +137,6 @@ function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, reset_du!(viscous_flux[dim], dg) end - # @threaded for i in eachindex(u) - # for dim in eachdim(mesh) - # viscous_flux[dim][i] = equations.diffusivity * u_grad[dim][i] - # end - # end - @unpack local_viscous_flux_threaded = parabolic_cache local_u_values_threaded = parabolic_cache.local_flux_values_threaded From 695e77923f2189ace008e15000f4dca707ad9ebf Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 23:56:39 -0500 Subject: [PATCH 39/86] renaming variables --- src/solvers/dgmulti/dg_parabolic.jl | 55 ++++++++++++++--------------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 4bab291aef5..2c9927a5b3c 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -10,16 +10,19 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, grad_u_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) viscous_flux = similar.(u_grad) - flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) + div_flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) local_viscous_flux_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] - local_flux_values_threaded = [similar(flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] + local_flux_values_threaded = [similar(div_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] - return (; u_transformed, u_grad, viscous_flux, u_face_values, grad_u_face_values, flux_face_values, + return (; u_transformed, u_grad, viscous_flux, + weak_differentiation_matrices, + u_face_values, grad_u_face_values, div_flux_face_values, local_viscous_flux_threaded, local_flux_values_threaded) end -# Transform variables prior to taking the gradient. Defaults to doing nothing. +# Transform solution variables prior to taking the gradient +# (e.g., conservative to primitive variables). Defaults to doing nothing. # TODO: can we avoid copying data? function transform_variables!(u_transformed, u, equations) @threaded for i in eachindex(u) @@ -33,7 +36,7 @@ function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::Abs apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) end -function calc_gradient_surface_integral(u_grad, u, flux_face_values, +function calc_gradient_surface_integral(u_grad, u, div_flux_face_values, mesh, equations::AbstractParabolicEquations, dg::DGMulti, cache, parabolic_cache) @unpack local_flux_values_threaded = parabolic_cache @@ -42,7 +45,7 @@ function calc_gradient_surface_integral(u_grad, u, flux_face_values, for dim in eachdim(mesh) for i in eachindex(local_flux_values) # compute [u] * (nx, ny, nz) - local_flux_values[i] = flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] + local_flux_values[i] = div_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] end apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(u_grad[dim], :, e), local_flux_values) end @@ -70,24 +73,20 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, mesh, equations, dg.surface_integral, dg, cache) # compute fluxes at interfaces - @unpack u_face_values, flux_face_values = cache + @unpack u_face_values, div_flux_face_values = parabolic_cache @unpack mapM, mapP, Jf = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] - # compute flux if node is not a boundary node - if idM != idP - uP = u_face_values[idP] - # {u} - u = 0.5 * (uP - uM) - flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? - end + uP = u_face_values[idP] + div_flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? end - calc_boundary_flux!(flux_face_values, u, t, Gradient(), boundary_conditions, + calc_boundary_flux!(div_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, mesh, equations, dg, cache, parabolic_cache) # compute surface contributions - calc_gradient_surface_integral(u_grad, u, flux_face_values, + calc_gradient_surface_integral(u_grad, u, div_flux_face_values, mesh, equations, dg, cache, parabolic_cache) for dim in eachdim(mesh) @@ -109,6 +108,7 @@ end # "lispy tuple programming" instead of for loop for type stability function calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions, mesh, equations, dg::DGMulti, cache, parabolic_cache) + # peel off first boundary condition calc_single_boundary_flux!(flux, u, t, operator_type, first(boundary_conditions), first(keys(boundary_conditions)), mesh, equations, dg, cache, parabolic_cache) @@ -187,37 +187,34 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh end # interpolates from solution coefficients to face quadrature points - @unpack grad_u_face_values = parabolic_cache + viscous_flux_face_values = parabolic_cache.grad_u_face_values for dim in eachdim(mesh) - prolong2interfaces!(grad_u_face_values[dim], viscous_flux[dim], mesh, equations, + prolong2interfaces!(viscous_flux_face_values[dim], viscous_flux[dim], mesh, equations, dg.surface_integral, dg, cache) end # compute fluxes at interfaces - @unpack grad_u_face_values, flux_face_values = parabolic_cache + @unpack div_flux_face_values = parabolic_cache @unpack mapM, mapP, nxyzJ = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) idM, idP = mapM[face_node_index], mapP[face_node_index] # compute f(u, ∇u) ⋅ n - flux_face_value = zero(eltype(flux_face_values)) + flux_face_value = zero(eltype(div_flux_face_values)) for dim in eachdim(mesh) - uM = grad_u_face_values[dim][idM] - # compute flux if node is not a boundary node - if idM != idP - uP = grad_u_face_values[dim][idP] + uM = viscous_flux_face_values[dim][idM] + uP = viscous_flux_face_values[dim][idP] # TODO: use strong/weak formulation? - flux_face_value = flux_face_value + 0.5 * (uP - uM) * nxyzJ[dim][face_node_index] - end + flux_face_value = flux_face_value + 0.5 * (uP - uM) * nxyzJ[dim][face_node_index] end - flux_face_values[idM] = flux_face_value + div_flux_face_values[idM] = flux_face_value end - calc_boundary_flux!(flux_face_values, u, t, Divergence(), boundary_conditions, - mesh, equations, dg, cache, parabolic_cache) + calc_boundary_flux!(div_flux_face_values, viscous_flux_face_values, t, Divergence(), + boundary_conditions, mesh, equations, dg, cache, parabolic_cache) # surface contributions - apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, flux_face_values) + apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, div_flux_face_values) invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0) end From e556cde6a9fb24c29f3180ed9bcb963cc28ab3b5 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 23:57:07 -0500 Subject: [PATCH 40/86] BC prepping --- src/solvers/dgmulti/dg_parabolic.jl | 33 ++++++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 2c9927a5b3c..6846008089b 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -123,9 +123,36 @@ calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{( mesh, equations, dg::DGMulti, cache, parabolic_cache) = nothing # TODO: finish -function calc_single_boundary_flux!(flux, u, t, operator_type, boundary_condition, - boundary_key, mesh, equations, dg::DGMulti, - cache, parabolic_cache) +function calc_single_boundary_flux!(div_flux_face_values, u_face_values, t, operator_type, boundary_condition, boundary_key, + mesh, equations, dg::DGMulti{NDIMS}, cache, parabolic_cache) where {NDIMS} + rd = dg.basis + md = mesh.md + + # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total). + # mesh.boundary_faces indexes into the columns of these face-reshaped arrays. + num_pts_per_face = rd.Nfq ÷ rd.Nfaces + num_faces_total = rd.Nfaces * md.num_elements + + # This function was originally defined as + # `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`. + # This results in allocations due to https://github.com/JuliaLang/julia/issues/36313. + # To avoid allocations, we use Tim Holy's suggestion: + # https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300. + reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ()) + + # loop through boundary faces, which correspond to columns of reshaped `u_face_values` + @unpack xyzf, nxyzJ, Jf = md + xyzf, nxyzJ, Jf = reshape_by_face.(xyzf), reshape_by_face.(nxyzJ), reshape_by_face(Jf) + for f in mesh.boundary_faces[boundary_key] + for i in Base.OneTo(num_pts_per_face) + face_normal = SVector{NDIMS}(getindex.(nxyzJ, i, f)) / Jf[i,f] + face_coordinates = SVector{NDIMS}(getindex.(xyzf, i, f)) + # div_flux_face_values[i,f] = boundary_condition(u_face_values[i,f], + # face_normal, face_coordinates, t, + # surface_flux, equations) * Jf[i,f] + end + end + return nothing end From 3edb18ce2b9fb9542a6d5416f4b3e51234d3b82b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 22 May 2022 23:57:26 -0500 Subject: [PATCH 41/86] cleanup of elixir --- examples/dgmulti_2d/elixir_advection_diffusion.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index 557a4f0d671..d66311d8fa3 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -11,11 +11,10 @@ parabolic_equations = ScalarDiffusion2D(1e-3) initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) initial_condition = initial_condition_zero -# example where we tag two separate boundary segments of the mesh +# tag two separate boundary segments inflow(x, tol=50*eps()) = abs(x[1]+1) inflow, :outflow => rest_of_boundary) - mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16), is_on_boundary=is_on_boundary) # define inviscid boundary conditions @@ -24,12 +23,13 @@ boundary_conditions = (; :inflow => boundary_condition_inflow) # define viscous boundary conditions boundary_condition_wall = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) -parabolic_boundary_conditions = (; boundary_condition_wall) +parabolic_boundary_conditions = (; :inflow => boundary_condition_inflow, + :outflow => boundary_condition_wall) semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg; - boundary_conditions)#, parabolic_boundary_conditions) + boundary_conditions, parabolic_boundary_conditions) -tspan = (0.0, .25) +tspan = (0.0, .5) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() From 455974a2cd9a24ad2f745b70019a6f243b857063 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 23 May 2022 00:28:22 -0500 Subject: [PATCH 42/86] more name fixing --- src/solvers/dgmulti/dg_parabolic.jl | 42 ++++++++++++++--------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 6846008089b..f5c7e061cac 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -2,22 +2,22 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, dg::DGMultiWeakForm, RealT, uEltype) nvars = nvariables(equations) - # u_parabolic stores "transformed" variables for computing the gradient + # u_transformed stores "transformed" variables for computing the gradient @unpack md = mesh u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) u_grad = ntuple(_ -> similar(u_transformed), ndims(mesh)) + viscous_flux = similar.(u_grad) + u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) + scalar_flux_face_values = similar(u_face_values) grad_u_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) - viscous_flux = similar.(u_grad) - div_flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) - local_viscous_flux_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] - local_flux_values_threaded = [similar(div_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] + local_flux_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] return (; u_transformed, u_grad, viscous_flux, - weak_differentiation_matrices, - u_face_values, grad_u_face_values, div_flux_face_values, + # weak_differentiation_matrices, + u_face_values, grad_u_face_values, scalar_flux_face_values, local_viscous_flux_threaded, local_flux_values_threaded) end @@ -36,7 +36,7 @@ function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::Abs apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) end -function calc_gradient_surface_integral(u_grad, u, div_flux_face_values, +function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, mesh, equations::AbstractParabolicEquations, dg::DGMulti, cache, parabolic_cache) @unpack local_flux_values_threaded = parabolic_cache @@ -45,7 +45,7 @@ function calc_gradient_surface_integral(u_grad, u, div_flux_face_values, for dim in eachdim(mesh) for i in eachindex(local_flux_values) # compute [u] * (nx, ny, nz) - local_flux_values[i] = div_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] + local_flux_values[i] = scalar_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] end apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(u_grad[dim], :, e), local_flux_values) end @@ -73,20 +73,20 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, mesh, equations, dg.surface_integral, dg, cache) # compute fluxes at interfaces - @unpack u_face_values, div_flux_face_values = parabolic_cache + @unpack u_face_values, scalar_flux_face_values = parabolic_cache @unpack mapM, mapP, Jf = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] uP = u_face_values[idP] - div_flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? + scalar_flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? end - calc_boundary_flux!(div_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, + calc_boundary_flux!(scalar_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, mesh, equations, dg, cache, parabolic_cache) # compute surface contributions - calc_gradient_surface_integral(u_grad, u, div_flux_face_values, + calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, mesh, equations, dg, cache, parabolic_cache) for dim in eachdim(mesh) @@ -123,7 +123,7 @@ calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{( mesh, equations, dg::DGMulti, cache, parabolic_cache) = nothing # TODO: finish -function calc_single_boundary_flux!(div_flux_face_values, u_face_values, t, operator_type, boundary_condition, boundary_key, +function calc_single_boundary_flux!(scalar_flux_face_values, u_face_values, t, operator_type, boundary_condition, boundary_key, mesh, equations, dg::DGMulti{NDIMS}, cache, parabolic_cache) where {NDIMS} rd = dg.basis md = mesh.md @@ -147,7 +147,7 @@ function calc_single_boundary_flux!(div_flux_face_values, u_face_values, t, oper for i in Base.OneTo(num_pts_per_face) face_normal = SVector{NDIMS}(getindex.(nxyzJ, i, f)) / Jf[i,f] face_coordinates = SVector{NDIMS}(getindex.(xyzf, i, f)) - # div_flux_face_values[i,f] = boundary_condition(u_face_values[i,f], + # scalar_flux_face_values[i,f] = boundary_condition(u_face_values[i,f], # face_normal, face_coordinates, t, # surface_flux, equations) * Jf[i,f] end @@ -221,27 +221,27 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh end # compute fluxes at interfaces - @unpack div_flux_face_values = parabolic_cache + @unpack scalar_flux_face_values = parabolic_cache @unpack mapM, mapP, nxyzJ = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) idM, idP = mapM[face_node_index], mapP[face_node_index] # compute f(u, ∇u) ⋅ n - flux_face_value = zero(eltype(div_flux_face_values)) + flux_face_value = zero(eltype(scalar_flux_face_values)) for dim in eachdim(mesh) uM = viscous_flux_face_values[dim][idM] uP = viscous_flux_face_values[dim][idP] - # TODO: use strong/weak formulation? + # TODO: use strong/weak formulation? flux_face_value = flux_face_value + 0.5 * (uP - uM) * nxyzJ[dim][face_node_index] end - div_flux_face_values[idM] = flux_face_value + scalar_flux_face_values[idM] = flux_face_value end - calc_boundary_flux!(div_flux_face_values, viscous_flux_face_values, t, Divergence(), + calc_boundary_flux!(scalar_flux_face_values, viscous_flux_face_values, t, Divergence(), boundary_conditions, mesh, equations, dg, cache, parabolic_cache) # surface contributions - apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, div_flux_face_values) + apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, scalar_flux_face_values) invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0) end From fd2489ebef4d35998e5c5d0a5516a3f4660e3b6b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 23 May 2022 00:33:27 -0500 Subject: [PATCH 43/86] Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Hendrik Ranocha --- .../semidiscretization_hyperbolic_parabolic.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 8cfe1ea688c..440e78008ac 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -11,7 +11,6 @@ A struct containing everything needed to describe a spatial semidiscretization of a mixed hyperbolic-parabolic conservation law. """ - struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations, InitialCondition, BoundaryConditions, ParabolicBoundaryConditions, SourceTerms, Solver, Cache, ParabolicCache} <: AbstractSemidiscretization From 4ef7835c09ab7101c290eacc1f5b355b007362f6 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 23 May 2022 14:08:23 -0500 Subject: [PATCH 44/86] adding Neumann BCs --- src/Trixi.jl | 2 +- src/equations/scalar_diffusion_2d.jl | 33 +++++++++++++++++++++++++++- 2 files changed, 33 insertions(+), 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index c435dc20dfc..cd3eae91c2b 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -133,7 +133,7 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D -export ScalarDiffusion2D +export ScalarDiffusion2D, BoundaryConditionNeumann export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, diff --git a/src/equations/scalar_diffusion_2d.jl b/src/equations/scalar_diffusion_2d.jl index 26cef6692b2..34a257f6b51 100644 --- a/src/equations/scalar_diffusion_2d.jl +++ b/src/equations/scalar_diffusion_2d.jl @@ -2,9 +2,40 @@ struct ScalarDiffusion2D{T} <: AbstractParabolicEquations{2, 1} diffusivity::T end -# no orientation here since the flux is vector-valued +# no orientation specified since the flux is vector-valued function flux(u, grad_u, equations::ScalarDiffusion2D) dudx, dudy = grad_u return equations.diffusivity * dudx, equations.diffusivity * dudy end +# Dirichlet-type boundary condition for use with a parabolic solver in weak form +@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, x, t, + operator_type::Gradient, + equations::ScalarDiffusion2D) + return boundary_condition.boundary_value_function(x, t, equations) +end + +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, normal::AbstractVector, x, t, + operator_type::Divergence, + equations::ScalarDiffusion2D) + return flux_inner +end + +# `boundary_value_function` should have signature `boundary_value_function(x, t, equations)` +# For Neumann BCs, this corresponds to kappa * (∇u ⋅ n) +struct BoundaryConditionNeumann{B} + boundary_value_function::B +end + +# Neumann type boundary condition for use with a parabolic solver in weak form +@inline function (boundary_condition::BoundaryConditionNeumann)(u_inner, normal::AbstractVector, x, t, + operator_type::Gradient, + equations::ScalarDiffusion2D) + return u_inner +end + +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, x, t, + operator_type::Divergence, + equations::ScalarDiffusion2D) + return boundary_condition.boundary_value_function(x, t, equations) +end \ No newline at end of file From 48a715a74943a98f0641768e8137d846d0484315 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 23 May 2022 14:10:08 -0500 Subject: [PATCH 45/86] switching to weak form weak form more weak form weak form --- src/solvers/dgmulti/dg_parabolic.jl | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index f5c7e061cac..6738c396dc6 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -2,6 +2,9 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, dg::DGMultiWeakForm, RealT, uEltype) nvars = nvariables(equations) + @unpack M, Drst = dg.basis + weak_differentiation_matrices = map(A -> -M \ (A' * M), Drst) + # u_transformed stores "transformed" variables for computing the gradient @unpack md = mesh u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) @@ -16,7 +19,7 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, local_flux_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] return (; u_transformed, u_grad, viscous_flux, - # weak_differentiation_matrices, + weak_differentiation_matrices, u_face_values, grad_u_face_values, scalar_flux_face_values, local_viscous_flux_threaded, local_flux_values_threaded) end @@ -56,6 +59,8 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, equations::AbstractParabolicEquations, boundary_conditions, dg::DGMulti, cache, parabolic_cache) + @unpack weak_differentiation_matrices = parabolic_cache + for dim in 1:length(u_grad) reset_du!(u_grad[dim], dg) end @@ -63,26 +68,26 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, # compute volume contributions to gradients @threaded for e in eachelement(mesh, dg) for i in eachdim(mesh), j in eachdim(mesh) - dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine - apply_to_each_field(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), + dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # TODO: assumes mesh is affine + apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj), view(u_grad[i], :, e), view(u, :, e)) end end - prolong2interfaces!(cache.u_face_values, u, - mesh, equations, dg.surface_integral, dg, cache) + @unpack u_face_values = parabolic_cache + prolong2interfaces!(u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) # compute fluxes at interfaces - @unpack u_face_values, scalar_flux_face_values = parabolic_cache + @unpack scalar_flux_face_values = parabolic_cache @unpack mapM, mapP, Jf = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] uP = u_face_values[idP] - scalar_flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation? + scalar_flux_face_values[idM] = 0.5 * (uP + uM) # TODO: use strong/weak formulation for curved meshes? end - calc_boundary_flux!(scalar_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, + calc_boundary_flux!(scalar_flux_face_values, nothing, t, Gradient(), boundary_conditions, mesh, equations, dg, cache, parabolic_cache) # compute surface contributions @@ -202,13 +207,15 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh equations::AbstractParabolicEquations, boundary_conditions, dg::DGMulti, cache, parabolic_cache) + @unpack weak_differentiation_matrices = parabolic_cache + reset_du!(du, dg) # compute volume contributions to divergence @threaded for e in eachelement(mesh, dg) for i in eachdim(mesh), j in eachdim(mesh) dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine - apply_to_each_field(mul_by_accum!(dg.basis.Drst[j], dxidxhatj), + apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj), view(du, :, e), view(viscous_flux[i], :, e)) end end @@ -232,12 +239,12 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh uM = viscous_flux_face_values[dim][idM] uP = viscous_flux_face_values[dim][idP] # TODO: use strong/weak formulation? - flux_face_value = flux_face_value + 0.5 * (uP - uM) * nxyzJ[dim][face_node_index] + flux_face_value = flux_face_value + 0.5 * (uP + uM) * nxyzJ[dim][face_node_index] end scalar_flux_face_values[idM] = flux_face_value end - calc_boundary_flux!(scalar_flux_face_values, viscous_flux_face_values, t, Divergence(), + calc_boundary_flux!(scalar_flux_face_values, nothing, t, Divergence(), boundary_conditions, mesh, equations, dg, cache, parabolic_cache) # surface contributions From 45d2141f4315a88ed7c778a85636cd456b721a8a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 23 May 2022 14:10:42 -0500 Subject: [PATCH 46/86] fixing cache variable names --- src/solvers/dgmulti/dg_parabolic.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 6738c396dc6..7b11c9b02c8 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -15,13 +15,14 @@ function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, scalar_flux_face_values = similar(u_face_values) grad_u_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) + local_u_values_threaded = [similar(u_transformed, dg.basis.Nq) for _ in 1:Threads.nthreads()] local_viscous_flux_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] - local_flux_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] + local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] return (; u_transformed, u_grad, viscous_flux, weak_differentiation_matrices, u_face_values, grad_u_face_values, scalar_flux_face_values, - local_viscous_flux_threaded, local_flux_values_threaded) + local_u_values_threaded, local_viscous_flux_threaded, local_flux_face_values_threaded) end # Transform solution variables prior to taking the gradient @@ -42,12 +43,12 @@ end function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, mesh, equations::AbstractParabolicEquations, dg::DGMulti, cache, parabolic_cache) - @unpack local_flux_values_threaded = parabolic_cache + @unpack local_flux_face_values_threaded = parabolic_cache @threaded for e in eachelement(mesh, dg) - local_flux_values = local_flux_values_threaded[Threads.threadid()] + local_flux_values = local_flux_face_values_threaded[Threads.threadid()] for dim in eachdim(mesh) for i in eachindex(local_flux_values) - # compute [u] * (nx, ny, nz) + # compute flux * (nx, ny, nz) local_flux_values[i] = scalar_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] end apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(u_grad[dim], :, e), local_flux_values) @@ -169,8 +170,7 @@ function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, reset_du!(viscous_flux[dim], dg) end - @unpack local_viscous_flux_threaded = parabolic_cache - local_u_values_threaded = parabolic_cache.local_flux_values_threaded + @unpack local_viscous_flux_threaded, local_u_values_threaded = parabolic_cache @threaded for e in eachelement(mesh, dg) From bbba3bac628cfe21d65ce2aff0d0742b2cdb8853 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 23 May 2022 14:10:53 -0500 Subject: [PATCH 47/86] simplified parabolic boundary flux treatement --- src/solvers/dgmulti/dg_parabolic.jl | 38 +++++++++++++---------------- 1 file changed, 17 insertions(+), 21 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 7b11c9b02c8..8aa5b23a576 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -128,37 +128,33 @@ end calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{(),Tuple{}}, mesh, equations, dg::DGMulti, cache, parabolic_cache) = nothing -# TODO: finish -function calc_single_boundary_flux!(scalar_flux_face_values, u_face_values, t, operator_type, boundary_condition, boundary_key, +# TODO: decide if we want to use the input `u_face_values` +function calc_single_boundary_flux!(flux_face_values, u_face_values, t, + operator_type, boundary_condition, boundary_key, mesh, equations, dg::DGMulti{NDIMS}, cache, parabolic_cache) where {NDIMS} rd = dg.basis md = mesh.md - # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total). - # mesh.boundary_faces indexes into the columns of these face-reshaped arrays. num_pts_per_face = rd.Nfq ÷ rd.Nfaces - num_faces_total = rd.Nfaces * md.num_elements - - # This function was originally defined as - # `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`. - # This results in allocations due to https://github.com/JuliaLang/julia/issues/36313. - # To avoid allocations, we use Tim Holy's suggestion: - # https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300. - reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ()) - - # loop through boundary faces, which correspond to columns of reshaped `u_face_values` @unpack xyzf, nxyzJ, Jf = md - xyzf, nxyzJ, Jf = reshape_by_face.(xyzf), reshape_by_face.(nxyzJ), reshape_by_face(Jf) for f in mesh.boundary_faces[boundary_key] for i in Base.OneTo(num_pts_per_face) - face_normal = SVector{NDIMS}(getindex.(nxyzJ, i, f)) / Jf[i,f] - face_coordinates = SVector{NDIMS}(getindex.(xyzf, i, f)) - # scalar_flux_face_values[i,f] = boundary_condition(u_face_values[i,f], - # face_normal, face_coordinates, t, - # surface_flux, equations) * Jf[i,f] + + # reverse engineer element + face node indices (avoids reshaping arrays) + e = ((f-1) ÷ rd.Nfaces) + 1 + fid = i + ((f-1) % rd.Nfaces) * num_pts_per_face + + face_normal = SVector{NDIMS}(getindex.(nxyzJ, fid, e)) / Jf[fid,e] + face_coordinates = SVector{NDIMS}(getindex.(xyzf, fid, e)) + + # for both the gradient and the divergence, the boundary flux is scalar valued. + # for the gradient, it is the solution; for divergence, it is the normal flux. + u_boundary = boundary_condition(flux_face_values[fid,e], + face_normal, face_coordinates, t, + operator_type, equations) + flux_face_values[fid,e] = u_boundary end end - return nothing end From 8a11d57b77df86feb791b75a9196cc6325ea4cde Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 23 May 2022 14:11:02 -0500 Subject: [PATCH 48/86] updated elixir with different types of BCs --- .../dgmulti_2d/elixir_advection_diffusion.jl | 27 ++++++++++++------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index d66311d8fa3..cc9f34fd369 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -5,8 +5,8 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), volume_integral = VolumeIntegralWeakForm()) -equations = LinearScalarAdvectionEquation2D(1.0, 1.0) -parabolic_equations = ScalarDiffusion2D(1e-3) +equations = LinearScalarAdvectionEquation2D(1.5, 1.0) +parabolic_equations = ScalarDiffusion2D(1e-2) initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) initial_condition = initial_condition_zero @@ -14,22 +14,31 @@ initial_condition = initial_condition_zero # tag two separate boundary segments inflow(x, tol=50*eps()) = abs(x[1]+1) inflow, :outflow => rest_of_boundary) +left(x, tol=50*eps()) = abs(x[1]+1) inflow, :outflow => rest_of_boundary, + :left => left, :right => right, :top => top, :bottom => bottom) mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16), is_on_boundary=is_on_boundary) +# BC types +boundary_condition_one = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)) +boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) +boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0)) + # define inviscid boundary conditions -boundary_condition_inflow = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)) -boundary_conditions = (; :inflow => boundary_condition_inflow) +boundary_conditions = (; :left => boundary_condition_one, :bottom => boundary_condition_zero) # define viscous boundary conditions -boundary_condition_wall = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) -parabolic_boundary_conditions = (; :inflow => boundary_condition_inflow, - :outflow => boundary_condition_wall) +parabolic_boundary_conditions = (; :left => boundary_condition_one, + :bottom => boundary_condition_zero, + :top => boundary_condition_zero, + :right => boundary_condition_neumann_zero) semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg; boundary_conditions, parabolic_boundary_conditions) -tspan = (0.0, .5) +tspan = (0.0, 1.5) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() From 478e831bcab05d6a726cdd08f8ecf9081e8f05a0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 24 May 2022 12:27:21 -0500 Subject: [PATCH 49/86] Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper --- .../semidiscretization_hyperbolic_parabolic.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 440e78008ac..e73231648a0 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -54,12 +54,13 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations end """ - SemidiscretizationHyperbolicParabolic(mesh, equations, initial_condition, solver; - source_terms=nothing, - boundary_conditions=boundary_condition_periodic, - RealT=real(solver), - uEltype=RealT, - initial_cache=NamedTuple()) + SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, solver; + source_terms=nothing, + boundary_conditions=boundary_condition_periodic, + boundary_conditions_parabolic=boundary_condition_periodic, + RealT=real(solver), + uEltype=RealT, + initial_cache=NamedTuple()) Construct a semidiscretization of a hyperbolic PDE. """ From cffc1391557a9e3213e67877cbd0cab94ea279b9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 24 May 2022 12:29:40 -0500 Subject: [PATCH 50/86] Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper --- .../semidiscretization_hyperbolic_parabolic.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index e73231648a0..30d5762ded5 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -97,15 +97,15 @@ end # semantics we want to use here. In particular, it allows us to re-use mutable parts, # e.g. `remake(semi).mesh === semi.mesh`. function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.solver), - mesh=semi.mesh, - equations=semi.equations, - parabolic_equations=semi.parabolic_equations, - initial_condition=semi.initial_condition, - solver=semi.solver, - source_terms=semi.source_terms, - boundary_conditions=semi.boundary_conditions, - parabolic_boundary_conditions=semi.parabolic_boundary_conditions - ) + mesh=semi.mesh, + equations=semi.equations, + parabolic_equations=semi.parabolic_equations, + initial_condition=semi.initial_condition, + solver=semi.solver, + source_terms=semi.source_terms, + boundary_conditions=semi.boundary_conditions, + parabolic_boundary_conditions=semi.parabolic_boundary_conditions + ) # TODO: Which parts do we want to `remake`? At least the solver needs some # special care if shock-capturing volume integrals are used (because of # the indicators and their own caches...). From 924a5b9ca411962ae9905ea53319b0bcbe12dc02 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 24 May 2022 12:38:15 -0500 Subject: [PATCH 51/86] Update src/equations/equations.jl Co-authored-by: Michael Schlottke-Lakemper --- src/equations/equations.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index c43bb88fa9d..31c15159a3a 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -330,6 +330,9 @@ abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEqu include("acoustic_perturbation_2d.jl") abstract type AbstractParabolicEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end + +# Linear scalar diffusion for use in linear scalar advection-diffusion problems +abstract type AbstractScalarDiffusionEquations{NDIMS, NVARS} <: AbstractParabolicEquations{NDIMS, NVARS} end include("scalar_diffusion_2d.jl") end # @muladd From f7a02842462aa613ba7eadd23280de42ab4c7a5f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 12:31:14 -0500 Subject: [PATCH 52/86] rhs! -> rhs_parabolic! --- .../semidiscretization_hyperbolic_parabolic.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 30d5762ded5..c9849d7d93d 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -223,9 +223,9 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() - @trixi_timeit timer() "parabolic rhs!" rhs!(du, u, t, mesh, parabolic_equations, initial_condition, - parabolic_boundary_conditions, source_terms, - solver, cache, parabolic_cache) + @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, parabolic_equations, initial_condition, + parabolic_boundary_conditions, source_terms, + solver, cache, parabolic_cache) runtime = time_ns() - time_start put!(semi.performance_counter, runtime) From 3a799f590d74b9d00b5f51720ab3ecc2d42d20db Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 12:31:28 -0500 Subject: [PATCH 53/86] move Gradient/Divergence types --- src/equations/equations.jl | 4 ++++ src/solvers/dgmulti/dg_parabolic.jl | 4 ---- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 31c15159a3a..ca0f4b2561a 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -331,6 +331,10 @@ include("acoustic_perturbation_2d.jl") abstract type AbstractParabolicEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +# operator types used for dispatch on parabolic boundary fluxes +struct Gradient end +struct Divergence end + # Linear scalar diffusion for use in linear scalar advection-diffusion problems abstract type AbstractScalarDiffusionEquations{NDIMS, NVARS} <: AbstractParabolicEquations{NDIMS, NVARS} end include("scalar_diffusion_2d.jl") diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 8aa5b23a576..b18e0f5fb38 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -100,10 +100,6 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, end end -# operator types used for dispatch on boundary fluxes -struct Gradient end -struct Divergence end - # do nothing for periodic domains function calc_boundary_flux!(flux, u, t, operator_type, ::BoundaryConditionPeriodic, mesh, equations::AbstractParabolicEquations, dg::DGMulti, From a9570944a0e7b6753fc87cebfab39a533b78cb61 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 12:31:37 -0500 Subject: [PATCH 54/86] add "no boundary condition" BC --- src/equations/equations.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ca0f4b2561a..f67363119d4 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -175,6 +175,10 @@ end return flux end +# Imposing no boundary condition just evaluates the flux at the inner state. +@inline function no_boundary_condition(flux, other_args...) + return flux +end # set sensible default values that may be overwritten by specific equations """ From 2e5146c52b7b47fc4294ccda63ae7de4e085e565 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 12:32:32 -0500 Subject: [PATCH 55/86] rename ScalarDiffusion -> LaplaceDiffusion r r --- .../dgmulti_2d/elixir_advection_diffusion.jl | 21 +++++------ .../elixir_advection_diffusion_periodic.jl | 2 +- src/Trixi.jl | 4 ++- src/equations/scalar_diffusion_2d.jl | 35 +++++++++++-------- test/test_dgmulti_2d.jl | 2 +- 5 files changed, 36 insertions(+), 28 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index cc9f34fd369..7b10a3e19d7 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -6,31 +6,32 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.5, 1.0) -parabolic_equations = ScalarDiffusion2D(1e-2) +parabolic_equations = LaplaceDiffusion2D(5e-2) initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) initial_condition = initial_condition_zero -# tag two separate boundary segments -inflow(x, tol=50*eps()) = abs(x[1]+1) inflow, :outflow => rest_of_boundary, - :left => left, :right => right, :top => top, :bottom => bottom) -mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16), is_on_boundary=is_on_boundary) +top(x, tol=50*eps()) = abs(x[2]-1) left, :right => right, :top => top, :bottom => bottom) +mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary) # BC types -boundary_condition_one = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)) +boundary_condition_left = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + .1 * x[2])) boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0)) # define inviscid boundary conditions -boundary_conditions = (; :left => boundary_condition_one, :bottom => boundary_condition_zero) +boundary_conditions = (; :left => boundary_condition_left, + :bottom => boundary_condition_zero, + :top => no_boundary_condition, + :right => no_boundary_condition) # define viscous boundary conditions -parabolic_boundary_conditions = (; :left => boundary_condition_one, +parabolic_boundary_conditions = (; :left => boundary_condition_left, :bottom => boundary_condition_zero, :top => boundary_condition_zero, :right => boundary_condition_neumann_zero) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl index 565f6c9d4cb..a435acb8209 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl @@ -6,7 +6,7 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) -parabolic_equations = ScalarDiffusion2D(1e-2) +parabolic_equations = LaplaceDiffusion2D(1e-2) function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) return SVector(exp(-100 * (x[1]^2 + x[2]^2))) diff --git a/src/Trixi.jl b/src/Trixi.jl index cd3eae91c2b..a2cf650c207 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -133,7 +133,9 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D -export ScalarDiffusion2D, BoundaryConditionNeumann +export no_boundary_condition + +export LaplaceDiffusion2D, BoundaryConditionNeumann export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, diff --git a/src/equations/scalar_diffusion_2d.jl b/src/equations/scalar_diffusion_2d.jl index 34a257f6b51..7e0a417d304 100644 --- a/src/equations/scalar_diffusion_2d.jl +++ b/src/equations/scalar_diffusion_2d.jl @@ -1,23 +1,29 @@ -struct ScalarDiffusion2D{T} <: AbstractParabolicEquations{2, 1} +@doc raw""" + LaplaceDiffusion2D{T} <: AbstractParabolicEquations{2, 1} + +`LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` +applied to each solution component. +""" +struct LaplaceDiffusion2D{T} <: AbstractScalarDiffusionEquations{2, 1} diffusivity::T end # no orientation specified since the flux is vector-valued -function flux(u, grad_u, equations::ScalarDiffusion2D) +function flux(u, grad_u, equations::LaplaceDiffusion2D) dudx, dudy = grad_u return equations.diffusivity * dudx, equations.diffusivity * dudy end # Dirichlet-type boundary condition for use with a parabolic solver in weak form -@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, x, t, - operator_type::Gradient, - equations::ScalarDiffusion2D) +@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::LaplaceDiffusion2D) return boundary_condition.boundary_value_function(x, t, equations) end -@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, normal::AbstractVector, x, t, - operator_type::Divergence, - equations::ScalarDiffusion2D) +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::LaplaceDiffusion2D) return flux_inner end @@ -27,15 +33,14 @@ struct BoundaryConditionNeumann{B} boundary_value_function::B end -# Neumann type boundary condition for use with a parabolic solver in weak form -@inline function (boundary_condition::BoundaryConditionNeumann)(u_inner, normal::AbstractVector, x, t, - operator_type::Gradient, - equations::ScalarDiffusion2D) +@inline function (boundary_condition::BoundaryConditionNeumann)(u_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations::LaplaceDiffusion2D) return u_inner end -@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, x, t, - operator_type::Divergence, - equations::ScalarDiffusion2D) +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations::LaplaceDiffusion2D) return boundary_condition.boundary_value_function(x, t, equations) end \ No newline at end of file diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 32f24fef774..c7eb8997d8a 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -25,7 +25,7 @@ isdir(outdir) && rm(outdir, recursive=true) initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2]) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) - equations_parabolic = Trixi.ScalarDiffusion2D(1.0) + equations_parabolic = Trixi.LaplaceDiffusion2D(1.0) semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) ode = semidiscretize(semi, (0.0, 0.01)) From 2c452935a12f421072e48d87525128f8b7328e60 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 12:32:52 -0500 Subject: [PATCH 56/86] rhs -> rhs_parabolic! again --- src/solvers/dgmulti/dg_parabolic.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index b18e0f5fb38..b984fa5915d 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -251,9 +251,9 @@ end # 2. compute f(u, grad(u)) # 3. compute div(u) # boundary conditions will be applied to both grad(u) and div(u). -function rhs!(du, u, t, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEquations, - initial_condition, boundary_conditions, source_terms, - dg::DGMulti, cache, parabolic_cache) +function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEquations, + initial_condition, boundary_conditions, source_terms, + dg::DGMulti, cache, parabolic_cache) reset_du!(du, dg) From 7159af3e57d37e6d3606d5882be1f2bd90bdd50d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 12:33:13 -0500 Subject: [PATCH 57/86] comments --- src/solvers/dgmulti/dg_parabolic.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index b984fa5915d..e14066e3f47 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -124,7 +124,7 @@ end calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{(),Tuple{}}, mesh, equations, dg::DGMulti, cache, parabolic_cache) = nothing -# TODO: decide if we want to use the input `u_face_values` +# TODO: DGMulti. Decide if we want to use the input `u_face_values` (currently unused) function calc_single_boundary_flux!(flux_face_values, u_face_values, t, operator_type, boundary_condition, boundary_key, mesh, equations, dg::DGMulti{NDIMS}, cache, parabolic_cache) where {NDIMS} @@ -236,6 +236,7 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh scalar_flux_face_values[idM] = flux_face_value end + # TODO: decide what to pass in calc_boundary_flux!(scalar_flux_face_values, nothing, t, Divergence(), boundary_conditions, mesh, equations, dg, cache, parabolic_cache) From d2cccd0d93e97c37f9797d6142f5dffe80fd9484 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 12:35:01 -0500 Subject: [PATCH 58/86] moving tests around --- test/test_dgmulti_2d.jl | 40 ---------------------------- test/test_parabolic.jl | 59 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 40 deletions(-) create mode 100644 test/test_parabolic.jl diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index c7eb8997d8a..90dd082c3ed 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -13,46 +13,6 @@ isdir(outdir) && rm(outdir, recursive=true) @testset "DGMulti 2D" begin - @trixi_testset "SemidiscretizationHyperbolicParabolic" begin - using Test - - dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), - surface_integral = SurfaceIntegralWeakForm(flux_central), - volume_integral = VolumeIntegralWeakForm()) - mesh = DGMultiMesh(dg, cells_per_dimension=(2, 2)) - - # test with polynomial initial condition x^2 * y - initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2]) - - equations = LinearScalarAdvectionEquation2D(1.0, 1.0) - equations_parabolic = Trixi.LaplaceDiffusion2D(1.0) - - semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) - ode = semidiscretize(semi, (0.0, 0.01)) - - @unpack cache, parabolic_cache, parabolic_equations = semi - @unpack u_grad = parabolic_cache - for dim in 1:length(u_grad) - fill!(u_grad[dim], zero(eltype(u_grad[dim]))) - end - - # pass in BoundaryConditionPeriodic() to fake a "do nothing" BC - Trixi.calc_gradient!(u_grad, ode.u0, mesh, parabolic_equations, Trixi.BoundaryConditionPeriodic(), dg, cache, parabolic_cache) - @unpack x, y = mesh.md - @test norm(getindex.(u_grad[1], 1) - 2 * x .* y) < 1e3 * eps() - @test norm(getindex.(u_grad[2], 1) - x.^2) < 1e3 * eps() - - u_flux = similar.(u_grad) - Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, parabolic_equations, dg, cache, parabolic_cache) - @test u_flux[1] ≈ u_grad[1] - @test u_flux[2] ≈ u_grad[2] - - du = similar(ode.u0) - # pass in BoundaryConditionPeriodic() to fake a "do nothing" BC - Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, parabolic_equations, Trixi.BoundaryConditionPeriodic(), dg, cache, parabolic_cache) - @test norm(getindex.(du, 1) - 2 * y) < 1e3 * eps() - end - @trixi_testset "elixir_euler_weakform.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), cells_per_dimension = (4, 4), diff --git a/test/test_parabolic.jl b/test/test_parabolic.jl new file mode 100644 index 00000000000..e6bbfeffad1 --- /dev/null +++ b/test/test_parabolic.jl @@ -0,0 +1,59 @@ +module TestExamplesDGMulti2D + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d") + +# Start with a clean environment: remove Trixi output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive=true) + +@testset "SemidiscretizationHyperbolicParabolic" begin + + @trixi_testset "DGMulti 2D" begin + + dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_central), + volume_integral = VolumeIntegralWeakForm()) + mesh = DGMultiMesh(dg, cells_per_dimension=(2, 2)) + + # test with polynomial initial condition x^2 * y + # test if we recover the exact second derivative + initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2]) + + equations = LinearScalarAdvectionEquation2D(1.0, 1.0) + equations_parabolic = Trixi.LaplaceDiffusion2D(1.0) + + semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) + ode = semidiscretize(semi, (0.0, 0.01)) + + @unpack cache, parabolic_cache, parabolic_equations = semi + @unpack u_grad = parabolic_cache + for dim in 1:length(u_grad) + fill!(u_grad[dim], zero(eltype(u_grad[dim]))) + end + + Trixi.calc_gradient!(u_grad, ode.u0, mesh, parabolic_equations, no_boundary_condition, dg, cache, parabolic_cache) + @unpack x, y = mesh.md + @test norm(getindex.(u_grad[1], 1) - 2 * x .* y) < 1e3 * eps() + @test norm(getindex.(u_grad[2], 1) - x.^2) < 1e3 * eps() + + u_flux = similar.(u_grad) + Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, parabolic_equations, dg, cache, parabolic_cache) + @test u_flux[1] ≈ u_grad[1] + @test u_flux[2] ≈ u_grad[2] + + du = similar(ode.u0) + Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, parabolic_equations, no_boundary_condition, dg, cache, parabolic_cache) + @test norm(getindex.(du, 1) - 2 * y) < 1e3 * eps() + end + +end + +# Clean up afterwards: delete Trixi output directory +@test_nowarn isdir(outdir) && rm(outdir, recursive=true) + +end # module From b73eb7436b2a14df516d3df8a340084da926f9c0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 13:29:29 -0500 Subject: [PATCH 59/86] scalarDiffusion -> LaplaceDiffusion --- src/equations/equations.jl | 4 ++-- .../{scalar_diffusion_2d.jl => laplace_diffusion_2d.jl} | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) rename src/equations/{scalar_diffusion_2d.jl => laplace_diffusion_2d.jl} (96%) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index f67363119d4..8c3f7d1ff01 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -340,7 +340,7 @@ struct Gradient end struct Divergence end # Linear scalar diffusion for use in linear scalar advection-diffusion problems -abstract type AbstractScalarDiffusionEquations{NDIMS, NVARS} <: AbstractParabolicEquations{NDIMS, NVARS} end -include("scalar_diffusion_2d.jl") +abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractParabolicEquations{NDIMS, NVARS} end +include("laplace_diffusion_2d.jl") end # @muladd diff --git a/src/equations/scalar_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl similarity index 96% rename from src/equations/scalar_diffusion_2d.jl rename to src/equations/laplace_diffusion_2d.jl index 7e0a417d304..28b88bcf8f4 100644 --- a/src/equations/scalar_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -4,7 +4,7 @@ `LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` applied to each solution component. """ -struct LaplaceDiffusion2D{T} <: AbstractScalarDiffusionEquations{2, 1} +struct LaplaceDiffusion2D{T} <: AbstractLaplaceDiffusionEquations{2, 1} diffusivity::T end From 363f8056649101da43d5d0245921b310261e7419 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 13:32:40 -0500 Subject: [PATCH 60/86] tuple-based constructor --- ...semidiscretization_hyperbolic_parabolic.jl | 28 ++++++++++++++++--- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index c9849d7d93d..5ff4935b2c2 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -54,16 +54,36 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations end """ - SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, solver; + SemidiscretizationHyperbolicParabolic(mesh, both_equations, initial_condition, solver; source_terms=nothing, - boundary_conditions=boundary_condition_periodic, - boundary_conditions_parabolic=boundary_condition_periodic, + both_boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), RealT=real(solver), uEltype=RealT, - initial_cache=NamedTuple()) + both_initial_caches=(NamedTuple(), NamedTuple())) Construct a semidiscretization of a hyperbolic PDE. """ +function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, + initial_condition, solver; + source_terms=nothing, + boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic), + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT=real(solver), uEltype=RealT, + initial_caches=(NamedTuple(), NamedTuple())) + + hyperbolic_equations, parabolic_equations = equations + hyperbolic_boundary_conditions, parabolic_boundary_conditions = boundary_conditions + initial_hyperbolic_cache, initial_parabolic_cache = initial_caches + + return SemidiscretizationHyperbolicParabolic(mesh, hyperbolic_equations, parabolic_equations, + initial_condition, solver; source_terms, + boundary_conditions=hyperbolic_boundary_conditions, + parabolic_boundary_conditions=parabolic_boundary_conditions, + RealT, uEltype, initial_cache=initial_hyperbolic_cache, + initial_parabolic_cache=initial_parabolic_cache) +end + function SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, solver; source_terms=nothing, From 3d30e54673dcc24b0c189b041353a9a5f2ae06e0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 25 May 2022 13:32:45 -0500 Subject: [PATCH 61/86] updating examples --- examples/dgmulti_2d/elixir_advection_diffusion.jl | 4 ++-- examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index 7b10a3e19d7..d7b081f5031 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -36,8 +36,8 @@ parabolic_boundary_conditions = (; :left => boundary_condition_left, :top => boundary_condition_zero, :right => boundary_condition_neumann_zero) -semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg; - boundary_conditions, parabolic_boundary_conditions) +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, parabolic_equations), initial_condition, dg; + boundary_conditions=(boundary_conditions, parabolic_boundary_conditions)) tspan = (0.0, 1.5) ode = semidiscretize(semi, tspan) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl index a435acb8209..fd4a653959f 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl @@ -14,7 +14,8 @@ end initial_condition = initial_condition_sharp_gaussian mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) -semi = SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, initial_condition, dg) +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, parabolic_equations), + initial_condition, dg) tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) From 29f81e5cace2cfb1a8863332626ffa3b4dbe06c2 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 26 May 2022 01:36:34 -0500 Subject: [PATCH 62/86] add ScalarPlotData2D docstring (#1145) Co-authored-by: Jesse Chan --- src/visualization/types.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/visualization/types.jl b/src/visualization/types.jl index 4b4ab9cfeaf..392f120c461 100644 --- a/src/visualization/types.jl +++ b/src/visualization/types.jl @@ -406,6 +406,12 @@ struct ScalarData{T} data::T end +""" + ScalarPlotData2D(u, semi::AbstractSemidiscretization; kwargs...) + +Returns an `PlotData2DTriangulated` object which is used to visualize a single scalar field. +`u` should be an array whose entries correspond to values of the scalar field at nodal points. +""" ScalarPlotData2D(u, semi::AbstractSemidiscretization; kwargs...) = ScalarPlotData2D(u, mesh_equations_solver_cache(semi)...; kwargs...) From 19910cbb259cd8f43378c1f18de9ebdacac11d64 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 04:07:08 -0500 Subject: [PATCH 63/86] Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha --- .../dgmulti_2d/elixir_advection_diffusion.jl | 19 +++++++++---------- src/equations/laplace_diffusion_2d.jl | 4 ++-- ...semidiscretization_hyperbolic_parabolic.jl | 2 +- src/solvers/dgmulti/dg.jl | 2 +- src/solvers/dgmulti/dg_parabolic.jl | 2 +- test/test_parabolic.jl | 10 +++++----- 6 files changed, 19 insertions(+), 20 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index d7b081f5031..e58cfa36b64 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -1,4 +1,3 @@ - using Trixi, OrdinaryDiffEq dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), @@ -12,15 +11,15 @@ initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVect initial_condition = initial_condition_zero # tag different boundary segments -left(x, tol=50*eps()) = abs(x[1]+1) left, :right => right, :top => top, :bottom => bottom) mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary) # BC types -boundary_condition_left = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + .1 * x[2])) +boundary_condition_left = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.1 * x[2])) boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0)) boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0)) @@ -49,8 +48,8 @@ callbacks = CallbackSet(summary_callback, alive_callback) ############################################################################### # run the simulation -tol = 1e-6 -tsave = LinRange(tspan..., 10) -sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, - saveat=tsave, callback=callbacks) +time_int_tol = 1e-6 +saveat = LinRange(tspan..., 10) +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, save_everystep=false, + saveat=saveat, callback=callbacks) summary_callback() # print the timer summary diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 28b88bcf8f4..99d6e8b6a52 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -1,8 +1,8 @@ @doc raw""" - LaplaceDiffusion2D{T} <: AbstractParabolicEquations{2, 1} + LaplaceDiffusion2D(diffusivity) `LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` -applied to each solution component. +with diffusivity ``\kappa`` applied to each solution component. """ struct LaplaceDiffusion2D{T} <: AbstractLaplaceDiffusionEquations{2, 1} diffusivity::T diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 5ff4935b2c2..c76b1f9f9e3 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -61,7 +61,7 @@ end uEltype=RealT, both_initial_caches=(NamedTuple(), NamedTuple())) -Construct a semidiscretization of a hyperbolic PDE. +Construct a semidiscretization of a hyperbolic-parabolic PDE. """ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, initial_condition, solver; diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index f4bbcc4d6e1..6f2a835dbd4 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -350,7 +350,7 @@ end # Todo: DGMulti. Specialize for modal DG on curved meshes using WADG # inverts Jacobian and scales by -1.0 -function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1.0) +function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1) @threaded for i in each_dof_global(mesh, dg, cache) du[i] *= scaling * cache.invJ[i] end diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index e14066e3f47..276d6543547 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -62,7 +62,7 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, @unpack weak_differentiation_matrices = parabolic_cache - for dim in 1:length(u_grad) + for dim in eachindex(u_grad) reset_du!(u_grad[dim], dg) end diff --git a/test/test_parabolic.jl b/test/test_parabolic.jl index e6bbfeffad1..90aef5facf3 100644 --- a/test/test_parabolic.jl +++ b/test/test_parabolic.jl @@ -25,21 +25,21 @@ isdir(outdir) && rm(outdir, recursive=true) initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2]) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) - equations_parabolic = Trixi.LaplaceDiffusion2D(1.0) + equations_parabolic = LaplaceDiffusion2D(1.0) semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) ode = semidiscretize(semi, (0.0, 0.01)) @unpack cache, parabolic_cache, parabolic_equations = semi @unpack u_grad = parabolic_cache - for dim in 1:length(u_grad) + for dim in eachindex(u_grad) fill!(u_grad[dim], zero(eltype(u_grad[dim]))) end Trixi.calc_gradient!(u_grad, ode.u0, mesh, parabolic_equations, no_boundary_condition, dg, cache, parabolic_cache) @unpack x, y = mesh.md - @test norm(getindex.(u_grad[1], 1) - 2 * x .* y) < 1e3 * eps() - @test norm(getindex.(u_grad[2], 1) - x.^2) < 1e3 * eps() + @test getindex.(u_grad[1], 1) ≈ 2 * x .* y + @test getindex.(u_grad[2], 1) ≈ x.^2 u_flux = similar.(u_grad) Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, parabolic_equations, dg, cache, parabolic_cache) @@ -48,7 +48,7 @@ isdir(outdir) && rm(outdir, recursive=true) du = similar(ode.u0) Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, parabolic_equations, no_boundary_condition, dg, cache, parabolic_cache) - @test norm(getindex.(du, 1) - 2 * y) < 1e3 * eps() + @test getindex.(du, 1) ≈ 2 * y end end From 3ab31623aeb595391d43e4e5f49360f3555dfb6c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 04:50:12 -0500 Subject: [PATCH 64/86] renaming for consistency parabolic_equations -> equations_parabolic no_boundary_condition -> boundary_condition_do_nothing --- .../dgmulti_2d/elixir_advection_diffusion.jl | 8 ++--- .../elixir_advection_diffusion_periodic.jl | 4 +-- src/Trixi.jl | 9 +++-- src/equations/equations.jl | 2 +- ...semidiscretization_hyperbolic_parabolic.jl | 36 +++++++++---------- src/solvers/dgmulti/dg_parabolic.jl | 10 +++--- test/test_parabolic.jl | 8 ++--- 7 files changed, 38 insertions(+), 39 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index e58cfa36b64..d67edb3579f 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -5,7 +5,7 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.5, 1.0) -parabolic_equations = LaplaceDiffusion2D(5e-2) +equations_parabolic = LaplaceDiffusion2D(5e-2) initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) initial_condition = initial_condition_zero @@ -26,8 +26,8 @@ boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> # define inviscid boundary conditions boundary_conditions = (; :left => boundary_condition_left, :bottom => boundary_condition_zero, - :top => no_boundary_condition, - :right => no_boundary_condition) + :top => boundary_condition_do_nothing, + :right => boundary_condition_do_nothing) # define viscous boundary conditions parabolic_boundary_conditions = (; :left => boundary_condition_left, @@ -35,7 +35,7 @@ parabolic_boundary_conditions = (; :left => boundary_condition_left, :top => boundary_condition_zero, :right => boundary_condition_neumann_zero) -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, parabolic_equations), initial_condition, dg; +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; boundary_conditions=(boundary_conditions, parabolic_boundary_conditions)) tspan = (0.0, 1.5) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl index fd4a653959f..59f2c5f02b2 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl @@ -6,7 +6,7 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) -parabolic_equations = LaplaceDiffusion2D(1e-2) +equations_parabolic = LaplaceDiffusion2D(1e-2) function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) return SVector(exp(-100 * (x[1]^2 + x[2]^2))) @@ -14,7 +14,7 @@ end initial_condition = initial_condition_sharp_gaussian mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, parabolic_equations), +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg) tspan = (0.0, 2.0) diff --git a/src/Trixi.jl b/src/Trixi.jl index a2cf650c207..5ee6e7afb95 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -130,13 +130,10 @@ export AcousticPerturbationEquations2D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D, LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D, InviscidBurgersEquation1D, + LaplaceDiffusion2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D -export no_boundary_condition - -export LaplaceDiffusion2D, BoundaryConditionNeumann - export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, @@ -157,8 +154,10 @@ export initial_condition_constant, initial_condition_density_wave, initial_condition_weak_blast_wave -export boundary_condition_periodic, +export boundary_condition_do_nothing, + boundary_condition_periodic, BoundaryConditionDirichlet, + BoundaryConditionNeumann, boundary_condition_noslip_wall, boundary_condition_slip_wall, boundary_condition_wall diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8c3f7d1ff01..e1a530dc40d 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -176,7 +176,7 @@ end end # Imposing no boundary condition just evaluates the flux at the inner state. -@inline function no_boundary_condition(flux, other_args...) +@inline function boundary_condition_do_nothing(flux, other_args...) return flux end diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index c76b1f9f9e3..99d637563ff 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -18,7 +18,7 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations mesh::Mesh equations::Equations - parabolic_equations::ParabolicEquations + equations_parabolic::ParabolicEquations # This guy is a bit messy since we abuse it as some kind of "exact solution" # although this doesn't really exist... @@ -38,16 +38,16 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations performance_counter::PerformanceCounter function SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations, InitialCondition, BoundaryConditions, ParabolicBoundaryConditions, SourceTerms, Solver, Cache, ParabolicCache}( - mesh::Mesh, equations::Equations, parabolic_equations::ParabolicEquations, initial_condition::InitialCondition, + mesh::Mesh, equations::Equations, equations_parabolic::ParabolicEquations, initial_condition::InitialCondition, boundary_conditions::BoundaryConditions, parabolic_boundary_conditions::ParabolicBoundaryConditions, source_terms::SourceTerms, solver::Solver, cache::Cache, parabolic_cache::ParabolicCache) where {Mesh, Equations, ParabolicEquations, InitialCondition, BoundaryConditions, ParabolicBoundaryConditions, SourceTerms, Solver, Cache, ParabolicCache} @assert ndims(mesh) == ndims(equations) - # Todo: assert nvariables(equations)==nvariables(parabolic_equations) + # Todo: assert nvariables(equations)==nvariables(equations_parabolic) performance_counter = PerformanceCounter() - new(mesh, equations, parabolic_equations, initial_condition, + new(mesh, equations, equations_parabolic, initial_condition, boundary_conditions, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache, performance_counter) end @@ -72,11 +72,11 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, RealT=real(solver), uEltype=RealT, initial_caches=(NamedTuple(), NamedTuple())) - hyperbolic_equations, parabolic_equations = equations + hyperbolic_equations, equations_parabolic = equations hyperbolic_boundary_conditions, parabolic_boundary_conditions = boundary_conditions initial_hyperbolic_cache, initial_parabolic_cache = initial_caches - return SemidiscretizationHyperbolicParabolic(mesh, hyperbolic_equations, parabolic_equations, + return SemidiscretizationHyperbolicParabolic(mesh, hyperbolic_equations, equations_parabolic, initial_condition, solver; source_terms, boundary_conditions=hyperbolic_boundary_conditions, parabolic_boundary_conditions=parabolic_boundary_conditions, @@ -84,7 +84,7 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, initial_parabolic_cache=initial_parabolic_cache) end -function SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equations, +function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, solver; source_terms=nothing, boundary_conditions=boundary_condition_periodic, @@ -99,13 +99,13 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, parabolic_equati _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) _parabolic_boundary_conditions = digest_boundary_conditions(parabolic_boundary_conditions, mesh, solver, cache) - parabolic_cache = (; create_cache(mesh, parabolic_equations, solver, RealT, uEltype)..., + parabolic_cache = (; create_cache(mesh, equations_parabolic, solver, RealT, uEltype)..., initial_parabolic_cache...) - SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(parabolic_equations), + SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), typeof(initial_condition), typeof(_boundary_conditions), typeof(_parabolic_boundary_conditions), typeof(source_terms), typeof(solver), typeof(cache), typeof(parabolic_cache)}( - mesh, equations, parabolic_equations, initial_condition, + mesh, equations, equations_parabolic, initial_condition, _boundary_conditions, _parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache) end @@ -119,7 +119,7 @@ end function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.solver), mesh=semi.mesh, equations=semi.equations, - parabolic_equations=semi.parabolic_equations, + equations_parabolic=semi.equations_parabolic, initial_condition=semi.initial_condition, solver=semi.solver, source_terms=semi.source_terms, @@ -130,7 +130,7 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s # special care if shock-capturing volume integrals are used (because of # the indicators and their own caches...). SemidiscretizationHyperbolicParabolic( - mesh, equations, parabolic_equations, initial_condition, solver; source_terms, boundary_conditions, parabolic_boundary_conditions, uEltype) + mesh, equations, equations_parabolic, initial_condition, solver; source_terms, boundary_conditions, parabolic_boundary_conditions, uEltype) end function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) @@ -139,7 +139,7 @@ function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) print(io, "SemidiscretizationHyperbolicParabolic(") print(io, semi.mesh) print(io, ", ", semi.equations) - print(io, ", ", semi.parabolic_equations) + print(io, ", ", semi.equations_parabolic) print(io, ", ", semi.initial_condition) print(io, ", ", semi.boundary_conditions) print(io, ", ", semi.parabolic_boundary_conditions) @@ -163,7 +163,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli summary_line(io, "#spatial dimensions", ndims(semi.equations)) summary_line(io, "mesh", semi.mesh) summary_line(io, "equations", semi.equations |> typeof |> nameof) - summary_line(io, "parabolic_equations", semi.parabolic_equations |> typeof |> nameof) + summary_line(io, "equations_parabolic", semi.equations_parabolic |> typeof |> nameof) summary_line(io, "initial condition", semi.initial_condition) # print_boundary_conditions(io, semi) @@ -236,14 +236,14 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) end function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) - @unpack mesh, parabolic_equations, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache = semi + @unpack mesh, equations_parabolic, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache = semi - u = wrap_array(u_ode, mesh, parabolic_equations, solver, parabolic_cache) - du = wrap_array(du_ode, mesh, parabolic_equations, solver, parabolic_cache) + u = wrap_array(u_ode, mesh, equations_parabolic, solver, parabolic_cache) + du = wrap_array(du_ode, mesh, equations_parabolic, solver, parabolic_cache) # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() - @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, parabolic_equations, initial_condition, + @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, equations_parabolic, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache) runtime = time_ns() - time_start diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 276d6543547..93edba06da9 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -252,22 +252,22 @@ end # 2. compute f(u, grad(u)) # 3. compute div(u) # boundary conditions will be applied to both grad(u) and div(u). -function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, parabolic_equations::AbstractParabolicEquations, +function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::AbstractParabolicEquations, initial_condition, boundary_conditions, source_terms, dg::DGMulti, cache, parabolic_cache) reset_du!(du, dg) @unpack u_transformed, u_grad, viscous_flux = parabolic_cache - transform_variables!(u_transformed, u, parabolic_equations) + transform_variables!(u_transformed, u, equations_parabolic) - calc_gradient!(u_grad, u_transformed, t, mesh, parabolic_equations, + calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, boundary_conditions, dg, cache, parabolic_cache) calc_viscous_fluxes!(viscous_flux, u_transformed, u_grad, - mesh, parabolic_equations, dg, cache, parabolic_cache) + mesh, equations_parabolic, dg, cache, parabolic_cache) - calc_divergence!(du, u_transformed, t, viscous_flux, mesh, parabolic_equations, + calc_divergence!(du, u_transformed, t, viscous_flux, mesh, equations_parabolic, boundary_conditions, dg, cache, parabolic_cache) return nothing diff --git a/test/test_parabolic.jl b/test/test_parabolic.jl index 90aef5facf3..1e95612151e 100644 --- a/test/test_parabolic.jl +++ b/test/test_parabolic.jl @@ -30,24 +30,24 @@ isdir(outdir) && rm(outdir, recursive=true) semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) ode = semidiscretize(semi, (0.0, 0.01)) - @unpack cache, parabolic_cache, parabolic_equations = semi + @unpack cache, parabolic_cache, equations_parabolic = semi @unpack u_grad = parabolic_cache for dim in eachindex(u_grad) fill!(u_grad[dim], zero(eltype(u_grad[dim]))) end - Trixi.calc_gradient!(u_grad, ode.u0, mesh, parabolic_equations, no_boundary_condition, dg, cache, parabolic_cache) + Trixi.calc_gradient!(u_grad, ode.u0, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, parabolic_cache) @unpack x, y = mesh.md @test getindex.(u_grad[1], 1) ≈ 2 * x .* y @test getindex.(u_grad[2], 1) ≈ x.^2 u_flux = similar.(u_grad) - Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, parabolic_equations, dg, cache, parabolic_cache) + Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, equations_parabolic, dg, cache, parabolic_cache) @test u_flux[1] ≈ u_grad[1] @test u_flux[2] ≈ u_grad[2] du = similar(ode.u0) - Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, parabolic_equations, no_boundary_condition, dg, cache, parabolic_cache) + Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, parabolic_cache) @test getindex.(du, 1) ≈ 2 * y end From a61648e2666bad5b940033167a4b903e90deccd9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 05:40:53 -0500 Subject: [PATCH 65/86] move Neumann BCs into equations.jl --- src/equations/equations.jl | 19 +++++++++++++++++-- src/equations/laplace_diffusion_2d.jl | 5 ----- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index e1a530dc40d..23470162618 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -175,6 +175,21 @@ end return flux end +""" + BoundaryConditionNeumann(boundary_value_function) + +Similar to `BoundaryConditionDirichlet`, but creates a Neumann boundary condition for parabolic +equations that uses the function `boundary_value_function` to specify the values at the boundary. +The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as +```julia +boundary_value_function(x, t, equations) +``` +where `x` specifies the coordinates, `t` is the current time, and `equation` is the corresponding system of equations. +""" +struct BoundaryConditionNeumann{B} + boundary_value_function::B +end + # Imposing no boundary condition just evaluates the flux at the inner state. @inline function boundary_condition_do_nothing(flux, other_args...) return flux @@ -333,14 +348,14 @@ include("lattice_boltzmann_3d.jl") abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("acoustic_perturbation_2d.jl") -abstract type AbstractParabolicEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +abstract type AbstractEquationsParabolic{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end # operator types used for dispatch on parabolic boundary fluxes struct Gradient end struct Divergence end # Linear scalar diffusion for use in linear scalar advection-diffusion problems -abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractParabolicEquations{NDIMS, NVARS} end +abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("laplace_diffusion_2d.jl") end # @muladd diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 99d6e8b6a52..3f3756ad62f 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -27,11 +27,6 @@ end return flux_inner end -# `boundary_value_function` should have signature `boundary_value_function(x, t, equations)` -# For Neumann BCs, this corresponds to kappa * (∇u ⋅ n) -struct BoundaryConditionNeumann{B} - boundary_value_function::B -end @inline function (boundary_condition::BoundaryConditionNeumann)(u_inner, normal::AbstractVector, x, t, operator_type::Gradient, From b4d7f6021abfe9ea548c8afea4d0c14a2f84e54f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 05:41:26 -0500 Subject: [PATCH 66/86] renaming, e.g., ParabolicEquations -> EquationsParabolic, etc --- ...semidiscretization_hyperbolic_parabolic.jl | 22 +++++++++---------- src/solvers/dgmulti/dg_parabolic.jl | 16 +++++++------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 99d637563ff..97a2628e110 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -11,36 +11,36 @@ A struct containing everything needed to describe a spatial semidiscretization of a mixed hyperbolic-parabolic conservation law. """ -struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations, InitialCondition, - BoundaryConditions, ParabolicBoundaryConditions, - SourceTerms, Solver, Cache, ParabolicCache} <: AbstractSemidiscretization +struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, + BoundaryConditions, BoundaryConditionsParabolic, + SourceTerms, Solver, Cache, CacheParabolic} <: AbstractSemidiscretization mesh::Mesh equations::Equations - equations_parabolic::ParabolicEquations + equations_parabolic::EquationsParabolic # 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 - parabolic_boundary_conditions::ParabolicBoundaryConditions + parabolic_boundary_conditions::BoundaryConditionsParabolic source_terms::SourceTerms solver::Solver - # TODO: do we want to introduce `parabolic_solver` for future specialization? + # TODO: introduce `solver_parabolic` for future specialization. cache::Cache - parabolic_cache::ParabolicCache + parabolic_cache::CacheParabolic performance_counter::PerformanceCounter - function SemidiscretizationHyperbolicParabolic{Mesh, Equations, ParabolicEquations, InitialCondition, BoundaryConditions, ParabolicBoundaryConditions, SourceTerms, Solver, Cache, ParabolicCache}( - mesh::Mesh, equations::Equations, equations_parabolic::ParabolicEquations, initial_condition::InitialCondition, - boundary_conditions::BoundaryConditions, parabolic_boundary_conditions::ParabolicBoundaryConditions, - source_terms::SourceTerms, solver::Solver, cache::Cache, parabolic_cache::ParabolicCache) where {Mesh, Equations, ParabolicEquations, InitialCondition, BoundaryConditions, ParabolicBoundaryConditions, SourceTerms, Solver, Cache, ParabolicCache} + function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic}( + mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition, + boundary_conditions::BoundaryConditions, parabolic_boundary_conditions::BoundaryConditionsParabolic, + source_terms::SourceTerms, solver::Solver, cache::Cache, parabolic_cache::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic} @assert ndims(mesh) == ndims(equations) # Todo: assert nvariables(equations)==nvariables(equations_parabolic) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 93edba06da9..11d6230d7ab 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -1,4 +1,4 @@ -function create_cache(mesh::DGMultiMesh, equations::AbstractParabolicEquations, +function create_cache(mesh::DGMultiMesh, equations::AbstractEquationsParabolic, dg::DGMultiWeakForm, RealT, uEltype) nvars = nvariables(equations) @@ -35,13 +35,13 @@ function transform_variables!(u_transformed, u, equations) end # interpolates from solution coefficients to face quadrature points -function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::AbstractParabolicEquations, +function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, surface_integral, dg::DGMulti, cache) apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) end function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, - mesh, equations::AbstractParabolicEquations, + mesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, parabolic_cache) @unpack local_flux_face_values_threaded = parabolic_cache @threaded for e in eachelement(mesh, dg) @@ -57,7 +57,7 @@ function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, end function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, - equations::AbstractParabolicEquations, + equations::AbstractEquationsParabolic, boundary_conditions, dg::DGMulti, cache, parabolic_cache) @unpack weak_differentiation_matrices = parabolic_cache @@ -102,7 +102,7 @@ end # do nothing for periodic domains function calc_boundary_flux!(flux, u, t, operator_type, ::BoundaryConditionPeriodic, - mesh, equations::AbstractParabolicEquations, dg::DGMulti, + mesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, parabolic_cache) return nothing end @@ -155,7 +155,7 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, end function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, - equations::AbstractParabolicEquations, + equations::AbstractEquationsParabolic, dg::DGMulti, cache, parabolic_cache) for dim in eachdim(mesh) @@ -196,7 +196,7 @@ function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, end function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh, - equations::AbstractParabolicEquations, + equations::AbstractEquationsParabolic, boundary_conditions, dg::DGMulti, cache, parabolic_cache) @unpack weak_differentiation_matrices = parabolic_cache @@ -252,7 +252,7 @@ end # 2. compute f(u, grad(u)) # 3. compute div(u) # boundary conditions will be applied to both grad(u) and div(u). -function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::AbstractParabolicEquations, +function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions, source_terms, dg::DGMulti, cache, parabolic_cache) From a49635ef8ba15cf47a056300b01fe484a3fe6952 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 05:41:33 -0500 Subject: [PATCH 67/86] rename test module --- test/test_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parabolic.jl b/test/test_parabolic.jl index 1e95612151e..4f2c4dea476 100644 --- a/test/test_parabolic.jl +++ b/test/test_parabolic.jl @@ -1,4 +1,4 @@ -module TestExamplesDGMulti2D +module TestExamplesParabolic2D using Test using Trixi From 3573cae71db58f74669fdb28fd01b2ee737d895d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 05:42:07 -0500 Subject: [PATCH 68/86] generalize LaplaceDiffusion2D to any number of variables --- src/equations/laplace_diffusion_2d.jl | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 3f3756ad62f..d9ed13ee11e 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -1,17 +1,29 @@ @doc raw""" - LaplaceDiffusion2D(diffusivity) + LaplaceDiffusion2D(diffusivity, equations) `LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` -with diffusivity ``\kappa`` applied to each solution component. +with diffusivity ``\kappa`` applied to each solution component defined by `equations`. + +If `equations` is not specified, a scalar `LaplaceDiffusion2D` is returned. """ -struct LaplaceDiffusion2D{T} <: AbstractLaplaceDiffusionEquations{2, 1} +struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N} diffusivity::T + equations::E +end + +function LaplaceDiffusion2D(diffusivity, equations=nothing) + if equations isa AbstractEquations + nvars = nvariables(equations) + else # assume scalar diffusion if no equations are specified + nvars = 1 + end + return LaplaceDiffusion2D{typeof(equations), nvars, typeof(diffusivity)}(diffusivity, equations) end # no orientation specified since the flux is vector-valued function flux(u, grad_u, equations::LaplaceDiffusion2D) dudx, dudy = grad_u - return equations.diffusivity * dudx, equations.diffusivity * dudy + return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy) end # Dirichlet-type boundary condition for use with a parabolic solver in weak form From 8bb3a73dd1738b0da17e4b846c1edfc1dcf149f9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 05:47:37 -0500 Subject: [PATCH 69/86] removing saveat times and standardizing `tol` names --- examples/dgmulti_2d/elixir_advection_diffusion.jl | 5 ++--- examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl | 6 ++++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index d67edb3579f..a68866a1859 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -49,7 +49,6 @@ callbacks = CallbackSet(summary_callback, alive_callback) # run the simulation time_int_tol = 1e-6 -saveat = LinRange(tspan..., 10) -sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, save_everystep=false, - saveat=saveat, callback=callbacks) +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) summary_callback() # print the timer summary diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl index 59f2c5f02b2..4f6d9fab580 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl @@ -29,6 +29,8 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) ############################################################################### # run the simulation -tol = 1e-6 -sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, callback=callbacks) +time_int_tol = 1e-6 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) + summary_callback() # print the timer summary From 26bbabaf3c259fc158732278250bf6a2588b4ddc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 05:53:26 -0500 Subject: [PATCH 70/86] parabolic_cache -> cache_parabolic cache --- ...semidiscretization_hyperbolic_parabolic.jl | 28 +++++----- src/solvers/dgmulti/dg_parabolic.jl | 54 +++++++++---------- test/test_parabolic.jl | 10 ++-- 3 files changed, 46 insertions(+), 46 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 97a2628e110..bf222bd7830 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -33,14 +33,14 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic # TODO: introduce `solver_parabolic` for future specialization. cache::Cache - parabolic_cache::CacheParabolic + cache_parabolic::CacheParabolic performance_counter::PerformanceCounter function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic}( mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition, boundary_conditions::BoundaryConditions, parabolic_boundary_conditions::BoundaryConditionsParabolic, - source_terms::SourceTerms, solver::Solver, cache::Cache, parabolic_cache::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic} + source_terms::SourceTerms, solver::Solver, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic} @assert ndims(mesh) == ndims(equations) # Todo: assert nvariables(equations)==nvariables(equations_parabolic) @@ -49,7 +49,7 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic new(mesh, equations, equations_parabolic, initial_condition, boundary_conditions, parabolic_boundary_conditions, - source_terms, solver, cache, parabolic_cache, performance_counter) + source_terms, solver, cache, cache_parabolic, performance_counter) end end @@ -74,14 +74,14 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, hyperbolic_equations, equations_parabolic = equations hyperbolic_boundary_conditions, parabolic_boundary_conditions = boundary_conditions - initial_hyperbolic_cache, initial_parabolic_cache = initial_caches + initial_hyperbolic_cache, initial_cache_parabolic = initial_caches return SemidiscretizationHyperbolicParabolic(mesh, hyperbolic_equations, equations_parabolic, initial_condition, solver; source_terms, boundary_conditions=hyperbolic_boundary_conditions, parabolic_boundary_conditions=parabolic_boundary_conditions, RealT, uEltype, initial_cache=initial_hyperbolic_cache, - initial_parabolic_cache=initial_parabolic_cache) + initial_cache_parabolic=initial_cache_parabolic) end function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, @@ -93,21 +93,21 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo # while `uEltype` is used as element type of solutions etc. RealT=real(solver), uEltype=RealT, initial_cache=NamedTuple(), - initial_parabolic_cache=NamedTuple()) + initial_cache_parabolic=NamedTuple()) cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...) _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) _parabolic_boundary_conditions = digest_boundary_conditions(parabolic_boundary_conditions, mesh, solver, cache) - parabolic_cache = (; create_cache(mesh, equations_parabolic, solver, RealT, uEltype)..., - initial_parabolic_cache...) + cache_parabolic = (; create_cache(mesh, equations_parabolic, solver, RealT, uEltype)..., + initial_cache_parabolic...) SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), typeof(initial_condition), typeof(_boundary_conditions), typeof(_parabolic_boundary_conditions), - typeof(source_terms), typeof(solver), typeof(cache), typeof(parabolic_cache)}( + typeof(source_terms), typeof(solver), typeof(cache), typeof(cache_parabolic)}( mesh, equations, equations_parabolic, initial_condition, _boundary_conditions, _parabolic_boundary_conditions, source_terms, - solver, cache, parabolic_cache) + solver, cache, cache_parabolic) end @@ -236,16 +236,16 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) end function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) - @unpack mesh, equations_parabolic, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, parabolic_cache = semi + @unpack mesh, equations_parabolic, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, cache_parabolic = semi - u = wrap_array(u_ode, mesh, equations_parabolic, solver, parabolic_cache) - du = wrap_array(du_ode, mesh, equations_parabolic, solver, parabolic_cache) + u = wrap_array(u_ode, mesh, equations_parabolic, solver, cache_parabolic) + du = wrap_array(du_ode, mesh, equations_parabolic, solver, cache_parabolic) # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, equations_parabolic, initial_condition, parabolic_boundary_conditions, source_terms, - solver, cache, parabolic_cache) + solver, cache, cache_parabolic) runtime = time_ns() - time_start put!(semi.performance_counter, runtime) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 11d6230d7ab..d098236b536 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -42,8 +42,8 @@ end function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, mesh, equations::AbstractEquationsParabolic, - dg::DGMulti, cache, parabolic_cache) - @unpack local_flux_face_values_threaded = parabolic_cache + dg::DGMulti, cache, cache_parabolic) + @unpack local_flux_face_values_threaded = cache_parabolic @threaded for e in eachelement(mesh, dg) local_flux_values = local_flux_face_values_threaded[Threads.threadid()] for dim in eachdim(mesh) @@ -58,9 +58,9 @@ end function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, - boundary_conditions, dg::DGMulti, cache, parabolic_cache) + boundary_conditions, dg::DGMulti, cache, cache_parabolic) - @unpack weak_differentiation_matrices = parabolic_cache + @unpack weak_differentiation_matrices = cache_parabolic for dim in eachindex(u_grad) reset_du!(u_grad[dim], dg) @@ -75,11 +75,11 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, end end - @unpack u_face_values = parabolic_cache + @unpack u_face_values = cache_parabolic prolong2interfaces!(u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) # compute fluxes at interfaces - @unpack scalar_flux_face_values = parabolic_cache + @unpack scalar_flux_face_values = cache_parabolic @unpack mapM, mapP, Jf = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] @@ -89,11 +89,11 @@ function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh, end calc_boundary_flux!(scalar_flux_face_values, nothing, t, Gradient(), boundary_conditions, - mesh, equations, dg, cache, parabolic_cache) + mesh, equations, dg, cache, cache_parabolic) # compute surface contributions calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values, - mesh, equations, dg, cache, parabolic_cache) + mesh, equations, dg, cache, cache_parabolic) for dim in eachdim(mesh) invert_jacobian!(u_grad[dim], mesh, equations, dg, cache; scaling=1.0) @@ -103,31 +103,31 @@ end # do nothing for periodic domains function calc_boundary_flux!(flux, u, t, operator_type, ::BoundaryConditionPeriodic, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, - cache, parabolic_cache) + cache, cache_parabolic) return nothing end # "lispy tuple programming" instead of for loop for type stability function calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions, - mesh, equations, dg::DGMulti, cache, parabolic_cache) + mesh, equations, dg::DGMulti, cache, cache_parabolic) # peel off first boundary condition calc_single_boundary_flux!(flux, u, t, operator_type, first(boundary_conditions), first(keys(boundary_conditions)), - mesh, equations, dg, cache, parabolic_cache) + mesh, equations, dg, cache, cache_parabolic) # recurse on the remainder of the boundary conditions calc_boundary_flux!(flux, u, t, operator_type, Base.tail(boundary_conditions), - mesh, equations, dg, cache, parabolic_cache) + mesh, equations, dg, cache, cache_parabolic) end # terminate recursion calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{(),Tuple{}}, - mesh, equations, dg::DGMulti, cache, parabolic_cache) = nothing + mesh, equations, dg::DGMulti, cache, cache_parabolic) = nothing # TODO: DGMulti. Decide if we want to use the input `u_face_values` (currently unused) function calc_single_boundary_flux!(flux_face_values, u_face_values, t, operator_type, boundary_condition, boundary_key, - mesh, equations, dg::DGMulti{NDIMS}, cache, parabolic_cache) where {NDIMS} + mesh, equations, dg::DGMulti{NDIMS}, cache, cache_parabolic) where {NDIMS} rd = dg.basis md = mesh.md @@ -156,13 +156,13 @@ end function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, - dg::DGMulti, cache, parabolic_cache) + dg::DGMulti, cache, cache_parabolic) for dim in eachdim(mesh) reset_du!(viscous_flux[dim], dg) end - @unpack local_viscous_flux_threaded, local_u_values_threaded = parabolic_cache + @unpack local_viscous_flux_threaded, local_u_values_threaded = cache_parabolic @threaded for e in eachelement(mesh, dg) @@ -197,9 +197,9 @@ end function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, - boundary_conditions, dg::DGMulti, cache, parabolic_cache) + boundary_conditions, dg::DGMulti, cache, cache_parabolic) - @unpack weak_differentiation_matrices = parabolic_cache + @unpack weak_differentiation_matrices = cache_parabolic reset_du!(du, dg) @@ -213,16 +213,16 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh end # interpolates from solution coefficients to face quadrature points - viscous_flux_face_values = parabolic_cache.grad_u_face_values + viscous_flux_face_values = cache_parabolic.grad_u_face_values for dim in eachdim(mesh) prolong2interfaces!(viscous_flux_face_values[dim], viscous_flux[dim], mesh, equations, dg.surface_integral, dg, cache) end # compute fluxes at interfaces - @unpack scalar_flux_face_values = parabolic_cache + @unpack scalar_flux_face_values = cache_parabolic @unpack mapM, mapP, nxyzJ = mesh.md - @threaded for face_node_index in each_face_node_global(mesh, dg, cache, parabolic_cache) + @threaded for face_node_index in each_face_node_global(mesh, dg, cache, cache_parabolic) idM, idP = mapM[face_node_index], mapP[face_node_index] # compute f(u, ∇u) ⋅ n @@ -238,7 +238,7 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh # TODO: decide what to pass in calc_boundary_flux!(scalar_flux_face_values, nothing, t, Divergence(), - boundary_conditions, mesh, equations, dg, cache, parabolic_cache) + boundary_conditions, mesh, equations, dg, cache, cache_parabolic) # surface contributions apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, scalar_flux_face_values) @@ -254,21 +254,21 @@ end # boundary conditions will be applied to both grad(u) and div(u). function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions, source_terms, - dg::DGMulti, cache, parabolic_cache) + dg::DGMulti, cache, cache_parabolic) reset_du!(du, dg) - @unpack u_transformed, u_grad, viscous_flux = parabolic_cache + @unpack u_transformed, u_grad, viscous_flux = cache_parabolic transform_variables!(u_transformed, u, equations_parabolic) calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic, - boundary_conditions, dg, cache, parabolic_cache) + boundary_conditions, dg, cache, cache_parabolic) calc_viscous_fluxes!(viscous_flux, u_transformed, u_grad, - mesh, equations_parabolic, dg, cache, parabolic_cache) + mesh, equations_parabolic, dg, cache, cache_parabolic) calc_divergence!(du, u_transformed, t, viscous_flux, mesh, equations_parabolic, - boundary_conditions, dg, cache, parabolic_cache) + boundary_conditions, dg, cache, cache_parabolic) return nothing diff --git a/test/test_parabolic.jl b/test/test_parabolic.jl index 4f2c4dea476..c06619b26c5 100644 --- a/test/test_parabolic.jl +++ b/test/test_parabolic.jl @@ -30,24 +30,24 @@ isdir(outdir) && rm(outdir, recursive=true) semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) ode = semidiscretize(semi, (0.0, 0.01)) - @unpack cache, parabolic_cache, equations_parabolic = semi - @unpack u_grad = parabolic_cache + @unpack cache, cache_parabolic, equations_parabolic = semi + @unpack u_grad = cache_parabolic for dim in eachindex(u_grad) fill!(u_grad[dim], zero(eltype(u_grad[dim]))) end - Trixi.calc_gradient!(u_grad, ode.u0, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, parabolic_cache) + Trixi.calc_gradient!(u_grad, ode.u0, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, cache_parabolic) @unpack x, y = mesh.md @test getindex.(u_grad[1], 1) ≈ 2 * x .* y @test getindex.(u_grad[2], 1) ≈ x.^2 u_flux = similar.(u_grad) - Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, equations_parabolic, dg, cache, parabolic_cache) + Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, equations_parabolic, dg, cache, cache_parabolic) @test u_flux[1] ≈ u_grad[1] @test u_flux[2] ≈ u_grad[2] du = similar(ode.u0) - Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, parabolic_cache) + Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, cache_parabolic) @test getindex.(du, 1) ≈ 2 * y end From 24eb1b31b3ab4d05aefc2eccdab518369440dc50 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:00:12 -0500 Subject: [PATCH 71/86] renaming `create_cache` -> `create_cache_parabolic` avoids conflicts when using dg::DGMultiFluxDiff --- .../semidiscretization_hyperbolic_parabolic.jl | 4 ++-- src/solvers/dgmulti/dg_parabolic.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index bf222bd7830..56013f9d16e 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -99,8 +99,8 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) _parabolic_boundary_conditions = digest_boundary_conditions(parabolic_boundary_conditions, mesh, solver, cache) - cache_parabolic = (; create_cache(mesh, equations_parabolic, solver, RealT, uEltype)..., - initial_cache_parabolic...) + cache_parabolic = (; create_cache_parabolic(mesh, equations_parabolic, solver, RealT, uEltype)..., + initial_cache_parabolic...) SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), typeof(initial_condition), typeof(_boundary_conditions), typeof(_parabolic_boundary_conditions), diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index d098236b536..f4c4799e021 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -1,5 +1,5 @@ -function create_cache(mesh::DGMultiMesh, equations::AbstractEquationsParabolic, - dg::DGMultiWeakForm, RealT, uEltype) +function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsParabolic, + dg::DGMulti, RealT, uEltype) nvars = nvariables(equations) @unpack M, Drst = dg.basis From 05c8b0cdb660c2456e038d9187c201cbfd82a793 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:05:47 -0500 Subject: [PATCH 72/86] Update examples/dgmulti_2d/elixir_advection_diffusion.jl Co-authored-by: Michael Schlottke-Lakemper --- examples/dgmulti_2d/elixir_advection_diffusion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index a68866a1859..c6ee7c13013 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -36,7 +36,7 @@ parabolic_boundary_conditions = (; :left => boundary_condition_left, :right => boundary_condition_neumann_zero) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; - boundary_conditions=(boundary_conditions, parabolic_boundary_conditions)) + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic)) tspan = (0.0, 1.5) ode = semidiscretize(semi, tspan) From 4f1d3e03fd6702e91a6572aabe2b2c688d6f7df0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:07:46 -0500 Subject: [PATCH 73/86] Update src/equations/equations.jl Co-authored-by: Michael Schlottke-Lakemper --- src/equations/equations.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 23470162618..7cb4b36bba6 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -191,8 +191,8 @@ struct BoundaryConditionNeumann{B} end # Imposing no boundary condition just evaluates the flux at the inner state. -@inline function boundary_condition_do_nothing(flux, other_args...) - return flux +@inline function boundary_condition_do_nothing(inner_flux_or_state, other_args...) + return inner_flux_or_state end # set sensible default values that may be overwritten by specific equations From d5e032c954b77590d75ee9702647b5d32151c9fd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:13:50 -0500 Subject: [PATCH 74/86] Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper --- .../semidiscretization_hyperbolic_parabolic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 56013f9d16e..b21064874f4 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -162,8 +162,8 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli summary_header(io, "SemidiscretizationHyperbolicParabolic") summary_line(io, "#spatial dimensions", ndims(semi.equations)) summary_line(io, "mesh", semi.mesh) - summary_line(io, "equations", semi.equations |> typeof |> nameof) - summary_line(io, "equations_parabolic", semi.equations_parabolic |> typeof |> nameof) + summary_line(io, "hyperbolic equations", semi.equations |> typeof |> nameof) + summary_line(io, "parabolic equations", semi.equations_parabolic |> typeof |> nameof) summary_line(io, "initial condition", semi.initial_condition) # print_boundary_conditions(io, semi) From 270f7ceab835bfd0ed54ba0ac609566fd6aeb3a3 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:15:16 -0500 Subject: [PATCH 75/86] parabolic_boundary_conditions -> boundary_conditions_parabolic similarly for hyperbolic d --- .../dgmulti_2d/elixir_advection_diffusion.jl | 4 +-- ...semidiscretization_hyperbolic_parabolic.jl | 34 +++++++++---------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index a68866a1859..61d5a684306 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -30,13 +30,13 @@ boundary_conditions = (; :left => boundary_condition_left, :right => boundary_condition_do_nothing) # define viscous boundary conditions -parabolic_boundary_conditions = (; :left => boundary_condition_left, +boundary_conditions_parabolic = (; :left => boundary_condition_left, :bottom => boundary_condition_zero, :top => boundary_condition_zero, :right => boundary_condition_neumann_zero) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; - boundary_conditions=(boundary_conditions, parabolic_boundary_conditions)) + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic)) tspan = (0.0, 1.5) ode = semidiscretize(semi, tspan) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 56013f9d16e..26bcd16e2f6 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -25,7 +25,7 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic initial_condition::InitialCondition boundary_conditions::BoundaryConditions - parabolic_boundary_conditions::BoundaryConditionsParabolic + boundary_conditions_parabolic::BoundaryConditionsParabolic source_terms::SourceTerms @@ -39,7 +39,7 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic}( mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition, - boundary_conditions::BoundaryConditions, parabolic_boundary_conditions::BoundaryConditionsParabolic, + boundary_conditions::BoundaryConditions, boundary_conditions_parabolic::BoundaryConditionsParabolic, source_terms::SourceTerms, solver::Solver, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic} @assert ndims(mesh) == ndims(equations) @@ -48,7 +48,7 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic performance_counter = PerformanceCounter() new(mesh, equations, equations_parabolic, initial_condition, - boundary_conditions, parabolic_boundary_conditions, + boundary_conditions, boundary_conditions_parabolic, source_terms, solver, cache, cache_parabolic, performance_counter) end end @@ -72,14 +72,14 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, RealT=real(solver), uEltype=RealT, initial_caches=(NamedTuple(), NamedTuple())) - hyperbolic_equations, equations_parabolic = equations - hyperbolic_boundary_conditions, parabolic_boundary_conditions = boundary_conditions + equations_hyperbolic, equations_parabolic = equations + boundary_conditions_hyperbolic, boundary_conditions_parabolic = boundary_conditions initial_hyperbolic_cache, initial_cache_parabolic = initial_caches - return SemidiscretizationHyperbolicParabolic(mesh, hyperbolic_equations, equations_parabolic, + return SemidiscretizationHyperbolicParabolic(mesh, equations_hyperbolic, equations_parabolic, initial_condition, solver; source_terms, - boundary_conditions=hyperbolic_boundary_conditions, - parabolic_boundary_conditions=parabolic_boundary_conditions, + boundary_conditions=boundary_conditions_hyperbolic, + boundary_conditions_parabolic=boundary_conditions_parabolic, RealT, uEltype, initial_cache=initial_hyperbolic_cache, initial_cache_parabolic=initial_cache_parabolic) end @@ -88,7 +88,7 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo initial_condition, solver; source_terms=nothing, boundary_conditions=boundary_condition_periodic, - parabolic_boundary_conditions=boundary_condition_periodic, + boundary_conditions_parabolic=boundary_condition_periodic, # `RealT` is used as real type for node locations etc. # while `uEltype` is used as element type of solutions etc. RealT=real(solver), uEltype=RealT, @@ -97,16 +97,16 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...) _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) - _parabolic_boundary_conditions = digest_boundary_conditions(parabolic_boundary_conditions, mesh, solver, cache) + _boundary_conditions_parabolic = digest_boundary_conditions(boundary_conditions_parabolic, mesh, solver, cache) cache_parabolic = (; create_cache_parabolic(mesh, equations_parabolic, solver, RealT, uEltype)..., initial_cache_parabolic...) SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic), - typeof(initial_condition), typeof(_boundary_conditions), typeof(_parabolic_boundary_conditions), + typeof(initial_condition), typeof(_boundary_conditions), typeof(_boundary_conditions_parabolic), typeof(source_terms), typeof(solver), typeof(cache), typeof(cache_parabolic)}( mesh, equations, equations_parabolic, initial_condition, - _boundary_conditions, _parabolic_boundary_conditions, source_terms, + _boundary_conditions, _boundary_conditions_parabolic, source_terms, solver, cache, cache_parabolic) end @@ -124,13 +124,13 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s solver=semi.solver, source_terms=semi.source_terms, boundary_conditions=semi.boundary_conditions, - parabolic_boundary_conditions=semi.parabolic_boundary_conditions + boundary_conditions_parabolic=semi.boundary_conditions_parabolic ) # TODO: Which parts do we want to `remake`? At least the solver needs some # special care if shock-capturing volume integrals are used (because of # the indicators and their own caches...). SemidiscretizationHyperbolicParabolic( - mesh, equations, equations_parabolic, initial_condition, solver; source_terms, boundary_conditions, parabolic_boundary_conditions, uEltype) + mesh, equations, equations_parabolic, initial_condition, solver; source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype) end function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) @@ -142,7 +142,7 @@ function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic) print(io, ", ", semi.equations_parabolic) print(io, ", ", semi.initial_condition) print(io, ", ", semi.boundary_conditions) - print(io, ", ", semi.parabolic_boundary_conditions) + print(io, ", ", semi.boundary_conditions_parabolic) print(io, ", ", semi.source_terms) print(io, ", ", semi.solver) print(io, ", cache(") @@ -236,7 +236,7 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) end function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) - @unpack mesh, equations_parabolic, initial_condition, parabolic_boundary_conditions, source_terms, solver, cache, cache_parabolic = semi + @unpack mesh, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, cache, cache_parabolic = semi u = wrap_array(u_ode, mesh, equations_parabolic, solver, cache_parabolic) du = wrap_array(du_ode, mesh, equations_parabolic, solver, cache_parabolic) @@ -244,7 +244,7 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol # TODO: Taal decide, do we need to pass the mesh? time_start = time_ns() @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, equations_parabolic, initial_condition, - parabolic_boundary_conditions, source_terms, + boundary_conditions_parabolic, source_terms, solver, cache, cache_parabolic) runtime = time_ns() - time_start put!(semi.performance_counter, runtime) From 822425ae37649f900e9ff5af999b197856fe9de3 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:16:44 -0500 Subject: [PATCH 76/86] Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper --- .../semidiscretization_hyperbolic_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index b21064874f4..12bc7fbe2bd 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -72,7 +72,7 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, RealT=real(solver), uEltype=RealT, initial_caches=(NamedTuple(), NamedTuple())) - hyperbolic_equations, equations_parabolic = equations + equations_hyperbolic, equations_parabolic = equations hyperbolic_boundary_conditions, parabolic_boundary_conditions = boundary_conditions initial_hyperbolic_cache, initial_cache_parabolic = initial_caches From 8986d5cd789922709930024282e4e51536dea5c7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:18:51 -0500 Subject: [PATCH 77/86] add equations as field to LaplaceDiffusion2D --- examples/dgmulti_2d/elixir_advection_diffusion.jl | 2 +- .../dgmulti_2d/elixir_advection_diffusion_periodic.jl | 2 +- src/equations/laplace_diffusion_2d.jl | 10 ++-------- 3 files changed, 4 insertions(+), 10 deletions(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index 61d5a684306..92ad1b94a7a 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -5,7 +5,7 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.5, 1.0) -equations_parabolic = LaplaceDiffusion2D(5e-2) +equations_parabolic = LaplaceDiffusion2D(5e-2, equations) initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0) initial_condition = initial_condition_zero diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl index 4f6d9fab580..ceeb01f71f1 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl @@ -6,7 +6,7 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial( volume_integral = VolumeIntegralWeakForm()) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) -equations_parabolic = LaplaceDiffusion2D(1e-2) +equations_parabolic = LaplaceDiffusion2D(1e-2, equations) function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) return SVector(exp(-100 * (x[1]^2 + x[2]^2))) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index d9ed13ee11e..b42137c99d3 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -11,14 +11,8 @@ struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N} equations::E end -function LaplaceDiffusion2D(diffusivity, equations=nothing) - if equations isa AbstractEquations - nvars = nvariables(equations) - else # assume scalar diffusion if no equations are specified - nvars = 1 - end - return LaplaceDiffusion2D{typeof(equations), nvars, typeof(diffusivity)}(diffusivity, equations) -end +LaplaceDiffusion2D(diffusivity, equations) = + LaplaceDiffusion2D{typeof(equations), nvariables(equations), typeof(diffusivity)}(diffusivity, equations) # no orientation specified since the flux is vector-valued function flux(u, grad_u, equations::LaplaceDiffusion2D) From 4c3a7edfafc1becc07bd74e396f33bd255632c91 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 06:48:04 -0500 Subject: [PATCH 78/86] moving Neumann BCs into equations.jl and adding default BC behavior --- src/equations/equations.jl | 32 ++++++++++++++++++++------- src/equations/laplace_diffusion_2d.jl | 15 +------------ 2 files changed, 25 insertions(+), 22 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 7cb4b36bba6..7cc6768d116 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -175,19 +175,39 @@ end return flux end +# operator types used for dispatch on parabolic boundary fluxes +struct Gradient end +struct Divergence end + +# Dirichlet-type boundary condition for use with parabolic equations. +@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector, + x, t, operator_type::Divergence, + equations) + return u_inner +end + + """ - BoundaryConditionNeumann(boundary_value_function) + BoundaryConditionNeumann(boundary_normal_flux_function) Similar to `BoundaryConditionDirichlet`, but creates a Neumann boundary condition for parabolic -equations that uses the function `boundary_value_function` to specify the values at the boundary. +equations that uses the function `boundary_normal_flux_function` to specify the values of the normal +flux at the boundary. The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as ```julia -boundary_value_function(x, t, equations) +boundary_normal_flux_function(x, t, equations) ``` where `x` specifies the coordinates, `t` is the current time, and `equation` is the corresponding system of equations. """ struct BoundaryConditionNeumann{B} - boundary_value_function::B + boundary_normal_flux_function::B +end + +# specify default behavior for a Neumann BC is to evaluate the flux at the inner state +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, + x, t, operator_type::Gradient, + equations) + return flux_inner end # Imposing no boundary condition just evaluates the flux at the inner state. @@ -350,10 +370,6 @@ include("acoustic_perturbation_2d.jl") abstract type AbstractEquationsParabolic{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end -# operator types used for dispatch on parabolic boundary fluxes -struct Gradient end -struct Divergence end - # Linear scalar diffusion for use in linear scalar advection-diffusion problems abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("laplace_diffusion_2d.jl") diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index b42137c99d3..ee8fa952b33 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -27,21 +27,8 @@ end return boundary_condition.boundary_value_function(x, t, equations) end -@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, normal::AbstractVector, - x, t, operator_type::Divergence, - equations::LaplaceDiffusion2D) - return flux_inner -end - - -@inline function (boundary_condition::BoundaryConditionNeumann)(u_inner, normal::AbstractVector, - x, t, operator_type::Gradient, - equations::LaplaceDiffusion2D) - return u_inner -end - @inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector, x, t, operator_type::Divergence, equations::LaplaceDiffusion2D) - return boundary_condition.boundary_value_function(x, t, equations) + return boundary_condition.boundary_normal_flux_function(x, t, equations) end \ No newline at end of file From a4a4b44b9ff09e4b28b69a145e3887f119083bbf Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 08:34:58 -0500 Subject: [PATCH 79/86] running parabolic tests, adding integration test --- test/runtests.jl | 4 ++++ test/{test_parabolic.jl => test_parabolic_2d.jl} | 10 +++++++++- 2 files changed, 13 insertions(+), 1 deletion(-) rename test/{test_parabolic.jl => test_parabolic_2d.jl} (86%) diff --git a/test/runtests.jl b/test/runtests.jl index e294d10746b..79a8e54e791 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -86,6 +86,10 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) include("test_dgmulti_3d.jl") end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "parabolic" + include("test_parabolic_2d.jl") + end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1" include("test_unit.jl") include("test_visualization.jl") diff --git a/test/test_parabolic.jl b/test/test_parabolic_2d.jl similarity index 86% rename from test/test_parabolic.jl rename to test/test_parabolic_2d.jl index c06619b26c5..56619d5e993 100644 --- a/test/test_parabolic.jl +++ b/test/test_parabolic_2d.jl @@ -13,7 +13,7 @@ isdir(outdir) && rm(outdir, recursive=true) @testset "SemidiscretizationHyperbolicParabolic" begin - @trixi_testset "DGMulti 2D" begin + @trixi_testset "DGMulti 2D rhs_parabolic!" begin dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), surface_integral = SurfaceIntegralWeakForm(flux_central), @@ -51,6 +51,14 @@ isdir(outdir) && rm(outdir, recursive=true) @test getindex.(du, 1) ≈ 2 * y end + @trixi_testset "elixir_advection_diffusion.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1) + l2 = [0.2485803335154642], + linf = [1.079606969242132] + ) + end + end # Clean up afterwards: delete Trixi output directory From 41b3480796e267733ca0cb08bc20f955f65dd0cd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 08:38:34 -0500 Subject: [PATCH 80/86] adding parabolic tests to ci --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f257da07a44..6aba503691e 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -69,6 +69,7 @@ jobs: - p4est_part1 - p4est_part2 - unstructured_dgmulti + - parabolic - paper_self_gravitating_gas_dynamics - misc_part1 - misc_part2 From f0beee87847ebc0285f42f47e7fd3c15008b081a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 09:26:06 -0500 Subject: [PATCH 81/86] missing comma --- test/test_parabolic_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 56619d5e993..5d6df53073e 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -53,7 +53,7 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_advection_diffusion.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion.jl"), - cells_per_dimension = (4, 4), tspan=(0.0, 0.1) + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.2485803335154642], linf = [1.079606969242132] ) From 696ebc4c4aeeca1dc10e8c579c1853513b194144 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 09:26:33 -0500 Subject: [PATCH 82/86] Update src/equations/laplace_diffusion_2d.jl Co-authored-by: Hendrik Ranocha --- src/equations/laplace_diffusion_2d.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index ee8fa952b33..774d93cf7b3 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -3,8 +3,6 @@ `LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` with diffusivity ``\kappa`` applied to each solution component defined by `equations`. - -If `equations` is not specified, a scalar `LaplaceDiffusion2D` is returned. """ struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N} diffusivity::T From 5bfa918651c024733591e658c3a47cb5c8f81b39 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 09:44:29 -0500 Subject: [PATCH 83/86] fix test --- examples/dgmulti_2d/elixir_advection_diffusion.jl | 2 ++ test/test_parabolic_2d.jl | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index 92ad1b94a7a..21effda986e 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -43,6 +43,8 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) callbacks = CallbackSet(summary_callback, alive_callback) ############################################################################### diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 5d6df53073e..7dbec5ae10e 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -25,7 +25,7 @@ isdir(outdir) && rm(outdir, recursive=true) initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2]) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) - equations_parabolic = LaplaceDiffusion2D(1.0) + equations_parabolic = LaplaceDiffusion2D(1.0, equations) semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) ode = semidiscretize(semi, (0.0, 0.01)) From 3bf4357cb8d8d66519749a643b4cf56acaf750f1 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 10:03:12 -0500 Subject: [PATCH 84/86] fix test typo --- test/test_parabolic_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 7dbec5ae10e..f4b7e72b6dc 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -27,7 +27,7 @@ isdir(outdir) && rm(outdir, recursive=true) equations = LinearScalarAdvectionEquation2D(1.0, 1.0) equations_parabolic = LaplaceDiffusion2D(1.0, equations) - semi = SemidiscretizationHyperbolicParabolic(mesh, equation, equation_parabolic, initial_condition, dg) + semi = SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, dg) ode = semidiscretize(semi, (0.0, 0.01)) @unpack cache, cache_parabolic, equations_parabolic = semi From 226795437c015fe35e5de28398c8efb7cf583449 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 10:26:40 -0500 Subject: [PATCH 85/86] fix tests (for real this time) --- test/test_parabolic_2d.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index f4b7e72b6dc..1e85fc34820 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -36,18 +36,24 @@ isdir(outdir) && rm(outdir, recursive=true) fill!(u_grad[dim], zero(eltype(u_grad[dim]))) end - Trixi.calc_gradient!(u_grad, ode.u0, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, cache_parabolic) + t = 0.0 + # pass in `boundary_condition_periodic` to fake "do-nothing" + # TODO: DGMulti. Make `boundary_condition_do_nothing` a callable struct like BoundaryConditionPeriodic + Trixi.calc_gradient!(u_grad, ode.u0, t, mesh, equations_parabolic, + boundary_condition_periodic, dg, cache, cache_parabolic) @unpack x, y = mesh.md @test getindex.(u_grad[1], 1) ≈ 2 * x .* y @test getindex.(u_grad[2], 1) ≈ x.^2 u_flux = similar.(u_grad) - Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, equations_parabolic, dg, cache, cache_parabolic) + Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, equations_parabolic, + dg, cache, cache_parabolic) @test u_flux[1] ≈ u_grad[1] @test u_flux[2] ≈ u_grad[2] du = similar(ode.u0) - Trixi.calc_divergence!(du, ode.u0, u_flux, mesh, equations_parabolic, boundary_condition_do_nothing, dg, cache, cache_parabolic) + Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic, + boundary_condition_periodic, dg, cache, cache_parabolic) @test getindex.(du, 1) ≈ 2 * y end From 63d67d77310cfad8572e74ff11edfb0fecb4158c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 27 May 2022 13:30:34 -0500 Subject: [PATCH 86/86] improve test coverage --- test/test_parabolic_2d.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 1e85fc34820..8da8468c8c0 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -28,7 +28,17 @@ isdir(outdir) && rm(outdir, recursive=true) equations_parabolic = LaplaceDiffusion2D(1.0, equations) semi = SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, dg) + @test_nowarn_debug show(stdout, semi) + @test_nowarn_debug show(stdout, MIME"text/plain"(), semi) + + @test nvariables(semi)==nvariables(equations) + @test Base.ndims(semi)==Base.ndims(mesh) + @test Base.real(semi)==Base.real(dg) + ode = semidiscretize(semi, (0.0, 0.01)) + u0 = similar(ode.u0) + Trixi.compute_coefficients!(u0, 0.0, semi) + @test u0 ≈ ode.u0 @unpack cache, cache_parabolic, equations_parabolic = semi @unpack u_grad = cache_parabolic