From e0c973d3966a846de0f311fb1608ac44178be8a0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 11 Jun 2021 15:39:06 +0200 Subject: [PATCH] First step towards FD methods and surface integrals/terms (#617) * introduce surface integrals This could be considered to be breaking since it changes a documented constructor of our `DG` type. On the other hand, the documentation of this constructor was just wrong (since it contained the basis twice but not the mortar). Hence, this is more like a bug fix, I think ;). * do not restrict volume_integral to AbstractVolumeIntegral for DGSEM constructor * use the central flux as default flux for surface integral/terms * remove temporary code * WIP: prepare 2D tree mesh for FD * update note on open issues * update compat bounds for SummationByPartsOperators * fix unit test * revisions based on code review * no need to use the SafeMode anymore * add NEWS * add tests * adapt comment on computation of reference_length [skip ci] * allow other surface terms/integral types for DGSEM * add additional TODO notes [skip ci] * adapt 1D/3D containers consistently * fix deinition of real * fix tests * add surface_integral to p4est methods * ignore warning importing deprecated binding Colors.RGB1 into PlotUtils --- NEWS.md | 2 + Project.toml | 2 + .../elixir_advection_extended.jl | 57 +++++ src/Trixi.jl | 11 +- src/auxiliary/precompile.jl | 24 +- src/callbacks_step/amr_dg1d.jl | 4 +- src/callbacks_step/amr_dg2d.jl | 2 +- src/callbacks_step/amr_dg3d.jl | 4 +- src/callbacks_step/analysis.jl | 4 +- src/callbacks_step/analysis_dg2d.jl | 2 +- src/equations/equations.jl | 3 + src/solvers/dg_common.jl | 7 +- src/solvers/dg_curved/dg.jl | 14 +- src/solvers/dg_curved/dg_1d.jl | 26 +-- src/solvers/dg_curved/dg_2d.jl | 35 +-- src/solvers/dg_curved/dg_3d.jl | 43 ++-- src/solvers/dg_p4est/dg_2d.jl | 28 ++- src/solvers/dg_tree/basis_lobatto_legendre.jl | 28 ++- src/solvers/dg_tree/containers_1d.jl | 123 ++++++----- src/solvers/dg_tree/containers_2d.jl | 209 +++++++++--------- src/solvers/dg_tree/containers_3d.jl | 189 ++++++++-------- src/solvers/dg_tree/dg.jl | 130 +++++++++-- src/solvers/dg_tree/dg_1d.jl | 43 ++-- src/solvers/dg_tree/dg_2d.jl | 81 ++++--- src/solvers/dg_tree/dg_2d_parallel.jl | 24 +- src/solvers/dg_tree/dg_3d.jl | 97 ++++---- src/solvers/dg_unstructured_quad/dg_2d.jl | 57 ++--- src/solvers/fdsbp_tree/fdsbp_2d.jl | 207 +++++++++++++++++ src/solvers/solvers.jl | 1 + test/test_examples_2d_part3.jl | 3 + test/test_tree_2d_fdsbp.jl | 20 ++ test/test_trixi.jl | 6 +- test/test_unit.jl | 8 +- 33 files changed, 989 insertions(+), 505 deletions(-) create mode 100644 examples/tree_2d_fdsbp/elixir_advection_extended.jl create mode 100644 src/solvers/fdsbp_tree/fdsbp_2d.jl create mode 100644 test/test_tree_2d_fdsbp.jl diff --git a/NEWS.md b/NEWS.md index 11afbc8ae5..9da9c8552d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -20,6 +20,8 @@ for human readability. - New structured, curvilinear, conforming mesh type `CurvedMesh` (experimental) - New unstructured, curvilinear, conforming mesh type `UnstructuredQuadMesh` in 2D (experimental) - New unstructured, curvilinear, adaptive (non-conforming) mesh type `P4estMesh` in 2D (experimental) +- Experimental support for finite difference (FD) summation-by-parts (SBP) methods via + [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl) #### Changed diff --git a/Project.toml b/Project.toml index 2e83688c01..2c4ceedbcb 100644 --- a/Project.toml +++ b/Project.toml @@ -22,6 +22,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b" +SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" @@ -42,6 +43,7 @@ Reexport = "1.0" Requires = "1.1" StaticArrays = "1.0" StrideArrays = "0.1.8" +SummationByPartsOperators = "0.5.4" TimerOutputs = "0.5" UnPack = "1.0" julia = "1.6" diff --git a/examples/tree_2d_fdsbp/elixir_advection_extended.jl b/examples/tree_2d_fdsbp/elixir_advection_extended.jl new file mode 100644 index 0000000000..7a16b35875 --- /dev/null +++ b/examples/tree_2d_fdsbp/elixir_advection_extended.jl @@ -0,0 +1,57 @@ +# TODO: FD +# !!! warning "Experimental feature" +# This is an experimental feature and may change in any future releases. + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advectionvelocity = (0.2, -0.3) +equations = LinearScalarAdvectionEquation2D(advectionvelocity) + +initial_condition = initial_condition_convergence_test + +D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(), + derivative_order=1, accuracy_order=4, + xmin=0.0, xmax=1.0, N=100) +solver = DG(D_SBP, nothing #= mortar =#, + SurfaceIntegralStrongForm(flux_lax_friedrichs), + VolumeIntegralStrongForm()) + +coordinates_min = (-1.0, -1.0) +coordinates_max = ( 1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=1, + n_cells_max=30_000, + periodicity=true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + extra_analysis_integrals=(energy_total,)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback) + + +############################################################################### +# run the simulation + +sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-8, reltol=1.0e-8, + save_everystep=false, callback=callbacks, maxiters=1e5) +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index 74269a6590..0874509a31 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -18,7 +18,7 @@ module Trixi # Include other packages that are used in Trixi # (standard library packages first, other packages next, all of them sorted alphabetically) -using LinearAlgebra: dot, mul!, norm, cross, normalize +using LinearAlgebra: diag, dot, mul!, norm, cross, normalize using Printf: @printf, @sprintf, println # import @reexport now to make it available for further imports/exports @@ -47,6 +47,11 @@ using TimerOutputs: TimerOutputs, @notimeit, TimerOutput, print_timer, reset_tim @reexport using UnPack: @unpack using UnPack: @pack! +# finite difference SBP operators +using SummationByPartsOperators: AbstractDerivativeOperator, DerivativeOperator, grid +import SummationByPartsOperators: integrate, left_boundary_weight, right_boundary_weight +@reexport using SummationByPartsOperators: + SummationByPartsOperators, derivative_operator # Define the entry points of our type hierarchy, e.g. # AbstractEquations, AbstractSemidiscretization etc. @@ -139,9 +144,11 @@ export TreeMesh, CurvedMesh, UnstructuredQuadMesh, P4estMesh export DG, DGSEM, LobattoLegendreBasis, - VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing, + VolumeIntegralWeakForm, VolumeIntegralStrongForm, + VolumeIntegralFluxDifferencing, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, + SurfaceIntegralWeakForm, SurfaceIntegralStrongForm, MortarL2 export nelements, nnodes, nvariables, diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 3f450ae358..0210baeee0 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -351,24 +351,24 @@ function _precompile_manual_() mortar_type = mortar_type_dgsem(RealT, nnodes_) # 1D, serial - @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype,1,polydeg}}) - @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype,1,polydeg}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{1,Trixi.SerialTree{1}},Trixi.ElementContainer1D{RealT,uEltype}}) # 2D, serial - @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}}) - @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}}) - @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg},mortar_type}) + @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.SerialTree{2}},Trixi.ElementContainer2D{RealT,uEltype},mortar_type}) # 2D, parallel - @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}}) - @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}}) - @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg},mortar_type}) - @assert Base.precompile(Tuple{typeof(Trixi.init_mpi_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype,1,polydeg}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype},mortar_type}) + @assert Base.precompile(Tuple{typeof(Trixi.init_mpi_interfaces),Array{Int,1},TreeMesh{2,Trixi.ParallelTree{2}},Trixi.ElementContainer2D{RealT,uEltype}}) # 3D, serial - @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype,1,polydeg}}) - @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype,1,polydeg}}) - @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype,1,polydeg},mortar_type}) + @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype}}) + @assert Base.precompile(Tuple{typeof(Trixi.init_mortars),Array{Int,1},TreeMesh{3,Trixi.SerialTree{3}},Trixi.ElementContainer3D{RealT,uEltype},mortar_type}) end # various `show` methods diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 15a9a1aa25..387fa63314 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -23,7 +23,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # re-initialize elements container @unpack elements = cache resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache) > old_n_elements resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) @@ -137,7 +137,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # re-initialize elements container @unpack elements = cache resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache) < old_n_elements resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 84f8aa45ca..46d21931fd 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -69,7 +69,7 @@ function reinitialize_containers!(mesh::TreeMesh{2}, equations, dg::DGSEM, cache # re-initialize elements container @unpack elements = cache resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) # re-initialize interfaces container @unpack interfaces = cache diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 056e53dccd..292893a18a 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -23,7 +23,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, # re-initialize elements container @unpack elements = cache resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache) > old_n_elements resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) @@ -177,7 +177,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, # re-initialize elements container @unpack elements = cache resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis.nodes) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache) < old_n_elements resize!(u_ode, nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index c306d965c6..2b0f227558 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -20,7 +20,7 @@ solution and integrated over the computational domain. See `Trixi.analyze`, `Trixi.pretty_form_utf`, `Trixi.pretty_form_ascii` for further information on how to create custom analysis quantities. """ -mutable struct AnalysisCallback{Analyzer<:SolutionAnalyzer, AnalysisIntegrals, InitialStateIntegrals, Cache} +mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegrals, Cache} start_time::Float64 interval::Int save_analysis::Bool @@ -178,7 +178,7 @@ function (analysis_callback::AnalysisCallback)(integrator) mpi_println() mpi_println("─"^100) # TODO: Taal refactor, polydeg is specific to DGSEM - mpi_println(" Simulation running '", get_name(equations), "' with polydeg = ", polydeg(solver)) + mpi_println(" Simulation running '", get_name(equations), "' with ", summary(solver)) mpi_println("─"^100) mpi_println(" #timesteps: " * @sprintf("% 14d", iter) * " " * diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index b985664b45..8ddd47942e 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -173,7 +173,7 @@ end function integrate(func::Func, u, mesh::Union{TreeMesh{2}, CurvedMesh{2}, UnstructuredQuadMesh, P4estMesh{2}}, - equations, dg::DGSEM, cache; normalize=true) where {Func} + equations, dg::DG, cache; normalize=true) where {Func} integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) return func(u_local, equations) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index f25e62c81d..bc165c928e 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -53,6 +53,9 @@ end @inline Base.ndims(::AbstractEquations{NDIMS}) where NDIMS = NDIMS +# equations act like scalars in broadcasting +Base.broadcastable(equations::AbstractEquations) = Ref(equations) + """ flux(u, orientation_or_normal, equations) diff --git a/src/solvers/dg_common.jl b/src/solvers/dg_common.jl index 7abf4e5e0e..8239802827 100644 --- a/src/solvers/dg_common.jl +++ b/src/solvers/dg_common.jl @@ -5,7 +5,7 @@ function allocate_coefficients(mesh::AbstractMesh, equations, dg::DG, cache) zeros(eltype(cache.elements), nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) end -@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DG, cache) +@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DGSEM, cache) @boundscheck begin @assert length(u_ode) == nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache) end @@ -46,6 +46,11 @@ end end end +# General fallback +@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DG, cache) + wrap_array_native(u_ode, mesh, equations, dg, cache) +end + # Like `wrap_array`, but guarantees to return a plain `Array`, which can be better # for interfacing with external C libraries (MPI, HDF5, visualization), # writing solution files etc. diff --git a/src/solvers/dg_curved/dg.jl b/src/solvers/dg_curved/dg.jl index a2a9bd1bef..78bead1b22 100644 --- a/src/solvers/dg_curved/dg.jl +++ b/src/solvers/dg_curved/dg.jl @@ -17,7 +17,8 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionPeriodic, - mesh::CurvedMesh, equations, dg::DG, cache, + mesh::CurvedMesh, equations, + surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @assert isperiodic(mesh, orientation) end @@ -25,22 +26,23 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition, - mesh::CurvedMesh, equations,dg::DG, cache, + mesh::CurvedMesh, equations, + surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements - @unpack surface_flux = dg + @unpack surface_flux = surface_integral u_inner = get_node_vars(u, equations, dg, node_indices..., element) x = get_node_coords(node_coordinates, equations, dg, node_indices..., element) - # If the mapping is orientation-reversing, the contravariant vectors' orientation - # is reversed as well. The normal vector must be oriented in the direction + # If the mapping is orientation-reversing, the contravariant vectors' orientation + # is reversed as well. The normal vector must be oriented in the direction # from `left_element` to `right_element`, or the numerical flux will be computed # incorrectly (downwind direction). sign_jacobian = sign(inverse_jacobian[node_indices..., element]) # Contravariant vector Ja^i is the normal vector - normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors, + normal = sign_jacobian * get_contravariant_vector(orientation, contravariant_vectors, node_indices..., element) # If the mapping is orientation-reversing, the normal vector will be reversed (see above). diff --git a/src/solvers/dg_curved/dg_1d.jl b/src/solvers/dg_curved/dg_1d.jl index cc88103512..25e491b988 100644 --- a/src/solvers/dg_curved/dg_1d.jl +++ b/src/solvers/dg_curved/dg_1d.jl @@ -13,15 +13,15 @@ function rhs!(du, u, t, # Calculate interface and boundary fluxes @timed timer() "interface flux" calc_interface_flux!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, u, t, boundary_conditions, mesh, equations, dg) + cache, u, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element @timed timer() "Jacobian" apply_jacobian!( @@ -36,8 +36,8 @@ end function calc_interface_flux!(cache, u, mesh::CurvedMesh{1}, - equations, dg::DG) - @unpack surface_flux = dg + equations, surface_integral, dg::DG) + @unpack surface_flux = surface_integral @threaded for element in eachelement(dg, cache) left_element = cache.elements.left_neighbors[1, element] @@ -46,11 +46,11 @@ function calc_interface_flux!(cache, u, mesh::CurvedMesh{1}, u_ll = get_node_vars(u, equations, dg, nnodes(dg), left_element) u_rr = get_node_vars(u, equations, dg, 1, element) - flux = surface_flux(u_ll, u_rr, 1, equations) + f1 = surface_flux(u_ll, u_rr, 1, equations) for v in eachvariable(equations) - cache.elements.surface_flux_values[v, 2, left_element] = flux[v] - cache.elements.surface_flux_values[v, 1, element] = flux[v] + cache.elements.surface_flux_values[v, 2, left_element] = f1[v] + cache.elements.surface_flux_values[v, 1, element] = f1[v] end end end @@ -61,21 +61,21 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::CurvedMesh{1}, equations, dg::DG) + mesh::CurvedMesh{1}, equations, surface_integral, dg::DG) @assert isperiodic(mesh) end function calc_boundary_flux!(cache, u, t, boundary_condition, - mesh::CurvedMesh{1}, equations, dg::DG) + mesh::CurvedMesh{1}, equations, surface_integral, dg::DG) calc_boundary_flux!(cache, u, t, (boundary_condition, boundary_condition), - mesh, equations, dg) + mesh, equations, surface_integral, dg) end function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple}, - mesh::CurvedMesh{1}, equations, dg::DG) - @unpack surface_flux = dg + mesh::CurvedMesh{1}, equations, surface_integral, dg::DG) + @unpack surface_flux = surface_integral @unpack surface_flux_values, node_coordinates = cache.elements orientation = 1 diff --git a/src/solvers/dg_curved/dg_2d.jl b/src/solvers/dg_curved/dg_2d.jl index bd195dfe5d..7fcf7cfb7a 100644 --- a/src/solvers/dg_curved/dg_2d.jl +++ b/src/solvers/dg_curved/dg_2d.jl @@ -13,15 +13,15 @@ function rhs!(du, u, t, # Calculate interface fluxes @timed timer() "interface flux" calc_interface_flux!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, u, t, boundary_conditions, mesh, equations, dg) + cache, u, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element @timed timer() "Jacobian" apply_jacobian!( @@ -78,7 +78,7 @@ end function calc_interface_flux!(cache, u, mesh::CurvedMesh{2}, - equations, dg::DG) + equations, surface_integral, dg::DG) @unpack elements = cache @threaded for element in eachelement(dg, cache) @@ -88,12 +88,14 @@ function calc_interface_flux!(cache, u, # Interfaces in x-direction (`orientation` = 1) calc_interface_flux!(elements.surface_flux_values, elements.left_neighbors[1, element], - element, 1, u, mesh, equations, dg, cache) + element, 1, u, mesh, equations, + surface_integral, dg, cache) # Interfaces in y-direction (`orientation` = 2) calc_interface_flux!(elements.surface_flux_values, elements.left_neighbors[2, element], - element, 2, u, mesh, equations, dg, cache) + element, 2, u, mesh, equations, + surface_integral, dg, cache) end return nothing @@ -103,13 +105,13 @@ end @inline function calc_interface_flux!(surface_flux_values, left_element, right_element, orientation, u, mesh::CurvedMesh{2}, equations, - dg::DG, cache) + surface_integral, dg::DG, cache) # This is slow for LSA, but for some reason faster for Euler (see #519) if left_element <= 0 # left_element = 0 at boundaries return nothing end - @unpack surface_flux = dg + @unpack surface_flux = surface_integral @unpack contravariant_vectors, inverse_jacobian = cache.elements right_direction = 2 * orientation @@ -157,23 +159,22 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::CurvedMesh{2}, equations, dg::DG) + mesh::CurvedMesh{2}, equations, surface_integral, dg::DG) @assert isperiodic(mesh) end function calc_boundary_flux!(cache, u, t, boundary_condition, - mesh::CurvedMesh{2}, equations, dg::DG) + mesh::CurvedMesh{2}, equations, surface_integral, dg::DG) calc_boundary_flux!(cache, u, t, (boundary_condition, boundary_condition, boundary_condition, boundary_condition), - mesh, equations, dg) + mesh, equations, surface_integral, dg) end function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple}, - mesh::CurvedMesh{2}, equations, dg::DG) - @unpack surface_flux = dg + mesh::CurvedMesh{2}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements linear_indices = LinearIndices(size(mesh)) @@ -185,7 +186,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for j in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (1, j), (j,), element) end @@ -196,7 +197,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for j in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (nnodes(dg), j), (j,), element) end end @@ -209,7 +210,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for i in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (i, 1), (i,), element) end @@ -220,7 +221,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for i in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (i, nnodes(dg)), (i,), element) end end diff --git a/src/solvers/dg_curved/dg_3d.jl b/src/solvers/dg_curved/dg_3d.jl index 343df4f4ca..21f71af593 100644 --- a/src/solvers/dg_curved/dg_3d.jl +++ b/src/solvers/dg_curved/dg_3d.jl @@ -13,15 +13,15 @@ function rhs!(du, u, t, # Calculate interface fluxes @timed timer() "interface flux" calc_interface_flux!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, u, t, boundary_conditions, mesh, equations, dg) + cache, u, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element @timed timer() "Jacobian" apply_jacobian!( @@ -209,9 +209,8 @@ end function calc_interface_flux!(cache, u, mesh::CurvedMesh{3}, - equations, dg::DG) + equations, surface_integral, dg::DG) @unpack elements = cache - @unpack surface_flux = dg @threaded for element in eachelement(dg, cache) # Interfaces in negative directions @@ -220,17 +219,20 @@ function calc_interface_flux!(cache, u, mesh::CurvedMesh{3}, # Interfaces in x-direction (`orientation` = 1) calc_interface_flux!(elements.surface_flux_values, elements.left_neighbors[1, element], - element, 1, u, mesh, equations, dg, cache) + element, 1, u, mesh, equations, + surface_integral, dg, cache) # Interfaces in y-direction (`orientation` = 2) calc_interface_flux!(elements.surface_flux_values, elements.left_neighbors[2, element], - element, 2, u, mesh, equations, dg, cache) + element, 2, u, mesh, equations, + surface_integral, dg, cache) # Interfaces in z-direction (`orientation` = 3) calc_interface_flux!(elements.surface_flux_values, elements.left_neighbors[3, element], - element, 3, u, mesh, equations, dg, cache) + element, 3, u, mesh, equations, + surface_integral, dg, cache) end return nothing @@ -240,13 +242,13 @@ end @inline function calc_interface_flux!(surface_flux_values, left_element, right_element, orientation, u, mesh::CurvedMesh{3}, equations, - dg::DG, cache) + surface_integral, dg::DG, cache) # This is slow for LSA, but for some reason faster for Euler (see #519) if left_element <= 0 # left_element = 0 at boundaries return surface_flux_values end - @unpack surface_flux = dg + @unpack surface_flux = surface_integral @unpack contravariant_vectors, inverse_jacobian = cache.elements right_direction = 2 * orientation @@ -304,23 +306,22 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::CurvedMesh{3}, equations, dg::DG) + mesh::CurvedMesh{3}, equations, surface_integral, dg::DG) @assert isperiodic(mesh) end function calc_boundary_flux!(cache, u, t, boundary_condition, - mesh::CurvedMesh{3}, equations, dg::DG) + mesh::CurvedMesh{3}, equations, surface_integral, dg::DG) calc_boundary_flux!(cache, u, t, (boundary_condition, boundary_condition, boundary_condition, boundary_condition, boundary_condition, boundary_condition), - mesh, equations, dg) + mesh, equations, surface_integral, dg) end function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple,Tuple}, - mesh::CurvedMesh{3}, equations, dg::DG) - @unpack surface_flux = dg + mesh::CurvedMesh{3}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements linear_indices = LinearIndices(size(mesh)) @@ -332,7 +333,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for k in eachnode(dg), j in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (1, j, k), (j, k), element) end @@ -343,7 +344,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for k in eachnode(dg), j in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 1, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (nnodes(dg), j, k), (j, k), element) end end @@ -356,7 +357,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for k in eachnode(dg), i in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (i, 1, k), (i, k), element) end @@ -367,7 +368,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for k in eachnode(dg), i in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 2, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (i, nnodes(dg), k), (i, k), element) end end @@ -380,7 +381,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for j in eachnode(dg), i in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 3, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (i, j, 1), (i, j), element) end @@ -391,7 +392,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::Union{NamedTuple, for j in eachnode(dg), i in eachnode(dg) calc_boundary_flux_by_direction!(surface_flux_values, u, t, 3, boundary_conditions[direction], - mesh, equations, dg, cache, + mesh, equations, surface_integral, dg, cache, direction, (i, j, nnodes(dg)), (i, j), element) end end diff --git a/src/solvers/dg_p4est/dg_2d.jl b/src/solvers/dg_p4est/dg_2d.jl index 625763176e..d4164f549e 100644 --- a/src/solvers/dg_p4est/dg_2d.jl +++ b/src/solvers/dg_p4est/dg_2d.jl @@ -15,7 +15,7 @@ end function prolong2interfaces!(cache, u, mesh::P4estMesh{2}, - equations, dg::DG) + equations, surface_integral, dg::DG) @unpack interfaces = cache size_ = (nnodes(dg), nnodes(dg)) @@ -47,8 +47,9 @@ end function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, nonconservative_terms::Val{false}, - equations, dg::DG, cache) - @unpack surface_flux = dg + equations, surface_integral, + dg::DG, cache) + @unpack surface_flux = surface_integral @unpack u, element_ids, node_indices = cache.interfaces size_ = (nnodes(dg), nnodes(dg)) @@ -93,7 +94,8 @@ end function prolong2boundaries!(cache, u, - mesh::P4estMesh{2}, equations, dg::DG) + mesh::P4estMesh{2}, equations, + surface_integral, dg::DG) @unpack boundaries = cache size_ = (nnodes(dg), nnodes(dg)) @@ -116,10 +118,11 @@ end function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, - mesh::P4estMesh, equations, dg::DG) + mesh::P4estMesh, equations, + surface_integral, dg::DG) @unpack boundaries = cache @unpack surface_flux_values, node_coordinates = cache.elements - @unpack surface_flux = dg + @unpack surface_flux = surface_integral size_ = (nnodes(dg), nnodes(dg)) @@ -161,7 +164,8 @@ end function prolong2mortars!(cache, u, mesh::P4estMesh{2}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DGSEM) @unpack element_ids, node_indices = cache.mortars size_ = (nnodes(dg), nnodes(dg)) @@ -208,10 +212,11 @@ end function calc_mortar_flux!(surface_flux_values, mesh::P4estMesh{2}, nonconservative_terms::Val{false}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DG, cache) + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) @unpack u, element_ids, node_indices = cache.mortars @unpack fstar_upper_threaded, fstar_lower_threaded = cache - @unpack surface_flux = dg + @unpack surface_flux = surface_integral size_ = (nnodes(dg), nnodes(dg)) @@ -305,8 +310,9 @@ end end -function calc_surface_integral!(du, mesh::P4estMesh, - equations, dg::DGSEM, cache) +function calc_surface_integral!(du, u, mesh::P4estMesh, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM, cache) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements diff --git a/src/solvers/dg_tree/basis_lobatto_legendre.jl b/src/solvers/dg_tree/basis_lobatto_legendre.jl index d9526625e0..6bd24c3ed4 100644 --- a/src/solvers/dg_tree/basis_lobatto_legendre.jl +++ b/src/solvers/dg_tree/basis_lobatto_legendre.jl @@ -84,6 +84,32 @@ end @inline polydeg(basis::LobattoLegendreBasis) = nnodes(basis) - 1 +@inline get_nodes(basis::LobattoLegendreBasis) = basis.nodes + + +""" + integrate(f, u, basis::LobattoLegendreBasis) + +Map the function `f` to the coefficients `u` and integrate with respect to the +quadrature rule given by `basis`. +""" +function integrate(f, u, basis::LobattoLegendreBasis) + @unpack weights = basis + + res = zero(f(first(u))) + for i in eachindex(u, weights) + res += f(u[i]) * weights[i] + end + return res +end + +# Return the first/last weight of the quadrature associated with `basis`. +# Since the mass matrix for nodal Lobatto-Legendre bases is diagonal, +# these weights are the only coefficients necessary for the scaling of +# surface terms/integrals in DGSEM. +left_boundary_weight(basis::LobattoLegendreBasis) = first(basis.weights) +right_boundary_weight(basis::LobattoLegendreBasis) = last(basis.weights) + struct LobattoLegendreMortarL2{RealT<:Real, NNODES, ForwardMatrix<:AbstractMatrix{RealT}, ReverseMatrix<:AbstractMatrix{RealT}} <: AbstractMortarL2{RealT} @@ -193,7 +219,7 @@ function SolutionAnalyzer(basis::LobattoLegendreBasis; analysis_polydeg=2*polyde # compute everything using `Float64` by default nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_) - vandermonde_ = polynomial_interpolation_matrix(basis.nodes, nodes_) + vandermonde_ = polynomial_interpolation_matrix(get_nodes(basis), nodes_) # type conversions to get the requested real type and enable possible # optimizations of runtime performance and latency diff --git a/src/solvers/dg_tree/containers_1d.jl b/src/solvers/dg_tree/containers_1d.jl index dd7de77cfd..113717f718 100644 --- a/src/solvers/dg_tree/containers_1d.jl +++ b/src/solvers/dg_tree/containers_1d.jl @@ -1,7 +1,6 @@ # Container data structure (structure-of-arrays style) for DG elements -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct ElementContainer1D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct ElementContainer1D{RealT<:Real, uEltype<:Real} <: AbstractContainer inverse_jacobian::Vector{RealT} # [elements] node_coordinates::Array{RealT, 3} # [orientation, i, elements] surface_flux_values::Array{uEltype, 3} # [variables, direction, elements] @@ -11,9 +10,9 @@ mutable struct ElementContainer1D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: _surface_flux_values::Vector{uEltype} end -nvariables(::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS -polydeg(::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::ElementContainer1D{RealT, uEltype}) where {RealT, uEltype} = uEltype +nvariables(elements::ElementContainer1D) = size(elements.surface_flux_values, 1) +nnodes(elements::ElementContainer1D) = size(elements.node_coordinates, 2) +Base.eltype(elements::ElementContainer1D) = eltype(elements.surface_flux_values) # Only one-dimensional `Array`s are `resize!`able in Julia. # Hence, we use `Vector`s as internal storage and `resize!` @@ -21,8 +20,8 @@ Base.eltype(::ElementContainer1D{RealT, uEltype}) where {RealT, uEltype} = uElty # `unsafe_wrap`ping multi-dimensional `Array`s around the # internal storage. function Base.resize!(elements::ElementContainer1D, capacity) - n_nodes = polydeg(elements) + 1 - nvars = nvariables(elements) + n_nodes = nnodes(elements) + n_variables = nvariables(elements) @unpack _node_coordinates, _surface_flux_values, inverse_jacobian, cell_ids = elements @@ -32,9 +31,9 @@ function Base.resize!(elements::ElementContainer1D, capacity) elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (1, n_nodes, capacity)) - resize!(_surface_flux_values, nvars * 2 * 1 * capacity) + resize!(_surface_flux_values, n_variables * 2 * 1 * capacity) elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values), - (nvars, 2 * 1, capacity)) + (n_variables, 2 * 1, capacity)) resize!(cell_ids, capacity) @@ -42,8 +41,7 @@ function Base.resize!(elements::ElementContainer1D, capacity) end -function ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function ElementContainer1D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real} nan_RealT = convert(RealT, NaN) nan_uEltype = convert(uEltype, NaN) @@ -54,13 +52,13 @@ function ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) w node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (1, n_nodes, capacity)) - _surface_flux_values = fill(nan_uEltype, NVARS * 2 * 1 * capacity) + _surface_flux_values = fill(nan_uEltype, n_variables * 2 * 1 * capacity) surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values), - (NVARS, 2 * 1, capacity)) + (n_variables, 2 * 1, capacity)) cell_ids = fill(typemin(Int), capacity) - return ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}( + return ElementContainer1D{RealT, uEltype}( inverse_jacobian, node_coordinates, surface_flux_values, cell_ids, _node_coordinates, _surface_flux_values) end @@ -70,21 +68,32 @@ end @inline nelements(elements::ElementContainer1D) = length(elements.cell_ids) # TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements)) @inline eachelement(elements::ElementContainer1D) = Base.OneTo(nelements(elements)) +@inline Base.real(elements::ElementContainer1D) = eltype(elements.node_coordinates) # Create element container and initialize element data function init_elements(cell_ids, mesh::TreeMesh1D, - equations::AbstractEquations{1, NVARS}, - basis::LobattoLegendreBasis{T, NNODES}, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real,NVARS, T, NNODES} + equations::AbstractEquations{1}, + basis, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real} # Initialize container n_elements = length(cell_ids) - elements = ElementContainer1D{RealT, uEltype, NVARS, NNODES-1}(n_elements) + elements = ElementContainer1D{RealT, uEltype}( + n_elements, nvariables(equations), nnodes(basis)) - init_elements!(elements, cell_ids, mesh, basis.nodes) + init_elements!(elements, cell_ids, mesh, basis) return elements end -function init_elements!(elements, cell_ids, mesh::TreeMesh1D, nodes) +function init_elements!(elements, cell_ids, mesh::TreeMesh1D, basis) + nodes = get_nodes(basis) + # Compute the length of the 1D reference interval by integrating + # the function with constant value unity on the corresponding + # element data type (using \circ) + reference_length = integrate(one ∘ eltype, nodes, basis) + # Compute the offset of the midpoint of the 1D reference interval + # (its difference from zero) + reference_offset = first(nodes) + reference_length / 2 + # Store cell ids elements.cell_ids .= cell_ids @@ -96,13 +105,17 @@ function init_elements!(elements, cell_ids, mesh::TreeMesh1D, nodes) # Get cell length dx = length_at_cell(mesh.tree, cell_id) - # Calculate inverse Jacobian as 1/(h/2) - elements.inverse_jacobian[element] = 2/dx + # Calculate inverse Jacobian + jacobian = dx / reference_length + elements.inverse_jacobian[element] = inv(jacobian) # Calculate node coordinates - for i in eachindex(nodes) + # Note that the `tree_coordinates` are the midpoints of the cells. + # Hence, we need to add an offset for `nodes` with a midpoint + # different from zero. + for i in eachnode(basis) elements.node_coordinates[1, i, element] = ( - mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i]) + mesh.tree.coordinates[1, cell_id] + jacobian * (nodes[i] - reference_offset)) end end @@ -112,8 +125,7 @@ end # Container data structure (structure-of-arrays style) for DG interfaces -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct InterfaceContainer1D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct InterfaceContainer1D{uEltype<:Real} <: AbstractContainer u::Array{uEltype, 3} # [leftright, variables, interfaces] neighbor_ids::Matrix{Int} # [leftright, interfaces] orientations::Vector{Int} # [interfaces] @@ -122,22 +134,20 @@ mutable struct InterfaceContainer1D{uEltype<:Real, NVARS, POLYDEG} <: AbstractCo _neighbor_ids::Vector{Int} end -nvariables(::InterfaceContainer1D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS -polydeg(::InterfaceContainer1D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::InterfaceContainer1D{uEltype}) where {uEltype} = uEltype +nvariables(interfaces::InterfaceContainer1D) = size(interfaces.u, 2) +Base.eltype(interfaces::InterfaceContainer1D) = eltype(interfaces.u) # See explanation of Base.resize! for the element container function Base.resize!(interfaces::InterfaceContainer1D, capacity) - n_nodes = polydeg(interfaces) + 1 - nvars = nvariables(interfaces) + n_variables = nvariables(interfaces) @unpack _u, _neighbor_ids, orientations = interfaces - resize!(_u, 2 * nvars * capacity) + resize!(_u, 2 * n_variables * capacity) interfaces.u = unsafe_wrap(Array, pointer(_u), - (2, nvars, capacity)) + (2, n_variables, capacity)) resize!(_neighbor_ids, 2 * capacity) - interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ), + interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, capacity)) resize!(orientations, capacity) @@ -146,22 +156,21 @@ function Base.resize!(interfaces::InterfaceContainer1D, capacity) end -function InterfaceContainer1D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function InterfaceContainer1D{uEltype}(capacity::Integer, n_variables, n_nodes) where {uEltype<:Real} nan = convert(uEltype, NaN) # Initialize fields with defaults - _u = fill(nan, 2 * NVARS * capacity) + _u = fill(nan, 2 * n_variables * capacity) u = unsafe_wrap(Array, pointer(_u), - (2, NVARS, capacity)) + (2, n_variables, capacity)) _neighbor_ids = fill(typemin(Int), 2 * capacity) - neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ), + neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, capacity)) orientations = fill(typemin(Int), capacity) - return InterfaceContainer1D{uEltype, NVARS, POLYDEG}( + return InterfaceContainer1D{uEltype}( u, neighbor_ids, orientations, _u, _neighbor_ids) end @@ -173,10 +182,11 @@ end # Create interface container and initialize interface data in `elements`. function init_interfaces(cell_ids, mesh::TreeMesh1D, - elements::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer1D) # Initialize container n_interfaces = count_required_interfaces(mesh, cell_ids) - interfaces = InterfaceContainer1D{uEltype, NVARS, POLYDEG}(n_interfaces) + interfaces = InterfaceContainer1D{eltype(elements)}( + n_interfaces, nvariables(elements), nnodes(elements)) # Connect elements with interfaces init_interfaces!(interfaces, elements, mesh) @@ -264,8 +274,7 @@ end # Container data structure (structure-of-arrays style) for DG boundaries -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct BoundaryContainer1D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct BoundaryContainer1D{RealT<:Real, uEltype<:Real} <: AbstractContainer u::Array{uEltype, 3} # [leftright, variables, boundaries] neighbor_ids::Vector{Int} # [boundaries] orientations::Vector{Int} # [boundaries] @@ -277,23 +286,21 @@ mutable struct BoundaryContainer1D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} < _node_coordinates::Vector{RealT} end -nvariables(::BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS -polydeg(::BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::BoundaryContainer1D{RealT, uEltype}) where {RealT, uEltype} = uEltype +nvariables(boundaries::BoundaryContainer1D) = size(boundaries.u, 2) +Base.eltype(boundaries::BoundaryContainer1D) = eltype(boundaries.u) # See explanation of Base.resize! for the element container function Base.resize!(boundaries::BoundaryContainer1D, capacity) - n_nodes = polydeg(boundaries) + 1 - nvars = nvariables(boundaries) + n_variables = nvariables(boundaries) @unpack _u, _node_coordinates, neighbor_ids, orientations, neighbor_sides = boundaries - resize!(_u, 2 * nvars * capacity) + resize!(_u, 2 * n_variables * capacity) boundaries.u = unsafe_wrap(Array, pointer(_u), - (2, nvars, capacity)) + (2, n_variables, capacity)) resize!(_node_coordinates, 1 * capacity) - boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates ), + boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (1, capacity)) resize!(neighbor_ids, capacity) @@ -306,15 +313,14 @@ function Base.resize!(boundaries::BoundaryContainer1D, capacity) end -function BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function BoundaryContainer1D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real} nan_RealT = convert(RealT, NaN) nan_uEltype = convert(uEltype, NaN) # Initialize fields with defaults - _u = fill(nan_uEltype, 2 * NVARS * capacity) + _u = fill(nan_uEltype, 2 * n_variables * capacity) u = unsafe_wrap(Array, pointer(_u), - (2, NVARS, capacity)) + (2, n_variables, capacity)) neighbor_ids = fill(typemin(Int), capacity) @@ -328,7 +334,7 @@ function BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) n_boundaries_per_direction = SVector(0, 0) - return BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}( + return BoundaryContainer1D{RealT, uEltype}( u, neighbor_ids, orientations, neighbor_sides, node_coordinates, n_boundaries_per_direction, _u, _node_coordinates) @@ -341,10 +347,11 @@ nboundaries(boundaries::BoundaryContainer1D) = length(boundaries.orientations) # Create boundaries container and initialize boundary data in `elements`. function init_boundaries(cell_ids, mesh::TreeMesh1D, - elements::ElementContainer1D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} + elements::ElementContainer1D) # Initialize container n_boundaries = count_required_boundaries(mesh, cell_ids) - boundaries = BoundaryContainer1D{RealT, uEltype, NVARS, POLYDEG}(n_boundaries) + boundaries = BoundaryContainer1D{real(elements), eltype(elements)}( + n_boundaries, nvariables(elements), nnodes(elements)) # Connect elements with boundaries init_boundaries!(boundaries, elements, mesh) diff --git a/src/solvers/dg_tree/containers_2d.jl b/src/solvers/dg_tree/containers_2d.jl index 88f0f48470..adf010bcc3 100644 --- a/src/solvers/dg_tree/containers_2d.jl +++ b/src/solvers/dg_tree/containers_2d.jl @@ -1,7 +1,6 @@ # Container data structure (structure-of-arrays style) for DG elements -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct ElementContainer2D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct ElementContainer2D{RealT<:Real, uEltype<:Real} <: AbstractContainer inverse_jacobian::Vector{RealT} # [elements] node_coordinates::Array{RealT, 4} # [orientation, i, j, elements] surface_flux_values::Array{uEltype, 4} # [variables, i, direction, elements] @@ -11,9 +10,9 @@ mutable struct ElementContainer2D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: _surface_flux_values::Vector{uEltype} end -nvariables(::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS -polydeg(::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::ElementContainer2D{RealT, uEltype}) where {RealT, uEltype} = uEltype +nvariables(elements::ElementContainer2D) = size(elements.surface_flux_values, 1) +nnodes(elements::ElementContainer2D) = size(elements.node_coordinates, 2) +Base.eltype(elements::ElementContainer2D) = eltype(elements.surface_flux_values) # Only one-dimensional `Array`s are `resize!`able in Julia. # Hence, we use `Vector`s as internal storage and `resize!` @@ -21,8 +20,8 @@ Base.eltype(::ElementContainer2D{RealT, uEltype}) where {RealT, uEltype} = uElty # `unsafe_wrap`ping multi-dimensional `Array`s around the # internal storage. function Base.resize!(elements::ElementContainer2D, capacity) - n_nodes = polydeg(elements) + 1 - nvars = nvariables(elements) + n_nodes = nnodes(elements) + n_variables = nvariables(elements) @unpack _node_coordinates, _surface_flux_values, inverse_jacobian, cell_ids = elements @@ -32,9 +31,9 @@ function Base.resize!(elements::ElementContainer2D, capacity) elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (2, n_nodes, n_nodes, capacity)) - resize!(_surface_flux_values, nvars * n_nodes * 2 * 2 * capacity) + resize!(_surface_flux_values, n_variables * n_nodes * 2 * 2 * capacity) elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values), - (nvars, n_nodes, 2 * 2, capacity)) + (n_variables, n_nodes, 2 * 2, capacity)) resize!(cell_ids, capacity) @@ -42,8 +41,7 @@ function Base.resize!(elements::ElementContainer2D, capacity) end -function ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function ElementContainer2D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real} nan_RealT = convert(RealT, NaN) nan_uEltype = convert(uEltype, NaN) @@ -54,14 +52,14 @@ function ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) w node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (2, n_nodes, n_nodes, capacity)) - _surface_flux_values = fill(nan_uEltype, NVARS * n_nodes * 2 * 2 * capacity) + _surface_flux_values = fill(nan_uEltype, n_variables * n_nodes * 2 * 2 * capacity) surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values), - (NVARS, n_nodes, 2 * 2, capacity)) + (n_variables, n_nodes, 2 * 2, capacity)) cell_ids = fill(typemin(Int), capacity) - return ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}( + return ElementContainer2D{RealT, uEltype}( inverse_jacobian, node_coordinates, surface_flux_values, cell_ids, _node_coordinates, _surface_flux_values) end @@ -71,22 +69,31 @@ end @inline nelements(elements::ElementContainer2D) = length(elements.cell_ids) # TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements)) @inline eachelement(elements::ElementContainer2D) = Base.OneTo(nelements(elements)) +@inline Base.real(elements::ElementContainer2D) = eltype(elements.node_coordinates) # Create element container and initialize element data function init_elements(cell_ids, mesh::TreeMesh2D, - equations::AbstractEquations{2, NVARS}, - basis::LobattoLegendreBasis{T, NNODES}, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real, NVARS, T, NNODES} + equations::AbstractEquations{2}, + basis, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real} # Initialize container n_elements = length(cell_ids) - elements = ElementContainer2D{RealT, uEltype, NVARS, NNODES-1}(n_elements) + elements = ElementContainer2D{RealT, uEltype}( + n_elements, nvariables(equations), nnodes(basis)) - init_elements!(elements, cell_ids, mesh, basis.nodes) + init_elements!(elements, cell_ids, mesh, basis) return elements end -function init_elements!(elements, cell_ids, mesh::TreeMesh2D, nodes) - n_nodes = length(nodes) +function init_elements!(elements, cell_ids, mesh::TreeMesh2D, basis) + nodes = get_nodes(basis) + # Compute the length of the 1D reference interval by integrating + # the function with constant value unity on the corresponding + # element data type (using \circ) + reference_length = integrate(one ∘ eltype, nodes, basis) + # Compute the offset of the midpoint of the 1D reference interval + # (its difference from zero) + reference_offset = first(nodes) + reference_length / 2 # Store cell ids elements.cell_ids .= cell_ids @@ -99,15 +106,19 @@ function init_elements!(elements, cell_ids, mesh::TreeMesh2D, nodes) # Get cell length dx = length_at_cell(mesh.tree, cell_id) - # Calculate inverse Jacobian as 1/(h/2) - elements.inverse_jacobian[element] = 2/dx + # Calculate inverse Jacobian + jacobian = dx / reference_length + elements.inverse_jacobian[element] = inv(jacobian) # Calculate node coordinates - for j in 1:n_nodes, i in 1:n_nodes + # Note that the `tree_coordinates` are the midpoints of the cells. + # Hence, we need to add an offset for `nodes` with a midpoint + # different from zero. + for j in eachnode(basis), i in eachnode(basis) elements.node_coordinates[1, i, j, element] = ( - mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i]) + mesh.tree.coordinates[1, cell_id] + jacobian * (nodes[i] - reference_offset)) elements.node_coordinates[2, i, j, element] = ( - mesh.tree.coordinates[2, cell_id] + dx/2 * nodes[j]) + mesh.tree.coordinates[2, cell_id] + jacobian * (nodes[j] - reference_offset)) end end @@ -117,32 +128,31 @@ end # Container data structure (structure-of-arrays style) for DG interfaces -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct InterfaceContainer2D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer - u::Array{uEltype, 4} # [leftright, variables, i, interfaces] - neighbor_ids::Matrix{Int} # [leftright, interfaces] - orientations::Vector{Int} # [interfaces] +mutable struct InterfaceContainer2D{uEltype<:Real} <: AbstractContainer + u::Array{uEltype, 4} # [leftright, variables, i, interfaces] + neighbor_ids::Array{Int, 2} # [leftright, interfaces] + orientations::Vector{Int} # [interfaces] # internal `resize!`able storage _u::Vector{uEltype} _neighbor_ids::Vector{Int} end -nvariables(::InterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS -polydeg(::InterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::InterfaceContainer2D{uEltype}) where {uEltype} = uEltype +nvariables(interfaces::InterfaceContainer2D) = size(interfaces.u, 2) +nnodes(interfaces::InterfaceContainer2D) = size(interfaces.u, 3) +Base.eltype(interfaces::InterfaceContainer2D) = eltype(interfaces.u) # See explanation of Base.resize! for the element container function Base.resize!(interfaces::InterfaceContainer2D, capacity) - n_nodes = polydeg(interfaces) + 1 - nvars = nvariables(interfaces) + n_nodes = nnodes(interfaces) + n_variables = nvariables(interfaces) @unpack _u, _neighbor_ids, orientations = interfaces - resize!(_u, 2 * nvars * n_nodes * capacity) + resize!(_u, 2 * n_variables * n_nodes * capacity) interfaces.u = unsafe_wrap(Array, pointer(_u), - (2, nvars, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) resize!(_neighbor_ids, 2 * capacity) - interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ), + interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, capacity)) resize!(orientations, capacity) @@ -151,23 +161,22 @@ function Base.resize!(interfaces::InterfaceContainer2D, capacity) end -function InterfaceContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function InterfaceContainer2D{uEltype}(capacity::Integer, n_variables, n_nodes) where {uEltype<:Real} nan = convert(uEltype, NaN) # Initialize fields with defaults - _u = fill(nan, 2 * NVARS * n_nodes * capacity) + _u = fill(nan, 2 * n_variables * n_nodes * capacity) u = unsafe_wrap(Array, pointer(_u), - (2, NVARS, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) _neighbor_ids = fill(typemin(Int), 2 * capacity) - neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ), + neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, capacity)) orientations = fill(typemin(Int), capacity) - return InterfaceContainer2D{uEltype, NVARS, POLYDEG}( + return InterfaceContainer2D{uEltype}( u, neighbor_ids, orientations, _u, _neighbor_ids) end @@ -179,10 +188,11 @@ end # Create interface container and initialize interface data in `elements`. function init_interfaces(cell_ids, mesh::TreeMesh2D, - elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer2D) # Initialize container n_interfaces = count_required_interfaces(mesh, cell_ids) - interfaces = InterfaceContainer2D{uEltype, NVARS, POLYDEG}(n_interfaces) + interfaces = InterfaceContainer2D{eltype(elements)}( + n_interfaces, nvariables(elements), nnodes(elements)) # Connect elements with interfaces init_interfaces!(interfaces, elements, mesh) @@ -281,8 +291,7 @@ end # Container data structure (structure-of-arrays style) for DG boundaries -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct BoundaryContainer2D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct BoundaryContainer2D{RealT<:Real, uEltype<:Real} <: AbstractContainer u::Array{uEltype, 4} # [leftright, variables, i, boundaries] neighbor_ids::Vector{Int} # [boundaries] orientations::Vector{Int} # [boundaries] @@ -294,20 +303,20 @@ mutable struct BoundaryContainer2D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} < _node_coordinates::Vector{RealT} end -nvariables(::BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS -polydeg(::BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::BoundaryContainer2D{RealT, uEltype}) where {RealT, uEltype} = uEltype +nvariables(boundaries::BoundaryContainer2D) = size(boundaries.u, 2) +nnodes(boundaries::BoundaryContainer2D) = size(boundaries.u, 3) +Base.eltype(boundaries::BoundaryContainer2D) = eltype(boundaries.u) # See explanation of Base.resize! for the element container function Base.resize!(boundaries::BoundaryContainer2D, capacity) - n_nodes = polydeg(boundaries) + 1 - nvars = nvariables(boundaries) + n_nodes = nnodes(boundaries) + n_variables = nvariables(boundaries) @unpack _u, _node_coordinates, neighbor_ids, orientations, neighbor_sides = boundaries - resize!(_u, 2 * nvars * n_nodes * capacity) + resize!(_u, 2 * n_variables * n_nodes * capacity) boundaries.u = unsafe_wrap(Array, pointer(_u), - (2, nvars, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) resize!(_node_coordinates, 2 * n_nodes * capacity) boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), @@ -323,15 +332,14 @@ function Base.resize!(boundaries::BoundaryContainer2D, capacity) end -function BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function BoundaryContainer2D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real} nan_RealT = convert(RealT, NaN) nan_uEltype = convert(uEltype, NaN) # Initialize fields with defaults - _u = fill(nan_uEltype, 2 * NVARS * n_nodes * capacity) + _u = fill(nan_uEltype, 2 * n_variables * n_nodes * capacity) u = unsafe_wrap(Array, pointer(_u), - (2, NVARS, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) neighbor_ids = fill(typemin(Int), capacity) @@ -345,7 +353,7 @@ function BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) n_boundaries_per_direction = SVector(0, 0, 0, 0) - return BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}( + return BoundaryContainer2D{RealT, uEltype}( u, neighbor_ids, orientations, neighbor_sides, node_coordinates, n_boundaries_per_direction, _u, _node_coordinates) @@ -358,10 +366,11 @@ end # Create boundaries container and initialize boundary data in `elements`. function init_boundaries(cell_ids, mesh::TreeMesh2D, - elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer2D) # Initialize container n_boundaries = count_required_boundaries(mesh, cell_ids) - boundaries = BoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(n_boundaries) + boundaries = BoundaryContainer2D{real(elements), eltype(elements)}( + n_boundaries, nvariables(elements), nnodes(elements)) # Connect elements with boundaries init_boundaries!(boundaries, elements, mesh) @@ -479,12 +488,10 @@ end # | | # lower = 1 | | # | | -# TODO: Taal refactor, remove NVARS, POLYDEG? -# TODO: Taal refactor, mutable struct or resize! for AMR? -mutable struct L2MortarContainer2D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct L2MortarContainer2D{uEltype<:Real} <: AbstractContainer u_upper::Array{uEltype, 4} # [leftright, variables, i, mortars] u_lower::Array{uEltype, 4} # [leftright, variables, i, mortars] - neighbor_ids::Matrix{Int} # [position, mortars] + neighbor_ids::Array{Int, 2} # [position, mortars] # Large sides: left -> 1, right -> 2 large_sides::Vector{Int} # [mortars] orientations::Vector{Int} # [mortars] @@ -494,27 +501,27 @@ mutable struct L2MortarContainer2D{uEltype<:Real, NVARS, POLYDEG} <: AbstractCon _neighbor_ids::Vector{Int} end -nvariables(::L2MortarContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS -polydeg(::L2MortarContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::L2MortarContainer2D{uEltype}) where {uEltype} = uEltype +nvariables(mortars::L2MortarContainer2D) = size(mortars.u_upper, 2) +nnodes(mortars::L2MortarContainer2D) = size(mortars.u_upper, 3) +Base.eltype(mortars::L2MortarContainer2D) = eltype(mortars.u_upper) # See explanation of Base.resize! for the element container function Base.resize!(mortars::L2MortarContainer2D, capacity) - n_nodes = polydeg(mortars) + 1 - nvars = nvariables(mortars) + n_nodes = nnodes(mortars) + n_variables = nvariables(mortars) @unpack _u_upper, _u_lower, _neighbor_ids, large_sides, orientations = mortars - resize!(_u_upper, 2 * nvars * n_nodes * capacity) + resize!(_u_upper, 2 * n_variables * n_nodes * capacity) mortars.u_upper = unsafe_wrap(Array, pointer(_u_upper), - (2, nvars, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) - resize!(_u_lower, 2 * nvars * n_nodes * capacity) + resize!(_u_lower, 2 * n_variables * n_nodes * capacity) mortars.u_lower = unsafe_wrap(Array, pointer(_u_lower), - (2, nvars, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) resize!(_neighbor_ids, 3 * capacity) - mortars.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ), + mortars.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (3, capacity)) resize!(large_sides, capacity) @@ -525,18 +532,17 @@ function Base.resize!(mortars::L2MortarContainer2D, capacity) end -function L2MortarContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function L2MortarContainer2D{uEltype}(capacity::Integer, n_variables, n_nodes) where {uEltype<:Real} nan = convert(uEltype, NaN) # Initialize fields with defaults - _u_upper = fill(nan, 2 * NVARS * n_nodes * capacity) + _u_upper = fill(nan, 2 * n_variables * n_nodes * capacity) u_upper = unsafe_wrap(Array, pointer(_u_upper), - (2, NVARS, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) - _u_lower = fill(nan, 2 * NVARS * n_nodes * capacity) + _u_lower = fill(nan, 2 * n_variables * n_nodes * capacity) u_lower = unsafe_wrap(Array, pointer(_u_lower), - (2, NVARS, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) _neighbor_ids = fill(typemin(Int), 3 * capacity) neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), @@ -546,7 +552,7 @@ function L2MortarContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) where { orientations = fill(typemin(Int), capacity) - return L2MortarContainer2D{uEltype, NVARS, POLYDEG}( + return L2MortarContainer2D{uEltype}( u_upper, u_lower, neighbor_ids, large_sides, orientations, _u_upper, _u_lower, _neighbor_ids) end @@ -576,11 +582,12 @@ end # Create mortar container and initialize mortar data in `elements`. function init_mortars(cell_ids, mesh::TreeMesh2D, - elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}, - ::LobattoLegendreMortarL2) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer2D, + ::LobattoLegendreMortarL2) # Initialize containers n_mortars = count_required_mortars(mesh, cell_ids) - mortars = L2MortarContainer2D{uEltype, NVARS, POLYDEG}(n_mortars) + mortars = L2MortarContainer2D{eltype(elements)}( + n_mortars, nvariables(elements), nnodes(elements)) # Connect elements with mortars init_mortars!(mortars, elements, mesh) @@ -686,7 +693,7 @@ end # Container data structure (structure-of-arrays style) for DG MPI interfaces -mutable struct MPIInterfaceContainer2D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct MPIInterfaceContainer2D{uEltype<:Real} <: AbstractContainer u::Array{uEltype, 4} # [leftright, variables, i, interfaces] local_element_ids::Vector{Int} # [interfaces] orientations::Vector{Int} # [interfaces] @@ -695,19 +702,19 @@ mutable struct MPIInterfaceContainer2D{uEltype<:Real, NVARS, POLYDEG} <: Abstrac _u::Vector{uEltype} end -nvariables(::MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS -polydeg(::MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::MPIInterfaceContainer2D{uEltype}) where {uEltype} = uEltype +nvariables(mpi_interfaces::MPIInterfaceContainer2D) = size(mpi_interfaces.u, 2) +nnodes(mpi_interfaces::MPIInterfaceContainer2D) = size(mpi_interfaces.u, 3) +Base.eltype(mpi_interfaces::MPIInterfaceContainer2D) = eltype(mpi_interfaces.u) # See explanation of Base.resize! for the element container function Base.resize!(mpi_interfaces::MPIInterfaceContainer2D, capacity) - n_nodes = polydeg(mpi_interfaces) + 1 - nvars = nvariables(mpi_interfaces) + n_nodes = nnodes(mpi_interfaces) + n_variables = nvariables(mpi_interfaces) @unpack _u, local_element_ids, orientations, remote_sides = mpi_interfaces - resize!(_u, 2 * nvars * n_nodes * capacity) + resize!(_u, 2 * n_variables * n_nodes * capacity) mpi_interfaces.u = unsafe_wrap(Array, pointer(_u), - (2, nvars, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) resize!(local_element_ids, capacity) @@ -719,14 +726,13 @@ function Base.resize!(mpi_interfaces::MPIInterfaceContainer2D, capacity) end -function MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function MPIInterfaceContainer2D{uEltype}(capacity::Integer, n_variables, n_nodes) where {uEltype<:Real} nan = convert(uEltype, NaN) # Initialize fields with defaults - _u = fill(nan, 2 * NVARS * n_nodes * capacity) + _u = fill(nan, 2 * n_variables * n_nodes * capacity) u = unsafe_wrap(Array, pointer(_u), - (2, NVARS, n_nodes, capacity)) + (2, n_variables, n_nodes, capacity)) local_element_ids = fill(typemin(Int), capacity) @@ -734,7 +740,7 @@ function MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}(capacity::Integer) whe remote_sides = fill(typemin(Int), capacity) - return MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}( + return MPIInterfaceContainer2D{uEltype}( u, local_element_ids, orientations, remote_sides, _u) end @@ -747,10 +753,11 @@ nmpiinterfaces(mpi_interfaces::MPIInterfaceContainer2D) = length(mpi_interfaces. # Create MPI interface container and initialize MPI interface data in `elements`. function init_mpi_interfaces(cell_ids, mesh::TreeMesh2D, - elements::ElementContainer2D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer2D) # Initialize container n_mpi_interfaces = count_required_mpi_interfaces(mesh, cell_ids) - mpi_interfaces = MPIInterfaceContainer2D{uEltype, NVARS, POLYDEG}(n_mpi_interfaces) + mpi_interfaces = MPIInterfaceContainer2D{eltype(elements)}( + n_mpi_interfaces, nvariables(elements), nnodes(elements)) # Connect elements with interfaces init_mpi_interfaces!(mpi_interfaces, elements, mesh) diff --git a/src/solvers/dg_tree/containers_3d.jl b/src/solvers/dg_tree/containers_3d.jl index 0fbb81bc2c..7c5dfaf437 100644 --- a/src/solvers/dg_tree/containers_3d.jl +++ b/src/solvers/dg_tree/containers_3d.jl @@ -1,7 +1,6 @@ # Container data structure (structure-of-arrays style) for DG elements -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct ElementContainer3D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct ElementContainer3D{RealT<:Real, uEltype<:Real} <: AbstractContainer inverse_jacobian::Vector{RealT} # [elements] node_coordinates::Array{RealT, 5} # [orientation, i, j, k, elements] surface_flux_values::Array{uEltype, 5} # [variables, i, j, direction, elements] @@ -11,9 +10,9 @@ mutable struct ElementContainer3D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: _surface_flux_values::Vector{uEltype} end -nvariables(::ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS -polydeg(::ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::ElementContainer3D{RealT, uEltype}) where {RealT, uEltype} = uEltype +nvariables(elements::ElementContainer3D) = size(elements.surface_flux_values, 1) +nnodes(elements::ElementContainer3D) = size(elements.node_coordinates, 2) +Base.eltype(elements::ElementContainer3D) = eltype(elements.surface_flux_values) # Only one-dimensional `Array`s are `resize!`able in Julia. # Hence, we use `Vector`s as internal storage and `resize!` @@ -21,8 +20,8 @@ Base.eltype(::ElementContainer3D{RealT, uEltype}) where {RealT, uEltype} = uElty # `unsafe_wrap`ping multi-dimensional `Array`s around the # internal storage. function Base.resize!(elements::ElementContainer3D, capacity) - n_nodes = polydeg(elements) + 1 - nvars = nvariables(elements) + n_nodes = nnodes(elements) + n_variables = nvariables(elements) @unpack _node_coordinates, _surface_flux_values, inverse_jacobian, cell_ids = elements @@ -32,9 +31,9 @@ function Base.resize!(elements::ElementContainer3D, capacity) elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (3, n_nodes, n_nodes, n_nodes, capacity)) - resize!(_surface_flux_values, nvars * n_nodes * n_nodes * 2 * 3 * capacity) + resize!(_surface_flux_values, n_variables * n_nodes * n_nodes * 2 * 3 * capacity) elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values), - (nvars, n_nodes, n_nodes, 2 * 3, capacity)) + (n_variables, n_nodes, n_nodes, 2 * 3, capacity)) resize!(cell_ids, capacity) @@ -42,8 +41,7 @@ function Base.resize!(elements::ElementContainer3D, capacity) end -function ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function ElementContainer3D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real} nan_RealT = convert(RealT, NaN) nan_uEltype = convert(uEltype, NaN) @@ -54,14 +52,14 @@ function ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) w node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (3, n_nodes, n_nodes, n_nodes, capacity)) - _surface_flux_values = fill(nan_uEltype, NVARS * n_nodes * n_nodes * 2 * 3 * capacity) + _surface_flux_values = fill(nan_uEltype, n_variables * n_nodes * n_nodes * 2 * 3 * capacity) surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values), - (NVARS, n_nodes, n_nodes, 2 * 3, capacity)) + (n_variables, n_nodes, n_nodes, 2 * 3, capacity)) cell_ids = fill(typemin(Int), capacity) - return ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}( + return ElementContainer3D{RealT, uEltype}( inverse_jacobian, node_coordinates, surface_flux_values, cell_ids, _node_coordinates, _surface_flux_values) end @@ -71,22 +69,31 @@ end nelements(elements::ElementContainer3D) = length(elements.cell_ids) # TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements)) @inline eachelement(elements::ElementContainer3D) = Base.OneTo(nelements(elements)) +@inline Base.real(elements::ElementContainer3D) = eltype(elements.node_coordinates) # Create element container and initialize element data function init_elements(cell_ids, mesh::TreeMesh3D, - equations::AbstractEquations{3, NVARS}, - basis::LobattoLegendreBasis{T, NNODES}, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real, NVARS, T, NNODES} + equations::AbstractEquations{3}, + basis, ::Type{RealT}, ::Type{uEltype}) where {RealT<:Real, uEltype<:Real} # Initialize container n_elements = length(cell_ids) - elements = ElementContainer3D{RealT, uEltype, NVARS, NNODES-1}(n_elements) + elements = ElementContainer3D{RealT, uEltype}( + n_elements, nvariables(equations), nnodes(basis)) - init_elements!(elements, cell_ids, mesh, basis.nodes) + init_elements!(elements, cell_ids, mesh, basis) return elements end -function init_elements!(elements, cell_ids, mesh::TreeMesh3D, nodes) - n_nodes = length(nodes) +function init_elements!(elements, cell_ids, mesh::TreeMesh3D, basis) + nodes = get_nodes(basis) + # Compute the length of the 1D reference interval by integrating + # the function with constant value unity on the corresponding + # element data type (using \circ) + reference_length = integrate(one ∘ eltype, nodes, basis) + # Compute the offset of the midpoint of the 1D reference interval + # (its difference from zero) + reference_offset = first(nodes) + reference_length / 2 # Store cell ids elements.cell_ids .= cell_ids @@ -99,17 +106,21 @@ function init_elements!(elements, cell_ids, mesh::TreeMesh3D, nodes) # Get cell length dx = length_at_cell(mesh.tree, cell_id) - # Calculate inverse Jacobian as 1/(h/2) - elements.inverse_jacobian[element] = 2/dx + # Calculate inverse Jacobian + jacobian = dx / reference_length + elements.inverse_jacobian[element] = inv(jacobian) # Calculate node coordinates - for k in 1:n_nodes, j in 1:n_nodes, i in 1:n_nodes + # Note that the `tree_coordinates` are the midpoints of the cells. + # Hence, we need to add an offset for `nodes` with a midpoint + # different from zero. + for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) elements.node_coordinates[1, i, j, k, element] = ( - mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i]) + mesh.tree.coordinates[1, cell_id] + jacobian * (nodes[i] - reference_offset)) elements.node_coordinates[2, i, j, k, element] = ( - mesh.tree.coordinates[2, cell_id] + dx/2 * nodes[j]) + mesh.tree.coordinates[2, cell_id] + jacobian * (nodes[j] - reference_offset)) elements.node_coordinates[3, i, j, k, element] = ( - mesh.tree.coordinates[3, cell_id] + dx/2 * nodes[k]) + mesh.tree.coordinates[3, cell_id] + jacobian * (nodes[k] - reference_offset)) end end @@ -119,8 +130,7 @@ end # Container data structure (structure-of-arrays style) for DG interfaces -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct InterfaceContainer3D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct InterfaceContainer3D{uEltype<:Real} <: AbstractContainer u::Array{uEltype, 5} # [leftright, variables, i, j, interfaces] neighbor_ids::Matrix{Int} # [leftright, interfaces] orientations::Vector{Int} # [interfaces] @@ -129,19 +139,19 @@ mutable struct InterfaceContainer3D{uEltype<:Real, NVARS, POLYDEG} <: AbstractCo _neighbor_ids::Vector{Int} end -nvariables(::InterfaceContainer3D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS -polydeg(::InterfaceContainer3D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::InterfaceContainer3D{uEltype}) where {uEltype} = uEltype +nvariables(interfaces::InterfaceContainer3D) = size(interfaces.u, 2) +nnodes(interfaces::InterfaceContainer3D) = size(interfaces.u, 3) +Base.eltype(interfaces::InterfaceContainer3D) = eltype(interfaces.u) # See explanation of Base.resize! for the element container function Base.resize!(interfaces::InterfaceContainer3D, capacity) - n_nodes = polydeg(interfaces) + 1 - nvars = nvariables(interfaces) + n_nodes = nnodes(interfaces) + n_variables = nvariables(interfaces) @unpack _u, _neighbor_ids, orientations = interfaces - resize!(_u, 2 * nvars * n_nodes * n_nodes * capacity) + resize!(_u, 2 * n_variables * n_nodes * n_nodes * capacity) interfaces.u = unsafe_wrap(Array, pointer(_u), - (2, nvars, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) resize!(_neighbor_ids, 2 * capacity) interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), @@ -153,23 +163,22 @@ function Base.resize!(interfaces::InterfaceContainer3D, capacity) end -function InterfaceContainer3D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function InterfaceContainer3D{uEltype}(capacity::Integer, n_variables, n_nodes) where {uEltype<:Real} nan = convert(uEltype, NaN) # Initialize fields with defaults - _u = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity) + _u = fill(nan, 2 * n_variables * n_nodes * n_nodes * capacity) u = unsafe_wrap(Array, pointer(_u), - (2, NVARS, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) _neighbor_ids = fill(typemin(Int), 2 * capacity) - neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ), + neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, capacity)) orientations = fill(typemin(Int), capacity) - return InterfaceContainer3D{uEltype, NVARS, POLYDEG}( + return InterfaceContainer3D{uEltype}( u, neighbor_ids, orientations, _u, _neighbor_ids) end @@ -181,10 +190,11 @@ ninterfaces(interfaces::InterfaceContainer3D) = length(interfaces.orientations) # Create interface container and initialize interface data in `elements`. function init_interfaces(cell_ids, mesh::TreeMesh3D, - elements::ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer3D) # Initialize container n_interfaces = count_required_interfaces(mesh, cell_ids) - interfaces = InterfaceContainer3D{uEltype, NVARS, POLYDEG}(n_interfaces) + interfaces = InterfaceContainer3D{eltype(elements)}( + n_interfaces, nvariables(elements), nnodes(elements)) # Connect elements with interfaces init_interfaces!(interfaces, elements, mesh) @@ -279,8 +289,7 @@ end # Container data structure (structure-of-arrays style) for DG boundaries -# TODO: Taal refactor, remove NVARS, POLYDEG? -mutable struct BoundaryContainer3D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct BoundaryContainer3D{RealT<:Real, uEltype<:Real} <: AbstractContainer u::Array{uEltype, 5} # [leftright, variables, i, j, boundaries] neighbor_ids::Vector{Int} # [boundaries] orientations::Vector{Int} # [boundaries] @@ -292,23 +301,23 @@ mutable struct BoundaryContainer3D{RealT<:Real, uEltype<:Real, NVARS, POLYDEG} < _node_coordinates::Vector{RealT} end -nvariables(::BoundaryContainer3D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = NVARS -polydeg(::BoundaryContainer3D{RealT, uEltype, NVARS, POLYDEG}) where {RealT, uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::BoundaryContainer3D{RealT, uEltype}) where {RealT, uEltype} = uEltype +nvariables(boundaries::BoundaryContainer3D) = size(boundaries.u, 2) +nnodes(boundaries::BoundaryContainer3D) = size(boundaries.u, 3) +Base.eltype(boundaries::BoundaryContainer3D) = eltype(boundaries.u) # See explanation of Base.resize! for the element container function Base.resize!(boundaries::BoundaryContainer3D, capacity) - n_nodes = polydeg(boundaries) + 1 - nvars = nvariables(boundaries) + n_nodes = nnodes(boundaries) + n_variables = nvariables(boundaries) @unpack _u, _node_coordinates, neighbor_ids, orientations, neighbor_sides = boundaries - resize!(_u, 2 * nvars * n_nodes * n_nodes * capacity) + resize!(_u, 2 * n_variables * n_nodes * n_nodes * capacity) boundaries.u = unsafe_wrap(Array, pointer(_u), - (2, nvars, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) resize!(_node_coordinates, 3 * n_nodes * n_nodes * capacity) - boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates ), + boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), (3, n_nodes, n_nodes, capacity)) resize!(neighbor_ids, capacity) @@ -321,15 +330,14 @@ function Base.resize!(boundaries::BoundaryContainer3D, capacity) end -function BoundaryContainer3D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function BoundaryContainer3D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real} nan_RealT = convert(RealT, NaN) nan_uEltype = convert(uEltype, NaN) # Initialize fields with defaults - _u = fill(nan_uEltype, 2 * NVARS * n_nodes * n_nodes * capacity) + _u = fill(nan_uEltype, 2 * n_variables * n_nodes * n_nodes * capacity) u = unsafe_wrap(Array, pointer(_u), - (2, NVARS, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) neighbor_ids = fill(typemin(Int), capacity) @@ -343,7 +351,7 @@ function BoundaryContainer3D{RealT, uEltype, NVARS, POLYDEG}(capacity::Integer) n_boundaries_per_direction = SVector(0, 0, 0, 0, 0, 0) - return BoundaryContainer3D{RealT, uEltype, NVARS, POLYDEG}( + return BoundaryContainer3D{RealT, uEltype}( u, neighbor_ids, orientations, neighbor_sides, node_coordinates, n_boundaries_per_direction, _u, _node_coordinates) @@ -356,10 +364,11 @@ nboundaries(boundaries::BoundaryContainer3D) = length(boundaries.orientations) # Create boundaries container and initialize boundary data in `elements`. function init_boundaries(cell_ids, mesh::TreeMesh3D, - elements::ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer3D) # Initialize container n_boundaries = count_required_boundaries(mesh, cell_ids) - boundaries = BoundaryContainer3D{RealT, uEltype, NVARS, POLYDEG}(n_boundaries) + boundaries = BoundaryContainer3D{real(elements), eltype(elements)}( + n_boundaries, nvariables(elements), nnodes(elements)) # Connect elements with boundaries init_boundaries!(boundaries, elements, mesh) @@ -495,12 +504,12 @@ end # # Left and right are used *both* for the numbering of the mortar faces *and* for the position of the # elements with respect to the axis orthogonal to the mortar. -mutable struct L2MortarContainer3D{uEltype<:Real, NVARS, POLYDEG} <: AbstractContainer +mutable struct L2MortarContainer3D{uEltype<:Real} <: AbstractContainer u_upper_left ::Array{uEltype, 5} # [leftright, variables, i, j, mortars] u_upper_right::Array{uEltype, 5} # [leftright, variables, i, j, mortars] u_lower_left ::Array{uEltype, 5} # [leftright, variables, i, j, mortars] u_lower_right::Array{uEltype, 5} # [leftright, variables, i, j, mortars] - neighbor_ids ::Matrix{Int} # [position, mortars] + neighbor_ids ::Array{Int, 2} # [position, mortars] # Large sides: left -> 1, right -> 2 large_sides ::Vector{Int} # [mortars] orientations::Vector{Int} # [mortars] @@ -512,35 +521,35 @@ mutable struct L2MortarContainer3D{uEltype<:Real, NVARS, POLYDEG} <: AbstractCon _neighbor_ids ::Vector{Int} end -nvariables(::L2MortarContainer3D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = NVARS -polydeg(::L2MortarContainer3D{uEltype, NVARS, POLYDEG}) where {uEltype, NVARS, POLYDEG} = POLYDEG -Base.eltype(::L2MortarContainer3D{uEltype}) where {uEltype} = uEltype +nvariables(mortars::L2MortarContainer3D) = size(mortars.u_upper_left, 2) +nnodes(mortars::L2MortarContainer3D) = size(mortars.u_upper_left, 3) +Base.eltype(mortars::L2MortarContainer3D) = eltype(mortars.u_upper_left) # See explanation of Base.resize! for the element container function Base.resize!(mortars::L2MortarContainer3D, capacity) - n_nodes = polydeg(mortars) + 1 - nvars = nvariables(mortars) + n_nodes = nnodes(mortars) + n_variables = nvariables(mortars) @unpack _u_upper_left, _u_upper_right, _u_lower_left, _u_lower_right, _neighbor_ids, large_sides, orientations = mortars - resize!(_u_upper_left, 2 * nvars * n_nodes * n_nodes * capacity) + resize!(_u_upper_left, 2 * n_variables * n_nodes * n_nodes * capacity) mortars.u_upper_left = unsafe_wrap(Array, pointer(_u_upper_left), - (2, nvars, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) - resize!(_u_upper_right, 2 * nvars * n_nodes * n_nodes * capacity) + resize!(_u_upper_right, 2 * n_variables * n_nodes * n_nodes * capacity) mortars.u_upper_right = unsafe_wrap(Array, pointer(_u_upper_right), - (2, nvars, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) - resize!(_u_lower_left, 2 * nvars * n_nodes * n_nodes * capacity) + resize!(_u_lower_left, 2 * n_variables * n_nodes * n_nodes * capacity) mortars.u_lower_left = unsafe_wrap(Array, pointer(_u_lower_left), - (2, nvars, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) - resize!(_u_lower_right, 2 * nvars * n_nodes * n_nodes * capacity) + resize!(_u_lower_right, 2 * n_variables * n_nodes * n_nodes * capacity) mortars.u_lower_right = unsafe_wrap(Array, pointer(_u_lower_right), - (2, nvars, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) resize!(_neighbor_ids, 5 * capacity) - mortars.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ), + mortars.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (5, capacity)) resize!(large_sides, capacity) @@ -551,26 +560,25 @@ function Base.resize!(mortars::L2MortarContainer3D, capacity) end -function L2MortarContainer3D{uEltype, NVARS, POLYDEG}(capacity::Integer) where {uEltype<:Real, NVARS, POLYDEG} - n_nodes = POLYDEG + 1 +function L2MortarContainer3D{uEltype}(capacity::Integer, n_variables, n_nodes) where {uEltype<:Real} nan = convert(uEltype, NaN) # Initialize fields with defaults - _u_upper_left = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity) + _u_upper_left = fill(nan, 2 * n_variables * n_nodes * n_nodes * capacity) u_upper_left = unsafe_wrap(Array, pointer(_u_upper_left), - (2, NVARS, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) - _u_upper_right = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity) + _u_upper_right = fill(nan, 2 * n_variables * n_nodes * n_nodes * capacity) u_upper_right = unsafe_wrap(Array, pointer(_u_upper_right), - (2, NVARS, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) - _u_lower_left = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity) + _u_lower_left = fill(nan, 2 * n_variables * n_nodes * n_nodes * capacity) u_lower_left = unsafe_wrap(Array, pointer(_u_lower_left), - (2, NVARS, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) - _u_lower_right = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity) + _u_lower_right = fill(nan, 2 * n_variables * n_nodes * n_nodes * capacity) u_lower_right = unsafe_wrap(Array, pointer(_u_lower_right), - (2, NVARS, n_nodes, n_nodes, capacity)) + (2, n_variables, n_nodes, n_nodes, capacity)) _neighbor_ids = fill(typemin(Int), 5 * capacity) neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), @@ -580,7 +588,7 @@ function L2MortarContainer3D{uEltype, NVARS, POLYDEG}(capacity::Integer) where { orientations = fill(typemin(Int), capacity) - return L2MortarContainer3D{uEltype, NVARS, POLYDEG}( + return L2MortarContainer3D{uEltype}( u_upper_left, u_upper_right, u_lower_left, u_lower_right, neighbor_ids, large_sides, orientations, @@ -620,11 +628,12 @@ end # Create mortar container and initialize mortar data in `elements`. function init_mortars(cell_ids, mesh::TreeMesh3D, - elements::ElementContainer3D{RealT, uEltype, NVARS, POLYDEG}, - mortar::LobattoLegendreMortarL2) where {RealT<:Real, uEltype<:Real, NVARS, POLYDEG} + elements::ElementContainer3D, + mortar::LobattoLegendreMortarL2) # Initialize containers n_mortars = count_required_mortars(mesh, cell_ids) - mortars = L2MortarContainer3D{uEltype, NVARS, POLYDEG}(n_mortars) + mortars = L2MortarContainer3D{eltype(elements)}( + n_mortars, nvariables(elements), nnodes(elements)) # Connect elements with mortars init_mortars!(mortars, elements, mesh) diff --git a/src/solvers/dg_tree/dg.jl b/src/solvers/dg_tree/dg.jl index 5fc4e3c98e..6675acf99a 100644 --- a/src/solvers/dg_tree/dg.jl +++ b/src/solvers/dg_tree/dg.jl @@ -4,6 +4,13 @@ abstract type AbstractVolumeIntegral end get_element_variables!(element_variables, u, mesh, equations, volume_integral::AbstractVolumeIntegral, dg, cache) = nothing +""" + VolumeIntegralStrongForm + +The classical strong form volume integral type for FD/DG methods. +""" +struct VolumeIntegralStrongForm <: AbstractVolumeIntegral end + """ VolumeIntegralWeakForm @@ -132,17 +139,78 @@ end +abstract type AbstractSurfaceIntegral end + """ - DG(basis, basis, surface_flux, volume_integral) + SurfaceIntegralWeakForm(surface_flux=flux_central) + +The classical weak form surface integral type for DG methods as explained in standard +textbooks such as +- Kopriva (2009) + Implementing Spectral Methods for Partial Differential Equations: + Algorithms for Scientists and Engineers + [doi: 10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5) + +See also [`VolumeIntegralWeakForm`](@ref). +""" +struct SurfaceIntegralWeakForm{SurfaceFlux} <: AbstractSurfaceIntegral + surface_flux::SurfaceFlux +end + +SurfaceIntegralWeakForm() = SurfaceIntegralWeakForm(flux_central) + +function Base.show(io::IO, ::MIME"text/plain", integral::SurfaceIntegralWeakForm) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + setup = [ + "surface flux" => integral.surface_flux + ] + summary_box(io, "SurfaceIntegralWeakForm", setup) + end +end + +""" + SurfaceIntegralStrongForm(surface_flux=flux_central) + +The classical strong form surface integral type for FD/DG methods. + +See also [`VolumeIntegralStrongForm`](@ref). +""" +struct SurfaceIntegralStrongForm{SurfaceFlux} <: AbstractSurfaceIntegral + surface_flux::SurfaceFlux +end + +SurfaceIntegralStrongForm() = SurfaceIntegralStrongForm(flux_central) + +function Base.show(io::IO, ::MIME"text/plain", integral::SurfaceIntegralStrongForm) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + setup = [ + "surface flux" => integral.surface_flux + ] + summary_box(io, "SurfaceIntegralStrongForm", setup) + end +end + + + +""" + DG(; basis, mortar, surface_integral, volume_integral) Create a discontinuous Galerkin method. If [`basis isa LobattoLegendreBasis`](@ref LobattoLegendreBasis), this creates a [`DGSEM`](@ref). """ -struct DG{Basis, Mortar, SurfaceFlux, VolumeIntegral} +struct DG{Basis, Mortar, SurfaceIntegral, VolumeIntegral} basis::Basis mortar::Mortar - surface_flux::SurfaceFlux + surface_integral::SurfaceIntegral volume_integral::VolumeIntegral end @@ -152,7 +220,7 @@ function Base.show(io::IO, dg::DG) print(io, "DG{", real(dg), "}(") print(io, dg.basis) print(io, ", ", dg.mortar) - print(io, ", ", dg.surface_flux) + print(io, ", ", dg.surface_integral) print(io, ", ", dg.volume_integral) print(io, ")") end @@ -164,10 +232,10 @@ function Base.show(io::IO, mime::MIME"text/plain", dg::DG) show(io, dg) else summary_header(io, "DG{" * string(real(dg)) * "}") - summary_line(io, "polynomial degree", polydeg(dg)) summary_line(io, "basis", dg.basis) summary_line(io, "mortar", dg.mortar) - summary_line(io, "surface flux", dg.surface_flux) + summary_line(io, "surface integral", dg.surface_integral |> typeof |> nameof) + show(increment_indent(io), mime, dg.surface_integral) summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof) if !(dg.volume_integral isa VolumeIntegralWeakForm) show(increment_indent(io), mime, dg.volume_integral) @@ -176,11 +244,9 @@ function Base.show(io::IO, mime::MIME"text/plain", dg::DG) end end -@inline Base.real(dg::DG) = real(dg.basis) +Base.summary(io::IO, dg::DG) = print(io, "DG(" * summary(dg.basis) * ")") -# TODO: Taal refactor, use case? -# Deprecate in favor of nnodes or order_of_accuracy? -@inline polydeg(dg::DG) = polydeg(dg.basis) +@inline Base.real(dg::DG) = real(dg.basis) @inline ndofs(mesh::TreeMesh, dg::DG, cache) = nelements(cache.elements) * nnodes(dg)^ndims(mesh) @@ -190,7 +256,7 @@ function get_element_variables!(element_variables, u, mesh, equations, dg::DG, c end -# TODO: Taal performance, 1:nnodes(dg) vs. Base.OneTo(nnodes(dg)) vs. SOneTo(nnodes(dg)) +# TODO: Taal performance, 1:nnodes(dg) vs. Base.OneTo(nnodes(dg)) vs. SOneTo(nnodes(dg)) for DGSEM @inline eachnode(dg::DG) = Base.OneTo(nnodes(dg)) @inline eachelement(dg::DG, cache) = Base.OneTo(nelements(dg, cache)) @inline eachinterface(dg::DG, cache) = Base.OneTo(ninterfaces(dg, cache)) @@ -264,35 +330,49 @@ include("l2projection.jl") include("basis_lobatto_legendre.jl") """ - DGSEM([RealT=Float64,] polydeg::Integer, - surface_flux=flux_central, - volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm(), - mortar=MortarL2(basis)) + DGSEM(; RealT=Float64, polydeg::Integer, + surface_flux=flux_central, + surface_integral=SurfaceIntegralWeakForm(surface_flux), + volume_integral=VolumeIntegralWeakForm(), + mortar=MortarL2(basis)) Create a discontinuous Galerkin spectral element method (DGSEM) using a [`LobattoLegendreBasis`](@ref) with polynomials of degree `polydeg`. """ -const DGSEM = DG{Basis, Mortar, SurfaceFlux, VolumeIntegral} where {Basis<:LobattoLegendreBasis, Mortar, SurfaceFlux, VolumeIntegral} +const DGSEM = DG{Basis} where {Basis<:LobattoLegendreBasis} +# TODO: Deprecated in v0.3 (no longer documented) function DGSEM(basis::LobattoLegendreBasis, surface_flux=flux_central, - volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm(), + volume_integral=VolumeIntegralWeakForm(), + mortar=MortarL2(basis)) + + surface_integral = SurfaceIntegralWeakForm(surface_flux) + return DG{typeof(basis), typeof(mortar), typeof(surface_integral), typeof(volume_integral)}( + basis, mortar, surface_integral, volume_integral) +end + +# TODO: Deprecated in v0.3 (no longer documented) +function DGSEM(basis::LobattoLegendreBasis, + surface_integral::AbstractSurfaceIntegral, + volume_integral=VolumeIntegralWeakForm(), mortar=MortarL2(basis)) - return DG{typeof(basis), typeof(mortar), typeof(surface_flux), typeof(volume_integral)}( - basis, mortar, surface_flux, volume_integral) + return DG{typeof(basis), typeof(mortar), typeof(surface_integral), typeof(volume_integral)}( + basis, mortar, surface_integral, volume_integral) end +# TODO: Deprecated in v0.3 (no longer documented) function DGSEM(RealT, polydeg::Integer, surface_flux=flux_central, - volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm(), + volume_integral=VolumeIntegralWeakForm(), mortar=MortarL2(LobattoLegendreBasis(RealT, polydeg))) basis = LobattoLegendreBasis(RealT, polydeg) return DGSEM(basis, surface_flux, volume_integral, mortar) end -DGSEM(polydeg, surface_flux=flux_central, volume_integral::AbstractVolumeIntegral=VolumeIntegralWeakForm()) = DGSEM(Float64, polydeg, surface_flux, volume_integral) +DGSEM(polydeg, surface_flux=flux_central, volume_integral=VolumeIntegralWeakForm()) = DGSEM(Float64, polydeg, surface_flux, volume_integral) # The constructor using only keyword arguments is convenient for elixirs since # it allows to modify the polynomial degree and other parameters via @@ -300,11 +380,16 @@ DGSEM(polydeg, surface_flux=flux_central, volume_integral::AbstractVolumeIntegra function DGSEM(; RealT=Float64, polydeg::Integer, surface_flux=flux_central, + surface_integral=SurfaceIntegralWeakForm(surface_flux), volume_integral=VolumeIntegralWeakForm()) basis = LobattoLegendreBasis(RealT, polydeg) - return DGSEM(basis, surface_flux, volume_integral) + return DGSEM(basis, surface_integral, volume_integral) end +@inline polydeg(dg::DGSEM) = polydeg(dg.basis) + +Base.summary(io::IO, dg::DGSEM) = print(io, "DGSEM(polydeg=$(polydeg(dg)))") + """ @@ -356,3 +441,4 @@ include("dg_2d_parallel.jl") # 3D DG implementation include("containers_3d.jl") include("dg_3d.jl") + diff --git a/src/solvers/dg_tree/dg_1d.jl b/src/solvers/dg_tree/dg_1d.jl index a37ccc10d1..f560be735a 100644 --- a/src/solvers/dg_tree/dg_1d.jl +++ b/src/solvers/dg_tree/dg_1d.jl @@ -95,25 +95,25 @@ function rhs!(du, u, t, # Prolong solution to interfaces @timed timer() "prolong2interfaces" prolong2interfaces!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate interface fluxes @timed timer() "interface flux" calc_interface_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg, cache) + dg.surface_integral, dg, cache) # Prolong solution to boundaries @timed timer() "prolong2boundaries" prolong2boundaries!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, mesh, equations, dg) + cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element @timed timer() "Jacobian" apply_jacobian!( @@ -283,7 +283,7 @@ end function prolong2interfaces!(cache, u, - mesh::TreeMesh{1}, equations, dg::DG) + mesh::TreeMesh{1}, equations, surface_integral, dg::DG) @unpack interfaces = cache @threaded for interface in eachinterface(dg, cache) @@ -303,8 +303,8 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, nonconservative_terms::Val{false}, equations, - dg::DG, cache) - @unpack surface_flux = dg + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @threaded for interface in eachinterface(dg, cache) @@ -337,9 +337,9 @@ end function prolong2boundaries!(cache, u, - mesh::TreeMesh{1}, equations, dg::DG) + mesh::TreeMesh{1}, equations, surface_integral, dg::DG) @unpack boundaries = cache - @unpack orientations, neighbor_sides = boundaries + @unpack neighbor_sides = boundaries @threaded for boundary in eachboundary(dg, cache) element = boundaries.neighbor_ids[boundary] @@ -362,13 +362,13 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::TreeMesh{1}, equations, dg::DG) + mesh::TreeMesh{1}, equations, surface_integral, dg::DG) @assert isempty(eachboundary(dg, cache)) end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition, - mesh::TreeMesh{1}, equations, dg::DG) + mesh::TreeMesh{1}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -379,13 +379,13 @@ function calc_boundary_flux!(cache, t, boundary_condition, # Calc boundary fluxes in each direction for direction in eachindex(firsts) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_condition, - equations, dg, cache, + equations, surface_integral, dg, cache, direction, firsts[direction], lasts[direction]) end end function calc_boundary_flux!(cache, t, boundary_conditions::Union{NamedTuple,Tuple}, - mesh::TreeMesh{1}, equations, dg::DG) + mesh::TreeMesh{1}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -395,15 +395,18 @@ function calc_boundary_flux!(cache, t, boundary_conditions::Union{NamedTuple,Tup # Calc boundary fluxes in each direction calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1], - equations, dg, cache, 1, firsts[1], lasts[1]) + equations, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2], - equations, dg, cache, 2, firsts[2], lasts[2]) + equations, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) end function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,3}, t, - boundary_condition, equations, dg::DG, cache, + boundary_condition, equations, + surface_integral, dg::DG, cache, direction, first_boundary, last_boundary) - @unpack surface_flux = dg + @unpack surface_flux = surface_integral @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries @threaded for boundary in first_boundary:last_boundary @@ -431,8 +434,8 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A end -function calc_surface_integral!(du, mesh::Union{TreeMesh{1}, CurvedMesh{1}}, - equations, dg::DGSEM, cache) +function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, CurvedMesh{1}}, + equations, surface_integral, dg::DGSEM, cache) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements diff --git a/src/solvers/dg_tree/dg_2d.jl b/src/solvers/dg_tree/dg_2d.jl index 393d7a0b06..952d21c5f5 100644 --- a/src/solvers/dg_tree/dg_2d.jl +++ b/src/solvers/dg_tree/dg_2d.jl @@ -127,35 +127,35 @@ function rhs!(du, u, t, # Prolong solution to interfaces @timed timer() "prolong2interfaces" prolong2interfaces!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate interface fluxes @timed timer() "interface flux" calc_interface_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg, cache) + dg.surface_integral, dg, cache) # Prolong solution to boundaries @timed timer() "prolong2boundaries" prolong2boundaries!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, mesh, equations, dg) + cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Prolong solution to mortars @timed timer() "prolong2mortars" prolong2mortars!( - cache, u, mesh, equations, dg.mortar, dg) + cache, u, mesh, equations, dg.mortar, dg.surface_integral, dg) # Calculate mortar fluxes @timed timer() "mortar flux" calc_mortar_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg.mortar, dg, cache) + dg.mortar, dg.surface_integral, dg, cache) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element @timed timer() "Jacobian" apply_jacobian!( @@ -565,7 +565,7 @@ Calculate the finite volume fluxes inside the elements (**with non-conservative end function prolong2interfaces!(cache, u, - mesh::TreeMesh{2}, equations, dg::DG) + mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @unpack interfaces = cache @unpack orientations = interfaces @@ -594,8 +594,8 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, nonconservative_terms::Val{false}, equations, - dg::DG, cache) - @unpack surface_flux = dg + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @threaded for interface in eachinterface(dg, cache) @@ -628,7 +628,7 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, nonconservative_terms::Val{true}, equations, - dg::DG, cache) + surface_integral, dg::DG, cache) @unpack u, neighbor_ids, orientations = cache.interfaces fstar_threaded = cache.fstar_upper_threaded noncons_diamond_primary_threaded = cache.noncons_diamond_upper_threaded @@ -641,7 +641,7 @@ function calc_interface_flux!(surface_flux_values, noncons_diamond_secondary = noncons_diamond_secondary_threaded[Threads.threadid()] # Calculate flux - calc_fstar!(fstar, equations, dg, u, interface, orientations[interface]) + calc_fstar!(fstar, equations, surface_integral, dg, u, interface, orientations[interface]) # Compute the nonconservative numerical "flux" along an interface # Done twice because left/right orientation matters så @@ -683,7 +683,7 @@ end function prolong2boundaries!(cache, u, - mesh::TreeMesh{2}, equations, dg::DG) + mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @unpack boundaries = cache @unpack orientations, neighbor_sides = boundaries @@ -723,13 +723,13 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::TreeMesh{2}, equations, dg::DG) + mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @assert isempty(eachboundary(dg, cache)) end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition, - mesh::TreeMesh{2}, equations, dg::DG) + mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -740,13 +740,13 @@ function calc_boundary_flux!(cache, t, boundary_condition, # Calc boundary fluxes in each direction for direction in eachindex(firsts) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_condition, - equations, dg, cache, + equations, surface_integral, dg, cache, direction, firsts[direction], lasts[direction]) end end function calc_boundary_flux!(cache, t, boundary_conditions::Union{NamedTuple,Tuple}, - mesh::TreeMesh{2}, equations, dg::DG) + mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -756,19 +756,24 @@ function calc_boundary_flux!(cache, t, boundary_conditions::Union{NamedTuple,Tup # Calc boundary fluxes in each direction calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1], - equations, dg, cache, 1, firsts[1], lasts[1]) + equations, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2], - equations, dg, cache, 2, firsts[2], lasts[2]) + equations, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[3], - equations, dg, cache, 3, firsts[3], lasts[3]) + equations, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[4], - equations, dg, cache, 4, firsts[4], lasts[4]) + equations, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) end function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t, - boundary_condition, equations, dg::DG, cache, + boundary_condition, equations, + surface_integral ,dg::DG, cache, direction, first_boundary, last_boundary) - @unpack surface_flux = dg + @unpack surface_flux = surface_integral @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries @threaded for boundary in first_boundary:last_boundary @@ -800,7 +805,7 @@ end function prolong2mortars!(cache, u, mesh::TreeMesh{2}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) + mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DGSEM) @threaded for mortar in eachmortar(dg, cache) @@ -887,8 +892,9 @@ end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{2}, nonconservative_terms::Val{false}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DG, cache) - @unpack neighbor_ids, u_lower, u_upper, orientations = cache.mortars + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack u_lower, u_upper, orientations = cache.mortars @unpack fstar_upper_threaded, fstar_lower_threaded = cache @threaded for mortar in eachmortar(dg, cache) @@ -898,8 +904,8 @@ function calc_mortar_flux!(surface_flux_values, # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper, equations, dg, u_upper, mortar, orientation) - calc_fstar!(fstar_lower, equations, dg, u_lower, mortar, orientation) + calc_fstar!(fstar_upper, equations, surface_integral, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_lower, equations, surface_integral, dg, u_lower, mortar, orientation) mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, @@ -912,8 +918,9 @@ end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{2}, nonconservative_terms::Val{true}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DG, cache) - @unpack neighbor_ids, u_lower, u_upper, orientations, large_sides = cache.mortars + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack u_lower, u_upper, orientations, large_sides = cache.mortars @unpack fstar_upper_threaded, fstar_lower_threaded, noncons_diamond_upper_threaded, noncons_diamond_lower_threaded = cache @@ -927,8 +934,8 @@ function calc_mortar_flux!(surface_flux_values, # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper, equations, dg, u_upper, mortar, orientation) - calc_fstar!(fstar_lower, equations, dg, u_lower, mortar, orientation) + calc_fstar!(fstar_upper, equations, surface_integral, dg, u_upper, mortar, orientation) + calc_fstar!(fstar_lower, equations, surface_integral, dg, u_lower, mortar, orientation) # Compute the nonconservative numerical terms along the upper and lower interface # Done twice because left/right orientation matters @@ -970,9 +977,10 @@ function calc_mortar_flux!(surface_flux_values, return nothing end -@inline function calc_fstar!(destination::AbstractArray{<:Any,2}, equations, dg::DGSEM, +@inline function calc_fstar!(destination::AbstractArray{<:Any,2}, equations, + surface_integral, dg::DGSEM, u_interfaces, interface, orientation) - @unpack surface_flux = dg + @unpack surface_flux = surface_integral for i in eachnode(dg) # Call pointwise two-point numerical flux function @@ -1060,8 +1068,9 @@ end end -function calc_surface_integral!(du, mesh::Union{TreeMesh{2}, CurvedMesh{2}}, - equations, dg::DGSEM, cache) +function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, CurvedMesh{2}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DG, cache) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements diff --git a/src/solvers/dg_tree/dg_2d_parallel.jl b/src/solvers/dg_tree/dg_2d_parallel.jl index a43d54af51..089fd8433f 100644 --- a/src/solvers/dg_tree/dg_2d_parallel.jl +++ b/src/solvers/dg_tree/dg_2d_parallel.jl @@ -269,7 +269,7 @@ function rhs!(du, u, t, # Prolong solution to MPI interfaces @timed timer() "prolong2mpiinterfaces" prolong2mpiinterfaces!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Start to send MPI data @timed timer() "start MPI send" start_mpi_send!( @@ -287,31 +287,31 @@ function rhs!(du, u, t, # Prolong solution to interfaces # TODO: Taal decide order of arguments, consistent vs. modified cache first? @timed timer() "prolong2interfaces" prolong2interfaces!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate interface fluxes @timed timer() "interface flux" calc_interface_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg, cache) + dg.surface_integral, dg, cache) # Prolong solution to boundaries @timed timer() "prolong2boundaries" prolong2boundaries!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, mesh, equations, dg) + cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Prolong solution to mortars @timed timer() "prolong2mortars" prolong2mortars!( - cache, u, mesh, equations, dg.mortar, dg) + cache, u, mesh, equations, dg.mortar, dg.surface_integral, dg) # Calculate mortar fluxes @timed timer() "mortar flux" calc_mortar_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg.mortar, dg, cache) + dg.mortar, dg.surface_integral, dg, cache) # Finish to receive MPI data @timed timer() "finish MPI receive" finish_mpi_receive!( @@ -321,11 +321,11 @@ function rhs!(du, u, t, @timed timer() "MPI interface flux" calc_mpi_interface_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg, cache) + dg.surface_integral, dg, cache) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element @timed timer() "Jacobian" apply_jacobian!( @@ -344,7 +344,7 @@ end function prolong2mpiinterfaces!(cache, u, mesh::ParallelTreeMesh{2}, - equations, dg::DG) + equations, surface_integral, dg::DG) @unpack mpi_interfaces = cache @threaded for interface in eachmpiinterface(dg, cache) @@ -380,8 +380,8 @@ end function calc_mpi_interface_flux!(surface_flux_values, mesh::ParallelTreeMesh{2}, nonconservative_terms::Val{false}, equations, - dg::DG, cache) - @unpack surface_flux = dg + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral @unpack u, local_element_ids, orientations, remote_sides = cache.mpi_interfaces @threaded for interface in eachmpiinterface(dg, cache) diff --git a/src/solvers/dg_tree/dg_3d.jl b/src/solvers/dg_tree/dg_3d.jl index 1b401376d7..4640fd4cc9 100644 --- a/src/solvers/dg_tree/dg_3d.jl +++ b/src/solvers/dg_tree/dg_3d.jl @@ -157,35 +157,35 @@ function rhs!(du, u, t, # Prolong solution to interfaces @timed timer() "prolong2interfaces" prolong2interfaces!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate interface fluxes @timed timer() "interface flux" calc_interface_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg, cache) + dg.surface_integral, dg, cache) # Prolong solution to boundaries @timed timer() "prolong2boundaries" prolong2boundaries!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, mesh, equations, dg) + cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Prolong solution to mortars @timed timer() "prolong2mortars" prolong2mortars!( - cache, u, mesh, equations, dg.mortar, dg) + cache, u, mesh, equations, dg.mortar, dg.surface_integral, dg) # Calculate mortar fluxes @timed timer() "mortar flux" calc_mortar_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg.mortar, dg, cache) + dg.mortar, dg.surface_integral, dg, cache) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element @timed timer() "Jacobian" apply_jacobian!( @@ -514,7 +514,7 @@ end function prolong2interfaces!(cache, u, - mesh::TreeMesh{3}, equations, dg::DG) + mesh::TreeMesh{3}, equations, surface_integral, dg::DG) @unpack interfaces = cache @unpack orientations = interfaces @@ -549,8 +549,8 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{3}, nonconservative_terms::Val{false}, equations, - dg::DG, cache) - @unpack surface_flux = dg + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral @unpack u, neighbor_ids, orientations = cache.interfaces @threaded for interface in eachinterface(dg, cache) @@ -582,7 +582,7 @@ end function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{3}, nonconservative_terms::Val{true}, equations, - dg::DG, cache) + surface_integral, dg::DG, cache) @unpack u, neighbor_ids, orientations = cache.interfaces fstar_threaded = cache.fstar_upper_left_threaded noncons_diamond_primary_threaded = cache.noncons_diamond_upper_left_threaded @@ -595,7 +595,7 @@ function calc_interface_flux!(surface_flux_values, noncons_diamond_secondary = noncons_diamond_secondary_threaded[Threads.threadid()] # Calculate flux - calc_fstar!(fstar, equations, dg, u, interface, orientations[interface]) + calc_fstar!(fstar, equations, surface_integral, dg, u, interface, orientations[interface]) # Compute the nonconservative numerical "flux" along an interface # Done twice because left/right orientation matters så @@ -638,7 +638,7 @@ end function prolong2boundaries!(cache, u, - mesh::TreeMesh{3}, equations, dg::DG) + mesh::TreeMesh{3}, equations, surface_integral, dg::DG) @unpack boundaries = cache @unpack orientations, neighbor_sides = boundaries @@ -691,13 +691,13 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::TreeMesh{3}, equations, dg::DG) + mesh::TreeMesh{3}, equations, surface_integral, dg::DG) @assert isempty(eachboundary(dg, cache)) end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition, - mesh::TreeMesh{3}, equations, dg::DG) + mesh::TreeMesh{3}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -708,13 +708,13 @@ function calc_boundary_flux!(cache, t, boundary_condition, # Calc boundary fluxes in each direction for direction in eachindex(firsts) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_condition, - equations, dg, cache, + equations, surface_integral, dg, cache, direction, firsts[direction], lasts[direction]) end end function calc_boundary_flux!(cache, t, boundary_conditions::Union{NamedTuple,Tuple}, - mesh::TreeMesh{3}, equations, dg::DG) + mesh::TreeMesh{3}, equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack n_boundaries_per_direction = cache.boundaries @@ -724,23 +724,30 @@ function calc_boundary_flux!(cache, t, boundary_conditions::Union{NamedTuple,Tup # Calc boundary fluxes in each direction calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1], - equations, dg, cache, 1, firsts[1], lasts[1]) + equations, surface_integral, dg, cache, + 1, firsts[1], lasts[1]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2], - equations, dg, cache, 2, firsts[2], lasts[2]) + equations, surface_integral, dg, cache, + 2, firsts[2], lasts[2]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[3], - equations, dg, cache, 3, firsts[3], lasts[3]) + equations, surface_integral, dg, cache, + 3, firsts[3], lasts[3]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[4], - equations, dg, cache, 4, firsts[4], lasts[4]) + equations, surface_integral, dg, cache, + 4, firsts[4], lasts[4]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[5], - equations, dg, cache, 5, firsts[5], lasts[5]) + equations, surface_integral, dg, cache, + 5, firsts[5], lasts[5]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[6], - equations, dg, cache, 6, firsts[6], lasts[6]) + equations, surface_integral, dg, cache, + 6, firsts[6], lasts[6]) end function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,5}, t, - boundary_condition, equations, dg::DG, cache, + boundary_condition, equations, + surface_integral, dg::DG, cache, direction, first_boundary, last_boundary) - @unpack surface_flux = dg + @unpack surface_flux = surface_integral @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries @threaded for boundary in first_boundary:last_boundary @@ -772,7 +779,8 @@ end function prolong2mortars!(cache, u, mesh::TreeMesh{3}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DGSEM) # temporary buffer for projections @unpack fstar_tmp1_threaded = cache @@ -902,9 +910,9 @@ end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{3}, nonconservative_terms::Val{false}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DG, cache) - @unpack (u_lower_left, u_lower_right, u_upper_left, u_upper_right, - neighbor_ids, orientations) = cache.mortars + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack u_lower_left, u_lower_right, u_upper_left, u_upper_right, orientations = cache.mortars @unpack (fstar_upper_left_threaded, fstar_upper_right_threaded, fstar_lower_left_threaded, fstar_lower_right_threaded, fstar_tmp1_threaded) = cache @@ -919,10 +927,10 @@ function calc_mortar_flux!(surface_flux_values, # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper_left, equations, dg, u_upper_left, mortar, orientation) - calc_fstar!(fstar_upper_right, equations, dg, u_upper_right, mortar, orientation) - calc_fstar!(fstar_lower_left, equations, dg, u_lower_left, mortar, orientation) - calc_fstar!(fstar_lower_right, equations, dg, u_lower_right, mortar, orientation) + calc_fstar!(fstar_upper_left, equations, surface_integral, dg, u_upper_left, mortar, orientation) + calc_fstar!(fstar_upper_right, equations, surface_integral, dg, u_upper_right, mortar, orientation) + calc_fstar!(fstar_lower_left, equations, surface_integral, dg, u_lower_left, mortar, orientation) + calc_fstar!(fstar_lower_right, equations, surface_integral, dg, u_lower_right, mortar, orientation) mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, mortar, @@ -937,9 +945,9 @@ end function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{3}, nonconservative_terms::Val{true}, equations, - mortar_l2::LobattoLegendreMortarL2, dg::DG, cache) - @unpack (u_lower_left, u_lower_right, u_upper_left, u_upper_right, - neighbor_ids, orientations) = cache.mortars + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack u_lower_left, u_lower_right, u_upper_left, u_upper_right, orientations = cache.mortars @unpack (fstar_upper_left_threaded, fstar_upper_right_threaded, fstar_lower_left_threaded, fstar_lower_right_threaded, noncons_diamond_upper_left_threaded, noncons_diamond_upper_right_threaded, @@ -961,10 +969,10 @@ function calc_mortar_flux!(surface_flux_values, # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper_left, equations, dg, u_upper_left, mortar, orientation) - calc_fstar!(fstar_upper_right, equations, dg, u_upper_right, mortar, orientation) - calc_fstar!(fstar_lower_left, equations, dg, u_lower_left, mortar, orientation) - calc_fstar!(fstar_lower_right, equations, dg, u_lower_right, mortar, orientation) + calc_fstar!(fstar_upper_left, equations, surface_integral, dg, u_upper_left, mortar, orientation) + calc_fstar!(fstar_upper_right, equations, surface_integral, dg, u_upper_right, mortar, orientation) + calc_fstar!(fstar_lower_left, equations, surface_integral, dg, u_lower_left, mortar, orientation) + calc_fstar!(fstar_lower_right, equations, surface_integral, dg, u_lower_right, mortar, orientation) # Compute the nonconservative numerical terms along the upper/lower interfaces # Done twice because left/right orientation matters @@ -1044,9 +1052,10 @@ function calc_mortar_flux!(surface_flux_values, return nothing end -@inline function calc_fstar!(destination::AbstractArray{<:Any,3}, equations, dg::DGSEM, +@inline function calc_fstar!(destination::AbstractArray{<:Any,3}, equations, + surface_integral, dg::DGSEM, u_interfaces, interface, orientation) - @unpack surface_flux = dg + @unpack surface_flux = surface_integral for j in eachnode(dg), i in eachnode(dg) # Call pointwise two-point numerical flux function @@ -1145,8 +1154,8 @@ end end -function calc_surface_integral!(du, mesh::Union{TreeMesh{3}, CurvedMesh{3}}, - equations, dg::DGSEM, cache) +function calc_surface_integral!(du, u, mesh::Union{TreeMesh{3}, CurvedMesh{3}}, + equations, surface_integral, dg::DGSEM, cache) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements diff --git a/src/solvers/dg_unstructured_quad/dg_2d.jl b/src/solvers/dg_unstructured_quad/dg_2d.jl index 4f832e4f81..625c4b210b 100644 --- a/src/solvers/dg_unstructured_quad/dg_2d.jl +++ b/src/solvers/dg_unstructured_quad/dg_2d.jl @@ -7,7 +7,7 @@ function create_cache(mesh::UnstructuredQuadMesh, equations, polydeg_ = polydeg(dg.basis) nvars = nvariables(equations) - elements = init_elements(RealT, uEltype, mesh, dg.basis.nodes, nvars, polydeg_) + elements = init_elements(RealT, uEltype, mesh, get_nodes(dg.basis), nvars, polydeg_) interfaces = init_interfaces(uEltype, mesh, nvars, polydeg_) @@ -40,25 +40,25 @@ function rhs!(du, u, t, # Prolong solution to interfaces @timed timer() "prolong2interfaces" prolong2interfaces!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate interface fluxes @timed timer() "interface flux" calc_interface_flux!( cache.elements.surface_flux_values, mesh, have_nonconservative_terms(equations), equations, - dg, cache) + dg.surface_integral, dg, cache) # Prolong solution to boundaries @timed timer() "prolong2boundaries" prolong2boundaries!( - cache, u, mesh, equations, dg) + cache, u, mesh, equations, dg.surface_integral, dg) # Calculate boundary fluxes @timed timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, mesh, equations, dg) + cache, t, boundary_conditions, mesh, equations, dg.surface_integral, dg) # Calculate surface integrals @timed timer() "surface integral" calc_surface_integral!( - du, mesh, equations, dg, cache) + du, u, mesh, equations, dg.surface_integral, dg, cache) # Apply Jacobian from mapping to reference element # Note! this routine is reused from dg_curved/dg_2d.jl @@ -159,7 +159,7 @@ end # Note! this routine is for quadrilateral elements with "right-handed" orientation function prolong2interfaces!(cache, u, mesh::UnstructuredQuadMesh, - equations, dg::DG) + equations, surface_integral, dg::DG) @unpack interfaces = cache @threaded for interface in eachinterface(dg, cache) @@ -214,8 +214,9 @@ end # quadrilateral mesh function calc_interface_flux!(surface_flux_values, mesh::UnstructuredQuadMesh, - nonconservative_terms::Val{false}, equations, dg::DG, cache) - @unpack surface_flux = dg + nonconservative_terms::Val{false}, equations, + surface_integral, dg::DG, cache) + @unpack surface_flux = surface_integral @unpack u, start_index, index_increment, element_ids, element_side_ids = cache.interfaces @unpack normal_directions = cache.elements @@ -266,7 +267,7 @@ end # move the approximate solution onto physical boundaries within a "right-handed" element function prolong2boundaries!(cache, u, mesh::UnstructuredQuadMesh, - equations, dg::DG) + equations, surface_integral, dg::DG) @unpack boundaries = cache @threaded for boundary in eachboundary(dg, cache) @@ -298,7 +299,8 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::Union{UnstructuredQuadMesh, P4estMesh}, equations, dg::DG) + mesh::Union{UnstructuredQuadMesh, P4estMesh}, + equations, surface_integral, dg::DG) @assert isempty(eachboundary(dg, cache)) end @@ -306,11 +308,11 @@ end # Function barrier for type stability function calc_boundary_flux!(cache, t, boundary_conditions, mesh::Union{UnstructuredQuadMesh, P4estMesh}, - equations, dg::DG) + equations, surface_integral, dg::DG) @unpack boundary_condition_types, boundary_indices = boundary_conditions calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, - mesh, equations, dg) + mesh, equations, surface_integral, dg) return nothing end @@ -320,7 +322,7 @@ end function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N,Any}, BC_indices::NTuple{N,Vector{Int}}, mesh::Union{UnstructuredQuadMesh, P4estMesh}, - equations, dg::DG) where {N} + equations, surface_integral, dg::DG) where {N} # Extract the boundary condition type and index vector boundary_condition = first(BCs) boundary_condition_indices = first(BC_indices) @@ -329,11 +331,13 @@ function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N,Any}, remaining_boundary_condition_indices = Base.tail(BC_indices) # process the first boundary condition type - calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, mesh, equations, dg) + calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, + mesh, equations, surface_integral, dg) # recursively call this method with the unprocessed boundary types calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions, - remaining_boundary_condition_indices, mesh, equations, dg) + remaining_boundary_condition_indices, + mesh, equations, surface_integral, dg) return nothing end @@ -341,13 +345,14 @@ end # terminate the type-stable iteration over tuples function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, mesh::Union{UnstructuredQuadMesh, P4estMesh}, - equations, dg::DG) - return nothing + equations, surface_integral, dg::DG) + nothing end function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, - mesh::UnstructuredQuadMesh, equations, dg::DG) + mesh::UnstructuredQuadMesh, equations, + surface_integral, dg::DG) @unpack surface_flux_values = cache.elements @unpack element_id, element_side_id = cache.boundaries @@ -361,7 +366,8 @@ function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, # calc boundary flux on the current boundary interface for node in eachnode(dg) - calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh, equations, dg, cache, + calc_boundary_flux!(surface_flux_values, t, boundary_condition, + mesh, equations, surface_integral, dg, cache, node, side, element, boundary) end end @@ -371,11 +377,12 @@ end # inlined version of the boundary flux calculation along a physical interface where the # boundary flux values are set according to a particular `boundary_condition` function @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::UnstructuredQuadMesh, equations, dg::DG, cache, + mesh::UnstructuredQuadMesh, equations, + surface_integral, dg::DG, cache, node_index, side_index, element_index, boundary_index) @unpack normal_directions = cache.elements @unpack u, node_coordinates = cache.boundaries - @unpack surface_flux = dg + @unpack surface_flux = surface_integral # pull the inner solution state from the boundary u values on the boundary element u_inner = get_node_vars(u, equations, dg, node_index, boundary_index) @@ -410,8 +417,8 @@ end # ----------------- ----------------- # 3 1 # Therefore, we require a different surface integral routine here despite their similar structure. -function calc_surface_integral!(du, mesh::UnstructuredQuadMesh, - equations, dg::DGSEM, cache) +function calc_surface_integral!(du, u, mesh::UnstructuredQuadMesh, + equations, surface_integral, dg::DGSEM, cache) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements @@ -448,7 +455,7 @@ function max_discrete_metric_identities(dg::DGSEM, cache) metric_id_dx = zeros(eltype(contravariant_vectors), nnodes(dg), nnodes(dg)) metric_id_dy = zeros(eltype(contravariant_vectors), nnodes(dg), nnodes(dg)) - max_metric_ids = zero(dg.basis.nodes[1]) + max_metric_ids = zero(eltype(contravariant_vectors)) for i in 1:ndims_, element in eachelement(dg, cache) # compute D*Ja_1^i + Ja_2^i*D^T diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl new file mode 100644 index 0000000000..2f4f9da18a --- /dev/null +++ b/src/solvers/fdsbp_tree/fdsbp_2d.jl @@ -0,0 +1,207 @@ + +# General interface methods for SummationByPartsOperators.jl and Trixi.jl +# TODO: FD. Move to another file +nnodes(D::AbstractDerivativeOperator) = size(D, 1) +eachnode(D::AbstractDerivativeOperator) = Base.OneTo(nnodes(D)) +get_nodes(D::AbstractDerivativeOperator) = grid(D) + + +# For dispatch +# TODO: FD. Move to another file +const FDSBP = DG{Basis} where {Basis<:AbstractDerivativeOperator} + + +# 2D containers +# TODO: FD. Move to another file +init_mortars(cell_ids, mesh, elements, mortar::Nothing) = nothing + + +# 2D caches +function create_cache(mesh::TreeMesh{2}, equations, + volume_integral::VolumeIntegralStrongForm, dg, uEltype) + + prototype = Array{SVector{nvariables(equations), uEltype}, ndims(mesh)}( + undef, ntuple(_ -> nnodes(dg), ndims(mesh))...) + f_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] + + return (; f_threaded,) +end + +create_cache(mesh, equations, mortar::Nothing, uEltype) = NamedTuple() +nmortars(mortar::Nothing) = 0 + + +# 2D RHS +function calc_volume_integral!(du, u, + mesh::TreeMesh{2}, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralStrongForm, + dg::FDSBP, cache) + D = dg.basis # SBP derivative operator + @unpack f_threaded = cache + + # SBP operators from SummationByPartsOperators.jl implement the basic interface + # of matrix-vector multiplication. Thus, we pass an "array of structures", + # packing all variables per node in an `SVector`. + if nvariables(equations) == 1 + # `reinterpret(reshape, ...)` removes the leading dimension only if more + # than one variable is used. + u_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, du), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + else + u_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u) + du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, du) + end + + # Use the tensor product structure to compute the discrete derivatives of + # the fluxes line-by-line and add them to `du` for each element. + @threaded for element in eachelement(dg, cache) + f_element = f_threaded[Threads.threadid()] + u_element = view(u_vectors, :, :, element) + + # x direction + @. f_element = flux(u_element, 1, equations) + for j in eachnode(dg) + mul!(view(du_vectors, :, j, element), D, view(f_element, :, j), + one(eltype(du)), one(eltype(du))) + end + + # y direction + @. f_element = flux(u_element, 2, equations) + for i in eachnode(dg) + mul!(view(du_vectors, i, :, element), D, view(f_element, i, :), + one(eltype(du)), one(eltype(du))) + end + end + + return nothing +end + + +function prolong2mortars!(cache, u, mesh, equations, mortar::Nothing, + surface_integral, dg::DG) + @assert isempty(eachmortar(dg, cache)) +end + +function calc_mortar_flux!(surface_flux_values, mesh, + nonconservative_terms, equations, + mortar::Nothing, + surface_integral, dg::DG, cache) + @assert isempty(eachmortar(dg, cache)) +end + + +function calc_surface_integral!(du, u, mesh::TreeMesh{2}, + equations, surface_integral::SurfaceIntegralStrongForm, + dg::DG, cache) + inv_weight_left = inv(left_boundary_weight(dg.basis)) + inv_weight_right = inv(right_boundary_weight(dg.basis)) + @unpack surface_flux_values = cache.elements + + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + # surface at -x + u_node = get_node_vars(u, equations, dg, 1, l, element) + f_node = flux(u_node, 1, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) + surface_contribution = -(f_num - f_node) * inv_weight_left + add_to_node_vars!(du, surface_contribution, equations, dg, 1, l, element) + + # surface at +x + u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) + f_node = flux(u_node, 1, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) + surface_contribution = +(f_num - f_node) * inv_weight_right + add_to_node_vars!(du, surface_contribution, equations, dg, nnodes(dg), l, element) + + # surface at -y + u_node = get_node_vars(u, equations, dg, l, 1, element) + f_node = flux(u_node, 2, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) + surface_contribution = -(f_num - f_node) * inv_weight_left + add_to_node_vars!(du, surface_contribution, equations, dg, l, 1, element) + + # surface at +y + u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) + f_node = flux(u_node, 2, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) + surface_contribution = +(f_num - f_node) * inv_weight_right + add_to_node_vars!(du, surface_contribution, equations, dg, l, nnodes(dg), element) + end + end + + return nothing +end + + +# AnalysisCallback +# TODO: FD. Move to another file +SolutionAnalyzer(D::AbstractDerivativeOperator) = D + +function integrate_via_indices(func::Func, u, + mesh::TreeMesh{2}, equations, + dg::FDSBP, cache, args...; normalize=true) where {Func} + # TODO: FD. This is rather inefficient right now and allocates... + weights = diag(SummationByPartsOperators.mass_matrix(dg.basis)) + + # Initialize integral with zeros of the right shape + integral = zero(func(u, 1, 1, 1, equations, dg, args...)) + + # Use quadrature to numerically integrate over entire domain + for element in eachelement(dg, cache) + volume_jacobian_ = volume_jacobian(element, mesh, cache) + for j in eachnode(dg), i in eachnode(dg) + integral += volume_jacobian_ * weights[i] * weights[j] * func(u, i, j, element, equations, dg, args...) + end + end + + # Normalize with total volume + if normalize + integral = integral / total_volume(mesh) + end + + return integral +end + +function calc_error_norms(func, u, t, analyzer, + mesh::TreeMesh{2}, equations, initial_condition, + dg::FDSBP, cache, cache_analysis) + # TODO: FD. This is rather inefficient right now and allocates... + weights = diag(SummationByPartsOperators.mass_matrix(dg.basis)) + @unpack node_coordinates = cache.elements + + # Set up data structures + l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations)) + linf_error = copy(l2_error) + + # Iterate over all elements for error calculations + for element in eachelement(dg, cache) + # Calculate errors at each node + volume_jacobian_ = volume_jacobian(element, mesh, cache) + + for j in eachnode(analyzer), i in eachnode(analyzer) + u_exact = initial_condition( + get_node_coords(node_coordinates, equations, dg, i, j, element), t, equations) + diff = func(u_exact, equations) - func( + get_node_vars(u, equations, dg, i, j, element), equations) + l2_error += diff.^2 * (weights[i] * weights[j] * volume_jacobian_) + linf_error = @. max(linf_error, abs(diff)) + end + end + + # For L2 error, divide by total volume + total_volume_ = total_volume(mesh) + l2_error = @. sqrt(l2_error / total_volume_) + + return l2_error, linf_error +end + + +# TODO: FD. Visualization +# We need to define better interfaces for all the plotting stuff. +# Right now, the easiest solution is to use scatter plots such as +# x = semi.cache.elements.node_coordinates[1, :, :, :] |> vec +# y = semi.cache.elements.node_coordinates[2, :, :, :] |> vec +# scatter(x, y, sol.u[end]) diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index 3b33822f7f..3d32435716 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -3,3 +3,4 @@ include("dg_curved/dg.jl") include("dg_unstructured_quad/dg.jl") include("dg_p4est/dg.jl") include("dg_common.jl") +include("fdsbp_tree/fdsbp_2d.jl") diff --git a/test/test_examples_2d_part3.jl b/test/test_examples_2d_part3.jl index 5d4bc8f1ae..70286dbdfe 100644 --- a/test/test_examples_2d_part3.jl +++ b/test/test_examples_2d_part3.jl @@ -27,6 +27,9 @@ isdir(outdir) && rm(outdir, recursive=true) # P4estMesh include("test_examples_2d_p4est.jl") + + # FDSBP methods on the TreeMesh + include("test_tree_2d_fdsbp.jl") end # Clean up afterwards: delete Trixi output directory diff --git a/test/test_tree_2d_fdsbp.jl b/test/test_tree_2d_fdsbp.jl new file mode 100644 index 0000000000..9a135dce68 --- /dev/null +++ b/test/test_tree_2d_fdsbp.jl @@ -0,0 +1,20 @@ +module TestTree2DFDSBP + +using Test +using Trixi + +include("test_trixi.jl") + +# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory +EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_fdsbp") + +@testset "Linear scalar advection" begin + @testset "elixir_advection_extended.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"), + l2 = [2.29497319e-06], + linf = [8.48027655e-06], + rtol = 5.0e-7) # These results change a little bit and depend on the CI system + end +end + +end # module diff --git a/test/test_trixi.jl b/test/test_trixi.jl index 1d40e4f217..8a4afa0631 100644 --- a/test/test_trixi.jl +++ b/test/test_trixi.jl @@ -116,7 +116,11 @@ macro test_nowarn_mod(expr) # We also ignore simple module redefinitions for convenience. Thus, we # check whether every line of `stderr_content` is of the form of a # module replacement warning. - @test occursin(r"^(WARNING: replacing module .+\.\n)*$", stderr_content) + # TODO: Upstream (PlotUtils). This should be removed again once the + # deprecated stuff is fixed upstream. + if stderr_content != "WARNING: importing deprecated binding Colors.RGB1 into PlotUtils.\nWARNING: importing deprecated binding Colors.RGB4 into PlotUtils.\n" + @test occursin(r"^(WARNING: replacing module .+\.\n)*$", stderr_content) + end ret finally rm(fname, force=true) diff --git a/test/test_unit.jl b/test/test_unit.jl index 538dd2ff37..fbee8cb423 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -193,7 +193,7 @@ Cassette.@context Ctx @testset "calc_jacobian_matrix" begin @testset "identity map" begin basis = LobattoLegendreBasis(5) - nodes = basis.nodes + nodes = Trixi.get_nodes(basis) jacobian_matrix = Array{Float64, 5}(undef, 2, 2, 6, 6, 1) node_coordinates = Array{Float64, 4}(undef, 2, 6, 6, 1) @@ -207,7 +207,7 @@ Cassette.@context Ctx @testset "maximum exact polydeg" begin basis = LobattoLegendreBasis(3) - nodes = basis.nodes + nodes = Trixi.get_nodes(basis) jacobian_matrix = Array{Float64, 5}(undef, 2, 2, 4, 4, 1) # f(x, y) = [x^3, xy^2] @@ -386,9 +386,9 @@ Cassette.@context Ctx end @testset "DG L2 mortar container debug output" begin - c2d = Trixi.L2MortarContainer2D{Float64, 1, 1}(1) + c2d = Trixi.L2MortarContainer2D{Float64}(1, 1, 1) @test isnothing(display(c2d)) - c3d = Trixi.L2MortarContainer3D{Float64, 1, 1}(1) + c3d = Trixi.L2MortarContainer3D{Float64}(1, 1, 1) @test isnothing(display(c3d)) end