From 30bdb234dc06b0a33e6b0a69373cc3ead44cfce1 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Fri, 8 Nov 2024 16:22:02 +0100 Subject: [PATCH 1/2] update formatter check to 1.0.60 to be consistent with Trixi.jl (#45) --- .github/workflows/FormatCheck.yml | 2 +- test/test_2d_moist_euler.jl | 28 ++++++++++++++-------------- test/test_spherical_advection.jl | 8 ++++---- 3 files changed, 19 insertions(+), 19 deletions(-) diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 8001533..f0eaf85 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -29,7 +29,7 @@ jobs: # TODO: Change the call below to # format(".") run: | - julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))' + julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))' julia -e 'using JuliaFormatter; format(".")' - name: Format check run: | diff --git a/test/test_2d_moist_euler.jl b/test/test_2d_moist_euler.jl index d669e96..646ea0e 100644 --- a/test/test_2d_moist_euler.jl +++ b/test/test_2d_moist_euler.jl @@ -16,7 +16,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.0006660124630171347, 0.008969786054960861, 0.0, - 0.0, + 0.0 ], linf=[ 1.0312042909910168e-5, @@ -24,7 +24,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.006392107590872236, 0.07612038028310053, 0.0, - 0.0, + 0.0 ], polydeg=3, tspan=(0.0, 0.1)) @@ -46,7 +46,7 @@ end 7.938812668709457, 4500.747616411578, 0.00015592413050814787, - 0.00014163475049532796, + 0.00014163475049532796 ], linf=[ 0.1427479052298466, @@ -54,7 +54,7 @@ end 91.56822550162855, 49528.383866247605, 0.0019364397602254623, - 0.0013259689889851285, + 0.0013259689889851285 ], polydeg=3, cells_per_dimension=(16, 16), @@ -77,7 +77,7 @@ end 0.0006974588377288118, 1.715668353329522, 8.831269143134121e-7, - 1.025579538944668e-6, + 1.025579538944668e-6 ], linf=[ 8.055695643149896e-5, @@ -85,7 +85,7 @@ end 0.005897639251024437, 19.24776030163048, 1.0043133039065386e-5, - 1.1439046776775402e-5, + 1.1439046776775402e-5 ], polydeg=3, cells_per_dimension=(16, 8), @@ -109,7 +109,7 @@ end 0.01172830034500581, 9.898560584459009, 0.0, - 0.0, + 0.0 ], linf=[ 0.0017602202683439927, @@ -117,7 +117,7 @@ end 0.5945327351674782, 489.89171406268724, 0.0, - 0.0, + 0.0 ], polydeg=3, cells_per_dimension=(10, 8), @@ -140,7 +140,7 @@ end 2.609465609739115e-5, 6.323484379066432e-5, 0.0, - 0.0, + 0.0 ], linf=[ 7.549984224430872e-5, @@ -148,7 +148,7 @@ end 0.00015964938767742964, 0.0005425860570698049, 0.0, - 0.0, + 0.0 ], polydeg=3, cells_per_dimension=(10, 8), @@ -171,7 +171,7 @@ end 0.015104690877128945, 0.5242098451814421, 5.474006801215573e-10, - 1.1103883907226752e-10, + 1.1103883907226752e-10 ], linf=[ 0.00013219484616722177, @@ -179,7 +179,7 @@ end 0.03789645369775574, 3.90888311638264, 3.938382289041286e-9, - 6.892033377287209e-10, + 6.892033377287209e-10 ], polydeg=3, cells_per_dimension=(10, 8), @@ -203,7 +203,7 @@ end 0.07460345535073129, 5.943049264963717, 4.471792794168725e-9, - 7.10320253652373e-10, + 7.10320253652373e-10 ], linf=[ 0.0007084183215528839, @@ -211,7 +211,7 @@ end 0.3697160082709827, 27.843155286573165, 2.1168438904322837e-8, - 3.691699932047233e-9, + 3.691699932047233e-9 ], polydeg=3, cells_per_dimension=(10, 8), diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index c9e9300..bb735a5 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -15,14 +15,14 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") 8.23414117e-03, 1.84210648e-03, 0.00000000e+00, - 6.44302430e-02, + 6.44302430e-02 ], linf=[ 9.48950488e-02, 9.64811952e-02, 1.37453400e-02, 0.00000000e+00, - 4.09322999e-01, + 4.09322999e-01 ], element_local_mapping=false) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -42,14 +42,14 @@ end 8.74316738e-03, 1.99380928e-03, 0.00000000e+00, - 6.45266612e-02, + 6.45266612e-02 ], linf=[ 1.01363241e-01, 1.03082705e-01, 1.41424723e-02, 0.00000000e+00, - 4.10471741e-01, + 4.10471741e-01 ], element_local_mapping=true) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 49ec206e95ab6cfa43105428418c7dc737acb42c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 8 Nov 2024 17:26:52 +0100 Subject: [PATCH 2/2] Add 3D Cartesian Shallow Water Equations (#36) Co-authored-by: Tristan Montoya --- Project.toml | 2 + ...xir_euler_spherical_advection_cartesian.jl | 49 +-- ...wwater_cubed_sphere_shell_EC_correction.jl | 125 ++++++ ...wwater_cubed_sphere_shell_EC_projection.jl | 120 ++++++ ...hallowwater_cubed_sphere_shell_standard.jl | 109 +++++ src/TrixiAtmo.jl | 9 +- src/callbacks_step/stepsize_dg2d.jl | 38 ++ src/equations/equations.jl | 1 + src/equations/shallow_water_3d.jl | 386 ++++++++++++++++++ src/semidiscretization/semidiscretization.jl | 20 + ...retization_hyperbolic_2d_manifold_in_3d.jl | 5 +- .../containers_2d_manifold_in_3d.jl | 161 +++++++- src/solvers/dgsem_p4est/dg.jl | 24 ++ .../dg_2d_manifold_in_3d_cartesian.jl | 104 ++++- src/solvers/solvers.jl | 2 +- test/runtests.jl | 4 + test/test_3d_shallow_water.jl | 91 +++++ test/test_spherical_advection.jl | 36 +- 18 files changed, 1218 insertions(+), 68 deletions(-) create mode 100644 examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl create mode 100644 examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl create mode 100644 examples/elixir_shallowwater_cubed_sphere_shell_standard.jl create mode 100644 src/callbacks_step/stepsize_dg2d.jl create mode 100644 src/equations/shallow_water_3d.jl create mode 100644 src/semidiscretization/semidiscretization.jl create mode 100644 test/test_3d_shallow_water.jl diff --git a/Project.toml b/Project.toml index b70d087..8b17e5e 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0-DEV" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" @@ -17,6 +18,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] LinearAlgebra = "1" MuladdMacro = "0.2.2" +LoopVectorization = "0.12.118" NLsolve = "4.5.1" Reexport = "1.0" Static = "0.8, 1" diff --git a/examples/elixir_euler_spherical_advection_cartesian.jl b/examples/elixir_euler_spherical_advection_cartesian.jl index b887955..6e82cb0 100644 --- a/examples/elixir_euler_spherical_advection_cartesian.jl +++ b/examples/elixir_euler_spherical_advection_cartesian.jl @@ -57,7 +57,8 @@ end # Source term function to transform the Euler equations into the linear advection equations with variable advection velocity function source_terms_convert_to_linear_advection(u, du, x, t, - equations::CompressibleEulerEquations3D) + equations::CompressibleEulerEquations3D, + normal_direction) v1 = u[2] / u[1] v2 = u[3] / u[1] v3 = u[4] / u[1] @@ -70,45 +71,17 @@ function source_terms_convert_to_linear_advection(u, du, x, t, return SVector(0.0, s2, s3, s4, s5) end -# Custom RHS that applies a source term that depends on du (used to convert the 3D Euler equations into the 3D linear advection -# equations with position-dependent advection speed) -function rhs_semi_custom!(du_ode, u_ode, semi, t) - # Compute standard Trixi RHS - Trixi.rhs!(du_ode, u_ode, semi, t) - - # Now apply the custom source term - Trixi.@trixi_timeit Trixi.timer() "custom source term" begin - @unpack solver, equations, cache = semi - @unpack node_coordinates = cache.elements - - # Wrap the solution and RHS - u = Trixi.wrap_array(u_ode, semi) - du = Trixi.wrap_array(du_ode, semi) - - Trixi.@threaded for element in eachelement(solver, cache) - for j in eachnode(solver), i in eachnode(solver) - u_local = Trixi.get_node_vars(u, equations, solver, i, j, element) - du_local = Trixi.get_node_vars(du, equations, solver, i, j, element) - x_local = Trixi.get_node_coords(node_coordinates, equations, solver, - i, j, element) - source = source_terms_convert_to_linear_advection(u_local, du_local, - x_local, t, equations) - Trixi.add_to_node_vars!(du, source, equations, solver, i, j, element) - end - end - end -end - initial_condition = initial_condition_advection_sphere element_local_mapping = false -mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, - initial_refinement_level = 0, - element_local_mapping = element_local_mapping) +mesh = P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, + initial_refinement_level = 0, + element_local_mapping = element_local_mapping) # A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convert_to_linear_advection) ############################################################################### # ODE solvers, callbacks etc. @@ -117,12 +90,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, pi) ode = semidiscretize(semi, tspan) -# Create custom discretization that runs with the custom RHS -ode_semi_custom = ODEProblem(rhs_semi_custom!, - ode.u0, - ode.tspan, - semi) - # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup # and resets the timers summary_callback = SummaryCallback() @@ -148,7 +115,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode_semi_custom, CarpenterKennedy2N54(williamson_condition = false), +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep = false, callback = callbacks); diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl new file mode 100644 index 0000000..9f12e2c --- /dev/null +++ b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl @@ -0,0 +1,125 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiAtmo +############################################################################### +# Entropy conservation for the spherical shallow water equations in Cartesian +# form obtained through an entropy correction term + +equations = ShallowWaterEquations3D(gravity_constant = 9.81) + +# Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux +polydeg = 3 +volume_flux = flux_fjordholm_etal #flux_wintermeyer_etal +surface_flux = flux_fjordholm_etal #flux_wintermeyer_etal #flux_lax_friedrichs +solver = DGSEM(polydeg = polydeg, + surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# 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::ShallowWaterEquations3D) + # Gaussian density + rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) + + # 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, 0), equations) +end + +# Source term function to apply the entropy correction term and the Lagrange multiplier to the semi-discretization. +# The vector contravariant_normal_vector contains the normal contravariant vectors scaled with the inverse Jacobian. +# The Lagrange multiplier is applied with the vector x, which contains the exact normal direction to the 2D manifold. +function source_terms_entropy_correction(u, du, x, t, + equations::ShallowWaterEquations3D, + contravariant_normal_vector) + + # Entropy correction term + s2 = -contravariant_normal_vector[1] * equations.gravity * u[1]^2 + s3 = -contravariant_normal_vector[2] * equations.gravity * u[1]^2 + s4 = -contravariant_normal_vector[3] * equations.gravity * u[1]^2 + + new_du = SVector(du[1], du[2] + s2, du[3] + s3, du[4] + s4, du[5]) + + # Apply Lagrange multipliers using the exact normal vector to the 2D manifold (x) + s = source_terms_lagrange_multiplier(u, new_du, x, t, equations, x) + + return SVector(s[1], s[2] + s2, s[3] + s3, s[4] + s4, s[5]) +end + +initial_condition = initial_condition_advection_sphere + +mesh = P4estMeshCubedSphere2D(5, 2.0, polydeg = polydeg, + initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + metric_terms = MetricTermsInvariantCurl(), + source_terms = source_terms_entropy_correction) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +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 = 10, + save_analysis = true, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight, energy_total)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +# 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, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl new file mode 100644 index 0000000..5b43e82 --- /dev/null +++ b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl @@ -0,0 +1,120 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiAtmo + +############################################################################### +# Entropy conservation for the spherical shallow water equations in Cartesian +# form obtained through the projection of the momentum onto the divergence-free +# tangential contravariant vectors + +equations = ShallowWaterEquations3D(gravity_constant = 9.81) + +# Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux +polydeg = 3 +volume_flux = flux_wintermeyer_etal # flux_fjordholm_etal +surface_flux = flux_wintermeyer_etal # flux_fjordholm_etal #flux_lax_friedrichs +solver = DGSEM(polydeg = polydeg, + surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# 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::ShallowWaterEquations3D) + # Gaussian density + rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) + + # 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, 0), equations) +end + +initial_condition = initial_condition_advection_sphere + +mesh = P4estMeshCubedSphere2D(5, 2.0, polydeg = polydeg, + initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + metric_terms = MetricTermsInvariantCurl(), + source_terms = source_terms_lagrange_multiplier) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +ode = semidiscretize(semi, tspan) + +# Clean the initial condition +for element in eachelement(solver, semi.cache) + for j in eachnode(solver), i in eachnode(solver) + u0 = Trixi.wrap_array(ode.u0, semi) + + contravariant_normal_vector = Trixi.get_contravariant_vector(3, + semi.cache.elements.contravariant_vectors, + i, j, element) + clean_solution_lagrange_multiplier!(u0[:, i, j, element], equations, + contravariant_normal_vector) + end +end + +# 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,), + extra_analysis_integrals = (waterheight, energy_total)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +# 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, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl new file mode 100644 index 0000000..bd03404 --- /dev/null +++ b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl @@ -0,0 +1,109 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiAtmo +############################################################################### +# Entropy consistency test for the spherical shallow water equations in Cartesian +# form using the standard DGSEM with LLF dissipation + +equations = ShallowWaterEquations3D(gravity_constant = 9.81) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs 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::ShallowWaterEquations3D) + # Gaussian density + rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) + + # 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, 0), equations) +end + +# Source term function to apply the "standard" Lagrange multiplier to the semi-discretization +# We call source_terms_lagrange_multiplier, but pass x as the normal direction, as this is the +# standard way to compute the Lagrange multipliers. +function source_terms_lagrange_multiplier_standard(u, du, x, t, + equations::ShallowWaterEquations3D, + contravariant_normal_vector) + return source_terms_lagrange_multiplier(u, du, x, t, equations, x) +end + +initial_condition = initial_condition_advection_sphere + +mesh = P4estMeshCubedSphere2D(5, 2.0, polydeg = polydeg, + initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_lagrange_multiplier_standard) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +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 = 10, + save_analysis = true, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight, energy_total)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +# 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, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index e6c4608..cf5f7f0 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -15,18 +15,25 @@ using StrideArrays: PtrArray using StaticArrayInterface: static_size using LinearAlgebra: norm using Reexport: @reexport +using LoopVectorization: @turbo @reexport using StaticArrays: SVector include("auxiliary/auxiliary.jl") include("equations/equations.jl") 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") -export CompressibleMoistEulerEquations2D +export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D export flux_chandrashekar, flux_LMARS +export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal, + lake_at_rest_error, source_terms_lagrange_multiplier, + clean_solution_lagrange_multiplier! +export P4estCubedSphere2D, MetricTermsCrossProduct, MetricTermsInvariantCurl export examples_dir end # module TrixiAtmo diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl new file mode 100644 index 0000000..822ffef --- /dev/null +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -0,0 +1,38 @@ +# 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) + # to avoid a division by zero if the speed vanishes everywhere, + # e.g. for steady-state linear advection + max_scaled_speed = nextfloat(zero(t)) + + @unpack contravariant_vectors, inverse_jacobian = cache.elements + + 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) + 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, + j, + element) + lambda1_transformed = abs(Ja11 * lambda1 + Ja12 * lambda2 + Ja13 * lambda3) + Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, i, + j, + element) + lambda2_transformed = abs(Ja21 * lambda1 + Ja22 * lambda2 + Ja23 * lambda3) + + inv_jacobian = abs(inverse_jacobian[i, j, element]) + + max_lambda1 = max(max_lambda1, lambda1_transformed * inv_jacobian) + max_lambda2 = max(max_lambda2, lambda2_transformed * inv_jacobian) + end + + max_scaled_speed = max(max_scaled_speed, max_lambda1 + max_lambda2) + end + + return 2 / (nnodes(dg) * max_scaled_speed) +end diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a833960..c8d0119 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -3,3 +3,4 @@ using Trixi: AbstractEquations abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("compressible_moist_euler_2d_lucas.jl") +include("shallow_water_3d.jl") diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl new file mode 100644 index 0000000..e7e8e19 --- /dev/null +++ b/src/equations/shallow_water_3d.jl @@ -0,0 +1,386 @@ +# 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 + +@doc raw""" + ShallowWaterEquations3D(; gravity, H0 = 0) + +Shallow water equations (SWE) in three space dimensions in conservation form (with constant bottom topography). +The equations are given by +```math +\begin{aligned} + \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h v_1) + + \frac{\partial}{\partial y}(h v_2) + \frac{\partial}{\partial z}(h v_3) &= 0 \\ + \frac{\partial}{\partial t}(h v_1) + \frac{\partial}{\partial x}\left(h v_1^2 + \frac{g}{2}h^2\right) + + \frac{\partial}{\partial y}(h v_1 v_2) + \frac{\partial}{\partial z}(h v_1 v_3) &= 0 \\ + \frac{\partial}{\partial t}(h v_2) + \frac{\partial}{\partial x}(h v_1 v_2) + + \frac{\partial}{\partial y}\left(h v_2^2 + \frac{g}{2}h^2\right) + \frac{\partial}{\partial z}(h v_2 v_3) &= 0 \\ + \frac{\partial}{\partial t}(h v_3) + \frac{\partial}{\partial x}(h v_1 v_3) + + \frac{\partial}{\partial y}(h v_2 v_3) + \frac{\partial}{\partial z}\left(h v_3^2 + \frac{g}{2}h^2\right) &= 0. +\end{aligned} +``` +The unknown quantities of the SWE are the water height ``h`` and the velocities ``\mathbf{v} = (v_1, v_2, v_3)^T``. +The gravitational constant is denoted by `g`. + +The 3D Shallow Water Equations (SWE) extend the 2D SWE to model shallow water flows on 2D manifolds embedded within 3D space. +To confine the flow to the 2D manifold, a source term incorporating a Lagrange multiplier is applied. +This term effectively removes momentum components that are normal to the manifold, ensuring the flow remains +constrained within the 2D surface. + +The additional quantity ``H_0`` is also available to store a reference value for the total water height that +is useful to set initial conditions or test the "lake-at-rest" well-balancedness. + +In addition to the unknowns, Trixi.jl currently stores the bottom topography values at the approximation points +despite being fixed in time. This is done for convenience of computing the bottom topography gradients +on the fly during the approximation as well as computing auxiliary quantities like the total water height ``H`` +or the entropy variables. +This affects the implementation and use of these equations in various ways: +* The flux values corresponding to the bottom topography must be zero. +* The bottom topography values must be included when defining initial conditions, boundary conditions or + source terms. +* [`AnalysisCallback`](https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#Trixi.AnalysisCallback) analyzes this variable. +* Trixi.jl's visualization tools will visualize the bottom topography by default. + +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 +""" +struct ShallowWaterEquations3D{RealT <: Real} <: + Trixi.AbstractShallowWaterEquations{3, 5} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height +end + +# Allow for flexibility to set the gravitational constant within an elixir depending on the +# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. +# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" +# well-balancedness test cases. +function ShallowWaterEquations3D(; gravity_constant, H0 = zero(gravity_constant)) + T = promote_type(typeof(gravity_constant), typeof(H0)) + ShallowWaterEquations3D(gravity_constant, H0) +end + +Trixi.have_nonconservative_terms(::ShallowWaterEquations3D) = False() # Deactivate non-conservative terms for the moment... +Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquations3D) = ("h", "h_v1", "h_v2", + "h_v3", "b") +# Note, we use the total water height, H = h + b, as the first primitive variable for easier +# visualization and setting initial conditions +Trixi.varnames(::typeof(cons2prim), ::ShallowWaterEquations3D) = ("H", "v1", "v2", "v3", + "b") + +# Calculate 1D flux for a single point +# Note, the bottom topography has no flux +@inline function Trixi.flux(u, orientation::Integer, equations::ShallowWaterEquations3D) + h, h_v1, h_v2, h_v3, _ = u + v1, v2, v3 = velocity(u, equations) + + p = 0.5 * equations.gravity * h^2 + if orientation == 1 + f1 = h_v1 + f2 = h_v1 * v1 + p + f3 = h_v1 * v2 + f4 = h_v1 * v3 + elseif orientation == 2 + f1 = h_v2 + f2 = h_v2 * v1 + f3 = h_v2 * v2 + p + f4 = h_v2 * v3 + else # orientation == 3 + f1 = h_v3 + f2 = h_v3 * v1 + f3 = h_v3 * v2 + f4 = h_v3 * v3 + p + end + return SVector(f1, f2, f3, f4, zero(eltype(u))) +end + +# Calculate 1D flux for a single point in the normal direction +# Note, this directional vector is not normalized and the bottom topography has no flux +@inline function Trixi.flux(u, normal_direction::AbstractVector, + equations::ShallowWaterEquations3D) + h = waterheight(u, equations) + v1, v2, v3 = velocity(u, equations) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + + v3 * normal_direction[3] + h_v_normal = h * v_normal + p = 0.5 * equations.gravity * h^2 + + f1 = h_v_normal + f2 = h_v_normal * v1 + p * normal_direction[1] + f3 = h_v_normal * v2 + p * normal_direction[2] + f4 = h_v_normal * v3 + p * normal_direction[3] + return SVector(f1, f2, f3, f4, zero(eltype(u))) +end + +""" + flux_wintermeyer_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterEquations2D) + +Total energy conservative (mathematical entropy for shallow water equations) split form. +When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. + +Further details are available in Theorem 1 of the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +""" +@inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterEquations3D) + # Unpack left and right state + h_ll, h_v1_ll, h_v2_ll, h_v3_ll, _ = u_ll + h_rr, h_v1_rr, h_v2_rr, h_v3_rr, _ = u_rr + + # Get the velocities on either side + v1_ll, v2_ll, v3_ll = velocity(u_ll, equations) + v1_rr, v2_rr, v3_rr = velocity(u_rr, equations) + + # Average each factor of products in flux + h_v1_avg = 0.5 * (h_v1_ll + h_v1_rr) + h_v2_avg = 0.5 * (h_v2_ll + h_v2_rr) + h_v3_avg = 0.5 * (h_v3_ll + h_v3_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (v3_ll + v3_rr) + p_avg = 0.5 * equations.gravity * h_ll * h_rr + + # Calculate fluxes depending on normal_direction + f1 = h_v1_avg * normal_direction[1] + h_v2_avg * normal_direction[2] + + h_v3_avg * normal_direction[3] + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = f1 * v3_avg + p_avg * normal_direction[3] + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_fjordholm_etal(u_ll, u_rr, orientation, + equations::ShallowWaterEquations1D) + +Total energy conservative (mathematical entropy for shallow water equations). When the bottom topography +is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. +For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). + +Details are available in Eq. (4.1) in the paper: +- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography + [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) +""" +@inline function Trixi.flux_fjordholm_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::ShallowWaterEquations3D) + # Unpack left and right state + h_ll = waterheight(u_ll, equations) + v1_ll, v2_ll, v3_ll = velocity(u_ll, equations) + h_rr = waterheight(u_rr, equations) + v1_rr, v2_rr, v3_rr = velocity(u_rr, equations) + + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + # Average each factor of products in flux + h_avg = 0.5 * (h_ll + h_rr) + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v3_avg = 0.5 * (v3_ll + v3_rr) + h2_avg = 0.5 * (h_ll^2 + h_rr^2) + p_avg = 0.5 * equations.gravity * h2_avg + v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr) + + # Calculate fluxes depending on normal_direction + f1 = h_avg * v_dot_n_avg + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + f4 = f1 * v3_avg + p_avg * normal_direction[3] + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + source_terms_lagrange_multiplier(u, du, x, t, + equations::ShallowWaterEquations3D, + normal_direction) + +Source term function to apply a Lagrange multiplier to the semi-discretization +in order to constrain the momentum to a 2D manifold. + +The vector normal_direction is perpendicular to the 2D manifold. By default, +this is the normal contravariant basis vector. +""" +function source_terms_lagrange_multiplier(u, du, x, t, + equations::ShallowWaterEquations3D, + normal_direction) + x_dot_div_f = (normal_direction[1] * du[2] + + normal_direction[2] * du[3] + + normal_direction[3] * du[4]) / + sum(normal_direction .^ 2) + + s2 = -normal_direction[1] * x_dot_div_f + s3 = -normal_direction[2] * x_dot_div_f + s4 = -normal_direction[3] * x_dot_div_f + + return SVector(0, s2, s3, s4, 0) +end + +""" + clean_solution_lagrange_multiplier!(u, equations::ShallowWaterEquations3D, normal_direction) + +Function to apply Lagrange multiplier discretely to the solution in order to constrain +the momentum to a 2D manifold. + +The vector normal_direction is perpendicular to the 2D manifold. By default, +this is the normal contravariant basis vector. +""" +function clean_solution_lagrange_multiplier!(u, equations::ShallowWaterEquations3D, + normal_direction) + x_dot_div_f = (normal_direction[1] * u[2] + + normal_direction[2] * u[3] + + normal_direction[3] * u[4]) / + sum(normal_direction .^ 2) + + u[2] -= normal_direction[1] * x_dot_div_f + u[3] -= normal_direction[2] * x_dot_div_f + u[4] -= normal_direction[3] * x_dot_div_f +end + +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::ShallowWaterEquations3D) + # Extract and compute the velocities in the normal direction + v1_ll, v2_ll, v3_ll = velocity(u_ll, equations) + v1_rr, v2_rr, v3_rr = velocity(u_rr, equations) + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + # 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) + + # 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) +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterEquations3D) + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + diss = -0.5 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) +end + +@inline function Trixi.max_abs_speeds(u, equations::ShallowWaterEquations3D) + h = waterheight(u, equations) + v1, v2, v3 = velocity(u, equations) + + c = sqrt(equations.gravity * h) + return abs(v1) + c, abs(v2) + c, abs(v3) + c +end + +# Helper function to extract the velocity vector from the conservative variables +@inline function velocity(u, equations::ShallowWaterEquations3D) + h, h_v1, h_v2, h_v3, _ = u + + v1 = h_v1 / h + v2 = h_v2 / h + v3 = h_v3 / h + return SVector(v1, v2, v3) +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::ShallowWaterEquations3D) + h, _, _, _, b = u + + H = h + b + v1, v2, v3 = velocity(u, equations) + return SVector(H, v1, v2, v3, b) +end + +# Convert conservative variables to entropy +# Note, only the first three are the entropy variables, the fourth entry still +# just carries the bottom topography values for convenience +@inline function Trixi.cons2entropy(u, equations::ShallowWaterEquations3D) + h, h_v1, h_v2, h_v3, b = u + + v1, v2, v3 = velocity(u, equations) + v_square = v1^2 + v2^2 + v3^2 + + w1 = equations.gravity * (h + b) - 0.5 * v_square + w2 = v1 + w3 = v2 + w4 = v3 + return SVector(w1, w2, w3, w4, b) +end + +# Convert entropy variables to conservative +@inline function Trixi.entropy2cons(w, equations::ShallowWaterEquations3D) + w1, w2, w3, w4, b = w + + h = (w1 + 0.5 * (w2^2 + w3^2 + w4^2)) / equations.gravity - b + h_v1 = h * w2 + h_v2 = h * w3 + h_v3 = h * w4 + return SVector(h, h_v1, h_v2, h_v3, b) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, equations::ShallowWaterEquations3D) + H, v1, v2, v3, b = prim + + h = H - b + h_v1 = h * v1 + h_v2 = h * v2 + h_v3 = h * v3 + return SVector(h, h_v1, h_v2, h_v3, b) +end + +@inline function waterheight(u, equations::ShallowWaterEquations3D) + return u[1] +end + +@inline function pressure(u, equations::ShallowWaterEquations3D) + h = waterheight(u, equations) + p = 0.5 * equations.gravity * h^2 + return p +end + +# Entropy function for the shallow water equations is the total energy +@inline function Trixi.entropy(cons, equations::ShallowWaterEquations3D) + energy_total(cons, equations) +end + +# Calculate total energy for a conservative state `cons` +@inline function energy_total(cons, equations::ShallowWaterEquations3D) + h, h_v1, h_v2, h_v3, b = cons + + e = (h_v1^2 + h_v2^2 + h_v3^2) / (2 * h) + 0.5 * equations.gravity * h^2 + + equations.gravity * h * b + return e +end + +# Calculate kinetic energy for a conservative state `cons` +@inline function energy_kinetic(u, equations::ShallowWaterEquations3D) + h, h_v1, h_v2, h_v3, _ = u + return (h_v1^2 + h_v2^2 + h_v3^2) / (2 * h) +end + +# Calculate potential energy for a conservative state `cons` +@inline function energy_internal(cons, equations::ShallowWaterEquations3D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) +end +end # @muladd diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl new file mode 100644 index 0000000..734d430 --- /dev/null +++ b/src/semidiscretization/semidiscretization.jl @@ -0,0 +1,20 @@ +""" + MetricTermsCrossProduct() + +Struct used for multiple dispatch on functions that compute the metric terms. +When the argument `metric_terms` is of type `MetricTermsCrossProduct`, the +contravariant vectors are computed using the cross-product form. +""" +struct MetricTermsCrossProduct end + +""" + MetricTermsInvariantCurl() + +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 +""" +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 af736e5..7421c99 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl @@ -12,8 +12,9 @@ function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{2}, # `RealT` is used as real type for node locations etc. # while `uEltype` is used as element type of solutions etc. RealT = real(solver), uEltype = RealT, - initial_cache = NamedTuple()) - cache = (; Trixi.create_cache(mesh, equations, solver, RealT, uEltype)..., + initial_cache = NamedTuple(), + metric_terms = MetricTermsCrossProduct()) + cache = (; Trixi.create_cache(mesh, equations, solver, RealT, metric_terms, uEltype)..., initial_cache...) _boundary_conditions = Trixi.digest_boundary_conditions(boundary_conditions, mesh, solver, diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl index 8aa6583..3199243 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl @@ -47,6 +47,7 @@ function Trixi.init_elements(mesh::Union{P4estMesh{2, 3, RealT}, T8codeMesh{2}}, equations::AbstractEquations{3}, basis, + metric_terms, ::Type{uEltype}) where {RealT <: Real, uEltype <: Real} nelements = Trixi.ncells(mesh) @@ -103,14 +104,15 @@ function Trixi.init_elements(mesh::Union{P4estMesh{2, 3, RealT}, _inverse_jacobian, _surface_flux_values) - init_elements_2d_manifold_in_3d!(elements, mesh, basis) + init_elements_2d_manifold_in_3d!(elements, mesh, basis, metric_terms) return elements end function init_elements_2d_manifold_in_3d!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) + basis::LobattoLegendreBasis, + metric_terms) (; node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian) = elements @@ -126,7 +128,8 @@ function init_elements_2d_manifold_in_3d!(elements, element, jacobian_matrix, node_coordinates, - basis) + basis, + metric_terms) # Compute the inverse Jacobian as the norm of the cross product of the covariant vectors for j in eachnode(basis), i in eachnode(basis) @@ -207,12 +210,20 @@ function calc_node_coordinates_2d_shell!(node_coordinates, return node_coordinates end +# Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, +# using eq (12) of : +# Giraldo, F. X. (2001). A spectral element shallow water model on spherical geodesic grids. +# International Journal for Numerical Methods in Fluids, 35(8), 869-901. +# https://doi.org/10.1002/1097-0363(20010430)35:8<869::AID-FLD116>3.0.CO;2-S +# This is equivalent to the cross-product form, but we end up with three contravariant vectors +# because there are three space dimensions. # This only works for a sphere function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArray{<:Any, 5}, element, jacobian_matrix, node_coordinates, - basis::LobattoLegendreBasis) + basis::LobattoLegendreBasis, + metric_terms::MetricTermsCrossProduct) @unpack derivative_matrix = basis for j in eachnode(basis), i in eachnode(basis) @@ -262,4 +273,146 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr return contravariant_vectors end + +# Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, +# using the invariant curl form for a cubed sphere. +# In the cubed sphere, the node coordinates are first mapped with polynomials in ξ and η, and then we add +# a linear ζ dependency, i.e.: +# xₗ(ξ, η, ζ) = ζ * xₗ(ξ, η) with ζ = 1 at the sphere surface. +# +# As a result, the covariant vectors with respect to ζ are xₗ_ζ = xₗ +# +# We compute the metric terms in two steps +# 1. Compute the 3D contravariant vectors, Jaⁱ=J*grad(ξⁱ), using the curl invariant form and xₗ_ζ = xₗ +# 2. Convert the 3D mapping Jacobian determinant J:=a_k.(a_i x a_j) to 2D J_2D=||a_i x a_j|| +## 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 + +function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArray{<:Any, + 5}, + element, + jacobian_matrix, node_coordinates, + basis::LobattoLegendreBasis, + metric_terms::MetricTermsInvariantCurl) + @unpack derivative_matrix = basis + + # 1. Compute the 3D contravariant vectors, Jaⁱₙ=J*grad(ξ), using the curl invariant form and xₗ_ζ = xₗ + # The general form is + # Jaⁱₙ = 0.5 * ( ∇ × (Xₘ ∇ Xₗ - Xₗ ∇ Xₘ) )ᵢ where (n, m, l) cyclic and ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ + + for n in 1:3 + # (n, m, l) cyclic + m = (n % 3) + 1 + l = ((n + 1) % 3) + 1 + + # Calculate Ja¹ₙ = 0.5 * [ (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_η - (Xₘ Xₗ_η - Xₗ Xₘ_η)_ζ ] + # For each of these, the first and second summand are computed in separate loops + # for performance reasons. + + # First summand 0.5 * (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_η + @turbo for j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to j-dimension to differentiate wrt η + result += 0.5 * derivative_matrix[j, ii] * + (node_coordinates[m, i, ii, element] * + node_coordinates[l, i, ii, element] - + node_coordinates[l, i, ii, element] * + node_coordinates[m, i, ii, element]) + end + + contravariant_vectors[n, 1, i, j, element] = result + end + + # Second summand -0.5 * (Xₘ Xₗ_η - Xₗ Xₘ_η)_ζ + @turbo for j in eachnode(basis), i in eachnode(basis) + # Due to the spherical-shell mapping, xₗ(ξ, η, ζ) = ζ * xₗ(ξ, η), we have + # 0.5 d/dζ (ζ^2 F(ξ, η)) = ζ F(ξ, η) + + result = (node_coordinates[m, i, j, element] * + jacobian_matrix[l, 2, i, j, element] - + node_coordinates[l, i, j, element] * + jacobian_matrix[m, 2, i, j, element]) + + contravariant_vectors[n, 1, i, j, element] -= result + end + + # Calculate Ja²ₙ = 0.5 * [ (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_ζ - (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_ξ ] + + # First summand 0.5 * (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_ζ + @turbo for j in eachnode(basis), i in eachnode(basis) + # Due to the spherical-shell mapping, xₗ(ξ, η, ζ) = ζ * xₗ(ξ, η), we have + # 0.5 d/dζ (ζ^2 F(ξ, η)) = ζ F(ξ, η) + + result = (node_coordinates[m, i, j, element] * + jacobian_matrix[l, 1, i, j, element] - + node_coordinates[l, i, j, element] * + jacobian_matrix[m, 1, i, j, element]) + + contravariant_vectors[n, 2, i, j, element] = result + end + + # Second summand -0.5 * (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_ξ + @turbo for j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to i-dimension to differentiate wrt ξ + result += 0.5 * derivative_matrix[i, ii] * + (node_coordinates[m, ii, j, element] * + node_coordinates[l, ii, j, element] - + node_coordinates[l, ii, j, element] * + node_coordinates[m, ii, j, element]) + end + + contravariant_vectors[n, 2, i, j, element] -= result + end + + # Calculate Ja³ₙ = 0.5 * [ (Xₘ Xₗ_η - Xₗ Xₘ_η)_ξ - (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_η ] + + # First summand 0.5 * (Xₘ Xₗ_η - Xₗ Xₘ_η)_ξ + @turbo for j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to i-dimension to differentiate wrt ξ + result += 0.5 * derivative_matrix[i, ii] * + (node_coordinates[m, ii, j, element] * + jacobian_matrix[l, 2, ii, j, element] - + node_coordinates[l, ii, j, element] * + jacobian_matrix[m, 2, ii, j, element]) + end + + contravariant_vectors[n, 3, i, j, element] = result + end + + # Second summand -0.5 * (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_η + @turbo for j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to j-dimension to differentiate wrt η + result += 0.5 * derivative_matrix[j, ii] * + (node_coordinates[m, i, ii, element] * + jacobian_matrix[l, 1, i, ii, element] - + node_coordinates[l, i, ii, element] * + jacobian_matrix[m, 1, i, ii, element]) + end + + contravariant_vectors[n, 3, i, j, element] -= result + end + end + + # 2. Convert the 3D mapping Jacobian determinant J:=a_k.(a_i x a_j) to 2D J_2D=||a_i x a_j|| + for j in eachnode(basis), i in eachnode(basis) + factor = 1 / norm(node_coordinates[:, i, j, element]) # = 1 / norm(jacobian_matrix[:, 3, i, j, element]) + for n in 1:3, m in 1:3 + contravariant_vectors[m, n, i, j, element] *= factor + end + end + + return contravariant_vectors +end end # @muladd diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 0d3e8e7..ea6f835 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -1,3 +1,27 @@ +# This method is called when a SemidiscretizationHyperbolic is constructed. +# It constructs the basic `cache` used throughout the simulation to compute +# the RHS etc. +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 + # in case someone has tampered with the `p4est` after creating the mesh + Trixi.balance!(mesh) + + elements = Trixi.init_elements(mesh, equations, dg.basis, metric_terms, uEltype) + interfaces = Trixi.init_interfaces(mesh, equations, dg.basis, elements) + boundaries = Trixi.init_boundaries(mesh, equations, dg.basis, elements) + mortars = Trixi.init_mortars(mesh, equations, dg.basis, elements) + + cache = (; elements, interfaces, boundaries, mortars) + + # Add specialized parts of the cache required to compute the volume integral etc. + cache = (; cache..., + Trixi.create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) + cache = (; cache..., Trixi.create_cache(mesh, equations, dg.mortar, uEltype)...) + + 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}() 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 dfb6b9d..4264c92 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,4 +1,77 @@ -# Weak-form kernel for 3D equations solved in 2D manifolds + +function Trixi.rhs!(du, u, t, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::AbstractEquations{3}, + boundary_conditions, source_terms::Source, + dg::DG, cache) where {Source} + + # Reset du + Trixi.@trixi_timeit Trixi.timer() "reset ∂u/∂t" Trixi.reset_du!(du, dg, cache) + + # Calculate volume integral + Trixi.@trixi_timeit Trixi.timer() "volume integral" begin + Trixi.calc_volume_integral!(du, u, mesh, + Trixi.have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) + end + + # Prolong solution to interfaces + Trixi.@trixi_timeit Trixi.timer() "prolong2interfaces" begin + Trixi.prolong2interfaces!(cache, u, mesh, equations, + dg.surface_integral, dg) + end + + # 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, + dg.surface_integral, dg, cache) + end + + # Prolong solution to boundaries + Trixi.@trixi_timeit Trixi.timer() "prolong2boundaries" begin + Trixi.prolong2boundaries!(cache, u, mesh, equations, + dg.surface_integral, dg) + end + + # Calculate boundary fluxes + Trixi.@trixi_timeit Trixi.timer() "boundary flux" begin + Trixi.calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations, + dg.surface_integral, dg) + end + + # Prolong solution to mortars + Trixi.@trixi_timeit Trixi.timer() "prolong2mortars" begin + Trixi.prolong2mortars!(cache, u, mesh, equations, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes + Trixi.@trixi_timeit Trixi.timer() "mortar flux" begin + Trixi.calc_mortar_flux!(cache.elements.surface_flux_values, mesh, + Trixi.have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Calculate surface integrals + Trixi.@trixi_timeit Trixi.timer() "surface integral" begin + Trixi.calc_surface_integral!(du, u, mesh, equations, + dg.surface_integral, dg, cache) + end + + # Apply Jacobian from mapping to reference element + Trixi.@trixi_timeit Trixi.timer() "Jacobian" Trixi.apply_jacobian!(du, mesh, equations, + dg, cache) + + # Calculate source terms + Trixi.@trixi_timeit Trixi.timer() "source terms" begin + calc_sources_2d_manifold_in_3d!(du, u, t, source_terms, equations, dg, cache) + end + + return nothing +end + +# Weak-form kernel for 3D equations solved on 2D manifolds @inline function Trixi.weak_form_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, @@ -43,3 +116,32 @@ return nothing end + +function calc_sources_2d_manifold_in_3d!(du, u, t, source_terms::Nothing, + equations::AbstractEquations{3}, dg::DG, cache) + return nothing +end + +function calc_sources_2d_manifold_in_3d!(du, u, t, source_terms, + equations::AbstractEquations{3}, dg::DG, cache) + @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements + + 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) + du_local = Trixi.get_node_vars(du, equations, dg, i, j, element) + x_local = Trixi.get_node_coords(node_coordinates, equations, dg, + i, j, element) + contravariant_normal_vector = Trixi.get_contravariant_vector(3, + contravariant_vectors, + i, j, + element) * + inverse_jacobian[i, j, element] + source = source_terms(u_local, du_local, x_local, t, equations, + contravariant_normal_vector) + Trixi.add_to_node_vars!(du, source, equations, dg, i, j, element) + end + end + + return nothing +end diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index b9f211d..d39e464 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -1,3 +1,3 @@ +include("dgsem_p4est/containers_2d_manifold_in_3d.jl") include("dgsem_p4est/dg.jl") include("dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl") -include("dgsem_p4est/containers_2d_manifold_in_3d.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 541cfd8..38d3e81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,4 +20,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection" include("test_spherical_advection.jl") end + + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "shallow_water_3d" + include("test_3d_shallow_water.jl") + end end diff --git a/test/test_3d_shallow_water.jl b/test/test_3d_shallow_water.jl new file mode 100644 index 0000000..0287a35 --- /dev/null +++ b/test/test_3d_shallow_water.jl @@ -0,0 +1,91 @@ +module TestSphericalAdvection + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") + +@trixiatmo_testset "elixir_shallowwater_cubed_sphere_shell_EC_correction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_cubed_sphere_shell_EC_correction.jl"), + l2=[ + 7.23800458e-02, + 9.98871590e-02, + 4.55606969e-02, + 3.17422064e-02, + 0.00000000e+00 + ], + linf=[ + 1.05686060e+00, + 1.04413842e+00, + 3.12374228e-01, + 3.19064636e-01, + 0.00000000e+00 + ]) + # 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_cubed_sphere_shell_EC_projection" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_cubed_sphere_shell_EC_projection.jl"), + l2=[ + 7.25131271e-02, + 1.01017589e-01, + 4.39697415e-02, + 3.08898940e-02, + 0.00000000e+00 + ], + linf=[ + 1.06007536e+00, + 1.05786719e+00, + 3.93826726e-01, + 2.34129714e-01, + 0.00000000e+00 + ]) + # 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_cubed_sphere_shell_standard" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_cubed_sphere_shell_standard.jl"), + l2=[ + 6.86661496e-02, + 9.43310712e-02, + 3.26202302e-02, + 2.02514293e-02, + 0.00000000e+00 + ], + linf=[ + 1.04610420e+00, + 1.05891219e+00, + 1.19809349e-01, + 5.54354032e-02, + 0.00000000e+00 + ]) + # 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 diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index bb735a5..6aa252e 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -11,18 +11,18 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_spherical_advection_cartesian.jl"), l2=[ - 8.44505073e-03, - 8.23414117e-03, - 1.84210648e-03, + 8.44498914e-03, + 8.23407970e-03, + 1.84210216e-03, 0.00000000e+00, - 6.44302430e-02 + 6.44302432e-02 ], linf=[ - 9.48950488e-02, - 9.64811952e-02, - 1.37453400e-02, + 9.48946345e-02, + 9.64807833e-02, + 1.37450721e-02, 0.00000000e+00, - 4.09322999e-01 + 0.40932295554303133 ], element_local_mapping=false) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -38,18 +38,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_spherical_advection_cartesian.jl"), l2=[ - 8.97501204e-03, - 8.74316738e-03, - 1.99380928e-03, - 0.00000000e+00, - 6.45266612e-02 + 0.00897499647958351, + 0.008743150662795171, + 0.001993813300529814, + 0.0, + 0.06452666141326718 ], linf=[ - 1.01363241e-01, - 1.03082705e-01, - 1.41424723e-02, - 0.00000000e+00, - 4.10471741e-01 + 0.10136306017376362, + 0.10308252591097666, + 0.014142319489653277, + 0.0, + 0.41047171697830986 ], element_local_mapping=true) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities)