diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl index 0e21dbdb..667297ff 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl @@ -46,9 +46,7 @@ mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = 3, #initial_refinement_level = 0, # Comment to use cells_per_dimension in the convergence test element_local_mapping = false) -# Convert initial condition given in terms of zonal and meridional velocity components to -# one given in terms of Cartesian momentum components -initial_condition_transformed = transform_to_cartesian(initial_condition, equations) +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, solver, diff --git a/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl b/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl index e7a75937..1820d83d 100644 --- a/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl +++ b/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl @@ -46,9 +46,7 @@ mesh = P4estMeshQuadIcosahedron2D(cells_per_dimension[1], EARTH_RADIUS, #initial_refinement_level = 0, polydeg = polydeg) -# Convert initial condition given in terms of zonal and meridional velocity components to -# one given in terms of Cartesian momentum components -initial_condition_transformed = transform_to_cartesian(initial_condition, equations) +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, solver, diff --git a/examples/elixir_spherical_advection_covariant.jl b/examples/elixir_spherical_advection_covariant_cubed_sphere.jl similarity index 83% rename from examples/elixir_spherical_advection_covariant.jl rename to examples/elixir_spherical_advection_covariant_cubed_sphere.jl index 257f13e5..8193213d 100644 --- a/examples/elixir_spherical_advection_covariant.jl +++ b/examples/elixir_spherical_advection_covariant_cubed_sphere.jl @@ -1,16 +1,18 @@ ############################################################################### # DGSEM for the linear advection equation on the cubed sphere ############################################################################### +# To run a convergence test, use +# convergence_test("../examples/elixir_spherical_advection_covariant_cubed_sphere.jl", 4, cells_per_dimension = (3,3)) using OrdinaryDiffEq, Trixi, TrixiAtmo ############################################################################### # Spatial discretization -cells_per_dimension = 5 +cells_per_dimension = (5, 5) initial_condition = initial_condition_gaussian -equations = CovariantLinearAdvectionEquation2D() +equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = GlobalCartesianCoordinates()) # Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, @@ -18,15 +20,12 @@ solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, # Create a 2D cubed sphere mesh the size of the Earth. For the covariant form to work # properly, we currently need polydeg to equal that of the solver, -# initial_refinement_level = 0, and element_local_mapping = true. -mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, +# initial_refinement_level = 0 (default), and element_local_mapping = true. +mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = Trixi.polydeg(solver), - initial_refinement_level = 0, element_local_mapping = true) -# Convert initial condition given in terms of zonal and meridional velocity components to -# one given in terms of contravariant velocity components -initial_condition_transformed = transform_to_contravariant(initial_condition, equations) +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, solver) @@ -48,7 +47,7 @@ analysis_callback = AnalysisCallback(semi, interval = 10, # The SaveSolutionCallback allows to save the solution to a file in regular intervals save_solution = SaveSolutionCallback(interval = 10, - solution_variables = contravariant2spherical) + solution_variables = contravariant2global) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl = 0.7) diff --git a/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl b/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl new file mode 100644 index 00000000..76999765 --- /dev/null +++ b/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl @@ -0,0 +1,66 @@ +############################################################################### +# DGSEM for the linear advection equation on the cubed sphere +############################################################################### +# To run a convergence test, use +# convergence_test("../examples/elixir_spherical_advection_covariant_quad_icosahedron.jl", 4, cells_per_dimension = (1,1)) + +using OrdinaryDiffEq, Trixi, TrixiAtmo + +############################################################################### +# Spatial discretization + +cells_per_dimension = (2, 2) +initial_condition = initial_condition_gaussian + +equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = GlobalCartesianCoordinates()) + +# Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +# Create a 2D cubed sphere mesh the size of the Earth. For the covariant form to work +# properly, we currently need polydeg to equal that of the solver, and +# initial_refinement_level = 0 (default) +mesh = P4estMeshQuadIcosahedron2D(cells_per_dimension[1], EARTH_RADIUS, + polydeg = Trixi.polydeg(solver)) + +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, solver) + +############################################################################### +# 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,)) + +# 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 = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index ba96306c..f7815ba3 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: norm, dot +using LinearAlgebra: norm, dot, det using Reexport: @reexport using LoopVectorization: @turbo using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace @@ -32,6 +32,7 @@ include("callbacks_step/callbacks_step.jl") export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, CovariantLinearAdvectionEquation2D +export GlobalCartesianCoordinates, GlobalSphericalCoordinates export flux_chandrashekar, flux_LMARS @@ -43,9 +44,9 @@ export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProdu MetricTermsInvariantCurl export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, SECONDS_PER_DAY -export spherical2contravariant, contravariant2spherical, spherical2cartesian, - transform_to_cartesian, transform_to_contravariant -export initial_condition_gaussian +export global2contravariant, contravariant2global, spherical2cartesian, + transform_initial_condition +export initial_condition_gaussian, initial_condition_gaussian_cartesian export examples_dir end # module TrixiAtmo diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index 2f4edbb8..96c9c29b 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -1,7 +1,44 @@ @muladd begin #! format: noindent -# Entropy time derivative which uses auxiliary variables +# For the covariant form, we want to integrate using the exact area element +# √G = (det(AᵀA))^(1/2), which is stored in cache.auxiliary_variables, not the approximate +# area element used in the Cartesian formulation, which stored in cache.elements. +function Trixi.integrate_via_indices(func::Func, u, + mesh::Union{StructuredMesh{2}, + StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, + equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, args...; + normalize = true) where {Func} + (; weights) = dg.basis + (; aux_node_vars) = cache.auxiliary_variables + + # Initialize integral with zeros of the right shape + integral = zero(func(u, 1, 1, 1, equations, dg, args...)) + total_volume = zero(real(mesh)) + + # Use quadrature to numerically integrate over entire domain + for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + sqrtG = area_element(aux_node, equations) + integral += weights[i] * weights[j] * sqrtG * + func(u, i, j, element, equations, dg, args...) + total_volume += weights[i] * weights[j] * sqrtG + end + end + + # Normalize with total volume + if normalize + integral = integral / total_volume + end + + return integral +end + +# Entropy time derivative for cons2entropy function which depends on auxiliary variables function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t, mesh::P4estMesh{2}, equations::AbstractCovariantEquations{2}, dg::DG, cache) @@ -28,7 +65,6 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, (; weights) = dg.basis (; node_coordinates) = cache.elements (; aux_node_vars) = cache.auxiliary_variables - # Set up data structures l2_error = zero(func(Trixi.get_node_vars(u, equations, dg, 1, 1, 1), equations)) linf_error = copy(l2_error) @@ -51,7 +87,7 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, u_numerical = Trixi.get_node_vars(u, equations, dg, i, j, element) diff = func(u_exact, equations) - func(u_numerical, equations) - # For the L2 error, integrate with respect to volume element stored in aux vars + # For the L2 error, integrate with respect to area element stored in aux vars sqrtG = area_element(aux_node, equations) l2_error += diff .^ 2 * (weights[i] * weights[j] * sqrtG) 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 4982812e..98079255 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 @@ -42,7 +42,9 @@ end function Trixi.save_solution_file(u, time, dt, timestep, mesh::P4estMesh{2}, - equations::AbstractEquations{3}, dg::DG, cache, + equations::Union{AbstractEquations{3}, + AbstractCovariantEquations{2}}, + dg::DG, cache, solution_callback, element_variables = Dict{Symbol, Any}(), node_variables = Dict{Symbol, Any}(); diff --git a/src/callbacks_step/save_solution_covariant.jl b/src/callbacks_step/save_solution_covariant.jl index d112146d..b39da3a3 100644 --- a/src/callbacks_step/save_solution_covariant.jl +++ b/src/callbacks_step/save_solution_covariant.jl @@ -1,13 +1,17 @@ +@muladd begin +#! format: noindent + # Convert to another set of variables using a solution_variables function that depends on # auxiliary variables function convert_variables(u, solution_variables, mesh::P4estMesh{2}, - equations, dg, cache) + equations::AbstractCovariantEquations{2}, dg, cache) (; aux_node_vars) = cache.auxiliary_variables # Extract the number of solution variables to be output # (may be different than the number of conservative variables) n_vars = length(solution_variables(Trixi.get_node_vars(u, equations, dg, 1, 1, 1), - get_node_aux_vars(aux_node_vars, equations, dg, 1, 1, + get_node_aux_vars(aux_node_vars, equations, dg, + 1, 1, 1), equations)) # Allocate storage for output data = Array{eltype(u)}(undef, n_vars, nnodes(dg), nnodes(dg), nelements(dg, cache)) @@ -25,29 +29,4 @@ function convert_variables(u, solution_variables, mesh::P4estMesh{2}, end return data end - -# Specialized save_solution_file method that supports a solution_variables function which -# depends on auxiliary variables. The conversion must be defined as solution_variables(u, -# aux_vars, equations), and an additional method must be defined as solution_variables(u, -# equations) = u, such that no conversion is done when auxiliary variables are not provided. -function Trixi.save_solution_file(u_ode, t, dt, iter, - semi::SemidiscretizationHyperbolic{<:Trixi.AbstractMesh, - <:AbstractCovariantEquations}, - solution_callback, - element_variables = Dict{Symbol, Any}(), - node_variables = Dict{Symbol, Any}(); - system = "") - mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(semi) - u = Trixi.wrap_array_native(u_ode, semi) - - # Perform the variable conversion at each node - data = convert_variables(u, solution_callback.solution_variables, mesh, equations, - solver, cache) - - # Call the existing Trixi.save_solution_file, which will use solution_variables(u, - # equations). Since the variables are already converted above, we therefore require - # solution_variables(u, equations) = u. - Trixi.save_solution_file(data, t, dt, iter, mesh, equations, solver, cache, - solution_callback, element_variables, - node_variables, system = system) -end +end # muladd diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index a94dc809..c38db84d 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -2,7 +2,8 @@ #! format: noindent @doc raw""" - CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3, 3} + CovariantLinearAdvectionEquation2D{GlobalCoordinateSystem} <: + AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} A variable-coefficient linear advection equation can be defined on a two-dimensional manifold $S \subset \mathbb{R}^3$ as @@ -14,7 +15,7 @@ as a system of equations in which the first variable is the scalar conserved qua and the second two are the contravariant components $v^1$ and $v^2$ used in the expansion with respect to the covariant basis vectors $\vec{a}_1$ and $\vec{a}_2$ as ```math -\vec{v} = v^1 \vec{a}_1 + v^2 \vec{a}_2. +\vec{v} = v^1 \vec{a}_1 + v^2 \vec{a}_2, ``` where $\vec{a}_1 = \partial \vec{x} / \partial \xi^1$ and $\vec{a}_2 = \partial \vec{x} / \partial \xi^2$ are the so-called covariant basis vectors, @@ -23,21 +24,33 @@ are spatially varying but assumed to be constant in time, so we do not apply any dissipation to such variables. The resulting system is then given on the reference element as ```math -\sqrt{G} \frac{\partial}{\partial t}\left[\begin{array}{c} h \\ v^1 \\ v^2 \end{array}\right] +\sqrt{G} \frac{\partial}{\partial t} +\left[\begin{array}{c} h \\ v^1 \\ v^2 \end{array}\right] + -\frac{\partial}{\partial \xi^1} \left[\begin{array}{c} \sqrt{G} h v^1 \\ 0 \\ 0 \end{array}\right] +\frac{\partial}{\partial \xi^1} +\left[\begin{array}{c} \sqrt{G} h v^1 \\ 0 \\ 0 \end{array}\right] + -\frac{\partial}{\partial \xi^2} \left[\begin{array}{c} \sqrt{G} h v^2 \\ 0 \\ 0 \end{array}\right] +\frac{\partial}{\partial \xi^2} +\left[\begin{array}{c} \sqrt{G} h v^2 \\ 0 \\ 0 \end{array}\right] = \left[\begin{array}{c} 0 \\ 0 \\ 0 \end{array}\right], ``` -where $G$ is the determinant of the covariant metric tensor expressed as a matrix with -entries $G_{ij} = \vec{a}_i \cdot \vec{a}_j$. Note that the variable advection velocity -components could also be stored as auxiliary variables, similarly to the geometric -information. +where $G$ is the determinant of the covariant metric tensor expressed as a 2 by 2 matrix +with entries $G_{ij} = \vec{a}_i \cdot \vec{a}_j$. Note that the variable advection +velocity components could alternatively be stored as auxiliary variables, similarly to the +geometric information. """ -struct CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3, 3} end +struct CovariantLinearAdvectionEquation2D{GlobalCoordinateSystem} <: + AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} + global_coordinate_system::GlobalCoordinateSystem + function CovariantLinearAdvectionEquation2D(; + global_coordinate_system = GlobalCartesianCoordinates()) + return new{typeof(global_coordinate_system)}(global_coordinate_system) + end +end +# The conservative variables are the scalar conserved quantity and two contravariant +# velocity components. function Trixi.varnames(::typeof(cons2cons), ::CovariantLinearAdvectionEquation2D) return ("scalar", "vcon1", "vcon2") end @@ -47,42 +60,44 @@ function velocity(u, ::CovariantLinearAdvectionEquation2D) return SVector(u[2], u[3]) end -# Convert contravariant velocity/momentum components to zonal and meridional components -@inline function contravariant2spherical(u, aux_vars, - equations::CovariantLinearAdvectionEquation2D) - vlon, vlat = basis_covariant(aux_vars, equations) * velocity(u, equations) - return SVector(u[1], vlon, vlat) +# Convert contravariant velocity/momentum components to the global coordinate system +@inline function contravariant2global(u, aux_vars, + equations::CovariantLinearAdvectionEquation2D) + vglo1, vglo2, vglo3 = basis_covariant(aux_vars, equations) * velocity(u, equations) + return SVector(u[1], vglo1, vglo2, vglo3) 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) - return SVector(u_spherical[1], vcon1, vcon2) +# Convert velocity/momentum components in the global coordinate system to contravariant +# components +@inline function global2contravariant(u, aux_vars, + equations::CovariantLinearAdvectionEquation2D) + vcon1, vcon2 = basis_contravariant(aux_vars, equations) * SVector(u[2], u[3], u[4]) + return SVector(u[1], vcon1, vcon2) end -function Trixi.varnames(::typeof(contravariant2spherical), +# Scalar conserved quantity and three global velocity components +function Trixi.varnames(::typeof(contravariant2global), ::CovariantLinearAdvectionEquation2D) - return ("scalar", "vlon", "vlat") + return ("scalar", "vglo1", "vglo2", "vglo3") end # We will define the "entropy variables" here to just be the scalar variable in the first -# slot, with zeros in the second and third positions +# slot, with zeros in all other positions @inline function Trixi.cons2entropy(u, aux_vars, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u)) - return SVector{3}(u[1], z, z) + return SVector(u[1], z, z) end -# Numerical flux plus dissipation for abstract covariant equations as a function of the +# Flux for abstract covariant equations as a function of the # state vector u, as well as the auxiliary variables aux_vars, which contain the geometric -# information. +# information required for the covariant form @inline function Trixi.flux(u, aux_vars, orientation::Integer, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u)) - return SVector(area_element(aux_vars, equations) * u[1] * u[orientation + 1], z, - z) + sqrtG = area_element(aux_vars, equations) + vcon = velocity(u, equations) + return SVector(sqrtG * u[1] * vcon[orientation], z, z) end # Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity @@ -92,18 +107,19 @@ end orientation_or_normal_direction, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u_ll)) + sqrtG = area_element(aux_vars_ll, equations) λ = dissipation.max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr, orientation_or_normal_direction, equations) - return -0.5f0 * area_element(aux_vars_ll, equations) * λ * - SVector(u_rr[1] - u_ll[1], z, z) + return -0.5f0 * sqrtG * λ * SVector(u_rr[1] - u_ll[1], z, z) end -# Maximum wave speed with respect to the a specific orientation +# Maximum contravariant wave speed with respect to specific basis vector @inline function Trixi.max_abs_speed_naive(u_ll, u_rr, aux_vars_ll, aux_vars_rr, orientation::Integer, equations::CovariantLinearAdvectionEquation2D) - return max(abs(velocity(u_ll, equations)[orientation]), - abs(velocity(u_rr, equations)[orientation])) + vcon_ll = velocity(u_ll, equations) # Contravariant components on left side + vcon_rr = velocity(u_rr, equations) # Contravariant components on right side + return max(abs(vcon_ll[orientation]), abs(vcon_rr[orientation])) end # Maximum wave speeds in each direction for CFL calculation diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ae6c74ff..8cc53932 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -6,6 +6,7 @@ using Trixi: AbstractEquations @doc raw""" AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT, + GlobalCoordinateSystem, NVARS} <: AbstractEquations{NDIMS, NVARS} Abstract type used to dispatch on systems of equations in covariant form, in which fluxes @@ -30,14 +31,34 @@ Some references on discontinuous Galerkin methods in covariant flux form are lis When using this equation type, functions which are evaluated pointwise, such as fluxes, 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 [`transform_to_contravariant`](@ref). +the geometric information needed for the covariant form. The type parameter +`GlobalCoordinateSystem` specifies the global coordinate system used to define the +covariant tangent basis, and may be either [`GlobalCartesianCoordinates`](@ref) or +[`GlobalSphericalCoordinates`](@ref). The `GlobalCoordinateSystem` type parameter also +specifies the coordinate system with respect to which the initial condition should be +prescribed. """ abstract type AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT, + GlobalCoordinateSystem, NVARS} <: AbstractEquations{NDIMS, NVARS} end +""" + GlobalCartesianCoordinates() + +Struct used for dispatch, specifying that the covariant tangent basis vectors should be +defined with respect to a global Cartesian coordinate system. +""" +struct GlobalCartesianCoordinates end + +""" + GlobalSphericalCoordinates() + +Struct used for dispatch, specifying that the covariant tangent basis vectors should be +defined with respect to a global spherical coordinate system. +""" +struct GlobalSphericalCoordinates end + """ have_aux_node_vars(equations) @@ -48,48 +69,118 @@ dispatching on the return type. """ @inline have_aux_node_vars(::AbstractEquations) = False() -# cons2cons method which takes in unused aux_vars variable -@inline Trixi.cons2cons(u, aux_vars, equations) = u +@doc raw""" + transform_initial_condition(initial_condition, equations) + +Takes in a function with the signature `initial_condition(x, t, equations)` which returns +an initial condition given in terms of global Cartesian or zonal/meridional velocity +components, and returns another function `initial_condition_transformed(x, t, equations)` +or `initial_condition_transformed(x, t, aux_vars, equations)` which returns the same +initial data, but transformed to the appropriate prognostic variables used internally by +the solver. For the covariant form, this involves a transformation of the global velocity +components to contravariant components using [`global2contravariant`](@ref) as well as a +conversion from primitive to conservative variables. For standard Cartesian formulations, +this simply involves a conversion from primitive to conservative variables. The intention +here is to have a set of test cases (for example, [`initial_condition_gaussian`](@ref)) for +which the initial condition is prescribed using a standardized set of primitive variables +in a global Cartesian coordinates, and transformed to the specific prognostic variables +required for a given model. +!!! note + When using the covariant formulation, the initial velocity components should be defined + in the coordinate system specified by the `GlobalCoordinateSystem` type parameter in + [`AbstractCovariantEquations`](@ref). +""" +function transform_initial_condition(initial_condition, ::AbstractCovariantEquations) + function initial_condition_transformed(x, t, aux_vars, equations) + return Trixi.prim2cons(global2contravariant(initial_condition(x, t, equations), + aux_vars, equations), aux_vars, + equations) + end + return initial_condition_transformed +end -# If no auxiliary variables are passed into the conversion to spherical coordinates, do not -# do any conversion. -@inline contravariant2spherical(u, equations) = u +# Version for standard Cartesian formulations +function transform_initial_condition(initial_condition, ::AbstractEquations) + function initial_condition_transformed(x, t, equations) + return Trixi.prim2cons(initial_condition(x, t, equations), equations) + end + return initial_condition_transformed +end -# For the covariant form, the auxiliary variables are the the NDIMS^2 entries of the -# covariant basis matrix +""" + contravariant2global(u, aux_vars, equations) + +Transform the vector `u` of solution variables with the momentum or velocity given in terms +of local contravariant components into the global coordinate system specified by the +`GlobalCoordinateSystem` type parameter in [`AbstractCovariantEquations`](@ref). `u` is a +vector type of the correct length `nvariables(equations)`. Notice the function doesn't +include any error checks for the purpose of efficiency, so please make sure your input is +correct. The inverse conversion is performed by [`global2contravariant`](@ref). +""" +function contravariant2global end + +""" + global2contravariant(u, aux_vars, equations) + +Transform the vector `u` of solution variables with momentum or velocity components +given with respect to the global coordinate system into local contravariant components. The +global coordinate system is specified by the `GlobalCoordinateSystem` type parameter in +[`AbstractCovariantEquations`](@ref). `u` is a vector type of the correct length +`nvariables(equations)`. Notice the function doesn't include any error checks for the +purpose of efficiency, so please make sure your input is correct. The inverse conversion is +performed by [`contravariant2global`](@ref). +""" +function global2contravariant end + +# cons2cons method which takes in unused aux_vars variable +@inline Trixi.cons2cons(u, aux_vars, ::AbstractEquations) = u +@inline Trixi.prim2cons(u, aux_vars, ::AbstractEquations) = u + +# For the covariant form, the auxiliary variables are the the NDIMS*NDIMS_AMBIENT entries +# of the exact covariant basis matrix, the NDIMS*NDIMS_AMBIENT entries of the exact +# contravariant basis matrix, and the area element. In the future, we will add other terms +# such as Christoffel symbols. @inline have_aux_node_vars(::AbstractCovariantEquations) = True() -@inline n_aux_node_vars(::AbstractCovariantEquations{NDIMS}) where {NDIMS} = NDIMS^2 +@inline n_aux_node_vars(::AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT}) where {NDIMS, +NDIMS_AMBIENT} = 2 * NDIMS * NDIMS_AMBIENT + 1 # Return auxiliary variable names for 2D covariant form @inline function auxvarnames(::AbstractCovariantEquations{2}) - return ("basis_covariant[1,1]", "basis_covariant[2,1]", - "basis_covariant[1,2]", "basis_covariant[2,2]") + return ("basis_covariant[1,1]", + "basis_covariant[2,1]", + "basis_covariant[3,1]", + "basis_covariant[1,2]", + "basis_covariant[2,2]", + "basis_covariant[3,2]", + "basis_contravariant[1,1]", + "basis_contravariant[2,1]", + "basis_contravariant[3,1]", + "basis_contravariant[1,2]", + "basis_contravariant[2,2]", + "basis_contravariant[3,2]", + "area_element") end # Extract the covariant basis vectors a_i from the auxiliary variables as a matrix A, -# where a_i = A[:, i] denotes the covariant tangent basis in a spherical/ellipsoidal -# coordinate system. +# where A[:, i] contains the components of the ith covariant tangent basis vector with +# respect to the global (Cartesian or spherical) coordinate system @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 area_element(aux_vars, ::AbstractCovariantEquations{2}) - return abs(aux_vars[1] * aux_vars[4] - aux_vars[2] * aux_vars[3]) + return SMatrix{3, 2}(aux_vars[1], aux_vars[2], aux_vars[3], + aux_vars[4], aux_vars[5], aux_vars[6]) end -@doc raw""" - transform_to_contravariant(initial_condition, equations) +# Extract the contravariant basis vectors a^i from the auxiliary variables as a matrix B, +# where B[i, :] contains the components of the ith contravariant tangent basis vector with +# respect to the global (Cartesian or spherical) coordinate system +@inline function basis_contravariant(aux_vars, ::AbstractCovariantEquations{2}) + return SMatrix{2, 3}(aux_vars[7], aux_vars[8], + aux_vars[9], aux_vars[10], + aux_vars[11], aux_vars[12]) +end -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, aux_vars, -equations)` which returns the same initial condition with the velocity or momentum vector -given in terms of contravariant components. -""" -function transform_to_contravariant(initial_condition, ::AbstractCovariantEquations) - function initial_condition_transformed(x, t, aux_vars, equations) - return spherical2contravariant(initial_condition(x, t), aux_vars, equations) - end - return initial_condition_transformed +# Extract the area element √G = (det(AᵀA))^(1/2) from the auxiliary variables +@inline function area_element(aux_vars, ::AbstractCovariantEquations{2}) + return aux_vars[13] end # Numerical flux plus dissipation for abstract covariant equations as a function of the diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index a3b68833..68652238 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -8,78 +8,103 @@ const EARTH_ROTATION_RATE = 7.292e-5 # rad/s const SECONDS_PER_DAY = 8.64e4 @doc raw""" - initial_condition_gaussian(x, t) + initial_condition_gaussian(x, t, equations) This Gaussian bell case is a smooth initial condition suitable for testing the convergence -of discretizations of the linear advection equation on a spherical domain of radius $a = 6. -37122 \times 10^3\ \mathrm{m}$, representing the surface of the Earth. The height field is -given in terms of the Cartesian coordinates $x$, $y$, $z$ as +of discretizations of the linear advection equation on a spherical domain of radius $6. +37122 \times 10^3\ \mathrm{m}$, representing the surface of the Earth. Denoting the +Euclidean norm as $\lVert \cdot \rVert$, the initial height field is given by ```math -h(x,y,z) = h_0 \exp\Big(-b_0 \big((x-x_0)^2 + (y-y_0)^2 + (z-z_0)^2\big) \Big), +h(\vec{x}) = h_0 \exp +\Big(-b_0 \big(\lVert \vec{x} - \vec{x}_0 \rVert / \lVert \vec{x} \rVert\big)^2 \Big), ``` -where $h_0 = 1 \times 10^3\ \mathrm{m}$ is the height of the bell, $b_0 = 5 / a$ is the -width parameter, and the Cartesian coordinates of the bell's centre at longitude -$\lambda_0 = 3\pi/2$ and latitude $\theta_0 = 0$ are given by +where $h_0 = 1 \times 10^3\ \mathrm{m}$ is the height of the bell, $b_0 = 5$ is the +width parameter, and $\vec{x}_0$ is the position of the centre of the bell, which is +initialized at a longitude of $3\pi/2$ and a latitude of zero. The velocity field +corresponds to a solid body rotation with a period of 12 days at an angle of +$\alpha = \pi/4$ from the polar axis. Denoting $\vec{\omega}$ as the corresponding angular +velocity vector, the velocity is therefore initialized as ```math -\begin{aligned} -x_0 &= a\cos\theta_0\cos\lambda_0,\\ -y_0 &= a\cos\theta_0\sin\lambda_0,\\ -z_0 &= a\sin\theta_0. -\end{aligned} +\vec{v}(\vec{x}) = \vec{\omega} \times \vec{x}. ``` -The velocity field corresponds to a solid body rotation with a period of 12 days, at an -angle of $\alpha = \pi/4$ from the polar axis. The zonal and meridional components of the -velocity field are then given in terms of the longitude $\lambda$ and latitude $\theta$ as -```math -\begin{aligned} -u(\lambda,\theta) &= V (\cos\theta\cos\alpha + \sin\theta\cos\lambda\sin\alpha),\\ -v(\lambda,\theta) &= -V \sin\lambda\sin\alpha, -\end{aligned} -``` -where we take $V = 2\pi a / 12\ \mathrm{days}$. This test case is adapted from Case 1 of -the test suite described in the following paper: +This problem is adapted from Case 1 of the test suite described in the following paper: - D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. [DOI: 10.1016/S0021-9991(05)80016-6](https://doi.org/10.1016/S0021-9991(05)80016-6) - -This function returns `SVector(h, vlon, vlat, b)`, where the first three entries are the -height, zonal velocity, and meridional velocity. The fourth entry, representing variable -bottom topography, is set to zero. The functions [`transform_to_contravariant`](@ref) and -[`transform_to_cartesian`](@ref) are available for converting to the prognostic variables -for Cartesian and covariant formulations. """ -@inline function initial_condition_gaussian(x, t) +@inline function initial_condition_gaussian(x, t, ::AbstractEquations) + RealT = eltype(x) + + a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) # radius of the sphere + omega = convert(RealT, 2π) / (12 * SECONDS_PER_DAY) # angular velocity + alpha = convert(RealT, π / 4) # angle of rotation + h_0 = 1000.0f0 # bump height in metres + b_0 = 5.0f0 # bump width parameter + lon_0, lat_0 = convert(RealT, 3π / 2), 0.0f0 # initial bump location + + # axis of rotation + axis = SVector(-cos(alpha), 0.0f0, sin(alpha)) + + # convert initial position to Cartesian coordinates + x_0 = SVector(a * cos(lat_0) * cos(lon_0), + a * cos(lat_0) * sin(lon_0), + a * sin(lat_0)) + + # apply rotation using Rodrigues' formula + axis_cross_x_0 = Trixi.cross(axis, x_0) + x_0 = x_0 + sin(omega * t) * axis_cross_x_0 + + (1 - cos(omega * t)) * Trixi.cross(axis, axis_cross_x_0) + + # compute Gaussian bump profile + h = h_0 * + exp(-b_0 * ((x[1] - x_0[1])^2 + (x[2] - x_0[2])^2 + (x[3] - x_0[3])^2) / (a^2)) + + # get Cartesian velocity components + vx, vy, vz = omega * Trixi.cross(axis, x) + + # Prescribe the rotated bell shape and Cartesian velocity components. + # The last variable is the bottom topography, which we set to zero + return SVector(h, vx, vy, vz, 0.0f0) +end + +# Version for spherical coordinates (note: the velocity is not well defined at the poles) +@inline function initial_condition_gaussian(x, t, + ::AbstractCovariantEquations{2, 3, + GlobalSphericalCoordinates}) RealT = eltype(x) - a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) #radius of the sphere - V = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY) # speed of rotation + a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) # radius of the sphere + omega = convert(RealT, 2π) / (12 * SECONDS_PER_DAY) # angular velocity alpha = convert(RealT, π / 4) # angle of rotation h_0 = 1000.0f0 # bump height in metres - b_0 = 5.0f0 / (a^2) # bump width + b_0 = 5.0f0 # bump width parameter lon_0, lat_0 = convert(RealT, 3π / 2), 0.0f0 # initial bump location + # axis of rotation + axis = SVector(-cos(alpha), 0.0f0, sin(alpha)) + # convert initial position to Cartesian coordinates x_0 = SVector(a * cos(lat_0) * cos(lon_0), a * cos(lat_0) * sin(lon_0), a * sin(lat_0)) + # apply rotation using Rodrigues' formula + axis_cross_x_0 = Trixi.cross(axis, x_0) + x_0 = x_0 + sin(omega * t) * axis_cross_x_0 + + (1 - cos(omega * t)) * Trixi.cross(axis, axis_cross_x_0) + # compute Gaussian bump profile - h = h_0 * exp(-b_0 * ((x[1] - x_0[1])^2 + (x[2] - x_0[2])^2 + (x[3] - x_0[3])^2)) + h = h_0 * + exp(-b_0 * ((x[1] - x_0[1])^2 + (x[2] - x_0[2])^2 + (x[3] - x_0[3])^2) / (a^2)) # get zonal and meridional components of the velocity - r_xy = x[1]^2 + x[2]^2 - if r_xy == 0.0 - lon = pi / 2 - else - lon = atan(x[2], x[1]) - end - - lat = asin(x[3] / a) - vlon = V * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha)) - vlat = -V * sin(lon) * sin(alpha) - - # the last variable is the bottom topography, which we set to zero - return SVector(h, vlon, vlat, 0.0) + lon, lat = atan(x[2], x[1]), asin(x[3] / a) + vlon = omega * a * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha)) + vlat = -omega * a * sin(lon) * sin(alpha) + + # Prescribe the rotated bell shape and spherical velocity components. + # The last variable is the bottom topography, which we set to zero + return SVector(h, vlon, vlat, 0.0f0, 0.0f0) end end # muladd diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl index 8340868c..930af2a0 100644 --- a/src/equations/shallow_water_3d.jl +++ b/src/equations/shallow_water_3d.jl @@ -380,47 +380,4 @@ end @inline function energy_internal(cons, equations::ShallowWaterEquations3D) return energy_total(cons, equations) - energy_kinetic(cons, equations) end - -# Transform zonal and meridional velocity/momentum components to Cartesian components -function spherical2cartesian(vlon, vlat, x) - # Co-latitude - colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) - - # Longitude - if sign(x[2]) == 0.0 - signy = 1.0 - else - signy = sign(x[2]) - end - r_xy = sqrt(x[1]^2 + x[2]^2) - if r_xy == 0.0 - lon = pi / 2 - else - lon = signy * acos(x[1] / r_xy) - end - - v1 = -cos(colat) * cos(lon) * vlat - sin(lon) * vlon - v2 = -cos(colat) * sin(lon) * vlat + cos(lon) * vlon - v3 = sin(colat) * vlat - - return SVector(v1, v2, v3) -end - -@doc raw""" - transform_to_cartesian(initial_condition, equations) - -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)` which returns the same initial condition -with the velocity or momentum vector given in terms of Cartesian components. -""" -function transform_to_cartesian(initial_condition, ::ShallowWaterEquations3D) - function initial_condition_transformed(x, t, equations) - h, vlon, vlat, b = initial_condition(x, t) - v1, v2, v3 = spherical2cartesian(vlon, vlat, x) - return Trixi.prim2cons(SVector(h, v1, v2, v3, b), equations) - end - return initial_condition_transformed -end end # @muladd diff --git a/src/meshes/p4est_cubed_sphere.jl b/src/meshes/p4est_cubed_sphere.jl index 7a5e563f..7b48b59a 100644 --- a/src/meshes/p4est_cubed_sphere.jl +++ b/src/meshes/p4est_cubed_sphere.jl @@ -43,7 +43,6 @@ The mesh will have no boundaries. form, and we require `initial_refinement_level = 0` for such cases. Furthermore, the calculation of the metric terms for the covariant form currently requires `polydeg` to be equal to the polynomial degree of the solver, and `element_local_mapping = true`. -!!! """ function P4estMeshCubedSphere2D(trees_per_face_dimension, radius; polydeg, RealT = Float64, diff --git a/src/meshes/p4est_icosahedron_quads.jl b/src/meshes/p4est_icosahedron_quads.jl index f15677c9..ed6e6f20 100644 --- a/src/meshes/p4est_icosahedron_quads.jl +++ b/src/meshes/p4est_icosahedron_quads.jl @@ -38,7 +38,6 @@ The mesh will have no boundaries. form, and we require `initial_refinement_level = 0` for such cases. Furthermore, the calculation of the metric terms for the covariant form currently requires `polydeg` to be equal to the polynomial degree of the solver. -!!! """ function P4estMeshQuadIcosahedron2D(trees_per_face_dimension, radius; polydeg, RealT = Float64, diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl index b222c57d..5f675095 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl @@ -125,13 +125,15 @@ end # This assumes a sphere, although the radius can be arbitrary, and general mesh topologies are supported. function init_elements_2d_manifold_in_3d!(elements, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::P4estMesh{2, 3}, basis::LobattoLegendreBasis, metric_terms) (; node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian) = elements - calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis) + # The standard calc_node_coordinates! can be used, since Trixi.jl now dispatches on + # P4estMesh{NDIMS, NDIMS_AMBIENT}, so it can be used here. + Trixi.calc_node_coordinates!(node_coordinates, mesh, basis) for element in 1:Trixi.ncells(mesh) @@ -161,69 +163,6 @@ function init_elements_2d_manifold_in_3d!(elements, return nothing end -# Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis -function calc_node_coordinates_2d_shell!(node_coordinates, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) - # Hanging nodes will cause holes in the mesh if its polydeg is higher - # than the polydeg of the solver. - @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" - @assert size(mesh.tree_node_coordinates, 1)==3 "Shell must use 3D node coordinates" - calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis.nodes) -end - -# Interpolate tree_node_coordinates to each quadrant at the specified nodes -function calc_node_coordinates_2d_shell!(node_coordinates, - mesh::P4estMesh{2}, - nodes::AbstractVector) - # We use `StrideArray`s here since these buffers are used in performance-critical - # places and the additional information passed to the compiler makes them faster - # than native `Array`s. - tmp1 = Trixi.StrideArray(undef, real(mesh), - Trixi.StaticInt(3), - Trixi.static_length(nodes), - Trixi.static_length(mesh.nodes)) - matrix1 = Trixi.StrideArray(undef, real(mesh), - Trixi.static_length(nodes), - Trixi.static_length(mesh.nodes)) - matrix2 = similar(matrix1) - baryweights_in = Trixi.barycentric_weights(mesh.nodes) - - # Macros from `p4est` - p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL - p4est_quadrant_len(l) = 1 << (Trixi.P4EST_MAXLEVEL - l) - - trees = Trixi.unsafe_wrap_sc(Trixi.p4est_tree_t, mesh.p4est.trees) - - for tree in eachindex(trees) - offset = trees[tree].quadrants_offset - quadrants = Trixi.unsafe_wrap_sc(Trixi.p4est_quadrant_t, trees[tree].quadrants) - - for i in eachindex(quadrants) - element = offset + i - quad = quadrants[i] - - quad_length = p4est_quadrant_len(quad.level) / p4est_root_len - - nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ - quad.x / p4est_root_len) .- 1 - nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ - quad.y / p4est_root_len) .- 1 - Trixi.polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, - baryweights_in) - Trixi.polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, - baryweights_in) - - Trixi.multiply_dimensionwise!(view(node_coordinates, :, :, :, element), - matrix1, matrix2, - view(mesh.tree_node_coordinates, :, :, :, - tree), tmp1) - end - end - - return node_coordinates -end - # Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, # using eq (12) of : # Giraldo, F. X. (2001). A spectral element shallow water model on spherical geodesic grids. 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 ed411c20..07d595d1 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 @@ -161,6 +161,10 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, equations::AbstractCovariantEquations{2, 3}, dg, elements) (; tree_node_coordinates) = mesh + (; aux_node_vars) = auxiliary_variables + + NDIMS = 2 + NDIMS_AMBIENT = 3 # Check that the degree of the mesh matches that of the solver @assert length(mesh.nodes) == nnodes(dg) @@ -182,33 +186,65 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, # Compute the auxiliary metric information at each node for j in eachnode(dg), i in eachnode(dg) - # compute the covariant basis matrix - basis_covariant = basis_covariant_guba_etal(v1, v2, v3, v4, - dg.basis.nodes[i], - dg.basis.nodes[j], radius) - # Set auxiliary node variables in the cache - auxiliary_variables.aux_node_vars[1:4, i, j, element] = SVector(basis_covariant) + + # Calculate the covariant basis in the desired global coordinate system + A = calc_basis_covariant(v1, v2, v3, v4, + dg.basis.nodes[i], dg.basis.nodes[j], + radius, equations.global_coordinate_system) + aux_node_vars[1:(NDIMS * NDIMS_AMBIENT), i, j, element] = SVector(A) + + G = A' * A # Covariant metric tensor (not used for now) + + # Contravariant basis + aux_node_vars[(NDIMS * NDIMS_AMBIENT + 1):(2 * NDIMS * NDIMS_AMBIENT), + i, j, element] = SVector(inv(G) * A') + + # Area element + aux_node_vars[n_aux_node_vars(equations), i, j, element] = sqrt(det(G)) end 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, v4, xi1, xi2, radius, + ::GlobalCartesianCoordinates) + + # Construct a bilinear mapping based on the four corner vertices + xe = 0.25f0 * ((1 - xi1) * (1 - xi2) * v1 + (1 + xi1) * (1 - xi2) * v2 + + (1 + xi1) * (1 + xi2) * v3 + (1 - xi1) * (1 + xi2) * v4) + + # Derivatives of bilinear map with respect to reference coordinates xi1, xi2 + dxedxi1 = 0.25f0 * + (-(1 - xi2) * v1 + (1 - xi2) * v2 + (1 + xi2) * v3 - (1 + xi2) * v4) + dxedxi2 = 0.25f0 * + (-(1 - xi1) * v1 - (1 + xi1) * v2 + (1 + xi1) * v3 + (1 - xi1) * v4) + + # 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 + # Analytically compute the transformation matrix A, such that G = AᵀA is the # 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. This formula is -# taken from Guba et al. (2014), and it is not specific to the cubed sphere nor -# is the resulting matrix singular at the poles -@inline function basis_covariant_guba_etal(v1, v2, v3, v4, xi1, xi2, radius) +# taken from Guba et al. (2014). +@inline function calc_basis_covariant(v1, v2, v3, v4, xi1, xi2, radius, + ::GlobalSphericalCoordinates) # Construct a bilinear mapping based on the four corner vertices - x_bilinear = 0.25f0 * - ((1 - xi1) * (1 - xi2) * v1 + (1 + xi1) * (1 - xi2) * v2 + - (1 + xi1) * (1 + xi2) * v3 + (1 - xi1) * (1 + xi2) * v4) + xe = 0.25f0 * ((1 - xi1) * (1 - xi2) * v1 + (1 + xi1) * (1 - xi2) * v2 + + (1 + xi1) * (1 + xi2) * v3 + (1 - xi1) * (1 + xi2) * v4) # Project the mapped local coordinates onto the sphere using a simple scaling - scaling_factor = radius / norm(x_bilinear) - x = scaling_factor * x_bilinear + scaling_factor = radius / norm(xe) + x = scaling_factor * xe # Convert Cartesian coordinates to longitude and latitude lon, lat = atan(x[2], x[1]), asin(x[3] / radius) @@ -226,12 +262,15 @@ end a33 = coslat # Compute the matrix A containing spherical components of the covariant basis - return 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) + 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) + + # Make zero component in the radial direction so the matrix has the right dimensions + return SMatrix{3, 2}(A[1, 1], A[2, 1], 0.0f0, A[1, 2], A[2, 2], 0.0f0) 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 8e926705..ec18e04c 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 @@ -139,9 +139,9 @@ end 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_spherical = contravariant2global(u_rr, aux_vars_rr, equations) + u_rr_transformed_to_ll = global2contravariant(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, @@ -153,9 +153,9 @@ end 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, - equations) + u_ll_spherical = contravariant2global(u_ll, aux_vars_ll, equations) + u_ll_transformed_to_rr = global2contravariant(u_ll_spherical, aux_vars_rr, + equations) if isodd(secondary_direction_index) flux_secondary = -surface_flux(u_ll_transformed_to_rr, u_rr, aux_vars_rr, aux_vars_rr, diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 72dd5e9a..e8f87fb6 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -39,17 +39,17 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_quad_icosahedron_shell_advection.jl"), l2=[ - 0.4570227714893776, - 11.807355540221533, - 4.311881740776439, - 11.807355540221321, + 0.45702277148735726, + 11.807355540175404, + 4.311881740745649, + 11.807355540181993, 0.0 ], linf=[ - 13.591965583197862, - 364.7641889538827, - 93.69731833991227, - 364.76418895387906, + 13.591965583200476, + 364.76418895396273, + 93.69731833993228, + 364.76418895397, 0.0 ]) # and small reference values @@ -67,17 +67,17 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_cubed_sphere_shell_advection.jl"), l2=[ - 0.893342967294085, - 22.84887991899238, - 9.75885058673732, - 22.84887991899168, + 0.8933429672952714, + 22.84887991902509, + 9.758850586757735, + 22.84887991902542, 0.0 ], linf=[ - 14.289456304679561, - 380.6958334082083, - 120.59259301648672, - 380.695833408201, + 14.289456304624764, + 380.6958334067349, + 120.59259301602742, + 380.69583340674217, 0.0 ], element_local_mapping=true) # and small reference values @@ -93,7 +93,7 @@ end @trixiatmo_testset "Spherical advection, covariant weak form, LLF surface flux" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_spherical_advection_covariant.jl"), + "elixir_spherical_advection_covariant_cubed_sphere.jl"), l2=[1.0007043506351705, 0.0, 0.0], linf=[14.235905681508598, 0.0, 0.0]) # Ensure that we do not have excessive memory allocations @@ -106,11 +106,27 @@ end end end +@trixiatmo_testset "Spherical advection, covariant weak form, LLF surface flux, global spherical coords" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_spherical_advection_covariant_cubed_sphere.jl"), + l2=[1.0007043506351705, 0.0, 0.0], + linf=[14.235905681508598, 0.0, 0.0], + global_coordinate_system=GlobalSphericalCoordinates()) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + # The covariant flux-differencing form should be equivalent to the weak form when the # arithmetic mean is used as the two-point flux @trixiatmo_testset "Spherical advection, covariant flux-differencing, central/LLF" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_spherical_advection_covariant.jl"), + "elixir_spherical_advection_covariant_cubed_sphere.jl"), l2=[1.0007043506351412, 0.0, 0.0], linf=[14.23590568150928, 0.0, 0.0], volume_integral=VolumeIntegralFluxDifferencing(flux_central)) @@ -127,7 +143,7 @@ end # Version with arithmetic mean used for both the volume and surface fluxes @trixiatmo_testset "Spherical advection, covariant flux-differencing, central/central" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_spherical_advection_covariant.jl"), + "elixir_spherical_advection_covariant_cubed_sphere.jl"), l2=[2.499889861385917, 0.0, 0.0], linf=[38.085244441156085, 0.0, 0.0], volume_integral=VolumeIntegralFluxDifferencing(flux_central), @@ -142,4 +158,19 @@ end end end +@trixiatmo_testset "Spherical advection on icosahedral grid, covariant weak form, LLF surface flux" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_spherical_advection_covariant_quad_icosahedron.jl"), + l2=[0.5183886767005157, 0.0, 0.0], + linf=[13.54834739856517, 0.0, 0.0]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + end # module