diff --git a/examples/elixir_shallowwater_covariant_geostrophic_balance.jl b/examples/elixir_shallowwater_covariant_geostrophic_balance.jl new file mode 100644 index 0000000..a1040d4 --- /dev/null +++ b/examples/elixir_shallowwater_covariant_geostrophic_balance.jl @@ -0,0 +1,81 @@ +############################################################################### +# DGSEM for the shallow water equations in covariant form on the cubed sphere +############################################################################### + +using OrdinaryDiffEq, Trixi, TrixiAtmo + +############################################################################### +# Parameters + +initial_condition = initial_condition_geostrophic_balance +polydeg = 3 +cells_per_dimension = 5 +n_saves = 10 +tspan = (0.0, 5.0 * SECONDS_PER_DAY) + +############################################################################### +# Spatial discretization + +mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg, + initial_refinement_level = 0, + element_local_mapping = true) + +equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, + EARTH_ROTATION_RATE, + global_coordinate_system = GlobalSphericalCoordinates()) + +# The covariant shallow water equations are treated as a nonconservative system in order to +# handle flux-differencing formulations of the covariant derivative. With +# VolumeIntegralWeakForm, there are actually no nonconservative terms, but we must still +# pass a no-op function "flux_nonconservative_weak_form" as the nonconservative surface flux +surface_flux = (flux_lax_friedrichs, flux_nonconservative_weak_form) + +# Create DG solver with polynomial degree = p +solver = DGSEM(polydeg = polydeg, volume_integral = VolumeIntegralWeakForm(), + surface_flux = surface_flux) + +# 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, solver, + source_terms = source_terms_weak_form) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0 to T +ode = semidiscretize(semi, tspan) + +# 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 = 200, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, + solution_variables = cons2cons) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.4) + +# 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 = 100.0, save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl new file mode 100644 index 0000000..881e621 --- /dev/null +++ b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl @@ -0,0 +1,81 @@ +############################################################################### +# Entropy-conservative DGSEM for the shallow water equations in covariant form +# on the cubed sphere +############################################################################### + +using OrdinaryDiffEq, Trixi, TrixiAtmo + +############################################################################### +# Parameters + +initial_condition = initial_condition_rossby_haurwitz +polydeg = 3 +cells_per_dimension = 5 +n_saves = 10 +tspan = (0.0, 7.0 * SECONDS_PER_DAY) + +############################################################################### +# Spatial discretization + +mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg, + initial_refinement_level = 0, + element_local_mapping = true) + +equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, + EARTH_ROTATION_RATE, + global_coordinate_system = GlobalCartesianCoordinates()) + +# Use entropy-conservative two-point fluxes for volume and surface terms +volume_flux = (flux_ec, flux_nonconservative_ec) +surface_flux = (flux_ec, flux_nonconservative_ec) + +# Create DG solver with polynomial degree = p +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# 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, solver, + source_terms = source_terms_ec) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0 to T +ode = semidiscretize(semi, tspan) + +# 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. Note that entropy should be conserved at the semi-discrete level. +analysis_callback = AnalysisCallback(semi, interval = 200, + save_analysis = true, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (entropy,)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, + solution_variables = cons2prim_and_vorticity) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.4) + +# 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 = 100.0, save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl index 667297f..727add3 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl @@ -46,6 +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) +# 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 diff --git a/examples/elixir_spherical_advection_covariant_cubed_sphere.jl b/examples/elixir_spherical_advection_covariant_cubed_sphere.jl index 8193213..1cc513e 100644 --- a/examples/elixir_spherical_advection_covariant_cubed_sphere.jl +++ b/examples/elixir_spherical_advection_covariant_cubed_sphere.jl @@ -25,6 +25,7 @@ mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = Trixi.polydeg(solver), element_local_mapping = true) +# 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 diff --git a/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl b/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl index 7699976..ddab0ab 100644 --- a/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl +++ b/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl @@ -1,5 +1,5 @@ ############################################################################### -# DGSEM for the linear advection equation on the cubed sphere +# DGSEM for the linear advection equation on a quadrilateral icosahedral grid ############################################################################### # To run a convergence test, use # convergence_test("../examples/elixir_spherical_advection_covariant_quad_icosahedron.jl", 4, cells_per_dimension = (1,1)) @@ -18,12 +18,13 @@ equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = Global 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 +# Create a 2D quadrilateral icosahedral 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)) +# 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 @@ -35,11 +36,12 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transform # 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 +# 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 +# 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,)) @@ -51,14 +53,16 @@ save_solution = SaveSolutionCallback(interval = 10, # 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 +# 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 +# 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); diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index f7815ba..f15b11a 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -31,22 +31,32 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") include("callbacks_step/callbacks_step.jl") export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, - CovariantLinearAdvectionEquation2D + CovariantLinearAdvectionEquation2D, CovariantShallowWaterEquations2D + export GlobalCartesianCoordinates, GlobalSphericalCoordinates export flux_chandrashekar, flux_LMARS +export flux_nonconservative_weak_form, flux_nonconservative_ec, source_terms_weak_form, + source_terms_ec + export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal, lake_at_rest_error, source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier! + export cons2prim_and_vorticity + export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct, MetricTermsInvariantCurl + export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, SECONDS_PER_DAY -export global2contravariant, contravariant2global, spherical2cartesian, + +export global2contravariant, contravariant2global, spherical2cartesian, cartesian2spherical, transform_initial_condition -export initial_condition_gaussian, initial_condition_gaussian_cartesian + +export initial_condition_gaussian, initial_condition_geostrophic_balance, + initial_condition_rossby_haurwitz export examples_dir end # module TrixiAtmo diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index 96c9c29..9b0c689 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -1,9 +1,28 @@ @muladd begin #! format: noindent +# When the equations are of type AbstractCovariantEquations, the functions which we would +# like to integrate depend on the solution as well as the auxiliary variables +function Trixi.integrate(func::Func, u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + equations::AbstractCovariantEquations{2}, dg::DG, + cache; normalize = true) where {Func} + (; aux_node_vars) = cache.auxiliary_variables + + Trixi.integrate_via_indices(u, mesh, equations, dg, cache; + normalize = normalize) do u, i, j, element, equations, + dg + u_local = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_local = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + return func(u_local, aux_local, equations) + end +end + # 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. +# J = (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}, @@ -23,10 +42,10 @@ function Trixi.integrate_via_indices(func::Func, u, 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 * + J = area_element(aux_node, equations) + integral += weights[i] * weights[j] * J * func(u, i, j, element, equations, dg, args...) - total_volume += weights[i] * weights[j] * sqrtG + total_volume += weights[i] * weights[j] * J end end @@ -65,6 +84,7 @@ 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) @@ -88,14 +108,14 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, diff = func(u_exact, equations) - func(u_numerical, equations) # 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) + J = area_element(aux_node, equations) + l2_error += diff .^ 2 * (weights[i] * weights[j] * J) # Compute Linf error as usual linf_error = @. max(linf_error, abs(diff)) # Increment total volume according to the volume element stored in aux vars - total_volume += weights[i] * weights[j] * sqrtG + total_volume += weights[i] * weights[j] * J end end @@ -104,4 +124,4 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, return l2_error, linf_error end -end # muladd +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 9807925..fabf353 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 @@ -3,7 +3,8 @@ # Convert to another set of variables using a solution_variables function function convert_variables(u, solution_variables, mesh::P4estMesh{2}, - equations::AbstractEquations{3}, dg, cache) + equations::AbstractEquations{3}, + dg, cache) (; contravariant_vectors) = cache.elements # Extract the number of solution variables to be output # (may be different than the number of conservative variables) @@ -129,7 +130,7 @@ end u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) - # Return th solution variables + # Return the solution variables return SVector(cons2prim(u_node, equations)..., relative_vorticity) end @@ -177,7 +178,9 @@ end end # Variable names for cons2prim_and_vorticity -Trixi.varnames(::typeof(cons2prim_and_vorticity), equations::ShallowWaterEquations3D) = (varnames(cons2prim, - equations)..., - "vorticity") +function Trixi.varnames(::typeof(cons2prim_and_vorticity), + equations::Union{ShallowWaterEquations3D, + AbstractCovariantEquations{2}}) + return (varnames(cons2prim, equations)..., "vorticity") end +end # @muladd diff --git a/src/callbacks_step/save_solution_covariant.jl b/src/callbacks_step/save_solution_covariant.jl index b39da3a..4cf8857 100644 --- a/src/callbacks_step/save_solution_covariant.jl +++ b/src/callbacks_step/save_solution_covariant.jl @@ -1,32 +1,73 @@ @muladd begin #! format: noindent -# Convert to another set of variables using a solution_variables function that depends on -# auxiliary variables +# Convert to another set of variables using a solution_variables function function convert_variables(u, solution_variables, mesh::P4estMesh{2}, 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, - 1), equations)) + n_vars = length(Trixi.varnames(solution_variables, equations)) + # Allocate storage for output data = Array{eltype(u)}(undef, n_vars, nnodes(dg), nnodes(dg), nelements(dg, cache)) # Loop over all nodes and convert variables, passing in auxiliary variables Trixi.@threaded for element in eachelement(dg, cache) for j in eachnode(dg), i in eachnode(dg) - u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) - aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) - u_node_converted = solution_variables(u_node, aux_node, equations) + if applicable(solution_variables, u, equations, dg, cache, i, j, element) + # The solution_variables function depends on the solution on other nodes + data_node = solution_variables(u, equations, dg, cache, i, j, element) + else + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, + element) + data_node = solution_variables(u_node, aux_node, equations) + end + for v in 1:n_vars - data[v, i, j, element] = u_node_converted[v] + data[v, i, j, element] = data_node[v] end end end return data end -end # muladd + +# Calculate the primitive variables and the relative vorticity at a given node +@inline function cons2prim_and_vorticity(u, equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, i, j, element) + (; aux_node_vars) = cache.auxiliary_variables + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + relative_vorticity = calc_vorticity_node(u, equations, dg, cache, i, j, element) + return SVector(cons2prim(u_node, aux_node, equations)..., relative_vorticity) +end + +# Calculate relative vorticity ζ = ∂₁v₂ - ∂₂v₁ for equations in covariant form +@inline function calc_vorticity_node(u, equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, i, j, element) + (; derivative_matrix) = dg.basis + (; aux_node_vars) = cache.auxiliary_variables + + dv2dxi1 = dv1dxi2 = zero(eltype(u)) + for ii in eachnode(dg) + u_node_ii = Trixi.get_node_vars(u, equations, dg, ii, j, element) + aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, element) + vcov = metric_covariant(aux_node_ii, equations) * + velocity_contravariant(u_node_ii, equations) + dv2dxi1 = dv2dxi1 + derivative_matrix[i, ii] * vcov[2] + end + + for jj in eachnode(dg) + u_node_jj = Trixi.get_node_vars(u, equations, dg, i, jj, element) + aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, element) + vcov = metric_covariant(aux_node_jj, equations) * + velocity_contravariant(u_node_jj, equations) + dv1dxi2 = dv1dxi2 + derivative_matrix[j, jj] * vcov[1] + end + + # compute the relative vorticity + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + return (dv2dxi1 - dv1dxi2) / area_element(aux_node, equations) +end +end # @muladd diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 45530be..0d9c888 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -69,4 +69,4 @@ function Trixi.max_dt(u, t, mesh::P4estMesh{2}, constant_speed::False, end return 2 / (nnodes(dg) * max_scaled_speed) end -end # muladd +end # @muladd diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index c38db84..af77ebe 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -24,21 +24,20 @@ 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} +J \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] +\left[\begin{array}{c} J 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] +\left[\begin{array}{c} J 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 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. +where $J = \lVert\vec{a}^1 \times \vec{a}^2 \rVert$ is the area element. Note that the +variable advection velocity components could alternatively be stored as auxiliary +variables, similarly to the geometric information. """ struct CovariantLinearAdvectionEquation2D{GlobalCoordinateSystem} <: AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} @@ -56,19 +55,19 @@ function Trixi.varnames(::typeof(cons2cons), ::CovariantLinearAdvectionEquation2 end # Convenience function to extract the velocity -function velocity(u, ::CovariantLinearAdvectionEquation2D) +function velocity_contravariant(u, ::CovariantLinearAdvectionEquation2D) return SVector(u[2], u[3]) end -# Convert contravariant velocity/momentum components to the global coordinate system +# Convert contravariant velocity 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) + vglo1, vglo2, vglo3 = basis_covariant(aux_vars, equations) * + velocity_contravariant(u, equations) return SVector(u[1], vglo1, vglo2, vglo3) end -# Convert velocity/momentum components in the global coordinate system to contravariant -# components +# Convert velocity 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]) @@ -89,15 +88,14 @@ end return SVector(u[1], z, z) end -# 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 required for the covariant form +# Flux as a function of the state vector u, as well as the auxiliary variables aux_vars, +# which contain the geometric information required for the covariant form @inline function Trixi.flux(u, aux_vars, orientation::Integer, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u)) - sqrtG = area_element(aux_vars, equations) - vcon = velocity(u, equations) - return SVector(sqrtG * u[1] * vcon[orientation], z, z) + J = area_element(aux_vars, equations) + vcon = velocity_contravariant(u, equations) + return SVector(J * u[1] * vcon[orientation], z, z) end # Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity @@ -107,24 +105,47 @@ end orientation_or_normal_direction, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u_ll)) - sqrtG = area_element(aux_vars_ll, equations) + J = 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 * sqrtG * λ * SVector(u_rr[1] - u_ll[1], z, z) + return -0.5f0 * J * λ * SVector(u_rr[1] - u_ll[1], z, z) end # 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) - vcon_ll = velocity(u_ll, equations) # Contravariant components on left side - vcon_rr = velocity(u_rr, equations) # Contravariant components on right side + vcon_ll = velocity_contravariant(u_ll, equations) # Contravariant components on left side + vcon_rr = velocity_contravariant(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 @inline function Trixi.max_abs_speeds(u, aux_vars, equations::CovariantLinearAdvectionEquation2D) - return abs.(velocity(u, equations)) + return abs.(velocity_contravariant(u, equations)) +end + +# If the initial velocity field is defined in Cartesian coordinates and the chosen global +# coordinate system is spherical, perform the appropriate conversion +@inline function cartesian2global(u, x, + ::CovariantLinearAdvectionEquation2D{GlobalSphericalCoordinates}) + vlon, vlat, vrad = cartesian2spherical(u[2], u[3], u[4], x) + return SVector(u[1], vlon, vlat, vrad) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is Cartesian, perform the appropriate conversion +@inline function spherical2global(u, x, + ::CovariantLinearAdvectionEquation2D{GlobalCartesianCoordinates}) + vx, vy, vz = spherical2cartesian(u[2], u[3], u[4], x) + return SVector(u[1], vx, vy, vz) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is spherical, do not convert +@inline function spherical2global(u, x, + ::CovariantLinearAdvectionEquation2D{GlobalSphericalCoordinates}) + return u end end # @muladd diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl new file mode 100644 index 0000000..1b5a010 --- /dev/null +++ b/src/equations/covariant_shallow_water.jl @@ -0,0 +1,349 @@ +@muladd begin +#! format: noindent + +@doc raw""" + CovariantShallowWaterEquations2D{GlobalCoordinateSystem} <: + AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} + +Denoting the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) by +$\nabla_b$ and summing over repeated indices, the shallow water equations can be expressed +on a two-dimensional surface in three-dimensional ambient space as +```math +\begin{aligned} +\partial_t h + \nabla_b (hv^b) &= 0,\\ +\partial_t (hv^a) + \nabla_b (hv^av^b) + gh G^{ab}\partial_b(h + b) +&= -fJ G^{ab}\varepsilon_{bc} hv^c, +\end{aligned} +``` +where $h$ is the fluid height, $v^a$ and $G^{ab}$ are the contravariant velocity and metric +tensor components, $g$ is the gravitational constant, $f$ is the Coriolis parameter, and +$J$ is the area element. Combining the inertial and gravitational terms in order to define +the momentum flux components +```math +\tau^{ab} = hv^a v^b + \frac{1}{2}G^{ab}gh^2, +``` +the covariant shallow water equations can be expressed as a system of conservation laws +with a source term: +```math +J \frac{\partial}{\partial t} +\left[\begin{array}{c} h \\ hv^1 \\ hv^2 \end{array}\right] ++ +\frac{\partial}{\partial \xi^1} +\left[\begin{array}{c} J h v^1 \\ J \tau^{11} \\ J \tau^{12} \end{array}\right] ++ +\frac{\partial}{\partial \xi^2} +\left[\begin{array}{c} J h v^2 \\ J \tau^{21} \\ J \tau^{22} \end{array}\right] += J \left[\begin{array}{c} 0 \\ +-\Gamma^1_{ab}\tau^{ab} - f J \big(G^{12}hv^1 - G^{11}hv^2\big) \\ +-\Gamma^2_{ab}\tau^{ab} - f J \big(G^{22}hv^1 - G^{21}hv^2\big) + \end{array}\right]. +``` +TrixiAtmo.jl implements standard weak-form DG methods for both the above system, as well as +entropy-stable split forms based on a novel flux-differencing discretization of the +covariant derivative. + +## References +- M. Baldauf (2020). Discontinuous Galerkin solver for the shallow-water equations in + covariant form on the sphere and the ellipsoid. Journal of Computational Physics + 410:109384. [DOI: 10.1016/j.jcp.2020.109384](https://doi.org/10.1016/j.jcp.2020.109384) +- L. Bao, R. D. Nair, and H. M. Tufo (2014). A mass and momentum flux-form high-order + discontinuous Galerkin shallow water model on the cubed-sphere. A mass and momentum + flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. + Journal of Computational Physics 271:224-243. + [DOI: 10.1016/j.jcp.2013.11.033](https://doi.org/10.1016/j.jcp.2013.11.033) + +!!! warning "Experimental implementation" + The use of entropy-stable split-form/flux-differencing formulations for covariant + equations is an experimental feature and may change in future releases. +""" +struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} <: + AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} + gravity::RealT # acceleration due to gravity + rotation_rate::RealT # rotation rate for Coriolis term + global_coordinate_system::GlobalCoordinateSystem + function CovariantShallowWaterEquations2D(gravity::RealT, + rotation_rate::RealT; + global_coordinate_system = GlobalCartesianCoordinates()) where {RealT <: + Real} + return new{typeof(global_coordinate_system), RealT}(gravity, rotation_rate) + end +end + +# Our implementation of flux-differencing formulation uses nonconservative terms, but the +# standard weak form does not. To handle both options, we have defined a dummy kernel for +# the nonconservative terms that does nothing when VolumeIntegralWeakForm is used with a +# nonconservative system. +Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True() + +# The conservative variables are the height and contravariant momentum components +function Trixi.varnames(::typeof(cons2cons), ::CovariantShallowWaterEquations2D) + return ("h", "h_vcon1", "h_vcon2") +end + +# The primitive variables are the height and contravariant velocity components +function Trixi.varnames(::typeof(cons2prim), ::CovariantShallowWaterEquations2D) + return ("h", "vcon1", "vcon2") +end + +# The change of variables contravariant2global converts the two local contravariant vector +# components u[2] and u[3] to the three global vector components specified by +# equations.global_coordinate_system (e.g. spherical or Cartesian). This transformation +# works for both primitive and conservative variables, although varnames refers +# specifically to transformations from conservative variables. +function Trixi.varnames(::typeof(contravariant2global), + ::CovariantShallowWaterEquations2D) + return ("h", "h_vglo1", "h_vglo2", "h_vglo3") +end + +# Convenience functions to extract physical variables from state vector +@inline waterheight(u, ::CovariantShallowWaterEquations2D) = u[1] +@inline velocity_contravariant(u, ::CovariantShallowWaterEquations2D) = SVector(u[2] / + u[1], + u[3] / + u[1]) +@inline momentum_contravariant(u, ::CovariantShallowWaterEquations2D) = SVector(u[2], + u[3]) + +@inline function Trixi.cons2prim(u, aux_vars, ::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + return SVector(h, h_vcon1 / h, h_vcon2 / h) +end + +@inline function Trixi.prim2cons(u, aux_vars, ::CovariantShallowWaterEquations2D) + h, vcon1, vcon2 = u + return SVector(h, h * vcon1, h * vcon2) +end + +@inline function Trixi.cons2entropy(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + h = waterheight(u, equations) + vcon = velocity_contravariant(u, equations) + vcov = metric_covariant(aux_vars, equations) * vcon + return SVector{3}(equations.gravity * h - 0.5f0 * dot(vcov, vcon), vcov[1], vcov[2]) +end + +# Convert contravariant momentum components to the global coordinate system +@inline function contravariant2global(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + h_vglo1, h_vglo2, h_vglo3 = basis_covariant(aux_vars, equations) * + momentum_contravariant(u, equations) + return SVector(waterheight(u, equations), h_vglo1, h_vglo2, h_vglo3) +end + +# Convert momentum components in the global coordinate system to contravariant components +@inline function global2contravariant(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + h_vcon1, h_vcon2 = basis_contravariant(aux_vars, equations) * + SVector(u[2], u[3], u[4]) + return SVector(u[1], h_vcon1, h_vcon2) +end + +# Entropy function (total energy per unit volume) +@inline function Trixi.entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D) + h = waterheight(u, equations) + vcon = velocity_contravariant(u, equations) + vcov = metric_covariant(aux_vars, equations) * vcon + return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2) +end + +# Flux as a function of the state vector u, as well as the auxiliary variables aux_vars, +# which contain the geometric information required for the covariant form +@inline function Trixi.flux(u, aux_vars, orientation::Integer, + equations::CovariantShallowWaterEquations2D) + # Geometric variables + Gcon = metric_contravariant(aux_vars, equations) + J = area_element(aux_vars, equations) + + # Physical variables + h = waterheight(u, equations) + h_vcon = momentum_contravariant(u, equations) + + # Compute and store the velocity and gravitational terms + vcon = h_vcon[orientation] / h + gravitational_term = 0.5f0 * equations.gravity * h^2 + + # Compute the momentum flux components in the desired orientation + momentum_flux_1 = h_vcon[1] * vcon + Gcon[1, orientation] * gravitational_term + momentum_flux_2 = h_vcon[2] * vcon + Gcon[2, orientation] * gravitational_term + + return SVector(J * h_vcon[orientation], J * momentum_flux_1, J * momentum_flux_2) +end + +# Symmetric part of entropy-conservative flux +@inline function Trixi.flux_ec(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation::Integer, + equations::CovariantShallowWaterEquations2D) + # Geometric variables + J_ll = area_element(aux_vars_ll, equations) + J_rr = area_element(aux_vars_rr, equations) + + # Physical variables + h_vcon_ll = momentum_contravariant(u_ll, equations) + h_vcon_rr = momentum_contravariant(u_rr, equations) + vcon_ll = velocity_contravariant(u_ll, equations) + vcon_rr = velocity_contravariant(u_rr, equations) + + # Scaled mass flux in conservative form + mass_flux = 0.5f0 * (J_ll * h_vcon_ll[orientation] + J_rr * h_vcon_rr[orientation]) + + # Half of scaled inertial flux in conservative form + momentum_flux = 0.25f0 * (J_ll * h_vcon_ll * vcon_ll[orientation] + + J_rr * h_vcon_rr * vcon_rr[orientation]) + + return SVector(mass_flux, momentum_flux[1], momentum_flux[2]) +end + +# Entropy-conservative flux with local Lax-Friedrichs dissipation +const flux_ec_llf = FluxPlusDissipation(flux_ec, + DissipationLocalLaxFriedrichs(max_abs_speed_naive)) + +# Non-symmetric part of entropy-conservative flux +@inline function flux_nonconservative_ec(u_ll, u_rr, aux_vars_ll, + aux_vars_rr, + orientation::Integer, + equations::CovariantShallowWaterEquations2D) + # Geometric variables + Gcov_ll = metric_covariant(aux_vars_ll, equations) + Gcov_rr = metric_covariant(aux_vars_rr, equations) + Gcon_ll = metric_contravariant(aux_vars_ll, equations) + J_ll = area_element(aux_vars_ll, equations) + J_rr = area_element(aux_vars_rr, equations) + + # Physical variables + h_vcon_ll = momentum_contravariant(u_ll, equations) + h_vcon_rr = momentum_contravariant(u_rr, equations) + vcov_ll = Gcov_ll * velocity_contravariant(u_ll, equations) + vcov_rr = Gcov_rr * velocity_contravariant(u_rr, equations) + + # Half of inertial term in non-conservative form + nonconservative_term_inertial = 0.5f0 * Gcon_ll * + (J_ll * h_vcon_ll[orientation] * vcov_rr + + J_rr * h_vcon_rr[orientation] * vcov_ll) + + # Gravity term in non-conservative form + nonconservative_term_gravitational = equations.gravity * J_ll * + Gcon_ll[:, orientation] * + waterheight(u_ll, equations) * + waterheight(u_rr, equations) + + nonconservative_term = nonconservative_term_inertial + + nonconservative_term_gravitational + + return SVector(zero(eltype(u_ll)), nonconservative_term[1], nonconservative_term[2]) +end + +# Geometric and Coriolis sources for a rotating sphere with VolumeIntegralWeakForm +@inline function source_terms_weak_form(u, x, t, aux_vars, + equations::CovariantShallowWaterEquations2D) + # Geometric variables + Gcon = metric_contravariant(aux_vars, equations) + Gamma1, Gamma2 = christoffel_symbols(aux_vars, equations) + J = area_element(aux_vars, equations) + + # Physical variables + h = waterheight(u, equations) + h_vcon = momentum_contravariant(u, equations) + v_con = velocity_contravariant(u, equations) + + # Doubly-contravariant flux tensor + momentum_flux = h_vcon * v_con' + 0.5f0 * equations.gravity * h^2 * Gcon + + # Coriolis parameter + f = 2 * equations.rotation_rate * x[3] / norm(x) # 2Ωsinθ + + # Geometric source term + s_geo = SVector(sum(Gamma1 .* momentum_flux), sum(Gamma2 .* momentum_flux)) + + # Combined source terms + source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon[1] - Gcon[1, 1] * h_vcon[2]) + source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon[1] - Gcon[2, 1] * h_vcon[2]) + + # Do not scale by Jacobian since apply_jacobian! is called before this + return SVector(zero(eltype(u)), -source_1, -source_2) +end + +# Geometric and Coriolis sources for a rotating sphere with VolumeIntegralFluxDifferencing +@inline function source_terms_ec(u, x, t, aux_vars, + equations::CovariantShallowWaterEquations2D) + # Geometric variables + Gcov = metric_covariant(aux_vars, equations) + Gcon = metric_contravariant(aux_vars, equations) + (Gamma1, Gamma2) = christoffel_symbols(aux_vars, equations) + J = area_element(aux_vars, equations) + + # Physical variables + h_vcon = momentum_contravariant(u, equations) + v_con = velocity_contravariant(u, equations) + + # Doubly-contravariant and mixed inertial flux tensors + h_vcon_vcon = h_vcon * v_con' + h_vcov_vcon = Gcov * h_vcon_vcon + + # Coriolis parameter + f = 2 * equations.rotation_rate * x[3] / norm(x) # 2Ωsinθ + + # Geometric source term + s_geo = 0.5f0 * (SVector(sum(Gamma1 .* h_vcon_vcon), sum(Gamma2 .* h_vcon_vcon)) - + Gcon * (Gamma1 * h_vcov_vcon[1, :] + Gamma2 * h_vcov_vcon[2, :])) + + # Combined source terms + source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon[1] - Gcon[1, 1] * h_vcon[2]) + source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon[1] - Gcon[2, 1] * h_vcon[2]) + + # Do not scale by Jacobian since apply_jacobian! is called before this + return SVector(zero(eltype(u)), -source_1, -source_2) +end + +# Maximum wave speed along the normal direction in reference space +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation, + equations::CovariantShallowWaterEquations2D) + # Geometric variables + Gcon_ll = metric_contravariant(aux_vars_ll, equations) + Gcon_rr = metric_contravariant(aux_vars_rr, equations) + + # Physical variables + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_ll, equations) + h_vcon_ll = momentum_contravariant(u_ll, equations) + h_vcon_rr = momentum_contravariant(u_rr, equations) + + return max(abs(h_vcon_ll[orientation] / h_ll) + + sqrt(Gcon_ll[orientation, orientation] * h_ll * equations.gravity), + abs(h_vcon_rr[orientation] / h_rr) + + sqrt(Gcon_rr[orientation, orientation] * h_rr * equations.gravity)) +end + +# Maximum wave speeds with respect to the covariant basis +@inline function Trixi.max_abs_speeds(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + vcon = velocity_contravariant(u, equations) + h = waterheight(u, equations) + Gcon = metric_contravariant(aux_vars, equations) + return abs(vcon[1]) + sqrt(Gcon[1, 1] * h * equations.gravity), + abs(vcon[2]) + sqrt(Gcon[2, 2] * h * equations.gravity) +end + +# If the initial velocity field is defined in Cartesian coordinates and the chosen global +# coordinate system is spherical, perform the appropriate conversion +@inline function cartesian2global(u, x, + ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) + h_vlon, h_vlat, h_vrad = cartesian2spherical(u[2], u[3], u[4], x) + return SVector(u[1], h_vlon, h_vlat, h_vrad) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is Cartesian, perform the appropriate conversion +@inline function spherical2global(u, x, + ::CovariantShallowWaterEquations2D{GlobalCartesianCoordinates}) + h_vx, h_vy, h_vz = spherical2cartesian(u[2], u[3], u[4], x) + return SVector(u[1], h_vx, h_vy, h_vz) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is spherical, do not convert +@inline function spherical2global(u, x, + ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) + return u +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8cc5393..3b394da 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -99,7 +99,7 @@ function transform_initial_condition(initial_condition, ::AbstractCovariantEquat return initial_condition_transformed end -# Version for standard Cartesian formulations +# Default 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) @@ -132,33 +132,71 @@ performed by [`contravariant2global`](@ref). """ function global2contravariant end -# cons2cons method which takes in unused aux_vars variable +# By default, the equations are assumed to be formulated in Cartesian coordinates. This +# function is specialized where needed. +function cartesian2global(u, x, equations::AbstractEquations) + return u +end + +# Default cons2cons and prim2cons methods which take 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. +# contravariant basis matrix, the exact area element, the upper-triangular covariant and +# contravariant metric tensor components, and the upper-triangular Christoffel symbols of +# the second kind @inline have_aux_node_vars(::AbstractCovariantEquations) = True() -@inline n_aux_node_vars(::AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT}) where {NDIMS, -NDIMS_AMBIENT} = 2 * NDIMS * NDIMS_AMBIENT + 1 + +# Add up the total number of auxiliary variables for equations in covariant form +@inline function n_aux_node_vars(::AbstractCovariantEquations{NDIMS, + NDIMS_AMBIENT}) where { + NDIMS, + NDIMS_AMBIENT + } + nvars_basis_covariant = NDIMS_AMBIENT * NDIMS + nvars_basis_contravariant = NDIMS * NDIMS_AMBIENT + nvars_area_element = 1 + nvars_metric_covariant = NDIMS * (NDIMS + 1) ÷ 2 + nvars_metric_contravariant = NDIMS * (NDIMS + 1) ÷ 2 + nvars_christoffel = NDIMS * NDIMS * (NDIMS + 1) ÷ 2 + + return nvars_basis_covariant + + nvars_basis_contravariant + + nvars_area_element + + nvars_metric_covariant + + nvars_metric_contravariant + + nvars_christoffel +end # Return auxiliary variable names for 2D covariant form @inline function auxvarnames(::AbstractCovariantEquations{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") + return ("basis_covariant[1,1]", # e₁ ⋅ a₁ + "basis_covariant[2,1]", # e₂ ⋅ a₁ + "basis_covariant[3,1]", # e₃ ⋅ a₁ + "basis_covariant[1,2]", # e₁ ⋅ a₂ + "basis_covariant[2,2]", # e₂ ⋅ a₂ + "basis_covariant[3,2]", # e₃ ⋅ a₂ + "basis_contravariant[1,1]", # e₁ ⋅ a¹ + "basis_contravariant[2,1]", # e₂ ⋅ a¹ + "basis_contravariant[3,1]", # e₃ ⋅ a¹ + "basis_contravariant[1,2]", # e₁ ⋅ a² + "basis_contravariant[2,2]", # e₂ ⋅ a² + "basis_contravariant[3,2]", # e₃ ⋅ a² + "area_element", # J = √(G₁₁G₂₂ - G₁₂G₂₁) = ||a₁ × a₂|| + "metric_covariant[1,1]", # G₁₁ + "metric_covariant[1,2]", # G₁₂ = G₂₁ + "metric_covariant[2,2]", # G₂₂ + "metric_contravariant[1,1]", # G¹¹ + "metric_contravariant[1,2]", # G¹² = G²¹ + "metric_contravariant[2,2]", # G²² + "christoffel_symbols[1][1,1]", # Γ¹₁₁ + "christoffel_symbols[1][1,2]", # Γ¹₁₂ = Γ¹₂₁ + "christoffel_symbols[2][2,2]", # Γ¹₂₂ + "christoffel_symbols[2][1,1]", # Γ²₁₁ + "christoffel_symbols[2][1,2]", # Γ²₁₂ = Γ²₂₁ + "christoffel_symbols[2][2,2]") # Γ²₂₂ end # Extract the covariant basis vectors a_i from the auxiliary variables as a matrix A, @@ -178,11 +216,29 @@ end aux_vars[11], aux_vars[12]) end -# Extract the area element √G = (det(AᵀA))^(1/2) from the auxiliary variables +# Extract the area element J = (det(AᵀA))^(1/2) from the auxiliary variables @inline function area_element(aux_vars, ::AbstractCovariantEquations{2}) return aux_vars[13] end +# Extract the covariant metric tensor components Gᵢⱼ from the auxiliary variables +@inline function metric_covariant(aux_vars, ::AbstractCovariantEquations{2}) + return SMatrix{2, 2}(aux_vars[14], aux_vars[15], + aux_vars[15], aux_vars[16]) +end + +# Extract the contravariant metric tensor components Gⁱʲ from the auxiliary variables +@inline function metric_contravariant(aux_vars, ::AbstractCovariantEquations{2}) + return SMatrix{2, 2}(aux_vars[17], aux_vars[18], + aux_vars[18], aux_vars[19]) +end + +# Extract the Christoffel symbols of the second kind Γⁱⱼₖ from the auxiliary variables +@inline function christoffel_symbols(aux_vars, ::AbstractCovariantEquations{2}) + return (SMatrix{2, 2}(aux_vars[20], aux_vars[21], aux_vars[21], aux_vars[22]), + SMatrix{2, 2}(aux_vars[23], aux_vars[24], aux_vars[24], aux_vars[25])) +end + # Numerical flux plus dissipation for abstract covariant equations as a function of the # state vectors u_ll and u_rr, as well as the auxiliary variable vectors aux_vars_ll and # aux_vars_rr, which contain the geometric information. We assume that u_ll and u_rr have @@ -212,11 +268,75 @@ end return 0.5f0 * (flux_ll + flux_rr) end +# Local Lax-Friedrichs dissipation for abstract covariant equations, where dissipation is +# applied to all conservative variables and the wave speed may depend on auxiliary variables +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, aux_vars_ll, + aux_vars_rr, + orientation_or_normal_direction, + equations::AbstractCovariantEquations) + λ = 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) * λ * (u_rr - u_ll) +end + +# Define the two-point nonconservative flux for weak form to be a no-op +@inline function flux_nonconservative_weak_form(u_ll::SVector{NVARS, RealT}, + u_rr::SVector{NVARS, RealT}, + aux_vars_ll, aux_vars_rr, + orientation_or_normal_direction, + equations::AbstractCovariantEquations{2, + NDIMS_AMBIENT, + GlobalCoordinateSystem, + NVARS}) where { + NDIMS_AMBIENT, + GlobalCoordinateSystem, + NVARS, + RealT + } + return zeros(SVector{NVARS, RealT}) +end + +# Convert a vector from a global spherical to Cartesian basis representation. A tangent +# vector will have vrad = 0. +@inline function spherical2cartesian(vlon, vlat, vrad, x) + # compute longitude and latitude + lon, lat = atan(x[2], x[1]), asin(x[3] / norm(x)) + + # compute trigonometric functions + sinlon, coslon = sincos(lon) + sinlat, coslat = sincos(lat) + + # return Cartesian components + vx = -sinlon * vlon - sinlat * coslon * vlat + coslat * coslon * vrad + vy = coslon * vlon - sinlat * sinlon * vlat + coslat * sinlon * vrad + vz = coslat * vlat + sinlat * vrad + return vx, vy, vz +end + +# Convert a vector from a global Cartesian to spherical basis representation. A tangent +# vector will have vrad = 0. +@inline function cartesian2spherical(vx, vy, vz, x) + # compute longitude and latitude + lon, lat = atan(x[2], x[1]), asin(x[3] / norm(x)) + + # compute trigonometric functions + sinlon, coslon = sincos(lon) + sinlat, coslat = sincos(lat) + + # return spherical components + vlon = -sinlon * vx + coslon * vy + vlat = -sinlat * coslon * vx - sinlat * sinlon * vy + coslat * vz + vrad = coslat * coslon * vx + coslat * sinlon * vy + sinlat * vz + + return vlon, vlat, vrad +end + abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end -include("reference_data.jl") include("covariant_advection.jl") +include("covariant_shallow_water.jl") include("compressible_moist_euler_2d_lucas.jl") include("shallow_water_3d.jl") +include("reference_data.jl") end # @muladd diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 6865223..1543d7c 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -1,8 +1,8 @@ @muladd begin #! format: noindent -# Physical constants -const EARTH_RADIUS = 6.37122e6 # m +# Physical constants in SI units (reference values from the Williamson et al. test suite) +const EARTH_RADIUS = 6.37122e6 # m const EARTH_GRAVITATIONAL_ACCELERATION = 9.80616 # m/s² const EARTH_ROTATION_RATE = 7.292e-5 # rad/s const SECONDS_PER_DAY = 8.64e4 @@ -11,7 +11,7 @@ const SECONDS_PER_DAY = 8.64e4 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 $6. +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. Denoting the Euclidean norm as $\lVert \cdot \rVert$, the initial height field is given by ```math @@ -33,9 +33,8 @@ This problem is adapted from Case 1 of the test suite described in the following 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) """ -@inline function initial_condition_gaussian(x, t, ::AbstractEquations) +@inline function initial_condition_gaussian(x, t, equations) 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 @@ -63,48 +62,119 @@ This problem is adapted from Case 1 of the test suite described in the following # 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) + # Convert primitive variables from Cartesian coordinates to the chosen global + # coordinate system, which depends on the equation type + return cartesian2global(SVector(h, vx, vy, vz, 0.0f0), x, equations) 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) +@doc raw""" + initial_condition_geostrophic_balance(x, t, equations) +Steady geostrophic balance for the spherical shallow water equations, corresponding to a +purely zonal velocity field given as a function of the latitude $\theta$ by +$v_\lambda(\theta) = v_0 \cos\theta$, where we define $v_0 = 2\pi a / (12 \ \mathrm{days})$ +in terms of the Earth's radius $a = 6.37122 \times 10^3\ \mathrm{m}$. The height field +then varies with the latitude as +```math +h(\theta) = \frac{1}{g} +\Big(gh_0 - \Big(a \Omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), +``` +where $gh_0 = 2.94 \times 10^4 \ \mathrm{m}^2/\mathrm{s}^2$, +$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and +$\Omega = 7.292 \times 10^{-5}\ \mathrm{s}^{-1}$. This problem corresponds to Case 2 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) +""" +@inline function initial_condition_geostrophic_balance(x, t, equations) + 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 + lat = asin(x[3] / a) - # 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)) + # compute zonal and meridional components of the velocity + v_0 = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY) + vlon, vlat = v_0 * cos(lat), zero(eltype(x)) - # 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 geopotential height + h = 1 / EARTH_GRAVITATIONAL_ACCELERATION * + (2.94f4 - (a * EARTH_ROTATION_RATE * v_0 + 0.5f0 * v_0^2) * (sin(lat))^2) - # 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)) + # Convert primitive variables from spherical coordinates to the chosen global + # coordinate system, which depends on the equation type + return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations) +end - # get zonal and meridional components of the velocity - 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) +@doc raw""" + initial_condition_rossby_haurwitz(x, t, equations) - # 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) +Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and +meridional velocity components are given, respectively, as functions of the longitude +$\lambda$ and latitude $\theta$ by +```math +\begin{aligned} +v_\lambda(\lambda,\theta) &= a \omega \cos \theta+a K \cos ^{R-1} \theta +\left(R \sin ^2 \theta-\cos ^2 \theta\right) \cos (R \lambda),\\ +v_\theta(\lambda,\theta) &= -a K R \cos ^{R-1} \theta \sin \theta \sin (R \lambda), +\end{aligned} +``` +where $\omega = K = 7.848 \times 10^{-6} \ \mathrm{s}^{-1}$ and $R = 4$ are given +constants, and $a = 6.37122 \times 10^3\ \mathrm{m}$ is the Earth's radius. Taking +$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\Omega = 7.292 \times 10^{-5} \ \mathrm{s}^{-1}$, +and $h_0 = 8000 \ \mathrm{m}$ and defining the functions +```math +\begin{aligned} +A(\theta) &= \frac{\omega}{2}(2 \Omega+\omega) \cos^2 \theta + +\frac{1}{4} K^2 \cos^{2 R} \theta\Big((R+1) \cos^2\theta +\left(2 R^2-R-2\right) - +\big(2 R^2 / \cos^2 \theta\big) \Big), \\ +B(\theta) &= \frac{2(\Omega+\omega) K}{(R+1)(R+2)} \cos ^R \theta\big((R^2+2 R+2) - +(R+1)^2 \cos^2 \theta\big), \\ +C(\theta) &= \frac{1}{4} K^2 \cos^{2 R} \theta\big((R+1) \cos^2 \theta-(R+2)\big), +\end{aligned} +``` +the initial height field is given by +```math +h(\lambda,\theta) = h_0 + +\frac{a^2}{g}\Big(A(\theta) + B(\theta)\cos(R\lambda) + C(\theta)\cos(2R\lambda) \Big). +``` +This problem corresponds to Case 6 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) +""" +@inline function initial_condition_rossby_haurwitz(x, t, equations) + RealT = eltype(x) + a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) + lon = atan(x[2], x[1]) + lat = asin(x[3] / a) + + h_0 = 8.0f3 + omega = 7.848f-6 + K = 7.848f-6 + R = 4.0f0 + + A = 0.5f0 * omega * (2 * EARTH_ROTATION_RATE + omega) * (cos(lat))^2 + + 0.25f0 * K^2 * (cos(lat))^(2 * R) * + ((R + 1) * (cos(lat))^2 + (2 * R^2 - R - 2) - 2 * R^2 / ((cos(lat))^2)) + B = 2 * (EARTH_ROTATION_RATE + omega) * K / ((R + 1) * (R + 2)) * (cos(lat))^R * + ((R^2 + 2R + 2) - (R + 1)^2 * (cos(lat))^2) + C = 0.25f0 * K^2 * (cos(lat))^(2 * R) * ((R + 1) * (cos(lat))^2 - (R + 2)) + + # compute geopotential height + h = h_0 + + (1 / EARTH_GRAVITATIONAL_ACCELERATION) * + (a^2 * A + a^2 * B * cos(R * lon) + a^2 * C * cos(2 * R * lon)) + + # compute zonal and meridional components of the velocity + vlon = a * omega * cos(lat) + + a * K * (cos(lat))^(R - 1) * (R * (sin(lat))^2 - + (cos(lat))^2) * cos(R * lon) + vlat = -a * K * R * (cos(lat))^(R - 1) * sin(lat) * sin(R * lon) + + # Convert primitive variables from spherical coordinates to the chosen global + # coordinate system, which depends on the equation type + return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations) end -end # muladd +end # @muladd diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl index 930af2a..c008b22 100644 --- a/src/equations/shallow_water_3d.jl +++ b/src/equations/shallow_water_3d.jl @@ -380,4 +380,13 @@ end @inline function energy_internal(cons, equations::ShallowWaterEquations3D) return energy_total(cons, equations) - energy_kinetic(cons, equations) end + +# Convert spherical velocity components to global Cartesian components +@inline function spherical2global(prim, x, ::ShallowWaterEquations3D) + h, vlat, vlon, vrad, b = prim # primitive variables + + v1, v2, v3 = spherical2cartesian(vlat, vlon, vrad, x) + + return SVector(h, v1, v2, v3, b) +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 07d595d..59ce7d9 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 @@ -163,9 +163,6 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, (; 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) @@ -187,22 +184,39 @@ 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) - # 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) + # Covariant basis in the desired global coordinate system as columns of a matrix + basis_covariant = calc_basis_covariant(v1, v2, v3, v4, + dg.basis.nodes[i], dg.basis.nodes[j], + radius, + equations.global_coordinate_system) + aux_node_vars[1:6, i, j, element] = SVector(basis_covariant) - G = A' * A # Covariant metric tensor (not used for now) + # Metric tensor + metric_covariant = basis_covariant' * basis_covariant + metric_contravariant = inv(basis_covariant' * basis_covariant) - # Contravariant basis - aux_node_vars[(NDIMS * NDIMS_AMBIENT + 1):(2 * NDIMS * NDIMS_AMBIENT), - i, j, element] = SVector(inv(G) * A') + # Contravariant basis vectors as rows of a matrix + basis_contravariant = metric_contravariant * basis_covariant' + aux_node_vars[7:12, i, j, element] = SVector(basis_contravariant) # Area element - aux_node_vars[n_aux_node_vars(equations), i, j, element] = sqrt(det(G)) + aux_node_vars[13, i, j, element] = sqrt(det(metric_covariant)) + + # Covariant metric tensor components + aux_node_vars[14:16, i, j, element] = SVector(metric_covariant[1, 1], + metric_covariant[1, 2], + metric_covariant[2, 2]) + + # Contravariant metric tensor components + aux_node_vars[17:19, i, j, element] = SVector(metric_contravariant[1, 1], + metric_contravariant[1, 2], + metric_contravariant[2, 2]) end + + # Christoffel symbols of the second kind (aux_node_vars[20:25, :, :, element]) + calc_christoffel_symbols!(aux_node_vars, mesh, equations, dg, element) end + return nothing end @@ -273,4 +287,67 @@ end # 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 + +# Calculate Christoffel symbols approximately using the collocation derivative +function calc_christoffel_symbols!(aux_node_vars, mesh::P4estMesh{2, 3}, + equations::AbstractCovariantEquations{2, 3}, dg, + element) + (; derivative_matrix) = dg.basis + + for j in eachnode(dg), i in eachnode(dg) + + # Numerically differentiate covariant metric components with respect to ξ¹ + dG11dxi1 = zero(eltype(aux_node_vars)) + dG12dxi1 = zero(eltype(aux_node_vars)) + dG22dxi1 = zero(eltype(aux_node_vars)) + for ii in eachnode(dg) + aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, + element) + Gcov_ii = metric_covariant(aux_node_ii, equations) + dG11dxi1 = dG11dxi1 + derivative_matrix[i, ii] * Gcov_ii[1, 1] + dG12dxi1 = dG12dxi1 + derivative_matrix[i, ii] * Gcov_ii[1, 2] + dG22dxi1 = dG22dxi1 + derivative_matrix[i, ii] * Gcov_ii[2, 2] + end + + # Numerically differentiate covariant metric components with respect to ξ² + dG11dxi2 = zero(eltype(aux_node_vars)) + dG12dxi2 = zero(eltype(aux_node_vars)) + dG22dxi2 = zero(eltype(aux_node_vars)) + for jj in eachnode(dg) + aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, + element) + Gcov_jj = metric_covariant(aux_node_jj, equations) + dG11dxi2 = dG11dxi2 + derivative_matrix[j, jj] * Gcov_jj[1, 1] + dG12dxi2 = dG12dxi2 + derivative_matrix[j, jj] * Gcov_jj[1, 2] + dG22dxi2 = dG22dxi2 + derivative_matrix[j, jj] * Gcov_jj[2, 2] + end + + # Compute Christoffel symbols of the first kind + christoffel_firstkind_1 = SMatrix{2, 2}(0.5f0 * dG11dxi1, + 0.5f0 * dG11dxi2, + 0.5f0 * dG11dxi2, + dG12dxi2 - 0.5f0 * dG22dxi1) + christoffel_firstkind_2 = SMatrix{2, 2}(dG12dxi1 - 0.5f0 * dG11dxi2, + 0.5f0 * dG22dxi1, + 0.5f0 * dG22dxi1, + 0.5f0 * dG22dxi2) + + # Raise indices to get Christoffel symbols of the second kind + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + Gcon = metric_contravariant(aux_node, equations) + aux_node_vars[20, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 1] + + Gcon[1, 2] * christoffel_firstkind_2[1, 1] + aux_node_vars[21, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 2] + + Gcon[1, 2] * christoffel_firstkind_2[1, 2] + aux_node_vars[22, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[2, 2] + + Gcon[1, 2] * christoffel_firstkind_2[2, 2] + + aux_node_vars[23, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[1, 1] + + Gcon[2, 2] * christoffel_firstkind_2[1, 1] + aux_node_vars[24, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[1, 2] + + Gcon[2, 2] * christoffel_firstkind_2[1, 2] + aux_node_vars[25, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[2, 2] + + Gcon[2, 2] * christoffel_firstkind_2[2, 2] + end +end +end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl index 1b494ef..c7028e5 100644 --- a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl @@ -152,4 +152,4 @@ function calc_sources_2d_manifold_in_3d!(du, u, t, source_terms, return nothing end -end # muladd +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 ec18e04..81fa008 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 @@ -113,8 +113,67 @@ end return nothing end -# Pointwise interface flux, transforming the contravariant prognostic variables into the -# local coordinate system +# Since the standard weak form for CovariantShallowWaterEquations2D has no nonconservative +# terms, but the equation type has have_nonconservative_terms = True(), we must define +# an appropriate kernel that just calls the standard weak form +@inline function Trixi.weak_form_kernel!(du, u, element, mesh::P4estMesh{2}, + nonconservative_terms::True, + equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, alpha = true) + Trixi.weak_form_kernel!(du, u, element, mesh, False(), equations, dg, cache, alpha) +end + +# Non-conservative flux differencing kernel which uses contravariant flux components, +# passing the geometric information contained in the auxiliary variables to the flux +# function +@inline function Trixi.flux_differencing_kernel!(du, u, element, mesh::P4estMesh{2}, + nonconservative_terms::True, + equations::AbstractCovariantEquations{2}, + volume_flux, dg::DGSEM, cache, + alpha = true) + (; derivative_split) = dg.basis + (; aux_node_vars) = cache.auxiliary_variables + symmetric_flux, nonconservative_flux = volume_flux + + # Apply the symmetric flux as usual + Trixi.flux_differencing_kernel!(du, u, element, mesh, False(), equations, + symmetric_flux, dg, cache, alpha) + + for j in eachnode(dg), i in eachnode(dg) + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + + # ξ¹ direction + integral_contribution = zero(u_node) + for ii in eachnode(dg) + u_node_ii = Trixi.get_node_vars(u, equations, dg, ii, j, element) + aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, + element) + flux1 = nonconservative_flux(u_node, u_node_ii, aux_node, aux_node_ii, 1, + equations) + integral_contribution = integral_contribution + + derivative_split[i, ii] * flux1 + end + + # ξ² direction + for jj in eachnode(dg) + u_node_jj = Trixi.get_node_vars(u, equations, dg, i, jj, element) + aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, + element) + flux2 = nonconservative_flux(u_node, u_node_jj, aux_node, aux_node_jj, 2, + equations) + integral_contribution = integral_contribution + + derivative_split[j, jj] * flux2 + end + + Trixi.multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution, + equations, dg, i, j, element) + end + + return nothing +end + +# Pointwise interface flux in local coordinates for problems without nonconservative terms @inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, nonconservative_terms::False, equations::AbstractCovariantEquations{2}, @@ -175,6 +234,103 @@ end end end +# Pointwise interface flux in local coordinates for problems with nonconservative terms +@inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, + nonconservative_terms::True, + equations::AbstractCovariantEquations{2}, + surface_integral, dg::DG, cache, + interface_index, normal_direction, + primary_node_index, + primary_direction_index, + primary_element_index, + secondary_node_index, + secondary_direction_index, + secondary_element_index) + (; u) = cache.interfaces + (; aux_surface_node_vars) = cache.auxiliary_variables + surface_flux, nonconservative_flux = surface_integral.surface_flux + + # Get surface values for solution and auxiliary variables + 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(aux_surface_node_vars, + equations, + dg, primary_node_index, + interface_index) + + # Compute flux in the primary element's coordinate system + u_rr_spherical = contravariant2global(u_rr, aux_vars_rr, equations) + u_rr_transformed_to_ll = global2contravariant(u_rr_spherical, aux_vars_ll, + equations) + primary_orientation = (primary_direction_index + 1) >> 1 + if isodd(primary_direction_index) + flux_primary = -(surface_flux(u_rr_transformed_to_ll, u_ll, aux_vars_ll, + aux_vars_ll, primary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_rr_transformed_to_ll, u_ll, + aux_vars_ll, aux_vars_ll, + primary_orientation, equations)) + else + flux_primary = surface_flux(u_ll, u_rr_transformed_to_ll, aux_vars_ll, + aux_vars_ll, primary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_ll, u_rr_transformed_to_ll, + aux_vars_ll, aux_vars_ll, + primary_orientation, equations) + end + + # Compute flux in the secondary element's coordinate system + u_ll_spherical = contravariant2global(u_ll, aux_vars_ll, equations) + u_ll_transformed_to_rr = global2contravariant(u_ll_spherical, aux_vars_rr, + equations) + secondary_orientation = (secondary_direction_index + 1) >> 1 + if isodd(secondary_direction_index) + flux_secondary = -(surface_flux(u_ll_transformed_to_rr, u_rr, aux_vars_rr, + aux_vars_rr, secondary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_ll_transformed_to_rr, u_rr, + aux_vars_rr, aux_vars_rr, + secondary_orientation, equations)) + else + flux_secondary = surface_flux(u_rr, u_ll_transformed_to_rr, aux_vars_rr, + aux_vars_rr, secondary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_rr, u_ll_transformed_to_rr, + aux_vars_rr, aux_vars_rr, + secondary_orientation, equations) + end + + # Update the surface flux values on both sides of the interface + for v in eachvariable(equations) + surface_flux_values[v, primary_node_index, primary_direction_index, + primary_element_index] = flux_primary[v] + surface_flux_values[v, secondary_node_index, secondary_direction_index, + secondary_element_index] = flux_secondary[v] + end +end + +# This function passes the auxiliary variables into the source term +function Trixi.calc_sources!(du, u, t, source_terms, + equations::AbstractCovariantEquations{2}, dg::DG, cache) + (; aux_node_vars) = cache.auxiliary_variables + + Trixi.@threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_local = Trixi.get_node_vars(u, equations, dg, i, j, element) + x_local = Trixi.get_node_coords(cache.elements.node_coordinates, equations, + dg, + i, j, element) + aux_local = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + du_local = source_terms(u_local, x_local, t, aux_local, equations) + Trixi.add_to_node_vars!(du, du_local, equations, dg, i, j, element) + end + end + + return nothing +end + +# Version of calc_sources! specialized for covariant formulation with no source term +function Trixi.calc_sources!(du, u, t, source_terms::Nothing, + equations::AbstractCovariantEquations{2}, dg::DG, cache) + return nothing +end + # Apply the exact Jacobian stored in auxiliary variables function Trixi.apply_jacobian!(du, mesh::P4estMesh{2}, equations::AbstractCovariantEquations{2}, diff --git a/test/runtests.jl b/test/runtests.jl index 38d3e81..8e65de2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,4 +24,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "shallow_water_3d" include("test_3d_shallow_water.jl") end + + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "shallow_water_2d_covariant" + include("test_2d_shallow_water_covariant.jl") + end end diff --git a/test/test_2d_shallow_water_covariant.jl b/test/test_2d_shallow_water_covariant.jl new file mode 100644 index 0000000..d5b1493 --- /dev/null +++ b/test/test_2d_shallow_water_covariant.jl @@ -0,0 +1,46 @@ +module TestSphericalAdvection + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") + +@trixiatmo_testset "elixir_shallowwater_covariant_geostrophic_balance" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_covariant_geostrophic_balance.jl"), + l2=[0.3065297703718947, + 0.00019984441967990253, + 8.7677772893388e-5], + linf=[ + 1.4786699241212773, + 0.0013754817113900697, + 0.0007564212354783176 + ], tspan=(0.0, 1.0 * SECONDS_PER_DAY)) + # 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 + +@trixiatmo_testset "elixir_shallowwater_covariant_rossby_haurwitz_EC" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_covariant_rossby_haurwitz_EC.jl"), + l2=[265.9858131601718, 0.1769083525032968, 0.2538332624036849], + linf=[579.0744773821989, 0.5054767269383715, 0.5628603103981238], + tspan=(0.0, 1.0 * SECONDS_PER_DAY)) + # 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