From 40800a6da262152c9c9e76689509b98d5077a137 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Mon, 25 Nov 2024 10:31:07 +0100 Subject: [PATCH] Covariant form as variable-coefficient problem using auxiliary variables (#47) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * linear advection on the sphere * formatter * very preliminary implementation of covariant form with p4est * removed untracked * removed allocations, now 100x faster * add DS_Store to gitignore * add DS_Store to gitignore * made equations actually 2D * integrated upstream changes * analysis callback for covariant form * cleanup * separate Andres's spherical shell containers from main Trixi * Added containers and 2D p4est mesh from my spherical shell implementation in Trixi.jl (cherry-picked from f45378e) Co-authored-by: Tristan Montoya * Added element container with PtrArray for performance and modified the structure * format * Fixed bug in the definition of init_elements and added nelements(). eltype() and ndims() functions for new struct * integrated Andres's new container type with covariant solver * added flux-differencing kernel for covariant form * integrated Andrés' PR with covariant solver * add tests for spherical advection in Cartesian coords * revert change to moist bubble case * revert change to moist bubble case * covariant weak form with tests * removed shallow water from this branch * reorganize a bit * hard-code dimension of node coordinates to avoid type instability * metric terms for shell of radius != 1 * test for covariant form is now on earth scale * fix name of cartesian test * moved i,j,element,cache into flux signatures for covariant form * added test for flux differencing covariant form * transformation between local coordinate systems now uses latitude-longitude as intermediate state, not Cartesian * format * make compatible with my PR trixi-framework/Trixi.jl#2068 * separate containers to store geometric information specific to covariant form * add element-local mapping from Guba et al. * add reference to Guba et al. to docstring * new mapping works with cartesian solver solver in covariant_form branch * exact metrics for element-local mapping * streamline tests and improve readability * rename default DGSEM mapping * didn't need new analysis cache * improve comments * refactor * remove DiffEqCallbacks for now * add better comments * improve comments * better documentation * add NDIMS_AMBIENT parameter to equations * fix docstring format * make link in docs work * fix link in docs * documentation for the covariant solver * added surface central flux and test * fix trixi_compat * specialize compute_coefficients! to dimension 2 * format * fix spelling * remove initial_condition parameter from rhs * no need to use Float32 for physical constant defaults * incorporate downstream changes * format * standardize test cases * move covariant form max_dt method to callbacks_step * Fix comments * generalize interface between equations and cache * relax relative tolerance on cartesian advection test) * improve docs * streamline initial condition definition and improve docs * tests for cartesian form now use standard and element-local mapping approaches * first attempt at aux vars * works with allocations * fix allocations * clean up docs * format and update comments * move auxvars into own cache field * improve comments * format comments * revert back to without basis matrix inverse * add warning for P4estMeshCubedSphere and rename volume->area element * rename function transforming initial condition * fix elixir and docs * poldeg(solver) -> Trixi.polydeg(solver) * reformat ugly line breaks * modularize covariant basis calculation for easier testing * remove duplicated cons2entropy function --------- Co-authored-by: Andrés Rueda-Ramírez --- Project.toml | 6 +- docs/Project.toml | 1 - docs/make.jl | 9 +- ...wwater_cubed_sphere_shell_EC_correction.jl | 8 +- ...wwater_cubed_sphere_shell_EC_projection.jl | 8 +- ...llowwater_cubed_sphere_shell_advection.jl} | 82 ++---- ...hallowwater_cubed_sphere_shell_standard.jl | 8 +- .../elixir_spherical_advection_covariant.jl | 68 +++++ src/TrixiAtmo.jl | 18 +- src/auxiliary/auxiliary.jl | 1 + src/callbacks_step/analysis_covariant.jl | 71 ++++++ src/callbacks_step/callbacks_step.jl | 2 + src/callbacks_step/stepsize_dg2d.jl | 40 ++- src/equations/covariant_advection.jl | 109 ++++++++ src/equations/equations.jl | 118 +++++++++ src/equations/reference_data.jl | 78 ++++++ src/equations/shallow_water_3d.jl | 60 ++++- src/meshes/p4est_cubed_sphere.jl | 18 +- src/semidiscretization/semidiscretization.jl | 9 +- ...retization_hyperbolic_2d_manifold_in_3d.jl | 4 +- ...containers_2d_manifold_in_3d_cartesian.jl} | 50 ++-- .../containers_2d_manifold_in_3d_covariant.jl | 237 ++++++++++++++++++ src/solvers/dgsem_p4est/dg.jl | 33 ++- .../dg_2d_manifold_in_3d_cartesian.jl | 23 +- .../dg_2d_manifold_in_3d_covariant.jl | 195 ++++++++++++++ src/solvers/solvers.jl | 2 - test/Project.toml | 2 +- test/test_spherical_advection.jl | 103 ++++++-- 28 files changed, 1197 insertions(+), 166 deletions(-) rename examples/{elixir_euler_spherical_advection_cartesian.jl => elixir_shallowwater_cubed_sphere_shell_advection.jl} (51%) create mode 100644 examples/elixir_spherical_advection_covariant.jl create mode 100644 src/callbacks_step/analysis_covariant.jl create mode 100644 src/callbacks_step/callbacks_step.jl create mode 100644 src/equations/covariant_advection.jl create mode 100644 src/equations/reference_data.jl rename src/solvers/dgsem_p4est/{containers_2d_manifold_in_3d.jl => containers_2d_manifold_in_3d_cartesian.jl} (93%) create mode 100644 src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl create mode 100644 src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl diff --git a/Project.toml b/Project.toml index 8b17e5e..1653fe1 100644 --- a/Project.toml +++ b/Project.toml @@ -5,8 +5,8 @@ version = "0.1.0-DEV" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" @@ -17,13 +17,13 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] LinearAlgebra = "1" -MuladdMacro = "0.2.2" LoopVectorization = "0.12.118" +MuladdMacro = "0.2.2" NLsolve = "4.5.1" Reexport = "1.0" Static = "0.8, 1" StaticArrayInterface = "1.5.1" StaticArrays = "1" StrideArrays = "0.1.28" -Trixi = "0.7, 0.8, 0.9" +Trixi = "0.9" julia = "1.9" diff --git a/docs/Project.toml b/docs/Project.toml index 3302727..342f407 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,7 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" -TrixiAtmo = "c9ed1054-d170-44a9-8ee2-d5566f5d1389" [compat] Documenter = "1.3" diff --git a/docs/make.jl b/docs/make.jl index 05ce510..6689a18 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,14 @@ -using TrixiAtmo using Documenter using DocumenterInterLinks +# Fix for https://github.com/trixi-framework/Trixi.jl/issues/668 +if (get(ENV, "CI", nothing) != "true") && + (get(ENV, "TRIXI_DOC_DEFAULT_ENVIRONMENT", nothing) != "true") + push!(LOAD_PATH, dirname(@__DIR__)) +end + +using TrixiAtmo + # Provide external links to the Trixi.jl docs (project root and inventory file) links = InterLinks("Trixi" => ("https://trixi-framework.github.io/Trixi.jl/stable/", "https://trixi-framework.github.io/Trixi.jl/stable/objects.inv")) diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl index 9f12e2c..956c7ac 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl @@ -45,12 +45,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio v0 = 1.0 # Velocity at the "equator" alpha = 0.0 #pi / 4 v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) - v_lat = -v0 * sin(phi) * sin(alpha) + vlat = -v0 * sin(phi) * sin(alpha) # Transform to Cartesian coordinate system - v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long - v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long - v3 = sin(colat) * v_lat + v1 = -cos(colat) * cos(phi) * vlat - sin(phi) * v_long + v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long + v3 = sin(colat) * vlat return prim2cons(SVector(rho, v1, v2, v3, 0), equations) end diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl index 5b43e82..21934c5 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl @@ -47,12 +47,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio v0 = 1.0 # Velocity at the "equator" alpha = 0.0 #pi / 4 v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) - v_lat = -v0 * sin(phi) * sin(alpha) + vlat = -v0 * sin(phi) * sin(alpha) # Transform to Cartesian coordinate system - v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long - v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long - v3 = sin(colat) * v_lat + v1 = -cos(colat) * cos(phi) * vlat - sin(phi) * v_long + v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long + v3 = sin(colat) * vlat return prim2cons(SVector(rho, v1, v2, v3, 0), equations) end diff --git a/examples/elixir_euler_spherical_advection_cartesian.jl b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl similarity index 51% rename from examples/elixir_euler_spherical_advection_cartesian.jl rename to examples/elixir_shallowwater_cubed_sphere_shell_advection.jl index 6e82cb0..0749c83 100644 --- a/examples/elixir_euler_spherical_advection_cartesian.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl @@ -6,58 +6,19 @@ using TrixiAtmo ############################################################################### # semidiscretization of the linear advection equation -# We use the Euler equations structure but modify the rhs! function to convert it to a -# variable-coefficient advection equation -equations = CompressibleEulerEquations3D(1.4) +cells_per_dimension = 5 +initial_condition = initial_condition_gaussian + +# We use the ShallowWaterEquations3D equations structure but modify the rhs! function to +# convert it to a variable-coefficient advection equation +equations = ShallowWaterEquations3D(gravity_constant = 0.0) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -polydeg = 3 -solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) - -# Initial condition for a Gaussian density profile with constant pressure -# and the velocity of a rotating solid body -function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEquations3D) - # Gaussian density - rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) - # Constant pressure - p = 1.0 - - # Spherical coordinates for the point x - if sign(x[2]) == 0.0 - signy = 1.0 - else - signy = sign(x[2]) - end - # Co-latitude - colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) - # Latitude (auxiliary variable) - lat = -colat + 0.5 * pi - # Longitude - r_xy = sqrt(x[1]^2 + x[2]^2) - if r_xy == 0.0 - phi = pi / 2 - else - phi = signy * acos(x[1] / r_xy) - end - - # Compute the velocity of a rotating solid body - # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system) - v0 = 1.0 # Velocity at the "equator" - alpha = 0.0 #pi / 4 - v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) - v_lat = -v0 * sin(phi) * sin(alpha) - - # Transform to Cartesian coordinate system - v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long - v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long - v3 = sin(colat) * v_lat - - return prim2cons(SVector(rho, v1, v2, v3, p), equations) -end +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) -# Source term function to transform the Euler equations into the linear advection equations with variable advection velocity +# Source term function to transform the Euler equations into a linear advection equation with variable advection velocity function source_terms_convert_to_linear_advection(u, du, x, t, - equations::CompressibleEulerEquations3D, + equations::ShallowWaterEquations3D, normal_direction) v1 = u[2] / u[1] v2 = u[3] / u[1] @@ -66,29 +27,27 @@ function source_terms_convert_to_linear_advection(u, du, x, t, s2 = du[1] * v1 - du[2] s3 = du[1] * v2 - du[3] s4 = du[1] * v3 - du[4] - s5 = 0.5 * (s2 * v1 + s3 * v2 + s4 * v3) - du[5] - return SVector(0.0, s2, s3, s4, s5) + return SVector(0.0, s2, s3, s4, 0.0) end -initial_condition = initial_condition_advection_sphere - -element_local_mapping = false - -mesh = P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, +# Create a 2D cubed sphere mesh the size of the Earth +mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = 3, initial_refinement_level = 0, - element_local_mapping = element_local_mapping) + 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) # A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver, source_terms = source_terms_convert_to_linear_advection) ############################################################################### # ODE solvers, callbacks etc. -# Create ODE problem with time span from 0.0 to π -tspan = (0.0, pi) -ode = semidiscretize(semi, tspan) +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 @@ -97,8 +56,7 @@ 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,), - extra_analysis_integrals = (Trixi.density,)) + extra_analysis_errors = (:conservation_error,)) # The SaveSolutionCallback allows to save the solution to a file in regular intervals save_solution = SaveSolutionCallback(interval = 10, diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl index bd03404..ba1c3bf 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl @@ -41,12 +41,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio v0 = 1.0 # Velocity at the "equator" alpha = 0.0 #pi / 4 v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) - v_lat = -v0 * sin(phi) * sin(alpha) + vlat = -v0 * sin(phi) * sin(alpha) # Transform to Cartesian coordinate system - v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long - v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long - v3 = sin(colat) * v_lat + v1 = -cos(colat) * cos(phi) * vlat - sin(phi) * v_long + v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long + v3 = sin(colat) * vlat return prim2cons(SVector(rho, v1, v2, v3, 0), equations) end diff --git a/examples/elixir_spherical_advection_covariant.jl b/examples/elixir_spherical_advection_covariant.jl new file mode 100644 index 0000000..6ad2115 --- /dev/null +++ b/examples/elixir_spherical_advection_covariant.jl @@ -0,0 +1,68 @@ +############################################################################### +# DGSEM for the linear advection equation on the cubed sphere +############################################################################### + +using OrdinaryDiffEq, Trixi, TrixiAtmo + +############################################################################### +# Spatial discretization + +cells_per_dimension = 5 +initial_condition = initial_condition_gaussian + +equations = CovariantLinearAdvectionEquation2D() + +# 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, +# initial_refinement_level = 0, and element_local_mapping = true. +mesh = P4estMeshCubedSphere2D(cells_per_dimension, 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) + +# 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 = cons2cons) + +# 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 cf5f7f0..c6004e6 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -8,15 +8,17 @@ See also: [trixi-framework/TrixiAtmo.jl](https://github.com/trixi-framework/Trix """ module TrixiAtmo +using Reexport: @reexport using Trixi using MuladdMacro: @muladd using Static: True, False using StrideArrays: PtrArray using StaticArrayInterface: static_size -using LinearAlgebra: norm +using LinearAlgebra: norm, dot using Reexport: @reexport using LoopVectorization: @turbo -@reexport using StaticArrays: SVector + +@reexport using StaticArrays: SVector, SMatrix include("auxiliary/auxiliary.jl") include("equations/equations.jl") @@ -24,9 +26,10 @@ include("meshes/meshes.jl") include("semidiscretization/semidiscretization.jl") include("solvers/solvers.jl") include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") -include("callbacks_step/stepsize_dg2d.jl") +include("callbacks_step/callbacks_step.jl") -export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D +export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, + CovariantLinearAdvectionEquation2D export flux_chandrashekar, flux_LMARS @@ -34,6 +37,11 @@ export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_int lake_at_rest_error, source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier! export P4estCubedSphere2D, MetricTermsCrossProduct, MetricTermsInvariantCurl -export examples_dir +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 examples_dir end # module TrixiAtmo diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 316a687..ae05020 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -7,6 +7,7 @@ read-only and should *not* be modified. To find out which files are available, use, e.g., `readdir`: # Examples + ```@example readdir(examples_dir()) ``` diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl new file mode 100644 index 0000000..3c8d3fb --- /dev/null +++ b/src/callbacks_step/analysis_covariant.jl @@ -0,0 +1,71 @@ +@muladd begin +#! format: noindent + +# Entropy time derivative which uses auxiliary variables +function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t, + mesh::P4estMesh{2}, + equations::AbstractCovariantEquations{2}, dg::DG, cache) + (; aux_node_vars) = cache.auxiliary_variables + + # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ + Trixi.integrate_via_indices(u, mesh, equations, dg, cache, + du) do u, i, j, element, equations, dg, du_node + # Get auxiliary variables, solution variables, and time derivative at given node + a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + du_node = Trixi.get_node_vars(du, equations, dg, i, j, element) + + # compute ∂S/∂u ⋅ ∂u/∂t, where the entropy variables ∂S/∂u depend on the solution + # and auxiliary variables + dot(cons2entropy(u_node, a_node, equations), du_node) + end +end + +# L2 and Linf error calculation for the covariant form +function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, + equations::AbstractCovariantEquations{2}, + initial_condition, dg::DGSEM, cache, cache_analysis) + (; 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) + total_volume = zero(real(mesh)) + + # Iterate over all elements for error calculations + for element in eachelement(dg, cache) + + # Calculate errors at each volume quadrature node + for j in eachnode(dg), i in eachnode(dg) + x_node = Trixi.get_node_coords(node_coordinates, equations, dg, i, j, + element) + + # Convert exact solution into contravariant components using geometric + # information stored in aux vars + a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + u_exact = initial_condition(x_node, t, a_node, equations) + + # Compute the difference as usual + 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 + sqrtG = area_element(a_node, equations) + l2_error += diff .^ 2 * (weights[i] * weights[j] * sqrtG) + + # 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 + end + end + + # For L2 error, divide by total volume + l2_error = @. sqrt(l2_error / total_volume) + + return l2_error, linf_error +end +end # muladd diff --git a/src/callbacks_step/callbacks_step.jl b/src/callbacks_step/callbacks_step.jl new file mode 100644 index 0000000..3a9168c --- /dev/null +++ b/src/callbacks_step/callbacks_step.jl @@ -0,0 +1,2 @@ +include("analysis_covariant.jl") +include("stepsize_dg2d.jl") diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 822ffef..874512f 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -1,8 +1,12 @@ +@muladd begin +#! format: noindent + # Specialization of max_dt function for 3D equations in 2D manifolds function Trixi.max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, - constant_speed::False, equations::AbstractEquations{3}, dg::DG, cache) + constant_speed::False, equations::AbstractEquations{3}, dg::DG, + cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection max_scaled_speed = nextfloat(zero(t)) @@ -16,11 +20,13 @@ function Trixi.max_dt(u, t, lambda1, lambda2, lambda3 = max_abs_speeds(u_node, equations) # Local speeds transformed to the reference element - Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, + Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, + i, j, element) lambda1_transformed = abs(Ja11 * lambda1 + Ja12 * lambda2 + Ja13 * lambda3) - Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, i, + Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, + i, j, element) lambda2_transformed = abs(Ja21 * lambda1 + Ja22 * lambda2 + Ja23 * lambda3) @@ -36,3 +42,31 @@ function Trixi.max_dt(u, t, return 2 / (nnodes(dg) * max_scaled_speed) end + +# Specialization of max_dt function for covariant formulation on 2D manifolds +function Trixi.max_dt(u, t, mesh::P4estMesh{2}, constant_speed::False, + equations::AbstractCovariantEquations{2}, + dg::DG, cache) + (; aux_node_vars) = cache.auxiliary_variables + + # to avoid a division by zero if the speed vanishes everywhere, + # e.g. for steady-state linear advection + max_scaled_speed = nextfloat(zero(t)) + + # Because the covariant form computes max_abs_speeds using the contravariant + # velocity components already, there is no need to transform them here + for element in eachelement(dg, cache) + max_lambda1 = max_lambda2 = zero(max_scaled_speed) + for j in eachnode(dg), i in eachnode(dg) + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + lambda1, lambda2 = Trixi.max_abs_speeds(u_node, a_node, equations) + max_lambda1 = max(max_lambda1, lambda1) + max_lambda2 = max(max_lambda2, lambda2) + end + + max_scaled_speed = max(max_scaled_speed, max_lambda1 + max_lambda2) + end + return 2 / (nnodes(dg) * max_scaled_speed) +end +end # muladd diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl new file mode 100644 index 0000000..e5c26d7 --- /dev/null +++ b/src/equations/covariant_advection.jl @@ -0,0 +1,109 @@ +@muladd begin +#! format: noindent + +@doc raw""" + CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3, 3} + +A variable-coefficient linear advection equation can be defined on a two-dimensional +manifold $S \subset \mathbb{R}^3$ as +```math +\partial_t h + \nabla_S \cdot (h \vec{v}) = 0, +``` +where $\nabla_S \cdot$ is the horizontal divergence operator on $S$. We treat this problem +as a system of equations in which the first variable is the scalar conserved quantity $h$, +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. +``` +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, +and $\xi^1$ and $\xi^2$ are the local reference space coordinates. The velocity components +are spatially varying but assumed to be constant in time, so we do not apply any flux or +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] ++ +\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] += +\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. +""" +struct CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3, 3} end + +function Trixi.varnames(::typeof(cons2cons), ::CovariantLinearAdvectionEquation2D) + return ("scalar", "vcon1", "vcon2") +end + +# Convenience function to extract the velocity +function velocity(u, ::CovariantLinearAdvectionEquation2D) + return SVector(u[2], u[3]) +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 +@inline function Trixi.cons2entropy(u, aux_vars, + equations::CovariantLinearAdvectionEquation2D) + z = zero(eltype(u)) + return SVector{3}(u[1], z, z) +end + +# Numerical flux plus dissipation 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. +@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) +end + +# Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity +# components, as they should remain unchanged in time +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, aux_vars_ll, + aux_vars_rr, + orientation_or_normal_direction, + equations::CovariantLinearAdvectionEquation2D) + z = zero(eltype(u_ll)) + λ = 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) +end + +# Maximum wave speed with respect to the a specific orientation +@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])) +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)) +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) +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) +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index c8d0119..690fede 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -1,6 +1,124 @@ +@muladd begin +#! format: noindent + using Trixi: AbstractEquations +@doc raw""" + AbstractCovariantEquations{NDIMS, + NDIMS_AMBIENT, + NVARS} <: AbstractEquations{NDIMS, NVARS} + +Abstract type used to dispatch on systems of equations in covariant form, in which fluxes +and prognostic variables are stored and computed in terms of their contravariant components +defining their expansions in terms of the local covariant tangent basis. The type parameter +`NDIMS` denotes the dimension of the manifold on which the equations are solved, while +`NDIMS_AMBIENT` is the dimension of the ambient space in which such a manifold is embedded. +Some references on discontinuous Galerkin methods in covariant flux form are listed below: + +- 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) +- M. Baldauf (2021). A horizontally explicit, vertically implicit (HEVI) discontinuous + Galerkin scheme for the 2-dimensional Euler and Navier-Stokes equations using + terrain-following coordinates. Journal of Computational Physics 446:110635. [DOI: 10.1016/ + j.jcp.2021.110635](https://doi.org/10.1016/j.jcp.2021.110635) +- 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) + +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). +""" +abstract type AbstractCovariantEquations{NDIMS, + NDIMS_AMBIENT, + NVARS} <: AbstractEquations{NDIMS, NVARS} end + +""" + have_aux_node_vars(equations) + +Trait function determining whether `equations` requires the use of auxiliary variables. +Classical conservation laws such as the [`CompressibleEulerEquations2D`](@ref) do not +require auxiliary variables. The return value will be `True()` or `False()` to allow +dispatching on the return type. +""" +@inline have_aux_node_vars(::AbstractEquations) = False() + +# For the covariant form, the auxiliary variables are the the NDIMS^2 entries of the +# covariant basis matrix +@inline have_aux_node_vars(::AbstractCovariantEquations) = True() +@inline n_aux_node_vars(::AbstractCovariantEquations{NDIMS}) where {NDIMS} = NDIMS^2 + +# 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]") +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. +@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]) +end + +@doc raw""" + transform_to_contravariant(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, 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 +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 +# been transformed into the same local coordinate system. +@inline function (numflux::Trixi.FluxPlusDissipation)(u_ll, u_rr, + aux_vars_ll, aux_vars_rr, + orientation_or_normal_direction, + equations::AbstractCovariantEquations) + + # The flux and dissipation need to be defined for the specific equation set + flux = numflux.numerical_flux(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation_or_normal_direction, equations) + diss = numflux.dissipation(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation_or_normal_direction, equations) + return flux + diss +end + +# Central flux 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 been transformed +# into the same local coordinate system. +@inline function Trixi.flux_central(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation_or_normal_direction, + equations::AbstractCovariantEquations) + flux_ll = Trixi.flux(u_ll, aux_vars_ll, orientation_or_normal_direction, equations) + flux_rr = Trixi.flux(u_rr, aux_vars_rr, orientation_or_normal_direction, equations) + return 0.5f0 * (flux_ll + flux_rr) +end + abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end + +include("reference_data.jl") +include("covariant_advection.jl") include("compressible_moist_euler_2d_lucas.jl") include("shallow_water_3d.jl") +end # @muladd diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl new file mode 100644 index 0000000..e7a59fb --- /dev/null +++ b/src/equations/reference_data.jl @@ -0,0 +1,78 @@ +@muladd begin +#! format: noindent + +# Physical constants +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 + +@doc raw""" + initial_condition_gaussian(x, t) + +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 +```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), +``` +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 +```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} +``` +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: +- 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) + RealT = eltype(x) + + a = EARTH_RADIUS # radius of the sphere in metres + V = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY) # speed of rotation + alpha = convert(RealT, π / 4) # angle of rotation + h_0 = 1000.0f0 # bump height in metres + b_0 = 5.0f0 / (a^2) # bump width + lon_0, lat_0 = convert(RealT, 3π / 2), 0.0f0 # initial bump location + + # 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 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)) + + # get zonal and meridional components of the velocity + lon, lat = atan(x[2], x[1]), 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) +end +end # muladd diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl index e7e8e19..8340868 100644 --- a/src/equations/shallow_water_3d.jl +++ b/src/equations/shallow_water_3d.jl @@ -1,7 +1,3 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin #! format: noindent @@ -46,9 +42,9 @@ This affects the implementation and use of these equations in various ways: References: - J. Coté (1988). "A Lagrange multiplier approach for the metric terms of semi-Lagrangian models on the sphere". - Quarterly Journal of the Royal Meteorological Society 114, 1347-1352. https://doi.org/10.1002/qj.49711448310 -- Giraldo (2001). "A spectral element shallow water model on spherical geodesic grids". - https://doi.org/10.1002/1097-0363(20010430)35:8%3C869::AID-FLD116%3E3.0.CO;2-S + Quarterly Journal of the Royal Meteorological Society 114, 1347-1352. [DOI: 10.1002/qj.49711448310](https://doi.org/10.1002/qj.49711448310) +- F. X. Giraldo (2001). "A spectral element shallow water model on spherical geodesic grids". + [DOI: 10.1002/1097-0363(20010430)35:8<869::AID-FLD116>3.0.CO;2-S](https://doi.org/10.1002/1097-0363(20010430)35:8%3C869::AID-FLD116%3E3.0.CO;2-S) """ struct ShallowWaterEquations3D{RealT <: Real} <: Trixi.AbstractShallowWaterEquations{3, 5} @@ -208,7 +204,7 @@ Details are available in Eq. (4.1) in the paper: end """ - source_terms_lagrange_multiplier(u, du, x, t, + source_terms_lagrange_multiplier(u, du, x, t, equations::ShallowWaterEquations3D, normal_direction) @@ -267,8 +263,9 @@ end # Compute the wave celerity on the left and right h_ll = waterheight(u_ll, equations) h_rr = waterheight(u_rr, equations) - c_ll = sqrt(equations.gravity * h_ll) - c_rr = sqrt(equations.gravity * h_rr) + + c_ll = sqrt(max(equations.gravity * h_ll, 0.0f0)) + c_rr = sqrt(max(equations.gravity * h_rr, 0.0f0)) # The normal velocities are already scaled by the norm return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) @@ -383,4 +380,47 @@ 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 ca4cd1e..7a5e563 100644 --- a/src/meshes/p4est_cubed_sphere.jl +++ b/src/meshes/p4est_cubed_sphere.jl @@ -31,11 +31,19 @@ The mesh will have no boundaries. Appendix A of [Guba et al. (2014)](https://doi.org/10.5194/gmd-7-2803-2014). If set to `true`, the four corner vertex positions for each element will be obtained through an equiangular gnomonic projection [(Ronchi et al. 1996)](https://doi.org/10.1006/jcph.1996. - 0047), and the tree node coordinates within the element will be obtained by first using - a bilinear mapping based on the four corner vertices, and then projecting the bilinearly - mapped nodes onto the spherical surface by normalizing the resulting Cartesian - coordinates and scaling by `radius`. If set to `false`, the equiangular gnomonic - projection will be used for all node coordinates. + 0047), and the tree node coordinates within the element (i.e. the field + `tree_node_coordinates`) will be obtained by first using a bilinear mapping based on the + four corner vertices, and then projecting the bilinearly mapped nodes onto the spherical + surface by normalizing the resulting Cartesian coordinates and scaling by `radius`. If + set to `false`, the equiangular gnomonic projection will be used for all tree node + coordinates. + +!!! warning + Adaptivity and MPI parallelization are not yet supported for equations in covariant + 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/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index 734d430..4069482 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -14,7 +14,12 @@ Struct used for multiple dispatch on functions that compute the metric terms. When the argument `metric_terms` is of type `MetricTermsInvariantCurl`, the contravariant vectors are computed using the invariant curl form. ## References -- Kopriva, D. A. (2006). Metric identities and the discontinuous spectral element method on curvilinear meshes. Journal of Scientific Computing 26, 301-327. https://doi.org/10.1007/s10915-005-9070-8 -- Vinokur, M. and Yee, H. C. (2001). Extension of efficient low dissipation high order schemes for 3-D curvilinear moving grids. In Caughey, D. A., and Hafez, M. M. (eds.), Frontiers of Computational Fluid Dynamics 2002, World Scientific, Singapore, pp. 129–164. https://doi.org/10.1142/9789812810793_0008 +- Kopriva, D. A. (2006). Metric identities and the discontinuous spectral element method on + curvilinear meshes. Journal of Scientific Computing 26, 301-327. + [DOI: 10.1007/s10915-005-9070-8](https://doi.org/10.1007/s10915-005-9070-8) +- Vinokur, M. and Yee, H. C. (2001). Extension of efficient low dissipation high order + schemes for 3-D curvilinear moving grids. In Caughey, D. A., and Hafez, M. M. (eds.), + Frontiers of Computational Fluid Dynamics 2002, World Scientific, Singapore, pp. 129–164. + [DOI: 10.1142/9789812810793_0008](https://doi.org/10.1142/9789812810793_0008) """ struct MetricTermsInvariantCurl end diff --git a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl index 7421c99..4944918 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl @@ -4,7 +4,9 @@ # AbstractEquations{3} for 3D equations. This change is necessary to support the Cartesian implementation # of a 2D manifold embedded in a 3D space. function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{2}, - equations::AbstractEquations{3}, + equations::Union{AbstractEquations{3}, + AbstractCovariantEquations{2, + 3}}, initial_condition, solver; source_terms = nothing, diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl similarity index 93% rename from src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl rename to src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl index 3199243..b222c57 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl @@ -1,12 +1,14 @@ @muladd begin #! format: noindent -# New p4est element container that allows the use of a PtrArray for the contravariant_vectors +# New p4est element container that allows the use of a PtrArray for the +# contravariant_vectors, and stores the auxiliary_variables as an additional field mutable struct P4estElementContainerPtrArray{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, NDIMSP2, NDIMSP3, ContravariantVectors <: - AbstractArray{RealT, NDIMSP3}} <: + AbstractArray{RealT, + NDIMSP3}} <: Trixi.AbstractContainer # Physical coordinates at each node node_coordinates::Array{RealT, NDIMSP2} # [orientation, node_i, node_j, node_k, element] @@ -29,7 +31,7 @@ mutable struct P4estElementContainerPtrArray{NDIMS, RealT <: Real, uEltype <: Re end @inline function Trixi.nelements(elements::P4estElementContainerPtrArray) - size(elements.node_coordinates, ndims(elements) + 2) + return size(elements.node_coordinates, ndims(elements) + 2) end @inline Base.ndims(::P4estElementContainerPtrArray{NDIMS}) where {NDIMS} = NDIMS @inline function Base.eltype(::P4estElementContainerPtrArray{NDIMS, @@ -38,14 +40,25 @@ end RealT, uEltype } - uEltype + return uEltype +end + +# Extract contravariant vector Ja^i (i = index) as SVector +# This function dispatches on the type of contravariant_vectors +static2val(::Trixi.StaticInt{N}) where {N} = Val{N}() +@inline function Trixi.get_contravariant_vector(index, contravariant_vectors::PtrArray, + indices...) + return SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), + static2val(static_size(contravariant_vectors, + Trixi.StaticInt(1))))) end # Create element container and initialize element data. # This function dispatches on the dimensions of the mesh and the equation (AbstractEquations{3}) function Trixi.init_elements(mesh::Union{P4estMesh{2, 3, RealT}, T8codeMesh{2}}, - equations::AbstractEquations{3}, + equations::Union{AbstractEquations{3}, + AbstractCovariantEquations{2, 3}}, basis, metric_terms, ::Type{uEltype}) where {RealT <: Real, uEltype <: Real} @@ -91,24 +104,26 @@ function Trixi.init_elements(mesh::Union{P4estMesh{2, 3, RealT}, ntuple(_ -> nnodes(basis), NDIMS - 1)..., NDIMS * 2, nelements)) + ContravariantVectors = typeof(contravariant_vectors) elements = P4estElementContainerPtrArray{NDIMS, RealT, uEltype, NDIMS + 1, NDIMS + 2, - NDIMS + 3, typeof(contravariant_vectors)}(node_coordinates, - jacobian_matrix, - contravariant_vectors, - inverse_jacobian, - surface_flux_values, - _node_coordinates, - _jacobian_matrix, - _contravariant_vectors, - _inverse_jacobian, - _surface_flux_values) + NDIMS + 3, + ContravariantVectors}(node_coordinates, + jacobian_matrix, + contravariant_vectors, + inverse_jacobian, + surface_flux_values, + _node_coordinates, + _jacobian_matrix, + _contravariant_vectors, + _inverse_jacobian, + _surface_flux_values) init_elements_2d_manifold_in_3d!(elements, mesh, basis, metric_terms) - return elements 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}}, basis::LobattoLegendreBasis, @@ -202,8 +217,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates, Trixi.multiply_dimensionwise!(view(node_coordinates, :, :, :, element), matrix1, matrix2, view(mesh.tree_node_coordinates, :, :, :, - tree), - tmp1) + tree), tmp1) end end 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 new file mode 100644 index 0000000..ed411c2 --- /dev/null +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl @@ -0,0 +1,237 @@ +@muladd begin +#! format: noindent + +# Container for storing values of auxiliary variables at volume/surface quadrature nodes +struct P4estAuxiliaryNodeVariableContainer{NDIMS, uEltype <: Real, NDIMSP2} + aux_node_vars::Array{uEltype, NDIMSP2} # [variable, i, j, k, element] + aux_surface_node_vars::Array{uEltype, NDIMSP2} # [variable, i, j, k, element] + + # internal `resize!`able storage + _aux_node_vars::Vector{uEltype} + _aux_surface_node_vars::Vector{uEltype} +end + +# Return the auxiliary variables at a given volume node index +@inline function get_node_aux_vars(aux_node_vars, equations, ::DG, indices...) + return SVector(ntuple(@inline(v->aux_node_vars[v, indices...]), + Val(n_aux_node_vars(equations)))) +end + +# Return the auxiliary variables at a given surface node index +@inline function get_surface_node_aux_vars(aux_surface_node_vars, equations, ::DG, + indices...) + aux_vars_ll = SVector(ntuple(@inline(v->aux_surface_node_vars[1, v, indices...]), + Val(n_aux_node_vars(equations)))) + aux_vars_rr = SVector(ntuple(@inline(v->aux_surface_node_vars[2, v, indices...]), + Val(n_aux_node_vars(equations)))) + return aux_vars_ll, aux_vars_rr +end + +# Create auxiliary node variable container and initialize auxiliary variables +function init_auxiliary_node_variables(mesh::Union{P4estMesh, T8codeMesh}, + equations, dg, elements, interfaces) + nelements = Trixi.ncells(mesh) + ninterfaces = Trixi.count_required_surfaces(mesh).interfaces + NDIMS = ndims(elements) + uEltype = eltype(elements) + + _aux_node_vars = Vector{uEltype}(undef, + n_aux_node_vars(equations) * + nnodes(dg)^NDIMS * nelements) + aux_node_vars = Trixi.unsafe_wrap(Array, pointer(_aux_node_vars), + (n_aux_node_vars(equations), + ntuple(_ -> nnodes(dg), NDIMS)..., + nelements)) + _aux_surface_node_vars = Vector{uEltype}(undef, + 2 * n_aux_node_vars(equations) * + nnodes(dg)^(NDIMS - 1) * + ninterfaces) + aux_surface_node_vars = Trixi.unsafe_wrap(Array, + pointer(_aux_surface_node_vars), + (2, n_aux_node_vars(equations), + ntuple(_ -> nnodes(dg), + NDIMS - 1)..., + ninterfaces)) + + auxiliary_variables = P4estAuxiliaryNodeVariableContainer{NDIMS, + uEltype, + NDIMS + 2}(aux_node_vars, + aux_surface_node_vars, + _aux_node_vars, + _aux_surface_node_vars) + + init_auxiliary_node_variables!(auxiliary_variables, mesh, equations, dg, elements) + init_auxiliary_surface_node_variables!(auxiliary_variables, mesh, equations, dg, + interfaces) + return auxiliary_variables +end + +# By default, the auxiliary surface node variables are just the volume node variables +# evaluated at the surface nodes, similarly to prolong2interfaces. This could, however, be +# specialized for other applications. +function init_auxiliary_surface_node_variables!(auxiliary_variables::P4estAuxiliaryNodeVariableContainer{2}, + mesh::Union{P4estMesh{2}, + T8codeMesh{2}}, equations, + dg::DG, interfaces) + (; aux_node_vars, aux_surface_node_vars) = auxiliary_variables + index_range = eachnode(dg) + + Trixi.@threaded for interface in axes(interfaces.node_indices, 2) + # Copy solution data from the primary element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + primary_element = interfaces.neighbor_ids[1, interface] + primary_indices = interfaces.node_indices[1, interface] + + i_primary_start, i_primary_step = Trixi.index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = Trixi.index_to_start_step_2d(primary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + for i in eachnode(dg) + for v in axes(aux_node_vars, 1) + aux_surface_node_vars[1, v, i, interface] = aux_node_vars[v, + i_primary, + j_primary, + primary_element] + end + i_primary += i_primary_step + j_primary += j_primary_step + end + + # Copy solution data from the secondary element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + secondary_element = interfaces.neighbor_ids[2, interface] + secondary_indices = interfaces.node_indices[2, interface] + + i_secondary_start, i_secondary_step = Trixi.index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = Trixi.index_to_start_step_2d(secondary_indices[2], + index_range) + + i_secondary = i_secondary_start + j_secondary = j_secondary_start + for i in eachnode(dg) + for v in axes(aux_node_vars, 1) + aux_surface_node_vars[2, v, i, interface] = aux_node_vars[v, + i_secondary, + j_secondary, + secondary_element] + end + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + return nothing +end + +# Get Cartesian node positions for the covariant form, dispatching on the dimension of the +# manifold as well as the ambient dimension +@inline function Trixi.get_node_coords(x, + ::AbstractCovariantEquations{NDIMS, + NDIMS_AMBIENT}, + ::DG, + indices...) where {NDIMS, NDIMS_AMBIENT} + return SVector(ntuple(@inline(idx->x[idx, indices...]), NDIMS_AMBIENT)) +end + +# Compute the auxiliary metric terms for the covariant form, assuming that the +# domain is a spherical shell. We do not make any assumptions on the mesh topology, but we +# require that the elements are constructed using the element-local mapping from the +# following paper: +# +# O. Guba, M. A. Taylor, P. A. Ullrich, J. R. Overfelt, and M. N. Levy (2014). The spectral +# element method (SEM) on variable-resolution grids: evaluating grid sensitivity and +# resolution-aware numerical viscosity. Geoscientific Model Development 7(6) 2803–2816. +# DOI: 10.5194/gmd-7-2803-2014 +# +# When creating a cubed sphere mesh using P4estMeshCubedSphere2D, this is enabled by +# passing the keyword argument "element_local_mapping = true" to the constructor and +# ensuring that the polynomial degree of the mesh is equal to that of the solver. +# +# Otherwise, the mesh's tree node coordinates will be interpolated to the solver's +# physical node coordinates, and this would introduce a polynomial approximation of the +# geometry, making the analytical metric terms computed here no longer correct. +function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, 3}, + equations::AbstractCovariantEquations{2, 3}, dg, + elements) + (; tree_node_coordinates) = mesh + + # Check that the degree of the mesh matches that of the solver + @assert length(mesh.nodes) == nnodes(dg) + + # The tree node coordinates are assumed to be on the spherical shell centred around the + # origin. Therefore, we can compute the radius once and use it throughout. + radius = norm(Trixi.get_node_coords(tree_node_coordinates, equations, dg, 1, 1, 1)) + + Trixi.@threaded for element in 1:Trixi.ncells(mesh) + # Extract the corner vertex positions from the P4estMesh + v1 = Trixi.get_node_coords(tree_node_coordinates, equations, dg, + 1, 1, element) + v2 = Trixi.get_node_coords(tree_node_coordinates, equations, dg, + nnodes(dg), 1, element) + v3 = Trixi.get_node_coords(tree_node_coordinates, equations, dg, + nnodes(dg), nnodes(dg), element) + v4 = Trixi.get_node_coords(tree_node_coordinates, equations, dg, + 1, nnodes(dg), element) + + # 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) + 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_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) + # 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) + + # Project the mapped local coordinates onto the sphere using a simple scaling + scaling_factor = radius / norm(x_bilinear) + x = scaling_factor * x_bilinear + + # Convert Cartesian coordinates to longitude and latitude + lon, lat = atan(x[2], x[1]), asin(x[3] / radius) + + # Compute trigonometric terms needed for analytical metrics + sinlon, coslon = sincos(lon) + sinlat, coslat = sincos(lat) + a11 = sinlon * sinlon * coslat * coslat + sinlat * sinlat + a12 = a21 = -sinlon * coslon * coslat * coslat + a13 = -coslon * sinlat * coslat + a22 = coslon * coslon * coslat * coslat + sinlat * sinlat + a23 = -sinlon * sinlat * coslat + a31 = -coslon * sinlat + a32 = -sinlon * sinlat + 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) +end +end # muladd diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index ea6f835..22331bf 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -1,6 +1,7 @@ # This method is called when a SemidiscretizationHyperbolic is constructed. # It constructs the basic `cache` used throughout the simulation to compute -# the RHS etc. +# the RHS etc, and is modified here to use custom metric terms as well as provide the +# option to # use auxiliary variables. function Trixi.create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::DG, ::Any, metric_terms, ::Type{uEltype}) where {uEltype <: Real} # Make sure to balance the `p4est` before creating any containers @@ -19,14 +20,30 @@ function Trixi.create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::D Trixi.create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) cache = (; cache..., Trixi.create_cache(mesh, equations, dg.mortar, uEltype)...) + # Add specialized parts of the cache for auxiliary node variables + cache = (; cache..., + create_cache_auxiliary(mesh, equations, + have_aux_node_vars(equations), + dg, elements, interfaces)...) return cache end -# Extract contravariant vector Ja^i (i = index) as SVector -# This function dispatches on the type of contravariant_vectors -static2val(::Trixi.StaticInt{N}) where {N} = Val{N}() -@inline function Trixi.get_contravariant_vector(index, contravariant_vectors::PtrArray, - indices...) - SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), - static2val(static_size(contravariant_vectors, Trixi.StaticInt(1))))) +# If there are auxiliary variables, initialize them +function create_cache_auxiliary(mesh, equations, have_aux_node_vars::True, dg, elements, + interfaces) + auxiliary_variables = init_auxiliary_node_variables(mesh, equations, dg, elements, + interfaces) + return (; auxiliary_variables) +end + +# Do nothing if there are no auxiliary variables +function create_cache_auxiliary(mesh, equations, have_aux_node_vars::False, dg, elements, + interfaces) + return NamedTuple() end + +include("containers_2d_manifold_in_3d_cartesian.jl") +include("containers_2d_manifold_in_3d_covariant.jl") + +include("dg_2d_manifold_in_3d_cartesian.jl") +include("dg_2d_manifold_in_3d_covariant.jl") 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 4264c92..7c1a81e 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 @@ -1,3 +1,5 @@ +@muladd begin +#! format: noindent function Trixi.rhs!(du, u, t, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, @@ -11,7 +13,8 @@ function Trixi.rhs!(du, u, t, # Calculate volume integral Trixi.@trixi_timeit Trixi.timer() "volume integral" begin Trixi.calc_volume_integral!(du, u, mesh, - Trixi.have_nonconservative_terms(equations), equations, + Trixi.have_nonconservative_terms(equations), + equations, dg.volume_integral, dg, cache) end @@ -24,7 +27,8 @@ function Trixi.rhs!(du, u, t, # Calculate interface fluxes Trixi.@trixi_timeit Trixi.timer() "interface flux" begin Trixi.calc_interface_flux!(cache.elements.surface_flux_values, mesh, - Trixi.have_nonconservative_terms(equations), equations, + Trixi.have_nonconservative_terms(equations), + equations, dg.surface_integral, dg, cache) end @@ -60,7 +64,8 @@ function Trixi.rhs!(du, u, t, end # Apply Jacobian from mapping to reference element - Trixi.@trixi_timeit Trixi.timer() "Jacobian" Trixi.apply_jacobian!(du, mesh, equations, + Trixi.@trixi_timeit Trixi.timer() "Jacobian" Trixi.apply_jacobian!(du, mesh, + equations, dg, cache) # Calculate source terms @@ -71,10 +76,11 @@ function Trixi.rhs!(du, u, t, return nothing end -# Weak-form kernel for 3D equations solved on 2D manifolds +# Weak-form kernel for 3D equations solved in 2D manifolds @inline function Trixi.weak_form_kernel!(du, u, element, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, + mesh::Union{StructuredMesh{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations::AbstractEquations{3}, @@ -93,7 +99,8 @@ end # Compute the contravariant flux by taking the scalar product of the # first contravariant vector Ja^1 and the flux vector - Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, j, + Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, + j, element) contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 for ii in eachnode(dg) @@ -104,7 +111,8 @@ end # Compute the contravariant flux by taking the scalar product of the # second contravariant vector Ja^2 and the flux vector - Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, i, j, + Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, i, + j, element) contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 for jj in eachnode(dg) @@ -145,3 +153,4 @@ function calc_sources_2d_manifold_in_3d!(du, u, t, source_terms, return nothing end +end # muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl new file mode 100644 index 0000000..f9f5316 --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl @@ -0,0 +1,195 @@ +@muladd begin +#! format: noindent + +# Compute coefficients for an initial condition that uses auxiliary variables +function Trixi.compute_coefficients!(u, func, t, mesh::P4estMesh{2}, + equations::AbstractCovariantEquations{2}, dg::DG, + cache) + (; aux_node_vars) = cache.auxiliary_variables + (; node_coordinates) = cache.elements + + Trixi.@threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + # Get physical Cartesian coordinates at node (i, j) + x_node = Trixi.get_node_coords(node_coordinates, equations, dg, i, j, + element) + + # Get auxiliary variables at node (i, j) + a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + + # Compute nodal values of the provided function + u_node = func(x_node, t, a_node, equations) + + # Set nodal variables in cache + Trixi.set_node_vars!(u, u_node, equations, dg, i, j, element) + end + end +end + +# Weak form kernel which uses contravariant flux components, passing the geometric +# information contained in the auxiliary variables to the flux function +@inline function Trixi.weak_form_kernel!(du, u, element, mesh::P4estMesh{2}, + nonconservative_terms::False, + equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, alpha = true) + (; derivative_dhat) = dg.basis + (; aux_node_vars) = cache.auxiliary_variables + + for j in eachnode(dg), i in eachnode(dg) + # Get solution variables and auxiliary variables at node (i, j) + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + + # Evaluate the contravariant flux components + flux1 = flux(u_node, a_node, 1, equations) + flux2 = flux(u_node, a_node, 2, equations) + + # Apply weak form derivative with respect to ξ¹ + for ii in eachnode(dg) + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], + flux1, equations, dg, ii, j, element) + end + + # Apply weak form derivative with respect to ξ² + for jj in eachnode(dg) + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], + flux2, equations, dg, i, jj, element) + end + end + + return nothing +end + +# 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::False, + equations::AbstractCovariantEquations{2}, + volume_flux, dg::DGSEM, cache, + alpha = true) + (; derivative_split) = dg.basis + (; aux_node_vars) = cache.auxiliary_variables + + for j in eachnode(dg), i in eachnode(dg) + # Get solution variables and auxiliary variables at node (i, j) + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + + # ξ¹ direction + for ii in (i + 1):nnodes(dg) + # Get solution variables and auxiliary variables at node (ii, j) + u_node_ii = Trixi.get_node_vars(u, equations, dg, ii, j, element) + a_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, element) + + # Evaluate contravariant component 1 of the variable-coefficient two-point flux + flux1 = volume_flux(u_node, u_node_ii, a_node, a_node_ii, 1, equations) + + # Multiply by entry of split derivative matrix and add to right-hand side + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[i, ii], flux1, + equations, dg, i, j, element) + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[ii, i], flux1, + equations, dg, ii, j, element) + end + + # ξ² direction + for jj in (j + 1):nnodes(dg) + # Get solution variables and auxiliary variables at node (i, jj) + u_node_jj = Trixi.get_node_vars(u, equations, dg, i, jj, element) + a_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, element) + + # Evaluate contravariant component 2 of the variable-coefficient two-point flux + flux2 = volume_flux(u_node, u_node_jj, a_node, a_node_jj, 2, equations) + + # Multiply by entry of split derivative matrix and add to right-hand side + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[j, jj], flux2, + equations, dg, i, j, element) + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[jj, j], flux2, + equations, dg, i, jj, element) + end + end + + return nothing +end + +# Pointwise interface flux, transforming the contravariant prognostic variables into the +# local coordinate system +@inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, + nonconservative_terms::False, + 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) = surface_integral + + # 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 = contravariant2spherical(u_rr, aux_vars_rr, equations) + u_rr_transformed_to_ll = spherical2contravariant(u_rr_spherical, aux_vars_ll, + equations) + if isodd(primary_direction_index) + flux_primary = -surface_flux(u_rr_transformed_to_ll, u_ll, + aux_vars_ll, aux_vars_ll, + (primary_direction_index + 1) >> 1, equations) + else + flux_primary = surface_flux(u_ll, u_rr_transformed_to_ll, + aux_vars_ll, aux_vars_ll, + (primary_direction_index + 1) >> 1, equations) + end + + # Compute flux in the secondary element's coordinate system + u_ll_spherical = contravariant2spherical(u_ll, aux_vars_ll, equations) + u_ll_transformed_to_rr = spherical2contravariant(u_ll_spherical, aux_vars_rr, + equations) + if isodd(secondary_direction_index) + flux_secondary = -surface_flux(u_ll_transformed_to_rr, u_rr, + aux_vars_rr, aux_vars_rr, + (secondary_direction_index + 1) >> 1, equations) + else + flux_secondary = surface_flux(u_rr, u_ll_transformed_to_rr, + aux_vars_rr, aux_vars_rr, + (secondary_direction_index + 1) >> 1, 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 + +# Apply the exact Jacobian stored in auxiliary variables +function Trixi.apply_jacobian!(du, mesh::P4estMesh{2}, + 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) + a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + factor = -1 / area_element(a_node, equations) + + for v in eachvariable(equations) + du[v, i, j, element] *= factor + end + end + end + + return nothing +end +end # @muladd diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index d39e464..cfcf2cf 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -1,3 +1 @@ -include("dgsem_p4est/containers_2d_manifold_in_3d.jl") include("dgsem_p4est/dg.jl") -include("dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl") diff --git a/test/Project.toml b/test/Project.toml index 67ebae1..d820deb 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,4 +8,4 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" NLsolve = "4.5.1" OrdinaryDiffEq = "6.53.2" Test = "1" -Trixi = "0.7, 0.8, 0.9" +Trixi = "0.9" diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 6aa252e..0d6b8e3 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -7,23 +7,24 @@ include("test_trixiatmo.jl") EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") -@trixiatmo_testset "Spherical advection, Cartesian form with standard mapping" begin +@trixiatmo_testset "Spherical advection, Cartesian weak form, LLF surface flux" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_spherical_advection_cartesian.jl"), + "elixir_shallowwater_cubed_sphere_shell_advection.jl"), l2=[ - 8.44498914e-03, - 8.23407970e-03, - 1.84210216e-03, - 0.00000000e+00, - 6.44302432e-02 + 0.7963221028128544, + 20.317489255709344, + 8.811382597221147, + 20.317512894705885, + 0.0 ], linf=[ - 9.48946345e-02, - 9.64807833e-02, - 1.37450721e-02, - 0.00000000e+00, - 0.40932295554303133 - ], element_local_mapping=false) + 10.872101730749478, + 289.65159632622454, + 95.16499199537975, + 289.65159632621726, + 0.0 + ]) + # and small reference values # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -34,23 +35,75 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") end end -@trixiatmo_testset "Spherical advection, Cartesian form with element-local mapping" begin +@trixiatmo_testset "Spherical advection, Cartesian weak form, element-local mapping" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_spherical_advection_cartesian.jl"), + "elixir_shallowwater_cubed_sphere_shell_advection.jl"), l2=[ - 0.00897499647958351, - 0.008743150662795171, - 0.001993813300529814, - 0.0, - 0.06452666141326718 + 0.8933384470346987, + 22.84334163690791, + 9.776600357597692, + 22.843351937152637, + 0.0 ], linf=[ - 0.10136306017376362, - 0.10308252591097666, - 0.014142319489653277, - 0.0, - 0.41047171697830986 + 14.289456304666146, + 380.6958334078372, + 120.59259301636666, + 380.6958334078299, + 0.0 ], element_local_mapping=true) + # and small reference values + # 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 "Spherical advection, covariant weak form, LLF surface flux" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_spherical_advection_covariant.jl"), + l2=[1.0007043506351705, 0.0, 0.0], + linf=[14.235905681508598, 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 + +# 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"), + l2=[1.0007043506351412, 0.0, 0.0], + linf=[14.23590568150928, 0.0, 0.0], + volume_integral=VolumeIntegralFluxDifferencing(flux_central)) + # 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 + +# 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"), + l2=[2.499889861385917, 0.0, 0.0], + linf=[38.085244441156085, 0.0, 0.0], + volume_integral=VolumeIntegralFluxDifferencing(flux_central), + surface_flux=flux_central) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let