diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index df937c2..a8f1263 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -33,7 +33,8 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, for j in eachnode(dg), i in eachnode(dg) x = Trixi.get_node_coords(node_coordinates, equations, dg, i, j, element) - aux_vars_node = get_node_aux_vars(cache.elements.auxiliary_variables, equations, + aux_vars_node = get_node_aux_vars(cache.elements.auxiliary_variables, + equations, dg, i, j, element) u_exact = initial_condition(x, t, aux_vars_node, equations) diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index 1a7a06a..4d956af 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -49,7 +49,8 @@ end @inline function Trixi.flux(u, aux_vars, orientation::Integer, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u)) - return SVector(volume_element(aux_vars, equations) * u[1] * u[orientation + 1], z, z) + return SVector(volume_element(aux_vars, equations) * u[1] * u[orientation + 1], z, + z) end # Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity @@ -74,7 +75,7 @@ end end # Maximum wave speeds in each direction for CFL calculation -@inline function Trixi.max_abs_speeds(u, aux_vars, +@inline function Trixi.max_abs_speeds(u, aux_vars, equations::CovariantLinearAdvectionEquation2D) return abs.(velocity(u, equations)) end @@ -89,7 +90,8 @@ end # Convert zonal and meridional velocity/momentum components to contravariant components @inline function spherical2contravariant(u_spherical, aux_vars, equations::CovariantLinearAdvectionEquation2D) - vcon1, vcon2 = basis_covariant(aux_vars, equations) \ velocity(u_spherical, equations) + vcon1, vcon2 = basis_covariant(aux_vars, equations) \ + velocity(u_spherical, equations) return SVector(u_spherical[1], vcon1, vcon2) end end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 0bd5cc0..b8dd36a 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -29,19 +29,20 @@ Some references on discontinuous Galerkin methods in covariant flux form are lis [DOI: 10.1016/j.jcp.2013.11.033](https://doi.org/10.1016/j.jcp.2013.11.033) When using this equation type, functions which are evaluated pointwise, such as fluxes, -source terms, and initial conditions take in the extra arguments `cache`, `node`, and -`element`, corresponding to the `cache` field of a `SemidiscretizationHyperbolic`, the node -index (for tensor-product elements, this should be a tuple of length `NDIMS`), and the -element index, respectively. To convert an initial condition given in terms of global -spherical velocity or momentum components to one given in terms of local contravariant -components, see [`spherical2contravariant`](@ref). +source terms, and initial conditions take in the extra argument `aux_vars`, which contains +the geometric information needed for the covariant form. To convert an initial condition +given in terms of global spherical velocity or momentum components to one given in terms of +local contravariant components, see [`spherical2contravariant`](@ref). """ abstract type AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT, NVARS} <: AbstractEquations{NDIMS, NVARS} end -@inline nauxvars(equations::AbstractCovariantEquations{NDIMS}) where {NDIMS} = NDIMS^2 + - 1 +# For the covariant form, the auxiliary variables are the the NDIMS^2 entries of the +# covariant basis matrix +@inline nauxvars(equations::AbstractCovariantEquations{NDIMS}) where {NDIMS} = NDIMS^2 + +# By default, there are no auxiliary variables @inline nauxvars(equations::AbstractEquations) = 0 @doc raw""" @@ -49,10 +50,8 @@ abstract type AbstractCovariantEquations{NDIMS, Takes in a function with the signature `initial_condition(x, t)` which returns an initial condition given in terms of zonal and meridional velocity or momentum components, and -returns another function with the signature -`initial_condition_transformed(x, t, equations, cache, node, element)` which returns -the same initial condition with the velocity or momentum vector given in terms of -contravariant components. +returns another function with the signature `initial_condition_transformed(x, t, aux_vars, +equations)` which returns the same initial condition with the velocity or momentum vector given in terms of contravariant components. """ function spherical2contravariant(initial_condition, ::AbstractCovariantEquations) function initial_condition_transformed(x, t, aux_vars, equations) @@ -87,8 +86,12 @@ end end # Extract geometric information from auxiliary variables -@inline volume_element(aux_vars, ::AbstractCovariantEquations{2}) = aux_vars[1] -@inline basis_covariant(aux_vars, ::AbstractCovariantEquations{2}) = SMatrix{2, 2}(aux_vars[2:5]) +@inline function basis_covariant(aux_vars, ::AbstractCovariantEquations{2}) + return SMatrix{2, 2}(aux_vars[1], aux_vars[2], aux_vars[3], aux_vars[4]) +end +@inline function volume_element(aux_vars, ::AbstractCovariantEquations{2}) + return abs(aux_vars[1] * aux_vars[4] - aux_vars[2] * aux_vars[3]) +end abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 6937058..954362f 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -73,6 +73,6 @@ for Cartesian and covariant formulations. vlat = -V * sin(lon) * sin(alpha) # the last variable is the bottom topography, which we set to zero - return SVector(h, vlon, vlat) + return SVector(h, vlon, vlat, 0.0) end end # muladd diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl index 1cfd34a..c70e0f7 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl @@ -1,36 +1,36 @@ +@muladd begin +#! format: noindent + # Compute the auxiliary metric terms for the covariant form, assuming that the # domain is a spherical shell. We do not make any assumptions on the mesh topology. function compute_auxiliary_variables!(elements, mesh::P4estMesh{2, 3}, - equations::AbstractCovariantEquations{2, 3}, basis) - (; inverse_jacobian, auxiliary_variables) = elements + equations::AbstractCovariantEquations{2, 3}, dg) - # The tree node coordinates are assumed to be on the spherical shell centred around the # origin. Therefore, we can compute the radius once and use it throughout. + # The tree node coordinates are assumed to be on the spherical shell centred around the + # origin. Therefore, we can compute the radius once and use it throughout. radius = sqrt(mesh.tree_node_coordinates[1, 1, 1, 1]^2 + mesh.tree_node_coordinates[2, 1, 1, 1]^2 + mesh.tree_node_coordinates[3, 1, 1, 1]^2) for element in 1:Trixi.ncells(mesh) - # For now, we will only use the corner nodes from the P4estMesh, and construct the - # mapping and its associated metric terms analytically without any polynomial - # approximation following the approach described in by Guba et al. (see - # https://doi.org/10.5194/gmd-7-2803-2014, Appendix A). If the option - # "element_local_mapping = true" is used when constructing the mesh using - # P4estMeshCubedSphere2D and the polynomial degree of the mesh is the same as that - # of the solver, the node positions are equal to those obtained using the standard - # calc_node_coordinates! function. + # Check that the degree of the mesh matches that of the solver + @assert length(mesh.nodes) == nnodes(dg) # Extract the corner vertex positions from the P4estMesh - nnodes = length(mesh.nodes) - v1 = SVector{3}(view(mesh.tree_node_coordinates, :, 1, 1, element)) - v2 = SVector{3}(view(mesh.tree_node_coordinates, :, nnodes, 1, element)) - v3 = SVector{3}(view(mesh.tree_node_coordinates, :, nnodes, nnodes, element)) - v4 = SVector{3}(view(mesh.tree_node_coordinates, :, 1, nnodes, element)) + v1 = Trixi.get_node_coords(mesh.tree_node_coordinates, equations, dg, + 1, 1, element) + v2 = Trixi.get_node_coords(mesh.tree_node_coordinates, equations, dg, + nnodes(dg), 1, element) + v3 = Trixi.get_node_coords(mesh.tree_node_coordinates, equations, dg, + nnodes(dg), nnodes(dg), element) + v4 = Trixi.get_node_coords(mesh.tree_node_coordinates, equations, dg, + 1, nnodes(dg), element) - for j in eachnode(basis), i in eachnode(basis) + for j in eachnode(dg), i in eachnode(dg) # get reference node coordinates - xi1, xi2 = basis.nodes[i], basis.nodes[j] + xi1, xi2 = dg.basis.nodes[i], dg.basis.nodes[j] # Construct a bilinear mapping based on the four corner vertices x_bilinear = 0.25f0 * @@ -61,16 +61,20 @@ function compute_auxiliary_variables!(elements, mesh::P4estMesh{2, 3}, # covariant metric tensor and a_i = A[1,i] * e_lon + A[2,i] * e_lat denotes # the covariant tangent basis, where e_lon and e_lat are the unit basis vectors # in the longitudinal and latitudinal directions, respectively. - A = 0.25f0 * scaling_factor * SMatrix{2, 3}(-sinlon, 0, coslon, 0, 0, 1) * - SMatrix{3, 3}(a11, a21, a31, a12, a22, a32, a13, a23, a33) * - SMatrix{3, 4}(v1[1], v1[2], v1[3], v2[1], v2[2], v2[3], - v3[1], v3[2], v3[3], v4[1], v4[2], v4[3]) * - SMatrix{4, 2}(-1 + xi2, 1 - xi2, 1 + xi2, -1 - xi2, - -1 + xi1, -1 - xi1, 1 + xi1, 1 - xi1) + basis_covariant = 0.25f0 * scaling_factor * + SMatrix{2, 3}(-sinlon, 0, coslon, 0, 0, 1) * + SMatrix{3, 3}(a11, a21, a31, a12, a22, a32, a13, a23, + a33) * + SMatrix{3, 4}(v1[1], v1[2], v1[3], v2[1], v2[2], v2[3], + v3[1], v3[2], v3[3], v4[1], v4[2], v4[3]) * + SMatrix{4, 2}(-1 + xi2, 1 - xi2, 1 + xi2, -1 - xi2, + -1 + xi1, -1 - xi1, 1 + xi1, 1 - xi1) # Set variables in the cache - auxiliary_variables[1, i, j, element] = 1 / inverse_jacobian[i, j, element] - auxiliary_variables[2:5, i, j, element] = SVector(A) + elements.auxiliary_variables[1, i, j, element] = basis_covariant[1, 1] + elements.auxiliary_variables[2, i, j, element] = basis_covariant[2, 1] + elements.auxiliary_variables[3, i, j, element] = basis_covariant[1, 2] + elements.auxiliary_variables[4, i, j, element] = basis_covariant[2, 2] end end @@ -87,7 +91,8 @@ end end # Return the auxiliary variables at a given volume node index -@inline function get_node_aux_vars(auxiliary_variables, equations, solver::DG, indices...) +@inline function get_node_aux_vars(auxiliary_variables, equations, solver::DG, + indices...) SVector(ntuple(@inline(v->auxiliary_variables[v, indices...]), Val(nauxvars(equations)))) end @@ -102,6 +107,7 @@ end return aux_vars_ll, aux_vars_rr end +# Interface container storing surface values of the auxiliary variables mutable struct P4estInterfaceContainerVariableCoefficient{NDIMS, uEltype <: Real, NDIMSP2} <: Trixi.AbstractContainer @@ -117,9 +123,8 @@ mutable struct P4estInterfaceContainerVariableCoefficient{NDIMS, uEltype <: Real _node_indices::Vector{NTuple{NDIMS, Symbol}} end -@inline function Trixi.ninterfaces(interfaces::P4estInterfaceContainerVariableCoefficient) - size(interfaces.neighbor_ids, 2) -end +@inline Trixi.ninterfaces(interfaces::P4estInterfaceContainerVariableCoefficient) = size(interfaces.neighbor_ids, + 2) @inline Base.ndims(::P4estInterfaceContainerVariableCoefficient{NDIMS}) where {NDIMS} = NDIMS # Create interface container and initialize interface data. @@ -201,62 +206,4 @@ end return interfaces end - -function prolong2interfaces_auxiliary!(cache, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, - dg::DG) - (; elements, interfaces) = cache - index_range = eachnode(dg) - - Trixi.@threaded for interface in Trixi.eachinterface(dg, cache) - # Copy solution data from the primary element using "delayed indexing" with - # a start value and a step size to get the correct face and orientation. - # Note that in the current implementation, the interface will be - # "aligned at the primary element", i.e., the index of the primary side - # will always run forwards. - primary_element = interfaces.neighbor_ids[1, interface] - primary_indices = interfaces.node_indices[1, interface] - - i_primary_start, i_primary_step = Trixi.index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = Trixi.index_to_start_step_2d(primary_indices[2], - index_range) - - i_primary = i_primary_start - j_primary = j_primary_start - for i in eachnode(dg) - for v in axes(elements.auxiliary_variables, 1) - interfaces.auxiliary_variables[1, v, i, interface] = elements.auxiliary_variables[v, - i_primary, - j_primary, - primary_element] - end - i_primary += i_primary_step - j_primary += j_primary_step - end - - # Copy solution data from the secondary element using "delayed indexing" with - # a start value and a step size to get the correct face and orientation. - secondary_element = interfaces.neighbor_ids[2, interface] - secondary_indices = interfaces.node_indices[2, interface] - - i_secondary_start, i_secondary_step = Trixi.index_to_start_step_2d(secondary_indices[1], - index_range) - j_secondary_start, j_secondary_step = Trixi.index_to_start_step_2d(secondary_indices[2], - index_range) - - i_secondary = i_secondary_start - j_secondary = j_secondary_start - for i in eachnode(dg) - for v in axes(elements.auxiliary_variables, 1) - interfaces.auxiliary_variables[2, v, i, interface] = elements.auxiliary_variables[v, - i_secondary, - j_secondary, - secondary_element] - end - i_secondary += i_secondary_step - j_secondary += j_secondary_step - end - end - - return nothing -end +end # muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl index 458a461..f54dc5f 100644 --- a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl @@ -63,14 +63,14 @@ function Trixi.rhs!(du, u, t, mesh::P4estMesh{2}, return nothing end -# Evaluate the initial condition in spherical vector components, then -# transform to contravariant components for use as prognostic variables +# Evaluate the initial condition in spherical vector components, then transform to +# contravariant components for use as prognostic variables. function Trixi.compute_coefficients!(u, func, t, mesh::P4estMesh{2}, equations::AbstractCovariantEquations{2}, dg::DG, cache) # Store the auxiliary variables at the volume quadrature nodes in cache.elements - compute_auxiliary_variables!(cache.elements, mesh, equations, dg.basis) + compute_auxiliary_variables!(cache.elements, mesh, equations, dg) # Copy the appropriate interface values of the auxiliary variables into cache.interfaces prolong2interfaces_auxiliary!(cache, mesh, dg) @@ -85,7 +85,6 @@ function Trixi.compute_coefficients!(u, func, t, mesh::P4estMesh{2}, Trixi.set_node_vars!(u, u_node, equations, dg, i, j, element) end end - end # Weak form kernel which uses contravariant flux components, passing the element container @@ -188,13 +187,13 @@ end u_ll, u_rr = Trixi.get_surface_node_vars(u, equations, dg, primary_node_index, interface_index) aux_vars_ll, aux_vars_rr = get_surface_node_aux_vars(auxiliary_variables, equations, - dg, primary_node_index, + dg, primary_node_index, interface_index) - + # Compute flux in the primary element's coordinate system u_rr_spherical = contravariant2spherical(u_rr, aux_vars_rr, equations) - u_rr_transformed_to_ll = spherical2contravariant(u_rr_spherical, aux_vars_ll, equations) - + u_rr_transformed_to_ll = spherical2contravariant(u_rr_spherical, aux_vars_ll, + equations) if isodd(primary_direction_index) flux_primary = -surface_flux(u_rr_transformed_to_ll, u_ll, aux_vars_ll, aux_vars_ll, @@ -204,10 +203,10 @@ end aux_vars_ll, aux_vars_ll, (primary_direction_index + 1) >> 1, equations) end - + # Compute flux in the secondary element's coordinate system u_ll_spherical = contravariant2spherical(u_ll, aux_vars_ll, equations) - u_ll_transformed_to_rr = spherical2contravariant(u_ll_spherical, aux_vars_rr, + u_ll_transformed_to_rr = spherical2contravariant(u_ll_spherical, aux_vars_rr, equations) if isodd(secondary_direction_index) flux_secondary = -surface_flux(u_ll_transformed_to_rr, u_rr, @@ -227,4 +226,83 @@ end secondary_element_index] = flux_secondary[v] end end + +# Apply the exact Jacobian stored in auxiliary variables +function Trixi.apply_jacobian!(du, mesh::P4estMesh{2}, + equations::AbstractCovariantEquations{2}, + dg::DG, cache) + Trixi.@threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + aux_vars_node = get_node_aux_vars(cache.elements.auxiliary_variables, + equations, dg, i, j, element) + factor = -1 / volume_element(aux_vars_node, equations) + + for v in eachvariable(equations) + du[v, i, j, element] *= factor + end + end + end + + return nothing +end + +# Prolong the auxiliary variables and store in cache.interfaces +function prolong2interfaces_auxiliary!(cache, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + dg::DG) + (; elements, interfaces) = cache + index_range = eachnode(dg) + + Trixi.@threaded for interface in Trixi.eachinterface(dg, cache) + # Copy solution data from the primary element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + primary_element = interfaces.neighbor_ids[1, interface] + primary_indices = interfaces.node_indices[1, interface] + + i_primary_start, i_primary_step = Trixi.index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = Trixi.index_to_start_step_2d(primary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + for i in eachnode(dg) + for v in axes(elements.auxiliary_variables, 1) + interfaces.auxiliary_variables[1, v, i, interface] = elements.auxiliary_variables[v, + i_primary, + j_primary, + primary_element] + end + i_primary += i_primary_step + j_primary += j_primary_step + end + + # Copy solution data from the secondary element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + secondary_element = interfaces.neighbor_ids[2, interface] + secondary_indices = interfaces.node_indices[2, interface] + + i_secondary_start, i_secondary_step = Trixi.index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = Trixi.index_to_start_step_2d(secondary_indices[2], + index_range) + + i_secondary = i_secondary_start + j_secondary = j_secondary_start + for i in eachnode(dg) + for v in axes(elements.auxiliary_variables, 1) + interfaces.auxiliary_variables[2, v, i, interface] = elements.auxiliary_variables[v, + i_secondary, + j_secondary, + secondary_element] + end + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + return nothing +end end # @muladd