diff --git a/Project.toml b/Project.toml index dbcff694..a52e073e 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +StartUpDG = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -29,6 +30,7 @@ NLsolve = "4.5.1" Printf = "1" QuadGK = "2" Reexport = "1.0" +StartUpDG = "1.2.0" Static = "0.8, 1" StaticArrayInterface = "1.5.1" StaticArrays = "1" diff --git a/examples/elixir_spherical_advection_covariant_prismed_icosahedron.jl b/examples/elixir_spherical_advection_covariant_prismed_icosahedron.jl new file mode 100644 index 00000000..f71a231c --- /dev/null +++ b/examples/elixir_spherical_advection_covariant_prismed_icosahedron.jl @@ -0,0 +1,77 @@ +############################################################################### +# DGSEM for the linear advection equation on a prismed icosahedral grid +############################################################################### +# To run a convergence test, use +# convergence_test("../examples/elixir_spherical_advection_covariant_prismed_icosahedron.jl", 4, initial_refinement_level = 1) + +using OrdinaryDiffEq, Trixi, TrixiAtmo + +############################################################################### +# Spatial discretization + +initial_condition = initial_condition_gaussian + +equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = GlobalCartesianCoordinates()) + +############################################################################### +# Build DG solver. + +tensor_polydeg = (2, 1) + +dg = DGMulti(element_type = Wedge(), + approximation_type = Polynomial(), + surface_flux = flux_central, + polydeg = tensor_polydeg) + +############################################################################### +# Build mesh. + +initial_refinement_level = 3 + +mesh = DGMultiMeshPrismIcosahedron(dg; + inner_radius = 0.999 * EARTH_RADIUS, + outer_radius = EARTH_RADIUS, + initial_refinement = initial_refinement_level) + +# Transform the initial condition to the proper set of conservative variables +initial_condition_transformed = transform_initial_condition(initial_condition, equations) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, dg) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0 to T +ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY)) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation +# setup and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the +# results +analysis_callback = AnalysisCallback(semi, interval = 10, + save_analysis = true, + extra_analysis_errors = (:conservation_error,), + uEltype = real(dg)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = contravariant2global) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.7) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE +# solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed +# callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, save_everystep = true, callback = callbacks) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index ed316d03..d0f7a7bb 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -15,7 +15,7 @@ using Printf: @sprintf using Static: True, False using StrideArrays: PtrArray using StaticArrayInterface: static_size -using LinearAlgebra: cross, norm, dot, det +using LinearAlgebra: Diagonal, cross, norm, dot, det using Reexport: @reexport using LoopVectorization: @turbo using QuadGK: quadgk @@ -35,6 +35,9 @@ using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace using Trixi: ln_mean, stolarsky_mean, inv_ln_mean +# DGMulti solvers +using StartUpDG: RefElemData, MeshData, AbstractElemShape + include("auxiliary/auxiliary.jl") include("equations/equations.jl") include("meshes/meshes.jl") @@ -69,9 +72,9 @@ export source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier! export cons2prim_and_vorticity, contravariant2global -export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct, - MetricTermsInvariantCurl, MetricTermsCovariantSphere, ChristoffelSymbolsAutodiff, - ChristoffelSymbolsCollocationDerivative +export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, DGMultiMeshPrismIcosahedron, + MetricTermsCrossProduct, MetricTermsInvariantCurl, MetricTermsCovariantSphere, + ChristoffelSymbolsAutodiff, ChristoffelSymbolsCollocationDerivative export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, SECONDS_PER_DAY diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index 60d64168..9b2b042a 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -77,6 +77,29 @@ function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t, end end +# Entropy time derivative for cons2entropy function which depends on auxiliary variables +function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t, + mesh::DGMultiMesh, equations::AbstractCovariantEquations, dg::DGMulti, cache) + rd = dg.basis + md = mesh.md + @unpack u_values, aux_quad_values = cache + + # interpolate u, du to quadrature points + du_values = similar(u_values) # TODO: DGMulti. Can we move this to the analysis cache somehow? + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), du_values, du) + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u) + + # compute ∫v(u) * du/dt = ∫dS/dt. We can directly compute v(u) instead of computing the entropy + # projection here, since the RHS will be projected to polynomials of degree N and testing with + # the L2 projection of v(u) would be equivalent to testing with v(u) due to the moment-preserving + # property of the L2 projection. + dS_dt = zero(eltype(first(du))) + for i in Base.OneTo(length(md.wJq)) + dS_dt += dot(cons2entropy(u_values[i], aux_quad_values[i], equations), du_values[i]) * md.wJq[i] + end + return dS_dt +end + # L2 and Linf error calculation for the covariant form function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, equations::AbstractCovariantEquations{2}, @@ -124,4 +147,27 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, return l2_error, linf_error end + +# L2 and Linf error calculation for the covariant form +function Trixi.calc_error_norms(func, u, t, analyzer, + mesh::DGMultiMesh{NDIMS}, equations::AbstractCovariantEquations, initial_condition, + dg::DGMulti{NDIMS}, cache, cache_analysis) where {NDIMS} + rd = dg.basis + md = mesh.md + @unpack u_values, aux_quad_values = cache + + # interpolate u to quadrature points + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u) + + component_l2_errors = zero(eltype(u_values)) + component_linf_errors = zero(eltype(u_values)) + for i in Trixi.each_quad_node_global(mesh, dg, cache) + u_exact = initial_condition(SVector(getindex.(md.xyzq, i)), t, aux_quad_values[i], equations) + error_at_node = func(u_values[i], equations) - func(u_exact, equations) + component_l2_errors += md.wJq[i] * error_at_node .^ 2 + component_linf_errors = max.(component_linf_errors, abs.(error_at_node)) + end + total_volume = sum(md.wJq) + return sqrt.(component_l2_errors ./ total_volume), component_linf_errors +end end # @muladd diff --git a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl index c55514e9..9003fa5b 100644 --- a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl +++ b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl @@ -118,6 +118,40 @@ function Trixi.save_solution_file(u, time, dt, timestep, return filename end +function Trixi.save_solution_file(u, time, dt, timestep, + mesh::DGMultiMesh, + equations::AbstractCovariantEquations, + dg::DG, cache, + solution_callback, + element_variables = Dict{Symbol, Any}(), + node_variables = Dict{Symbol, Any}(); + system = "") + @unpack output_directory, solution_variables = solution_callback + + # Filename based on current time step. + if isempty(system) + filename = joinpath(output_directory, @sprintf("solution_%09d.vtu", timestep)) + else + filename = joinpath(output_directory, + @sprintf("solution_%s_%09d.vtu", system, timestep)) + end + + # Convert to different set of variables if requested. + if solution_variables === cons2cons + data = u + n_vars = nvariables(equations) + else + data = convert_variables(u, solution_variables, mesh, equations, dg, cache) + # Find out variable count by looking at output from `solution_variables` function. + n_vars = length(data[1]) + end + # Create a dictionary with variable names and data + var_names = Trixi.varnames(solution_variables, equations) + var_dict = Dict(var_names[i] => getindex.(data, i) for i in 1:n_vars) + StartUpDG.export_to_vtk(dg.basis, mesh.md, var_dict, filename; + equi_dist_nodes = false) +end + # Calculate the primitive variables and the relative vorticity at a given node @inline function cons2prim_and_vorticity(u, normal_vector, mesh::P4estMesh{2}, diff --git a/src/callbacks_step/save_solution_covariant.jl b/src/callbacks_step/save_solution_covariant.jl index db2e9c1b..f17b33d6 100644 --- a/src/callbacks_step/save_solution_covariant.jl +++ b/src/callbacks_step/save_solution_covariant.jl @@ -33,6 +33,27 @@ function convert_variables(u, solution_variables, mesh::P4estMesh{2}, return data end +# Convert to another set of variables using a solution_variables function +function convert_variables(u, solution_variables, mesh::DGMultiMesh, + equations::AbstractCovariantEquations, dg, cache) + (; aux_values) = cache + # Extract the number of solution variables to be output + # (may be different than the number of conservative variables) + n_vars = length(Trixi.varnames(solution_variables, equations)) + + # Allocate storage for output, an Np x n_elements array of nvars arrays + # data = Array{eltype(u)}(undef, n_vars, (size(aux_node_vars)[2:end])...) + data = map(_ -> SVector{n_vars, eltype(u)}(zeros(n_vars)), u) + + # Loop over all nodes and convert variables, passing in auxiliary variables + for i in Trixi.each_dof_global(mesh, dg) + u_node = u[:, i] + aux_node = aux_values[i] + data[i] = solution_variables(u_node, aux_node, equations) + end + return data +end + # Calculate the primitive variables and relative vorticity at a given node. The velocity # components in the global coordinate system and the bottom topography are returned, such # that the outputs for CovariantShallowWaterEquations2D and ShallowWaterEquations3D are the diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 0d9c888d..ddad1542 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -69,4 +69,31 @@ function Trixi.max_dt(u, t, mesh::P4estMesh{2}, constant_speed::False, end return 2 / (nnodes(dg) * max_scaled_speed) end + +function Trixi.max_dt(u, t, mesh::DGMultiMesh, + constant_speed::False, equations::AbstractCovariantEquations{NDIMS}, + dg::DGMulti, cache) where {NDIMS} + @unpack md = mesh + rd = dg.basis + (; aux_values) = cache + + dt_min = Inf + for e in eachelement(mesh, dg, cache) + max_speeds = ntuple(_ -> nextfloat(zero(t)), NDIMS) + for i in 1:nnodes(dg) + u_node = u[i, e] + aux_node = aux_values[i, e] + detg = area_element(aux_node, equations) + lambda_i = max_abs_speeds(u_node, aux_node, equations) + max_speeds = max.(max_speeds, lambda_i) + end + dt_min = min(dt_min, 1 / sum(max_speeds)) + end + + # This mimics `max_dt` for `TreeMesh`, except that `nnodes(dg)` is replaced by + # `polydeg+1`. This is because `nnodes(dg)` returns the total number of + # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns + # the number of 1D nodes for `DGSEM` solvers. + return 2 * dt_min * Trixi.dt_polydeg_scaling(dg) +end end # @muladd diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index 84e998b5..225fa6ab 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -96,6 +96,15 @@ end return SVector(J * u[1] * vcon[orientation], z, z) end +@inline function Trixi.flux(u, aux_vars, normal_direction::AbstractVector, + equations::CovariantLinearAdvectionEquation2D) + z = zero(eltype(u)) + J = area_element(aux_vars, equations) + vcon = velocity_contravariant(u, equations) + a = dot(vcon, normal_direction) # velocity in normal direction + return SVector(J * u[1] * a, z, z) +end + # Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity # components, as they should remain unchanged in time @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, aux_vars_ll, @@ -118,6 +127,17 @@ end return max(abs(vcon_ll[orientation]), abs(vcon_rr[orientation])) end +@inline function Trixi.max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + normal_direction::AbstractVector, + equations::CovariantLinearAdvectionEquation2D) + vcon_ll = velocity_contravariant(u_ll, equations) # Contravariant components on left side + vcon_rr = velocity_contravariant(u_rr, equations) # Contravariant components on right side + # Calculate the velocity in the normal direction + a_ll = abs(dot(vcon_ll, normal_direction)) + a_rr = abs(dot(vcon_rr, normal_direction)) + return max(a_ll, a_rr) +end + # Maximum wave speeds in each direction for CFL calculation @inline function max_abs_speeds(u, aux_vars, equations::CovariantLinearAdvectionEquation2D) diff --git a/src/meshes/dgmulti_prismed_sphere.jl b/src/meshes/dgmulti_prismed_sphere.jl new file mode 100644 index 00000000..13c9e609 --- /dev/null +++ b/src/meshes/dgmulti_prismed_sphere.jl @@ -0,0 +1,115 @@ +# Construct a DGMultiMesh consisting of prismatic elements based on a refined icosahedron +# As of now this 3D mesh only serves to simulate 2D covariant equations on the sphere, with +# the spherical shell thus only consisting of a single layer of prismatic elements extruded +# from the triangular faces of the icosahedron. +function DGMultiMeshPrismIcosahedron(dg::DGMulti{NDIMS}; + inner_radius = 1.0 * EARTH_RADIUS, + outer_radius = 0.99 * EARTH_RADIUS, + initial_refinement = 3, + is_on_boundary = nothing) where {NDIMS} + + vertex_coordinates = calc_node_coordinates_icosahedron_vertices(outer_radius) + + Vxyz = ntuple(n -> vertex_coordinates[n, :], NDIMS) + + EToV = zeros(Int, 20, 3) + + for i in 1:size(EToV, 1) + EToV[i, :] = icosahedron_triangle_vertices_idx_map[i] + end + + for j in 1:initial_refinement + EToV_old = EToV + Vxyz_old = ntuple(n -> copy(Vxyz[n]), NDIMS) + old_to_new = Dict{Int, Int}() + edge_to_new = Dict{Tuple{Int, Int}, Int}() + EToV = zeros(Int, (size(EToV_old, 1) * 4, 3)) + + Vxyz = ntuple(n -> Vector{eltype(Vxyz[1])}(), NDIMS) + + for i in 1:size(EToV_old, 1) + idx_old = EToV_old[i, :] + + for k in idx_old + if !haskey(old_to_new, k) + old_to_new[k] = length(Vxyz[1]) + 1 + for n in 1:NDIMS + push!(Vxyz[n], Vxyz_old[n][k]) + end + end + end + + for k in idx_old, l in idx_old + if k < l + edge = (k, l) + if !haskey(edge_to_new, edge) + edge_to_new[edge] = length(Vxyz[1]) + 1 + vk = ntuple(n -> Vxyz_old[n][k], NDIMS) + vl = ntuple(n -> Vxyz_old[n][l], NDIMS) + midpoint = 0.5 .* (vk .+ vl) + midpoint = midpoint .* (outer_radius / norm(midpoint)) # Normalize to outer radius + for n in 1:NDIMS + push!(Vxyz[n], midpoint[n]) + end + end + end + end + id1 = old_to_new[idx_old[1]] + id2 = edge_to_new[Tuple(sort([idx_old[1], idx_old[2]]))] + id3 = old_to_new[idx_old[2]] + id4 = edge_to_new[Tuple(sort([idx_old[2], idx_old[3]]))] + id5 = old_to_new[idx_old[3]] + id6 = edge_to_new[Tuple(sort([idx_old[3], idx_old[1]]))] + + ids = [id1, id2, id3, id4, id5, id6] + + # Fill EToV for the 4 new triangles + for (sub_i, vertex_map) in enumerate(icosahedron_tri_vertices_idx_map) + EToV[(i - 1) * 4 + sub_i, :] = getindex.(Ref(ids), vertex_map) + end + + end + end + + num_vertices = size(Vxyz[1], 1) + vertices_inner = ntuple(n -> Vxyz[n] .* (inner_radius / outer_radius), NDIMS) + + Vxyz = ntuple(n -> vcat(Vxyz[n], vertices_inner[n]), NDIMS) + EToV = hcat(EToV .+ num_vertices, EToV) + + md = MeshData(Vxyz, EToV, dg.basis) + boundary_faces = StartUpDG.tag_boundary_faces(md, is_on_boundary) + return DGMultiMesh(dg, Trixi.GeometricTermsType(Trixi.Curved(), dg), md, boundary_faces) +end + +# We use a local numbering to obtain the triangle vertices of each triangular face +# +# Fig: Local quad vertex numbering for a triangular face (corner vertices of the triangular face in parenthesis) +# 5 (3) +# /\ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# 6/ 4\ +# /⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺\ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# 1/_________________\_/_________________3\ +# (1) 2 (2) + +# Index map for the vertices of each triangle on the triangular faces of the icosahedron (see Fig 4) +const icosahedron_tri_vertices_idx_map = ([1, 2, 6], # Tri 1 + [2, 3, 4], # Tri 2 + [4, 5, 6], # Tri 3 + [2, 4, 6]) # Tri 4 \ No newline at end of file diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index 8bd9c4b0..ab8d9d71 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -1,2 +1,3 @@ include("p4est_cubed_sphere.jl") include("p4est_icosahedron_quads.jl") +include("dgmulti_prismed_sphere.jl") \ No newline at end of file diff --git a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl index 885e5a3b..fd7e39f4 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl @@ -40,6 +40,38 @@ function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{2}, performance_counter) end +function Trixi.SemidiscretizationHyperbolic(mesh::DGMultiMesh, + equations::AbstractCovariantEquations, + initial_condition, + solver; + source_terms = nothing, + boundary_conditions = boundary_condition_periodic, + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT = real(solver), uEltype = RealT, + initial_cache = NamedTuple(), + metric_terms = MetricTermsCrossProduct(), + auxiliary_field = nothing) + cache = (; + Trixi.create_cache(mesh, equations, solver, RealT, metric_terms, + auxiliary_field, uEltype)..., initial_cache...) + _boundary_conditions = Trixi.digest_boundary_conditions(boundary_conditions, mesh, + solver, + cache) + + performance_counter = Trixi.PerformanceCounter() + + SemidiscretizationHyperbolic{typeof(mesh), typeof(equations), + typeof(initial_condition), + typeof(_boundary_conditions), typeof(source_terms), + typeof(solver), typeof(cache)}(mesh, equations, + initial_condition, + _boundary_conditions, + source_terms, solver, + cache, + performance_counter) +end + # Constructor for SemidiscretizationHyperbolic for the covariant form. Requires # compatibility between the mesh and equations (i.e. the same `NDIMS` and `NDIMS_AMBIENT`) # and sets the default metric terms to MetricTermsCovariantSphere. diff --git a/src/solvers/dgmulti/containers_manifold_covariant.jl b/src/solvers/dgmulti/containers_manifold_covariant.jl new file mode 100644 index 00000000..8be0c580 --- /dev/null +++ b/src/solvers/dgmulti/containers_manifold_covariant.jl @@ -0,0 +1,99 @@ +@muladd begin +#! format: noindnent + +function init_auxiliary_node_variables!(aux_values, mesh::DGMultiMesh, + equations::AbstractCovariantEquations{2, 3}, + dg::DGMulti{<:Any, <:Wedge}, + bottom_topography) + rd = dg.basis + (; V1) = rd + (; xyz) = mesh.md + md = mesh.md + n_aux = n_aux_node_vars(equations) + + # Compute the radius of the sphere from the first element's fourth vertex, such that we can use it + # throughout the computation. We assume that each Wedge element's last three corner vertices lie + # on the simulated sphere. + VX, VY, VZ = map(coords -> transpose(coords[:, 1]) / V1', xyz) + v_outer = getindex.([VX, VY, VZ], 4) + radius = norm(v_outer) + + for element in eachelement(mesh, dg) + # Compute corner vertices of the element + VX, VY, VZ = map(coords -> transpose(coords[:, element]) / V1', xyz) + v1, v2, v3 = map(i -> getindex.([VX, VY, VZ], i), 4:6) + + aux_node = Vector{eltype(aux_values[1, 1])}(undef, n_aux) + + # Compute the auxiliary metric information at each node + for i in 1:Trixi.nnodes(dg) + # Covariant basis in the desired global coordinate system as columns of a matrix + basis_covariant = calc_basis_covariant(v1, v2, v3, + rd.rst[1][i], rd.rst[2][i], + radius, + equations.global_coordinate_system) + + aux_node[1:6] = SVector(basis_covariant) + + + # Covariant metric tensor G := basis_covariant' * basis_covariant + metric_covariant = basis_covariant' * basis_covariant + + # Contravariant metric tensor inv(G) + metric_contravariant = inv(metric_covariant) + + # Contravariant basis vectors as rows of a matrix + basis_contravariant = metric_contravariant * basis_covariant' + + + aux_node[7:12] = SVector(basis_contravariant) + # Area element + aux_node[13] = sqrt(det(metric_covariant)) + + # Covariant metric tensor components + aux_node[14:16] = SVector(metric_covariant[1, 1], + metric_covariant[1, 2], + metric_covariant[2, 2]) + + # Contravariant metric tensor components + aux_node[17:19] = SVector(metric_contravariant[1, 1], + metric_contravariant[1, 2], + metric_contravariant[2, 2]) + # Bottom topography + if !isnothing(bottom_topography) + nothing + end + aux_values[i, element] = SVector{n_aux}(aux_node) + end + # TODO: implement Christoffel symbols + # Christoffel symbols of the second kind (aux_values[21:26, :, :, element]) + # calc_christoffel_symbols!(aux_values, mesh, equations, dg, element) + end + + return nothing +end + +# Analytically compute the transformation matrix A, such that G = AᵀA is the +# covariant metric tensor and a_i = A[1,i] * e_x + A[2,i] * e_y + A[3,i] * e_z denotes +# the covariant tangent basis, where e_x, e_y, and e_z are the Cartesian unit basis vectors. +@inline function calc_basis_covariant(v1, v2, v3, xi1, xi2, radius, + ::GlobalCartesianCoordinates) + + # Construct a bilinear mapping based on the four corner vertices + xe = 0.5f0 * (-(xi1 + xi2) * v1 + (1 + xi1) * v2 + + (1 + xi2) * v3) + # Derivatives of bilinear map with respect to reference coordinates xi1, xi2 + dxedxi1 = 0.5f0 * + (-v1 + v2) + dxedxi2 = 0.5f0 * + (-v1 + v3) + + # Use product/quotient rule on the projection + norm_xe = norm(xe) + dxdxi1 = radius / norm_xe * (dxedxi1 - dot(xe, dxedxi1) / norm_xe^2 * xe) + dxdxi2 = radius / norm_xe * (dxedxi2 - dot(xe, dxedxi2) / norm_xe^2 * xe) + + return SMatrix{3, 2}(dxdxi1[1], dxdxi1[2], dxdxi1[3], + dxdxi2[1], dxdxi2[2], dxdxi2[3]) +end +end # @muladd diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl new file mode 100644 index 00000000..87224fb7 --- /dev/null +++ b/src/solvers/dgmulti/dg.jl @@ -0,0 +1,64 @@ +# Constructs cache variables including auxiliary variables for covariant equations and DGMultiMeshes +function Trixi.create_cache(mesh::DGMultiMesh{NDIMS}, equations::AbstractCovariantEquations, dg::Trixi.DGMultiWeakForm, + metric_terms, auxiliary_field, RealT, + uEltype) where {NDIMS} + rd = dg.basis + md = mesh.md + + # volume quadrature weights, volume interpolation matrix, mass matrix, differentiation matrices + @unpack wq, Vq, M, Drst = rd + + # ∫f(u) * dv/dx_i = ∑_j (Vq*Drst[i])'*diagm(wq)*(rstxyzJ[i,j].*f(Vq*u)) + weak_differentiation_matrices = map(D -> -M \ ((Vq * D)' * Diagonal(wq)), Drst) + + nvars = nvariables(equations) + naux = n_aux_node_vars(equations) + + # storage for volume quadrature values, face quadrature values, flux values + u_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xq), dg) + u_face_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xf), dg) + flux_face_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xf), dg) + aux_values = Trixi.allocate_nested_array(uEltype, naux, size(md.x), dg) + aux_quad_values = Trixi.allocate_nested_array(uEltype, naux, size(md.xq), dg) + aux_face_values = Trixi.allocate_nested_array(uEltype, naux, size(md.xf), dg) + + if typeof(rd.approximation_type) <: + Union{SBP, Trixi.AbstractNonperiodicDerivativeOperator} + lift_scalings = rd.wf ./ rd.wq[rd.Fmask] # lift scalings for diag-norm SBP operators + else + lift_scalings = nothing + end + + # local storage for volume integral and source computations + local_values_threaded = [Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:Threads.nthreads()] + + # For curved meshes, we interpolate geometric terms from nodal points to quadrature points. + # For affine meshes, we just access one element of this interpolated data. + dxidxhatj = map(x -> rd.Vq * x, md.rstxyzJ) + + init_auxiliary_node_variables!(aux_values, mesh, equations, dg, auxiliary_field) + + # Interpolate auxiliary variables to quadrature and face points + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), aux_quad_values, aux_values) + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vf), aux_face_values, aux_values) + + + # interpolate J to quadrature points for weight-adjusted DG (WADG) + invJ = inv.(rd.Vq * md.J) + + # for scaling by curved geometric terms (not used by affine DGMultiMesh) + flux_threaded = [[Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:NDIMS] for _ in 1:Threads.nthreads()] + rotated_flux_threaded = [Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:Threads.nthreads()] + + cache = (; md, weak_differentiation_matrices, lift_scalings, invJ, dxidxhatj, + u_values, u_face_values, flux_face_values, + aux_values, aux_quad_values, aux_face_values, + local_values_threaded, flux_threaded, rotated_flux_threaded) + return cache +end + +include("containers_manifold_covariant.jl") +include("dg_manifold_covariant.jl") diff --git a/src/solvers/dgmulti/dg_manifold_covariant.jl b/src/solvers/dgmulti/dg_manifold_covariant.jl new file mode 100644 index 00000000..3b9196ce --- /dev/null +++ b/src/solvers/dgmulti/dg_manifold_covariant.jl @@ -0,0 +1,102 @@ +# Compute coefficients for an initial condition that uses auxiliary variables +function Trixi.compute_coefficients!(::Nothing, u, initial_condition, t, + mesh::DGMultiMesh, equations::AbstractCovariantEquations, + dg::DGMulti, cache) + md = mesh.md + rd = dg.basis + @unpack u_values, aux_quad_values = cache + # evaluate the initial condition at quadrature points + Trixi.@threaded for i in Trixi.each_quad_node_global(mesh, dg, cache) + x_node = SVector(getindex.(md.xyzq, i)) + aux_node = aux_quad_values[i] + u_values[i] = initial_condition(x_node, t, aux_node, equations) + end + + # multiplying by Pq computes the L2 projection + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Pq), u, u_values) +end + +# version for covariant equations on DGMultiMeshes +function Trixi.calc_volume_integral!(du, u, mesh::DGMultiMesh{NDIMS_AMBIENT, <:Trixi.NonAffine}, + have_nonconservative_terms::False, + equations::AbstractCovariantEquations{NDIMS}, + volume_integral::VolumeIntegralWeakForm, dg::DGMulti, + cache) where {NDIMS_AMBIENT, NDIMS} + rd = dg.basis + md = mesh.md + (; weak_differentiation_matrices, u_values, aux_quad_values, local_values_threaded) = cache + + Jq = rd.Vq * md.J + + # interpolate to quadrature points + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u) + + Trixi.@threaded for e in Trixi.eachelement(mesh, dg, cache) + flux_values = local_values_threaded[Threads.threadid()] + for i in 1:NDIMS + for j in Trixi.eachindex(flux_values) + u_node = u_values[j, e] + aux_node = aux_quad_values[j, e] + area_elem = area_element(aux_node, equations) + J_node = Jq[j, e] + # Rescale the flux, such that the volume integral becomes independent of the thickness + # of the spherical shell. We compute that thickness by taking the ratio of the element's + # volume and the area element. + flux_values[j] = flux(u_node, aux_node, i, equations) * (J_node / area_elem) + end + + Trixi.apply_to_each_field(Trixi.mul_by_accum!(weak_differentiation_matrices[i]), + view(du, :, e), flux_values) + end + end +end + +# Calculate flux at interface by passing auxiliary variables to the surface flux function +function Trixi.calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm, + mesh::DGMultiMesh, + have_nonconservative_terms::False, + equations::AbstractCovariantEquations{NDIMS}, + dg::DGMulti{NDIMS_AMBIENT, <:Wedge}) where {NDIMS_AMBIENT, NDIMS} + @unpack surface_flux = surface_integral + md = mesh.md + rd = dg.basis + @unpack mapM, mapP, Jf = md + @unpack nrstJ, Nfq = rd + @unpack u_face_values, flux_face_values, aux_face_values = cache + Trixi.@threaded for face_node_index in Trixi.each_face_node_global(mesh, dg, cache) + + # inner (idM -> minus) and outer (idP -> plus) indices + idM, idP = mapM[face_node_index], mapP[face_node_index] + uM = u_face_values[idM] + uP = u_face_values[idP] + auxM = aux_face_values[idM] + auxP = aux_face_values[idP] + # Transform uP to the same coordinate system as uM + uP_global = contravariant2global(uP, auxP, equations) + uP_transformed_to_M = global2contravariant(uP_global, auxM, equations) + + # Compute ref_index from face_node_index in order to access covariant normal vector + # in reference element coordinates. We only take the first NDIMS components, projecting + # the normal vector onto the (reference) tangent space of the manifold, possibly yielding a zero vector. + ref_index = mod(face_node_index - 1, Nfq) + 1 + normal = SVector(ntuple(k -> nrstJ[k][ref_index], NDIMS)) + normal = normal / norm(normal) + + # If the unprojected normal vector belonged one of the triangular faces of the reference wedge, + # the normalized projected normal vector consists of NaNs. In this case, we skip the flux computation, + if any(isnan, normal) + continue + end + + # Compute the norm of the normal vector transformed to physical coordinates + basis = basis_contravariant(auxM, equations) + norm_normal_transformed = norm(transpose(basis) * normal) + area_elem = area_element(auxM, equations) + + # Scale the fluxes by inverse norm of transformed normal vector and + # the area element to get the correct fluxes per unit area, + # making the surface integral independent of the thickness of the spherical shell. + factor = 1 / norm_normal_transformed * Jf[idM] / area_elem + flux_face_values[idM] = surface_flux(uM, uP_transformed_to_M, auxM, auxM, normal, equations) * factor + end +end diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index cfcf2cf9..9a852cd7 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -1 +1,2 @@ include("dgsem_p4est/dg.jl") +include("dgmulti/dg.jl")