From c63f51498e494fe6983d1493cd8349be18ff0bd9 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Thu, 12 Dec 2024 14:54:20 +0100 Subject: [PATCH 01/21] add additional aux vars --- src/equations/equations.jl | 63 ++++++++++++++----- .../containers_2d_manifold_in_3d_covariant.jl | 45 +++++++++---- 2 files changed, 78 insertions(+), 30 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8cc5393..8fa4bbc 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -138,27 +138,56 @@ function global2contravariant end # For the covariant form, the auxiliary variables are the the NDIMS*NDIMS_AMBIENT entries # of the exact covariant basis matrix, the NDIMS*NDIMS_AMBIENT entries of the exact -# contravariant basis matrix, and the area element. In the future, we will add other terms -# such as Christoffel symbols. +# contravariant basis matrix, the covariant and contravariant metric tensor components, the +# Christoffel symbols of the second kind, and the area element. @inline have_aux_node_vars(::AbstractCovariantEquations) = True() -@inline n_aux_node_vars(::AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT}) where {NDIMS, -NDIMS_AMBIENT} = 2 * NDIMS * NDIMS_AMBIENT + 1 +@inline function n_aux_node_vars(::AbstractCovariantEquations{NDIMS, + NDIMS_AMBIENT}) where { + NDIMS, + NDIMS_AMBIENT + } + nvars_basis_covariant = NDIMS_AMBIENT * NDIMS + nvars_basis_contravariant = NDIMS * NDIMS_AMBIENT + nvars_area_element = 1 + nvars_metric_covariant = NDIMS*(NDIMS+1) ÷ 2 + nvars_metric_contravariant = NDIMS*(NDIMS+1) ÷ 2 + nvars_christoffel = NDIMS * NDIMS*(NDIMS+1) ÷ 2 + + return nvars_basis_covariant + + nvars_basis_contravariant + + nvars_area_element + + nvars_metric_covariant + + nvars_metric_contravariant + + nvars_christoffel +end # Return auxiliary variable names for 2D covariant form @inline function auxvarnames(::AbstractCovariantEquations{2}) - return ("basis_covariant[1,1]", - "basis_covariant[2,1]", - "basis_covariant[3,1]", - "basis_covariant[1,2]", - "basis_covariant[2,2]", - "basis_covariant[3,2]", - "basis_contravariant[1,1]", - "basis_contravariant[2,1]", - "basis_contravariant[3,1]", - "basis_contravariant[1,2]", - "basis_contravariant[2,2]", - "basis_contravariant[3,2]", - "area_element") + return ("basis_covariant[1,1]", # e₁ ⋅ a₁ + "basis_covariant[2,1]", # e₂ ⋅ a₁ + "basis_covariant[3,1]", # e₃ ⋅ a₁ + "basis_covariant[1,2]", # e₁ ⋅ a₂ + "basis_covariant[2,2]", # e₂ ⋅ a₂ + "basis_covariant[3,2]", # e₃ ⋅ a₂ + "basis_contravariant[1,1]", # e₁ ⋅ a¹ + "basis_contravariant[2,1]", # e₂ ⋅ a¹ + "basis_contravariant[3,1]", # e₃ ⋅ a¹ + "basis_contravariant[1,2]", # e₁ ⋅ a² + "basis_contravariant[2,2]", # e₂ ⋅ a² + "basis_contravariant[3,2]", # e₃ ⋅ a² + "area_element", # √G = √(G₁₁G₂₂ - G₁₂G₂₁) = a₁ × a₂ + "metric_covariant[1,1]", # G₁₁ + "metric_covariant[1,2]", # G₁₂ = G₂₁ + "metric_covariant[2,2]", # G₂₂ + "metric_contravariant[1,1]", # G¹¹ + "metric_contravariant[1,2]", # G¹² = G²¹ + "metric_contravariant[2,2]", # G²² + "christoffel[1,1,1]", # Γ¹₁₁ + "christoffel[1,1,2]", # Γ¹₁₂ = Γ¹₂₁ + "christoffel[1,2,2]", # Γ¹₂₂ + "christoffel[2,1,1]", # Γ²₁₁ + "christoffel[2,1,2]", # Γ²₁₂ = Γ²₂₁ + "christoffel[2,2,2]") # Γ²₂₂ end # Extract the covariant basis vectors a_i from the auxiliary variables as a matrix A, diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl index 07d595d..174dbf5 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl @@ -186,21 +186,40 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, # Compute the auxiliary metric information at each node for j in eachnode(dg), i in eachnode(dg) - - # Calculate the covariant basis in the desired global coordinate system - A = calc_basis_covariant(v1, v2, v3, v4, - dg.basis.nodes[i], dg.basis.nodes[j], - radius, equations.global_coordinate_system) - aux_node_vars[1:(NDIMS * NDIMS_AMBIENT), i, j, element] = SVector(A) - - G = A' * A # Covariant metric tensor (not used for now) - - # Contravariant basis - aux_node_vars[(NDIMS * NDIMS_AMBIENT + 1):(2 * NDIMS * NDIMS_AMBIENT), - i, j, element] = SVector(inv(G) * A') + + + # Covariant basis in the desired global coordinate system as columns of a matrix + basis_covariant = calc_basis_covariant(v1, v2, v3, v4, + dg.basis.nodes[i], dg.basis.nodes[j], + radius, + equations.global_coordinate_system) + start_ind, end_ind = 1, length(basis_covariant) + aux_node_vars[start_ind:end_ind, i, j, element] = SVector(basis_covariant) + start_ind = end_ind + 1 + + # Metric tensor + metric_covariant = basis_covariant' * basis_covariant + metric_contravariant = inv(basis_covariant' * basis_covariant) + + # Contravariant basis vectors as rows of a matrix + basis_contravariant = metric_contravariant * basis_covariant' + end_ind = start_ind+length(basis_contravariant)-1 + aux_node_vars[start_ind:end_ind, i, j, element] = SVector(basis_contravariant) + start_ind = end_ind + 1 # Area element - aux_node_vars[n_aux_node_vars(equations), i, j, element] = sqrt(det(G)) + area_element = sqrt(det(metric_covariant)) + aux_node_vars[start_ind, i, j, element] = area_element + start_ind = start_ind + 1 + + # Covariant metric tensor components + end_ind = start_ind + 2 + aux_node_vars[start_ind:end_ind, i, j, element] = SVector(metric_covariant[1,1], metric_covariant[1,2], metric_covariant[2,2]) + start_ind = end_ind + 1 + + # Contravariant metric tensor components + end_ind = start_ind + 2 + aux_node_vars[start_ind:end_ind, i, j, element] = SVector(metric_contravariant[1,1], metric_contravariant[1,2], metric_contravariant[2,2]) end end return nothing From 32f49a9e3ff468d113f358ccdf4f1e8bbd8495e0 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Wed, 18 Dec 2024 23:59:13 +0100 Subject: [PATCH 02/21] compute christoffel symbols --- src/callbacks_step/analysis_covariant.jl | 3 +- src/equations/equations.jl | 43 ++++++--- .../containers_2d_manifold_in_3d_covariant.jl | 96 +++++++++++++++---- 3 files changed, 108 insertions(+), 34 deletions(-) diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index 96c9c29..1c5849e 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -3,7 +3,7 @@ # For the covariant form, we want to integrate using the exact area element # √G = (det(AᵀA))^(1/2), which is stored in cache.auxiliary_variables, not the approximate -# area element used in the Cartesian formulation, which stored in cache.elements. +# area element used in the Cartesian formulation, which stored in cache.elements function Trixi.integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, @@ -65,6 +65,7 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, (; weights) = dg.basis (; node_coordinates) = cache.elements (; aux_node_vars) = cache.auxiliary_variables + # Set up data structures l2_error = zero(func(Trixi.get_node_vars(u, equations, dg, 1, 1, 1), equations)) linf_error = copy(l2_error) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8fa4bbc..7c9c635 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -138,8 +138,9 @@ function global2contravariant end # For the covariant form, the auxiliary variables are the the NDIMS*NDIMS_AMBIENT entries # of the exact covariant basis matrix, the NDIMS*NDIMS_AMBIENT entries of the exact -# contravariant basis matrix, the covariant and contravariant metric tensor components, the -# Christoffel symbols of the second kind, and the area element. +# contravariant basis matrix, the exact area element, the upper-triangular covariant and +# contravariant metric tensor components, and the upper-triangular Christoffel symbols of +# the second kind @inline have_aux_node_vars(::AbstractCovariantEquations) = True() @inline function n_aux_node_vars(::AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT}) where { @@ -149,9 +150,9 @@ function global2contravariant end nvars_basis_covariant = NDIMS_AMBIENT * NDIMS nvars_basis_contravariant = NDIMS * NDIMS_AMBIENT nvars_area_element = 1 - nvars_metric_covariant = NDIMS*(NDIMS+1) ÷ 2 - nvars_metric_contravariant = NDIMS*(NDIMS+1) ÷ 2 - nvars_christoffel = NDIMS * NDIMS*(NDIMS+1) ÷ 2 + nvars_metric_covariant = NDIMS * (NDIMS + 1) ÷ 2 + nvars_metric_contravariant = NDIMS * (NDIMS + 1) ÷ 2 + nvars_christoffel = NDIMS * NDIMS * (NDIMS + 1) ÷ 2 return nvars_basis_covariant + nvars_basis_contravariant + @@ -175,19 +176,19 @@ end "basis_contravariant[1,2]", # e₁ ⋅ a² "basis_contravariant[2,2]", # e₂ ⋅ a² "basis_contravariant[3,2]", # e₃ ⋅ a² - "area_element", # √G = √(G₁₁G₂₂ - G₁₂G₂₁) = a₁ × a₂ + "area_element", # √G = √(G₁₁G₂₂ - G₁₂G₂₁) = ||a₁ × a₂|| "metric_covariant[1,1]", # G₁₁ "metric_covariant[1,2]", # G₁₂ = G₂₁ "metric_covariant[2,2]", # G₂₂ "metric_contravariant[1,1]", # G¹¹ "metric_contravariant[1,2]", # G¹² = G²¹ "metric_contravariant[2,2]", # G²² - "christoffel[1,1,1]", # Γ¹₁₁ - "christoffel[1,1,2]", # Γ¹₁₂ = Γ¹₂₁ - "christoffel[1,2,2]", # Γ¹₂₂ - "christoffel[2,1,1]", # Γ²₁₁ - "christoffel[2,1,2]", # Γ²₁₂ = Γ²₂₁ - "christoffel[2,2,2]") # Γ²₂₂ + "christoffel_symbols[1][1,1]", # Γ¹₁₁ + "christoffel_symbols[1][1,2]", # Γ¹₁₂ = Γ¹₂₁ + "christoffel_symbols[2][2,2]", # Γ¹₂₂ + "christoffel_symbols[2][1,1]", # Γ²₁₁ + "christoffel_symbols[2][1,2]", # Γ²₁₂ = Γ²₂₁ + "christoffel_symbols[2][2,2]") # Γ²₂₂ end # Extract the covariant basis vectors a_i from the auxiliary variables as a matrix A, @@ -212,6 +213,24 @@ end return aux_vars[13] end +# Extract the covariant metric tensor components Gᵢⱼ from the auxiliary variables +@inline function metric_covariant(aux_vars, ::AbstractCovariantEquations{2}) + return SMatrix{2, 2}(aux_vars[14], aux_vars[15], + aux_vars[15], aux_vars[16]) +end + +# Extract the contravariant metric tensor components Gⁱʲ from the auxiliary variables +@inline function metric_contravariant(aux_vars, ::AbstractCovariantEquations{2}) + return SMatrix{2, 2}(aux_vars[17], aux_vars[18], + aux_vars[18], aux_vars[19]) +end + +# Extract the Christoffel symbols of the second kind Γⁱⱼₖ from the auxiliary variables +@inline function christoffel_symbols(aux_vars, ::AbstractCovariantEquations{2}) + return (SMatrix{2, 2}(aux_vars[20], aux_vars[21], aux_vars[21], aux_vars[22]), + SMatrix{2, 2}(aux_vars[23], aux_vars[24], aux_vars[24], aux_vars[25])) +end + # Numerical flux plus dissipation for abstract covariant equations as a function of the # state vectors u_ll and u_rr, as well as the auxiliary variable vectors aux_vars_ll and # aux_vars_rr, which contain the geometric information. We assume that u_ll and u_rr have diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl index 174dbf5..508e443 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl @@ -163,9 +163,6 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, (; tree_node_coordinates) = mesh (; aux_node_vars) = auxiliary_variables - NDIMS = 2 - NDIMS_AMBIENT = 3 - # Check that the degree of the mesh matches that of the solver @assert length(mesh.nodes) == nnodes(dg) @@ -186,16 +183,13 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, # Compute the auxiliary metric information at each node for j in eachnode(dg), i in eachnode(dg) - # Covariant basis in the desired global coordinate system as columns of a matrix basis_covariant = calc_basis_covariant(v1, v2, v3, v4, - dg.basis.nodes[i], dg.basis.nodes[j], - radius, - equations.global_coordinate_system) - start_ind, end_ind = 1, length(basis_covariant) - aux_node_vars[start_ind:end_ind, i, j, element] = SVector(basis_covariant) - start_ind = end_ind + 1 + dg.basis.nodes[i], dg.basis.nodes[j], + radius, + equations.global_coordinate_system) + aux_node_vars[1:6, i, j, element] = SVector(basis_covariant) # Metric tensor metric_covariant = basis_covariant' * basis_covariant @@ -203,25 +197,26 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2, # Contravariant basis vectors as rows of a matrix basis_contravariant = metric_contravariant * basis_covariant' - end_ind = start_ind+length(basis_contravariant)-1 - aux_node_vars[start_ind:end_ind, i, j, element] = SVector(basis_contravariant) - start_ind = end_ind + 1 + aux_node_vars[7:12, i, j, element] = SVector(basis_contravariant) # Area element - area_element = sqrt(det(metric_covariant)) - aux_node_vars[start_ind, i, j, element] = area_element - start_ind = start_ind + 1 + aux_node_vars[13, i, j, element] = sqrt(det(metric_covariant)) # Covariant metric tensor components - end_ind = start_ind + 2 - aux_node_vars[start_ind:end_ind, i, j, element] = SVector(metric_covariant[1,1], metric_covariant[1,2], metric_covariant[2,2]) - start_ind = end_ind + 1 + aux_node_vars[14:16, i, j, element] = SVector(metric_covariant[1, 1], + metric_covariant[1, 2], + metric_covariant[2, 2]) # Contravariant metric tensor components - end_ind = start_ind + 2 - aux_node_vars[start_ind:end_ind, i, j, element] = SVector(metric_contravariant[1,1], metric_contravariant[1,2], metric_contravariant[2,2]) + aux_node_vars[17:19, i, j, element] = SVector(metric_contravariant[1, 1], + metric_contravariant[1, 2], + metric_contravariant[2, 2]) end + + # Christoffel symbols of the second kind (aux_node_vars[20:25, :, :, element]) + calc_christoffel_symbols!(aux_node_vars, mesh, equations, dg, element) end + return nothing end @@ -292,4 +287,63 @@ end # Make zero component in the radial direction so the matrix has the right dimensions return SMatrix{3, 2}(A[1, 1], A[2, 1], 0.0f0, A[1, 2], A[2, 2], 0.0f0) end + +# Calculate Christoffel symbols approximately using the collocation derivative +function calc_christoffel_symbols!(aux_node_vars, mesh::P4estMesh{2, 3}, + equations::AbstractCovariantEquations{2, 3}, dg, + element) + (; derivative_matrix) = dg.basis + + for j in eachnode(dg), i in eachnode(dg) + + # Numerically differentiate covariant metric components with respect to ξ¹ + dG11dxi1 = dG12dxi1 = dG22dxi1 = zero(eltype(aux_node_vars)) + for ii in eachnode(dg) + aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, + element) + Gcov_ii = metric_covariant(aux_node_ii, equations) + dG11dxi1 = dG11dxi1 + derivative_matrix[i, ii] * Gcov_ii[1, 1] + dG12dxi1 = dG12dxi1 + derivative_matrix[i, ii] * Gcov_ii[1, 2] + dG22dxi1 = dG22dxi1 + derivative_matrix[i, ii] * Gcov_ii[2, 2] + end + + # Numerically differentiate covariant metric components with respect to ξ² + dG11dxi2 = dG12dxi2 = dG22dxi2 = zero(eltype(aux_node_vars)) + for jj in eachnode(dg) + aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, + element) + Gcov_jj = metric_covariant(aux_node_jj, equations) + dG11dxi2 = dG11dxi2 + derivative_matrix[j, jj] * Gcov_jj[1, 1] + dG12dxi2 = dG12dxi2 + derivative_matrix[j, jj] * Gcov_jj[1, 2] + dG22dxi2 = dG22dxi2 + derivative_matrix[j, jj] * Gcov_jj[2, 2] + end + + # Compute Christoffel symbols of the first kind + christoffel_firstkind_1 = SMatrix{2, 2}(0.5f0 * dG11dxi1, + 0.5f0 * dG11dxi2, + 0.5f0 * dG11dxi2, + dG12dxi2 - 0.5f0 * dG22dxi1) + christoffel_firstkind_2 = SMatrix{2, 2}(dG12dxi1 - 0.5f0 * dG11dxi2, + 0.5f0 * dG22dxi1, + 0.5f0 * dG22dxi1, + 0.5f0 * dG22dxi2) + + # Raise indices to get Christoffel symbols of the second kind + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + Gcon = metric_covariant(aux_node, equations) + aux_node_vars[20, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 1] + + Gcon[1, 2] * christoffel_firstkind_2[1, 1] + aux_node_vars[21, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 2] + + Gcon[1, 2] * christoffel_firstkind_2[1, 2] + aux_node_vars[22, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[2, 2] + + Gcon[1, 2] * christoffel_firstkind_2[2, 2] + + aux_node_vars[23, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[1, 1] + + Gcon[2, 2] * christoffel_firstkind_2[1, 1] + aux_node_vars[24, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[1, 2] + + Gcon[2, 2] * christoffel_firstkind_2[1, 2] + aux_node_vars[25, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[2, 2] + + Gcon[2, 2] * christoffel_firstkind_2[2, 2] + end +end end # muladd From 668d63d07f9b9b8c53d010295764871a663542fb Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Thu, 19 Dec 2024 01:07:10 +0100 Subject: [PATCH 03/21] fix typo in christoffel calculation --- .../containers_2d_manifold_in_3d_covariant.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl index 508e443..9580f1a 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl @@ -297,7 +297,9 @@ function calc_christoffel_symbols!(aux_node_vars, mesh::P4estMesh{2, 3}, for j in eachnode(dg), i in eachnode(dg) # Numerically differentiate covariant metric components with respect to ξ¹ - dG11dxi1 = dG12dxi1 = dG22dxi1 = zero(eltype(aux_node_vars)) + dG11dxi1 = zero(eltype(aux_node_vars)) + dG12dxi1 = zero(eltype(aux_node_vars)) + dG22dxi1 = zero(eltype(aux_node_vars)) for ii in eachnode(dg) aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, element) @@ -308,7 +310,9 @@ function calc_christoffel_symbols!(aux_node_vars, mesh::P4estMesh{2, 3}, end # Numerically differentiate covariant metric components with respect to ξ² - dG11dxi2 = dG12dxi2 = dG22dxi2 = zero(eltype(aux_node_vars)) + dG11dxi2 = zero(eltype(aux_node_vars)) + dG12dxi2 = zero(eltype(aux_node_vars)) + dG22dxi2 = zero(eltype(aux_node_vars)) for jj in eachnode(dg) aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, element) @@ -330,7 +334,7 @@ function calc_christoffel_symbols!(aux_node_vars, mesh::P4estMesh{2, 3}, # Raise indices to get Christoffel symbols of the second kind aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) - Gcon = metric_covariant(aux_node, equations) + Gcon = metric_contravariant(aux_node, equations) aux_node_vars[20, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 1] + Gcon[1, 2] * christoffel_firstkind_2[1, 1] aux_node_vars[21, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 2] + From 9508cc240a4e61acf31ba7f77a115967bcc72a48 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Mon, 23 Dec 2024 17:36:41 -0500 Subject: [PATCH 04/21] covariant SWE weak form runs now --- .vscode/settings.json | 5 + ...xir_shallowwater_covariant_cubed_sphere.jl | 73 ++++ src/TrixiAtmo.jl | 15 +- src/equations/covariant_advection.jl | 11 +- src/equations/covariant_shallow_water.jl | 359 ++++++++++++++++++ src/equations/equations.jl | 27 ++ .../dg_2d_manifold_in_3d_covariant.jl | 110 +++++- 7 files changed, 588 insertions(+), 12 deletions(-) create mode 100644 .vscode/settings.json create mode 100644 examples/elixir_shallowwater_covariant_cubed_sphere.jl create mode 100644 src/equations/covariant_shallow_water.jl diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..fcfefc2 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "cSpell.enabledFileTypes": { + "julia": false + } +} \ No newline at end of file diff --git a/examples/elixir_shallowwater_covariant_cubed_sphere.jl b/examples/elixir_shallowwater_covariant_cubed_sphere.jl new file mode 100644 index 0000000..825b70c --- /dev/null +++ b/examples/elixir_shallowwater_covariant_cubed_sphere.jl @@ -0,0 +1,73 @@ +############################################################################### +# DGSEM for the shallow water equations on the cubed sphere +############################################################################### + +using OrdinaryDiffEq, Trixi, TrixiAtmo + +############################################################################### +# Parameters + +initial_condition = initial_condition_geostrophic_balance +polydeg = 3 +cells_per_dimension = 5 +n_saves = 10 + +############################################################################### +# Spatial discretization + +tspan = (0.0, 1.0 * SECONDS_PER_DAY) + +mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg, + initial_refinement_level = 0, + element_local_mapping = true) + +equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, + EARTH_ROTATION_RATE, + global_coordinate_system = GlobalSphericalCoordinates()) + +# Create DG solver with polynomial degree = p +solver = DGSEM(polydeg = polydeg, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_weak_form), + volume_integral = VolumeIntegralWeakForm()) + +initial_condition_transformed = transform_initial_condition(initial_condition, equations) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver, + source_terms = source_terms_weak_form) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0 to T +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 50, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, + solution_variables = cons2cons) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.4) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 100.0, save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index f7815ba..4dac71b 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -31,22 +31,29 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") include("callbacks_step/callbacks_step.jl") export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, - CovariantLinearAdvectionEquation2D + CovariantLinearAdvectionEquation2D, CovariantShallowWaterEquations2D export GlobalCartesianCoordinates, GlobalSphericalCoordinates -export flux_chandrashekar, flux_LMARS +export flux_chandrashekar, flux_LMARS, flux_split_covariant, flux_nonconservative_weak_form, + flux_nonconservative_split_covariant export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal, - lake_at_rest_error, source_terms_lagrange_multiplier, + lake_at_rest_error, source_terms_lagrange_multiplier, source_terms_weak_form, clean_solution_lagrange_multiplier! + export cons2prim_and_vorticity + export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct, MetricTermsInvariantCurl + export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, SECONDS_PER_DAY + export global2contravariant, contravariant2global, spherical2cartesian, transform_initial_condition -export initial_condition_gaussian, initial_condition_gaussian_cartesian + +export initial_condition_gaussian, initial_condition_gaussian_cartesian, + initial_condition_geostrophic_balance export examples_dir end # module TrixiAtmo diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index c38db84..8181eff 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -60,15 +60,14 @@ function velocity(u, ::CovariantLinearAdvectionEquation2D) return SVector(u[2], u[3]) end -# Convert contravariant velocity/momentum components to the global coordinate system +# Convert contravariant velocity components to the global coordinate system @inline function contravariant2global(u, aux_vars, equations::CovariantLinearAdvectionEquation2D) vglo1, vglo2, vglo3 = basis_covariant(aux_vars, equations) * velocity(u, equations) return SVector(u[1], vglo1, vglo2, vglo3) end -# Convert velocity/momentum components in the global coordinate system to contravariant -# components +# Convert velocity components in the global coordinate system to contravariant components @inline function global2contravariant(u, aux_vars, equations::CovariantLinearAdvectionEquation2D) vcon1, vcon2 = basis_contravariant(aux_vars, equations) * SVector(u[2], u[3], u[4]) @@ -89,9 +88,9 @@ end return SVector(u[1], z, z) end -# Flux for abstract covariant equations as a function of the -# state vector u, as well as the auxiliary variables aux_vars, which contain the geometric -# information required for the covariant form +# Flux for abstract covariant equations as a function of the state vector u, as well as the +# auxiliary variables aux_vars, which contain the geometric information required for the +# covariant form @inline function Trixi.flux(u, aux_vars, orientation::Integer, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u)) diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl new file mode 100644 index 0000000..4b46d5a --- /dev/null +++ b/src/equations/covariant_shallow_water.jl @@ -0,0 +1,359 @@ +@muladd begin +#! format: noindent + +struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} <: + AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} + gravity::RealT + rotation_rate::RealT + global_coordinate_system::GlobalCoordinateSystem + function CovariantShallowWaterEquations2D(gravity::RealT, + rotation_rate::RealT; + global_coordinate_system = GlobalCartesianCoordinates()) where {RealT <: + Real} + return new{typeof(global_coordinate_system), RealT}(gravity, rotation_rate) + end +end + +# The conservative variables are the height and contravariant momentum components +function Trixi.varnames(::typeof(cons2cons), ::CovariantShallowWaterEquations2D) + return ("h", "h_vcon1", "h_vcon2") +end + +# Convenience function to extract the velocity +function velocity(u, ::CovariantShallowWaterEquations2D) + return SVector(u[2] / u[1], u[3] / u[1]) +end + +# Convenience function to extract the momentum +function momentum(u, ::CovariantShallowWaterEquations2D) + return SVector(u[2], u[3]) +end + +# Our implementation of flux-differencing formulation uses nonconservative terms, but the +# standard weak form does not. To handle both options, we have defined a dummy kernel for +# the nonconservative terms that does nothing when VolumeIntegralWeakForm is used with a +# nonconservative system. +Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True() + +# Entropy function (total energy per unit volume) +@inline function Trixi.entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + Gcov = metric_covariant(aux_vars, equations) + vcon = SVector(h_vcon1 / h, h_vcon2 / h) + vcov = Gcov * vcon + return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2) +end + +# Entropy variables (partial derivatives of entropy with respect to conservative variables) +@inline function Trixi.cons2entropy(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + + Gcov = metric_covariant(aux_vars, equations) + vcon = SVector(h_vcon1 / h, h_vcon2 / h) + vcov = Gcov * vcon + w1 = equations.gravity * h - 0.5f0 * dot(vcov, vcon) + return SVector{3}(w1, vcov[1], vcov[2]) +end + +# Height and three global momentum components +function Trixi.varnames(::typeof(contravariant2global), + ::CovariantShallowWaterEquations2D) + return ("h", "h_vglo1", "h_vglo2", "h_vglo3") +end + +# Convert contravariant momentum components to the global coordinate system +@inline function contravariant2global(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + vglo1, vglo2, vglo3 = basis_covariant(aux_vars, equations) * momentum(u, equations) + return SVector(u[1], vglo1, vglo2, vglo3) +end + +# Convert momentum components in the global coordinate system to contravariant components +@inline function global2contravariant(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + vcon1, vcon2 = basis_contravariant(aux_vars, equations) * SVector(u[2], u[3], u[4]) + return SVector(u[1], vcon1, vcon2) +end + +# The flux for the covariant form takes in the element container and node/element indices +# in order to give the flux access to the geometric information +@inline function Trixi.flux(u, aux_vars, orientation::Integer, + equations::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + + J = area_element(aux_vars, equations) + Gcon = metric_contravariant(aux_vars, equations) + gravitational_term = 0.5f0 * equations.gravity * h^2 + h_vcon = SVector(h_vcon1, h_vcon2) + vcon_orientation = h_vcon[orientation] / h + + momentum_flux_1 = h_vcon1 * vcon_orientation + + Gcon[1, orientation] * gravitational_term + momentum_flux_2 = h_vcon2 * vcon_orientation + + Gcon[2, orientation] * gravitational_term + + return SVector(J * h_vcon[orientation], J * momentum_flux_1, J * momentum_flux_2) +end + +# Symmetric part of entropy-conservative flux +@inline function flux_split_covariant(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation::Integer, + equations::CovariantShallowWaterEquations2D) + h_ll, h_vcon1_ll, h_vcon2_ll = u_ll + h_rr, h_vcon1_rr, h_vcon2_rr = u_rr + + J_ll = area_element(aux_vars_ll, equations) + J_rr = area_element(aux_vars_rr, equations) + + h_vcon_ll = SVector(h_vcon1_ll, h_vcon2_ll) + h_vcon_rr = SVector(h_vcon1_rr, h_vcon2_rr) + + # Scaled mass flux in conservative form + mass_flux = 0.5f0 * (J_ll * h_vcon_ll[orientation] + J_rr * h_vcon_rr[orientation]) + + # Half of scaled inertial flux in conservative form + momentum_flux = 0.25f0 * (J_ll * h_vcon_ll * h_vcon_ll[orientation] / h_ll + + J_rr * h_vcon_rr * h_vcon_rr[orientation] / h_rr) + + return SVector(mass_flux, momentum_flux[1], momentum_flux[2]) +end + +# Split-covariant flux with local Lax-Friedrichs dissipation +const flux_split_covariant_lax_friedrichs = FluxPlusDissipation(flux_split_covariant, + DissipationLocalLaxFriedrichs(max_abs_speed_naive)) + +# Non-symmetric part of entropy-conservative flux +@inline function flux_nonconservative_split_covariant(u_ll, u_rr, aux_vars_ll, + aux_vars_rr, + orientation::Integer, + equations::CovariantShallowWaterEquations2D) + h_ll, h_vcon1_ll, h_vcon2_ll = u_ll + h_rr, h_vcon1_rr, h_vcon2_rr = u_rr + + # Geometric variables + Gcov_ll = metric_covariant(aux_vars_ll, equations) + Gcov_rr = metric_covariant(aux_vars_rr, equations) + Gcon_ll = metric_contravariant(aux_vars_ll, equations) + J_ll = area_element(aux_vars_ll, equations) + J_rr = area_element(aux_vars_rr, equations) + + # Contravariant and covariant momentum and velocity components + h_vcon_ll = SVector(h_vcon1_ll, h_vcon2_ll) + h_vcon_rr = SVector(h_vcon1_rr, h_vcon2_rr) + vcov_ll = Gcov_ll * h_vcon_ll / h_ll + vcov_rr = Gcov_rr * h_vcon_rr / h_rr + + # Half of inertial term in non-conservative form + nonconservative_term_inertial = 0.5f0 * Gcon_ll * + (J_ll * h_vcon_ll[orientation] * vcov_rr + + J_rr * h_vcon_rr[orientation] * vcov_ll) + + # Gravity term in non-conservative form + nonconservative_term_gravitational = equations.gravity * + J_ll * Gcon_ll[:, orientation] * + h_ll * h_rr # the same as h_ll * (h_ll + h_rr) + + nonconservative_term = nonconservative_term_inertial + + nonconservative_term_gravitational + + return SVector(zero(eltype(u_ll)), nonconservative_term[1], nonconservative_term[2]) +end + +@inline function source_terms_weak_form(u, x, t, aux_vars, + equations::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + + # Geometric variables + Gcon = metric_contravariant(aux_vars, equations) + (Gamma1, Gamma2) = christoffel_symbols(aux_vars, equations) + J = area_element(aux_vars, equations) + + # Physical variables + h_vcon = SVector{2}(h_vcon1, h_vcon2) + v_con = SVector{2}(h_vcon1 / h, h_vcon2 / h) + + # Doubly-contravariant and mixed inertial flux tensors + T = h_vcon * v_con' + 0.5f0 * equations.gravity * h^2 * Gcon + + # Coriolis parameter + f = 2 * equations.rotation_rate * x[3] / norm(x) # 2Ωsinθ + + # Geometric source term + s_geo = SVector(sum(Gamma1 .* T), sum(Gamma2 .* T)) + + # Combined source terms + source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon1 - Gcon[1, 1] * h_vcon2) + source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon1 - Gcon[2, 1] * h_vcon2) + + # Do not scale by Jacobian since apply_jacobian! is called before this + return SVector(zero(eltype(u)), -source_1, -source_2) +end + +# Geometric and Coriolis source terms for a rotating sphere +@inline function source_terms_split_covariant(u, x, t, aux_vars, + equations::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + + # Geometric variables + Gcov = metric_covariant(aux_vars, equations) + Gcon = metric_contravariant(aux_vars, equations) + (Gamma1, Gamma2) = christoffel_symbols(aux_vars, equations) + J = area_element(aux_vars, equations) + + # Physical variables + h_vcon = SVector{2}(h_vcon1, h_vcon2) + v_con = SVector{2}(h_vcon1 / h, h_vcon2 / h) + + # Doubly-contravariant and mixed inertial flux tensors + h_vcon_vcon = h_vcon * v_con' + h_vcov_vcon = Gcov * h_vcon_vcon + + # Coriolis parameter + f = 2 * equations.rotation_rate * x[3] / norm(x) # 2Ωsinθ + + # Geometric source term + s_geo = 0.5f0 * (SVector(sum(Gamma1 .* h_vcon_vcon), sum(Gamma2 .* h_vcon_vcon)) - + Gcon * (Gamma1 * h_vcov_vcon[1, :] + Gamma2 * h_vcov_vcon[2, :])) + + # Combined source terms + source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon1 - Gcon[1, 1] * h_vcon2) + source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon1 - Gcon[2, 1] * h_vcon2) + + # Do not scale by Jacobian since apply_jacobian! is called before this + return SVector(zero(eltype(u)), -source_1, -source_2) +end + +# Maximum wave speed along the normal direction in reference space +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation, + equations::CovariantShallowWaterEquations2D) + h_ll, h_vcon1_ll, h_vcon2_ll = u_ll + h_rr, h_vcon1_rr, h_vcon2_rr = u_rr + + h_vcon_ll = SVector(h_vcon1_ll, h_vcon2_ll) + h_vcon_rr = SVector(h_vcon1_rr, h_vcon2_rr) + + Gcon_ll = metric_contravariant(aux_vars_ll, equations) + Gcon_rr = metric_contravariant(aux_vars_rr, equations) + + phi_ll = max(h_ll * equations.gravity, 0) + phi_rr = max(h_rr * equations.gravity, 0) + + return max(abs(h_vcon_ll[orientation] / h_ll) + + sqrt(Gcon_ll[orientation, orientation] * phi_ll), + abs(h_vcon_rr[orientation] / h_rr) + + sqrt(Gcon_rr[orientation, orientation] * phi_rr)) +end + +# Maximum wave speeds with respect to the covariant basis +@inline function Trixi.max_abs_speeds(u, aux_vars, + equations::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + Gcon = metric_contravariant(aux_vars, equations) + phi = max(h * equations.gravity, 0) + return abs(h_vcon1 / h) + sqrt(Gcon[1, 1] * phi), + abs(h_vcon2 / h) + sqrt(Gcon[2, 2] * phi) +end + +# Steady geostrophically balanzed zonal flow +function initial_condition_geostrophic_balance(x, t, + equations::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) + (; gravity, rotation_rate) = equations + + radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) + lat = asin(x[3] / radius) + + # compute zonal and meridional components of the velocity + V = convert(eltype(x), 2π) * radius / (12 * SECONDS_PER_DAY) + vlon, vlat = V * cos(lat), zero(eltype(x)) + + # compute geopotential height + h = 1 / gravity * + (2.94f4 - (radius * rotation_rate * V + 0.5f0 * V^2) * (sin(lat))^2) + + # convert to conservative variables + return SVector(h, h * vlon, h * vlat, zero(eltype(x))) +end + +# Rossby-Haurwitz wave +function initial_condition_rossby_haurwitz(x, t, + equations::CovariantShallowWaterEquations2D{<:GlobalSphericalCoordinates}) + (; gravity, rotation_rate) = equations + + radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) + lon = atan(x[2], x[1]) + lat = asin(x[3] / radius) + + h_0 = 8.0f3 + K = 7.848f-6 + R = 4.0f0 + + A = 0.5f0 * K * (2 * rotation_rate + K) * (cos(lat))^2 + + 0.25f0 * K^2 * (cos(lat))^(2 * R) * + ((R + 1) * (cos(lat))^2 + + (2 * R^2 - R - 2) - 2 * R^2 / ((cos(lat))^2)) + B = 2 * (rotation_rate + K) * K / ((R + 1) * (R + 2)) * (cos(lat))^R * + ((R^2 + 2R + 2) - (R + 1)^2 * (cos(lat))^2) + C = 0.25f0 * K^2 * (cos(lat))^(2 * R) * ((R + 1) * (cos(lat))^2 - (R + 2)) + + h = h_0 + + (1 / gravity) * + (radius^2 * A + radius^2 * B * cos(R * lon) + radius^2 * C * cos(2 * R * lon)) + + vlon = radius * K * cos(lat) + + radius * K * (cos(lat))^(R - 1) * (R * (sin(lat))^2 - (cos(lat))^2) * + cos(R * lon) + vlat = -radius * K * R * (cos(lat))^(R - 1) * sin(lat) * sin(R * lon) + + # convert to conservative variables + return SVector(h, h * vlon, h * vlat) +end + +@inline function galewsky_velocity(θ, u_0, θ_0, θ_1) + if (θ_0 < θ) && (θ < θ_1) + u = u_0 / exp(-4 / (θ_1 - θ_0)^2) * exp(1 / (θ - θ_0) * 1 / (θ - θ_1)) + else + u = zero(θ) + end + return u +end + +@inline function galewsky_integrand(θ, u_0, θ_0, θ_1, a, + equations::CovariantShallowWaterEquations2D) + (; rotation_rate) = equations + u = galewsky_velocity(θ, u_0, θ_0, θ_1) + return u * (2 * rotation_rate * sin(θ) + u * tan(θ) / a) +end + +function initial_condition_barotropic_instability(x, t, + equations::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) + (; gravity) = equations + realT = eltype(x) + radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) + lon = atan(x[2], x[1]) + lat = asin(x[3] / radius) + + # compute zonal and meridional velocity components + u_0 = 80.0f0 + lat_0 = convert(realT, π / 7) + lat_1 = convert(realT, π / 2) - lat_0 + vlon = galewsky_velocity(lat, u_0, lat_0, lat_1) + vlat = zero(eltype(x)) + + # numerically integrate (here we use the QuadGK package) to get height + galewsky_integral, _ = quadgk(latp -> galewsky_integrand(latp, u_0, lat_0, lat_1, + radius, + equations), -π / 2, lat) + h = 10158.0f0 - radius / gravity * galewsky_integral + + # add perturbation to initiate instability + α, β = convert(realT, 1 / 3), convert(realT, 1 / 15) + lat_2 = convert(realT, π / 4) + if (-π < lon) && (lon < π) + h = h + 120.0f0 * cos(lat) * exp(-((lon / α)^2)) * exp(-((lat_2 - lat) / β)^2) + end + # convert to conservative variables + return SVector(h, h * vlon, h * vlat) +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 7c9c635..41b4112 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -260,11 +260,38 @@ end return 0.5f0 * (flux_ll + flux_rr) end +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, aux_vars_ll, + aux_vars_rr, + orientation_or_normal_direction, + equations::AbstractCovariantEquations) + sqrtG = area_element(aux_vars_ll, equations) + λ = dissipation.max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + orientation_or_normal_direction, equations) + return -0.5f0 * sqrtG * λ * (u_rr - u_ll) +end + +# Dummy flux for weak form +@inline function flux_nonconservative_weak_form(u_ll::SVector{NVARS, RealT}, + u_rr::SVector{NVARS, RealT}, + aux_vars_ll, aux_vars_rr, + orientation_or_normal_direction, + equations::AbstractCovariantEquations{2, + NDIMS_AMBIENT, + GlobalCoordinateSystem, + NVARS}) where { + NDIMS_AMBIENT, + GlobalCoordinateSystem, + NVARS, + RealT + } + return zeros(SVector{NVARS, RealT}) +end abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("reference_data.jl") include("covariant_advection.jl") +include("covariant_shallow_water.jl") include("compressible_moist_euler_2d_lucas.jl") include("shallow_water_3d.jl") end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl index ec18e04..ccc1f77 100644 --- a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl @@ -113,8 +113,17 @@ end return nothing end -# Pointwise interface flux, transforming the contravariant prognostic variables into the -# local coordinate system +# Since the standard weak form for CovariantShallowWaterEquations2D has no nonconservative +# terms, but the equation type has have_nonconservative_terms = True(), we must define +# an appropriate kernel that just calls the standard weak form +@inline function Trixi.weak_form_kernel!(du, u, element, mesh::P4estMesh{2}, + nonconservative_terms::True, + equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, alpha = true) + Trixi.weak_form_kernel!(du, u, element, mesh, False(), equations, dg, cache, alpha) +end + +# Pointwise interface flux in local coordinates for problems without nonconservative terms @inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, nonconservative_terms::False, equations::AbstractCovariantEquations{2}, @@ -175,6 +184,103 @@ end end end +# Pointwise interface flux in local coordinates for problems with nonconservative terms +@inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, + nonconservative_terms::True, + equations::AbstractCovariantEquations{2}, + surface_integral, dg::DG, cache, + interface_index, normal_direction, + primary_node_index, + primary_direction_index, + primary_element_index, + secondary_node_index, + secondary_direction_index, + secondary_element_index) + (; u) = cache.interfaces + (; aux_surface_node_vars) = cache.auxiliary_variables + surface_flux, nonconservative_flux = surface_integral.surface_flux + + # Get surface values for solution and auxiliary variables + u_ll, u_rr = Trixi.get_surface_node_vars(u, equations, dg, primary_node_index, + interface_index) + aux_vars_ll, aux_vars_rr = get_surface_node_aux_vars(aux_surface_node_vars, + equations, + dg, primary_node_index, + interface_index) + + # Compute flux in the primary element's coordinate system + u_rr_spherical = contravariant2global(u_rr, aux_vars_rr, equations) + u_rr_transformed_to_ll = global2contravariant(u_rr_spherical, aux_vars_ll, + equations) + primary_orientation = (primary_direction_index + 1) >> 1 + if isodd(primary_direction_index) + flux_primary = -(surface_flux(u_rr_transformed_to_ll, u_ll, aux_vars_ll, + aux_vars_ll, primary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_rr_transformed_to_ll, u_ll, + aux_vars_ll, aux_vars_ll, + primary_orientation, equations)) + else + flux_primary = surface_flux(u_ll, u_rr_transformed_to_ll, aux_vars_ll, + aux_vars_ll, primary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_ll, u_rr_transformed_to_ll, + aux_vars_ll, aux_vars_ll, + primary_orientation, equations) + end + + # Compute flux in the secondary element's coordinate system + u_ll_spherical = contravariant2global(u_ll, aux_vars_ll, equations) + u_ll_transformed_to_rr = global2contravariant(u_ll_spherical, aux_vars_rr, + equations) + secondary_orientation = (secondary_direction_index + 1) >> 1 + if isodd(secondary_direction_index) + flux_secondary = -(surface_flux(u_ll_transformed_to_rr, u_rr, aux_vars_rr, + aux_vars_rr, secondary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_ll_transformed_to_rr, u_rr, + aux_vars_rr, aux_vars_rr, + secondary_orientation, equations)) + else + flux_secondary = surface_flux(u_rr, u_ll_transformed_to_rr, aux_vars_rr, + aux_vars_rr, secondary_orientation, equations) + + 0.5f0 * nonconservative_flux(u_rr, u_ll_transformed_to_rr, + aux_vars_rr, aux_vars_rr, + secondary_orientation, equations) + end + + # Update the surface flux values on both sides of the interface + for v in eachvariable(equations) + surface_flux_values[v, primary_node_index, primary_direction_index, + primary_element_index] = flux_primary[v] + surface_flux_values[v, secondary_node_index, secondary_direction_index, + secondary_element_index] = flux_secondary[v] + end +end + +# This function passes the auxiliary variables into the source term +function Trixi.calc_sources!(du, u, t, source_terms, + equations::AbstractCovariantEquations{2}, dg::DG, cache) + (; aux_node_vars) = cache.auxiliary_variables + + Trixi.@threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_local = Trixi.get_node_vars(u, equations, dg, i, j, element) + x_local = Trixi.get_node_coords(cache.elements.node_coordinates, equations, + dg, + i, j, element) + aux_local = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + du_local = source_terms(u_local, x_local, t, aux_local, equations) + Trixi.add_to_node_vars!(du, du_local, equations, dg, i, j, element) + end + end + + return nothing +end + +# Version of calc_sources! specialized for covariant formulation with no source term +function Trixi.calc_sources!(du, u, t, source_terms::Nothing, + equations::AbstractCovariantEquations{2}, dg::DG, cache) + return nothing +end + # Apply the exact Jacobian stored in auxiliary variables function Trixi.apply_jacobian!(du, mesh::P4estMesh{2}, equations::AbstractCovariantEquations{2}, From 1b5e975c238d43a4fbce319be2276d0995d91038 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Mon, 23 Dec 2024 17:36:59 -0500 Subject: [PATCH 05/21] covariant SWE weak form runs now --- .vscode/settings.json | 5 ----- 1 file changed, 5 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index fcfefc2..0000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,5 +0,0 @@ -{ - "cSpell.enabledFileTypes": { - "julia": false - } -} \ No newline at end of file From b620511a3e206450a8b9045e379734f3d1d8090a Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Mon, 23 Dec 2024 18:17:15 -0500 Subject: [PATCH 06/21] no longer pass surface_integral into prolong2mortars --- src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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 7c1a81e..1b494ef 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 @@ -46,8 +46,7 @@ function Trixi.rhs!(du, u, t, # Prolong solution to mortars Trixi.@trixi_timeit Trixi.timer() "prolong2mortars" begin - Trixi.prolong2mortars!(cache, u, mesh, equations, - dg.mortar, dg.surface_integral, dg) + Trixi.prolong2mortars!(cache, u, mesh, equations, dg.mortar, dg) end # Calculate mortar fluxes From b8410a0237df72a11112bf13cccb0b303cf0d607 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Mon, 23 Dec 2024 18:57:41 -0500 Subject: [PATCH 07/21] bump compat for Trixi from 0.9 to 0.9.9 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7acc0d0..c99dfe2 100644 --- a/Project.toml +++ b/Project.toml @@ -29,5 +29,5 @@ Static = "0.8, 1" StaticArrayInterface = "1.5.1" StaticArrays = "1" StrideArrays = "0.1.28" -Trixi = "0.9" +Trixi = "0.9.9" julia = "1.9" From 1e725bc9d69d12bace876002a4d906f2ac5f2532 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Wed, 25 Dec 2024 18:19:31 -0500 Subject: [PATCH 08/21] add test for spherical SWE --- test/runtests.jl | 4 ++++ test/test_2d_shallow_water_covariant.jl | 32 +++++++++++++++++++++++++ 2 files changed, 36 insertions(+) create mode 100644 test/test_2d_shallow_water_covariant.jl diff --git a/test/runtests.jl b/test/runtests.jl index 38d3e81..8e65de2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,4 +24,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "shallow_water_3d" include("test_3d_shallow_water.jl") end + + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "shallow_water_2d_covariant" + include("test_2d_shallow_water_covariant.jl") + end end diff --git a/test/test_2d_shallow_water_covariant.jl b/test/test_2d_shallow_water_covariant.jl new file mode 100644 index 0000000..d5e54d6 --- /dev/null +++ b/test/test_2d_shallow_water_covariant.jl @@ -0,0 +1,32 @@ +module TestSphericalAdvection + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") + +@trixiatmo_testset "spherical shallow water, covariant weak form" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_covariant_cubed_sphere.jl"), + l2=[ + 0.3065314463985936, + 0.00019984467582352902, + 8.767819502807507e-5 + ], + linf=[ + 1.4786544678290738, + 0.0013754600033514114, + 0.0007564014737144256 + ]) + # 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 From f741b925bdee5ca37b23bb2d6c60c339bc610ad5 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Mon, 30 Dec 2024 00:06:13 -0500 Subject: [PATCH 09/21] approach to automatically transform to desired coordinate system --- src/equations/covariant_advection.jl | 23 ++++++++++ src/equations/covariant_shallow_water.jl | 58 ++++++++++++++---------- src/equations/equations.jl | 47 ++++++++++++++++++- src/equations/reference_data.jl | 51 ++++++--------------- src/equations/shallow_water_3d.jl | 9 ++++ 5 files changed, 127 insertions(+), 61 deletions(-) diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index 8181eff..5e08fea 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -126,4 +126,27 @@ end equations::CovariantLinearAdvectionEquation2D) return abs.(velocity(u, equations)) end + +# If the initial velocity field is defined in Cartesian coordinates and the chosen global +# coordinate system is spherical, perform the appropriate conversion +@inline function cartesian2global(u, x, + ::CovariantLinearAdvectionEquation2D{GlobalSphericalCoordinates}) + vlon, vlat, vrad = cartesian2spherical(u[2], u[3], u[4], x) + return SVector(u[1], vlon, vlat, vrad) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is Cartesian, perform the appropriate conversion +@inline function spherical2global(u, x, + ::CovariantLinearAdvectionEquation2D{GlobalCartesianCoordinates}) + vx, vy, vz = spherical2cartesian(u[2], u[3], u[4], x) + return SVector(u[1], vx, vy, vz) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is spherical, do not convert +@inline function spherical2global(u, x, + ::CovariantLinearAdvectionEquation2D{GlobalSphericalCoordinates}) + return u +end end # @muladd diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index 4b46d5a..6e636b4 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -44,11 +44,20 @@ Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True() return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2) end +@inline function Trixi.cons2prim(u, ::CovariantShallowWaterEquations2D) + h, h_vcon1, h_vcon2 = u + return SVector(h, h_vcon1 / h, h_vcon2 / h) +end + +@inline function Trixi.prim2cons(u, ::CovariantShallowWaterEquations2D) + h, vcon1, vcon2 = u + return SVector(h, h * vcon1, h * vcon2) +end + # Entropy variables (partial derivatives of entropy with respect to conservative variables) @inline function Trixi.cons2entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D) h, h_vcon1, h_vcon2 = u - Gcov = metric_covariant(aux_vars, equations) vcon = SVector(h_vcon1 / h, h_vcon2 / h) vcov = Gcov * vcon @@ -56,6 +65,7 @@ end return SVector{3}(w1, vcov[1], vcov[2]) end + # Height and three global momentum components function Trixi.varnames(::typeof(contravariant2global), ::CovariantShallowWaterEquations2D) @@ -75,7 +85,6 @@ end vcon1, vcon2 = basis_contravariant(aux_vars, equations) * SVector(u[2], u[3], u[4]) return SVector(u[1], vcon1, vcon2) end - # The flux for the covariant form takes in the element container and node/element indices # in order to give the flux access to the geometric information @inline function Trixi.flux(u, aux_vars, orientation::Integer, @@ -173,7 +182,7 @@ end h_vcon = SVector{2}(h_vcon1, h_vcon2) v_con = SVector{2}(h_vcon1 / h, h_vcon2 / h) - # Doubly-contravariant and mixed inertial flux tensors + # Doubly-contravariant flux tensor T = h_vcon * v_con' + 0.5f0 * equations.gravity * h^2 * Gcon # Coriolis parameter @@ -256,26 +265,6 @@ end abs(h_vcon2 / h) + sqrt(Gcon[2, 2] * phi) end -# Steady geostrophically balanzed zonal flow -function initial_condition_geostrophic_balance(x, t, - equations::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) - (; gravity, rotation_rate) = equations - - radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - lat = asin(x[3] / radius) - - # compute zonal and meridional components of the velocity - V = convert(eltype(x), 2π) * radius / (12 * SECONDS_PER_DAY) - vlon, vlat = V * cos(lat), zero(eltype(x)) - - # compute geopotential height - h = 1 / gravity * - (2.94f4 - (radius * rotation_rate * V + 0.5f0 * V^2) * (sin(lat))^2) - - # convert to conservative variables - return SVector(h, h * vlon, h * vlat, zero(eltype(x))) -end - # Rossby-Haurwitz wave function initial_condition_rossby_haurwitz(x, t, equations::CovariantShallowWaterEquations2D{<:GlobalSphericalCoordinates}) @@ -356,4 +345,27 @@ function initial_condition_barotropic_instability(x, t, # convert to conservative variables return SVector(h, h * vlon, h * vlat) end + +# If the initial velocity field is defined in Cartesian coordinates and the chosen global +# coordinate system is spherical, perform the appropriate conversion +@inline function cartesian2global(u, x, + ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) + h_vlon, h_vlat, h_vrad = cartesian2spherical(u[2], u[3], u[4], x) +return SVector(u[1], h_vlon, h_vlat, h_vrad) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is Cartesian, perform the appropriate conversion +@inline function spherical2global(u, x, + ::CovariantShallowWaterEquations2D{GlobalCartesianCoordinates}) +h_vx, h_vy, h_vz = spherical2cartesian(u[2], u[3], u[4], x) +return SVector(u[1], h_vx, h_vy, h_vz) +end + +# If the initial velocity field is defined in spherical coordinates and the chosen global +# coordinate system is spherical, do not convert +@inline function spherical2global(u, x, + ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) +return u +end end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 41b4112..8130451 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -132,6 +132,12 @@ performed by [`contravariant2global`](@ref). """ function global2contravariant end +# By default, the equations are assumed to be formulated in Cartesian coordinates. This +# function is specialized where needed. +function cartesian2global(u, x, equations::AbstractEquations) + return u +end + # cons2cons method which takes in unused aux_vars variable @inline Trixi.cons2cons(u, aux_vars, ::AbstractEquations) = u @inline Trixi.prim2cons(u, aux_vars, ::AbstractEquations) = u @@ -270,7 +276,7 @@ end return -0.5f0 * sqrtG * λ * (u_rr - u_ll) end -# Dummy flux for weak form +# Dummy two-point nonconservative flux for weak form @inline function flux_nonconservative_weak_form(u_ll::SVector{NVARS, RealT}, u_rr::SVector{NVARS, RealT}, aux_vars_ll, aux_vars_rr, @@ -286,12 +292,49 @@ end } return zeros(SVector{NVARS, RealT}) end + +# Convert a vector from a global spherical to Cartesian basis representation, where we note +# that the radial component is not necessarily zero +@inline function spherical2cartesian(vlon, vlat, vrad, x) + # compute longitude and latitude + lon, lat = atan(x[2], x[1]), asin(x[3] / norm(x)) + + # compute trigonometric functions + sinlon, coslon = sincos(lon) + sinlat, coslat = sincos(lat) + + # return Cartesian components + vx = -sinlon * vlon - sinlat * coslon * vlat + coslat * coslon * vrad + vy = coslon * vlon - sinlat * sinlon * vlat + coslat * sinlon * vrad + vz = coslat * vlat + sinlat * vrad + + return vx, vy, vz +end + +# Convert a vector from a global Cartesian to spherical basis representation, where we note +# that the radial component is not necessarily zero +@inline function cartesian2spherical(vx, vy, vz, x) + # compute longitude and latitude + lon, lat = atan(x[2], x[1]), asin(x[3] / norm(x)) + + # compute trigonometric functions + sinlon, coslon = sincos(lon) + sinlat, coslat = sincos(lat) + + # return spherical components + vlon = -sinlon * vx + coslon * vy + vlat = -sinlat * coslon * vx - sinlat * sinlon * vy + coslat * vz + vrad = coslat * coslon * vx + coslat * sinlon * vy + sinlat * vz # zero for any tangent vector + + return vlon, vlat, vrad +end + abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end -include("reference_data.jl") include("covariant_advection.jl") include("covariant_shallow_water.jl") include("compressible_moist_euler_2d_lucas.jl") include("shallow_water_3d.jl") +include("reference_data.jl") end # @muladd diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 6865223..4c531f6 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -33,7 +33,7 @@ This problem is adapted from Case 1 of the test suite described in the following spherical geometry. Journal of Computational Physics, 102(1):211-224. [DOI: 10.1016/S0021-9991(05)80016-6](https://doi.org/10.1016/S0021-9991(05)80016-6) """ -@inline function initial_condition_gaussian(x, t, ::AbstractEquations) +@inline function initial_condition_gaussian(x, t, equations::AbstractEquations) RealT = eltype(x) a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) # radius of the sphere @@ -65,46 +65,25 @@ This problem is adapted from Case 1 of the test suite described in the following # Prescribe the rotated bell shape and Cartesian velocity components. # The last variable is the bottom topography, which we set to zero - return SVector(h, vx, vy, vz, 0.0f0) + return cartesian2global(SVector(h, vx, vy, vz, 0.0f0), x, equations) end -# Version for spherical coordinates (note: the velocity is not well defined at the poles) -@inline function initial_condition_gaussian(x, t, - ::AbstractCovariantEquations{2, 3, - GlobalSphericalCoordinates}) - RealT = eltype(x) - - a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) # radius of the sphere - omega = convert(RealT, 2π) / (12 * SECONDS_PER_DAY) # angular velocity - alpha = convert(RealT, π / 4) # angle of rotation - h_0 = 1000.0f0 # bump height in metres - b_0 = 5.0f0 # bump width parameter - lon_0, lat_0 = convert(RealT, 3π / 2), 0.0f0 # initial bump location - - # axis of rotation - axis = SVector(-cos(alpha), 0.0f0, sin(alpha)) - - # convert initial position to Cartesian coordinates - x_0 = SVector(a * cos(lat_0) * cos(lon_0), - a * cos(lat_0) * sin(lon_0), - a * sin(lat_0)) +function initial_condition_geostrophic_balance(x, t, + equations::CovariantShallowWaterEquations2D) + (; gravity, rotation_rate) = equations - # apply rotation using Rodrigues' formula - axis_cross_x_0 = Trixi.cross(axis, x_0) - x_0 = x_0 + sin(omega * t) * axis_cross_x_0 + - (1 - cos(omega * t)) * Trixi.cross(axis, axis_cross_x_0) + radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) + lat = asin(x[3] / radius) - # compute Gaussian bump profile - h = h_0 * - exp(-b_0 * ((x[1] - x_0[1])^2 + (x[2] - x_0[2])^2 + (x[3] - x_0[3])^2) / (a^2)) + # compute zonal and meridional components of the velocity + V = convert(eltype(x), 2π) * radius / (12 * SECONDS_PER_DAY) + vlon, vlat = V * cos(lat), zero(eltype(x)) - # get zonal and meridional components of the velocity - lon, lat = atan(x[2], x[1]), asin(x[3] / a) - vlon = omega * a * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha)) - vlat = -omega * a * sin(lon) * sin(alpha) + # compute geopotential height + h = 1 / gravity * + (2.94f4 - (radius * rotation_rate * V + 0.5f0 * V^2) * (sin(lat))^2) - # Prescribe the rotated bell shape and spherical velocity components. - # The last variable is the bottom topography, which we set to zero - return SVector(h, vlon, vlat, 0.0f0, 0.0f0) + # convert primitive variables from spherical coordinates to the chosen coordinate system + return spherical2global(SVector(h, vlon, vlat, zero(eltype(x))), x, equations) end end # muladd diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl index 930af2a..c008b22 100644 --- a/src/equations/shallow_water_3d.jl +++ b/src/equations/shallow_water_3d.jl @@ -380,4 +380,13 @@ end @inline function energy_internal(cons, equations::ShallowWaterEquations3D) return energy_total(cons, equations) - energy_kinetic(cons, equations) end + +# Convert spherical velocity components to global Cartesian components +@inline function spherical2global(prim, x, ::ShallowWaterEquations3D) + h, vlat, vlon, vrad, b = prim # primitive variables + + v1, v2, v3 = spherical2cartesian(vlat, vlon, vrad, x) + + return SVector(h, v1, v2, v3, b) +end end # @muladd From a90bfa09608701c6946eb03259d3df076706fe57 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Tue, 31 Dec 2024 00:30:31 -0500 Subject: [PATCH 10/21] added missing aux_vars to prim2cons --- Project.toml | 1 + src/TrixiAtmo.jl | 2 ++ src/equations/covariant_shallow_water.jl | 8 ++++---- src/equations/equations.jl | 23 +++++++++++++---------- 4 files changed, 20 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index c99dfe2..fb2eaf3 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.1.0-DEV" [deps] HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 4dac71b..7f3b433 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -20,6 +20,8 @@ using Reexport: @reexport using LoopVectorization: @turbo using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace +using Infiltrator + @reexport using StaticArrays: SVector, SMatrix include("auxiliary/auxiliary.jl") diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index 6e636b4..07dd9c2 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -44,12 +44,12 @@ Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True() return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2) end -@inline function Trixi.cons2prim(u, ::CovariantShallowWaterEquations2D) +@inline function Trixi.cons2prim(u, aux_vars, ::CovariantShallowWaterEquations2D) h, h_vcon1, h_vcon2 = u return SVector(h, h_vcon1 / h, h_vcon2 / h) end -@inline function Trixi.prim2cons(u, ::CovariantShallowWaterEquations2D) +@inline function Trixi.prim2cons(u, aux_vars, ::CovariantShallowWaterEquations2D) h, vcon1, vcon2 = u return SVector(h, h * vcon1, h * vcon2) end @@ -65,7 +65,6 @@ end return SVector{3}(w1, vcov[1], vcov[2]) end - # Height and three global momentum components function Trixi.varnames(::typeof(contravariant2global), ::CovariantShallowWaterEquations2D) @@ -169,6 +168,7 @@ const flux_split_covariant_lax_friedrichs = FluxPlusDissipation(flux_split_covar return SVector(zero(eltype(u_ll)), nonconservative_term[1], nonconservative_term[2]) end +# Geometric and Coriolis sources for a rotating sphere with VolumeIntegralWeakForm @inline function source_terms_weak_form(u, x, t, aux_vars, equations::CovariantShallowWaterEquations2D) h, h_vcon1, h_vcon2 = u @@ -199,7 +199,7 @@ end return SVector(zero(eltype(u)), -source_1, -source_2) end -# Geometric and Coriolis source terms for a rotating sphere +# Geometric and Coriolis sources for a rotating sphere with VolumeIntegralFluxDifferencing @inline function source_terms_split_covariant(u, x, t, aux_vars, equations::CovariantShallowWaterEquations2D) h, h_vcon1, h_vcon2 = u diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8130451..dcfde48 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -94,12 +94,12 @@ function transform_initial_condition(initial_condition, ::AbstractCovariantEquat function initial_condition_transformed(x, t, aux_vars, equations) return Trixi.prim2cons(global2contravariant(initial_condition(x, t, equations), aux_vars, equations), aux_vars, - equations) + equations) end return initial_condition_transformed end -# Version for standard Cartesian formulations +# Default version for standard Cartesian formulations function transform_initial_condition(initial_condition, ::AbstractEquations) function initial_condition_transformed(x, t, equations) return Trixi.prim2cons(initial_condition(x, t, equations), equations) @@ -138,7 +138,7 @@ function cartesian2global(u, x, equations::AbstractEquations) return u end -# cons2cons method which takes in unused aux_vars variable +# Default cons2cons and prim2cons methods which take in unused aux_vars variable @inline Trixi.cons2cons(u, aux_vars, ::AbstractEquations) = u @inline Trixi.prim2cons(u, aux_vars, ::AbstractEquations) = u @@ -148,6 +148,8 @@ end # contravariant metric tensor components, and the upper-triangular Christoffel symbols of # the second kind @inline have_aux_node_vars(::AbstractCovariantEquations) = True() + +# Add up the total number of auxiliary variables for equations in covariant form @inline function n_aux_node_vars(::AbstractCovariantEquations{NDIMS, NDIMS_AMBIENT}) where { NDIMS, @@ -266,6 +268,8 @@ end return 0.5f0 * (flux_ll + flux_rr) end +# Local Lax-Friedrichs dissipation for abstract covariant equations, where dissipation is +# applied to all conservative variables and the wave speed may depend on auxiliary variables @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, aux_vars_ll, aux_vars_rr, orientation_or_normal_direction, @@ -276,7 +280,7 @@ end return -0.5f0 * sqrtG * λ * (u_rr - u_ll) end -# Dummy two-point nonconservative flux for weak form +# Define the two-point nonconservative flux for weak form to be a no-op @inline function flux_nonconservative_weak_form(u_ll::SVector{NVARS, RealT}, u_rr::SVector{NVARS, RealT}, aux_vars_ll, aux_vars_rr, @@ -293,8 +297,8 @@ end return zeros(SVector{NVARS, RealT}) end -# Convert a vector from a global spherical to Cartesian basis representation, where we note -# that the radial component is not necessarily zero +# Convert a vector from a global spherical to Cartesian basis representation. A tangent +# vector will have vrad = 0. @inline function spherical2cartesian(vlon, vlat, vrad, x) # compute longitude and latitude lon, lat = atan(x[2], x[1]), asin(x[3] / norm(x)) @@ -307,12 +311,11 @@ end vx = -sinlon * vlon - sinlat * coslon * vlat + coslat * coslon * vrad vy = coslon * vlon - sinlat * sinlon * vlat + coslat * sinlon * vrad vz = coslat * vlat + sinlat * vrad - return vx, vy, vz end -# Convert a vector from a global Cartesian to spherical basis representation, where we note -# that the radial component is not necessarily zero +# Convert a vector from a global Cartesian to spherical basis representation. A tangent +# vector will have vrad = 0. @inline function cartesian2spherical(vx, vy, vz, x) # compute longitude and latitude lon, lat = atan(x[2], x[1]), asin(x[3] / norm(x)) @@ -324,7 +327,7 @@ end # return spherical components vlon = -sinlon * vx + coslon * vy vlat = -sinlat * coslon * vx - sinlat * sinlon * vy + coslat * vz - vrad = coslat * coslon * vx + coslat * sinlon * vy + sinlat * vz # zero for any tangent vector + vrad = coslat * coslon * vx + coslat * sinlon * vy + sinlat * vz return vlon, vlat, vrad end From b9c3591e4b507cc8798704867e643ab203d76e49 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Tue, 31 Dec 2024 00:42:35 -0500 Subject: [PATCH 11/21] both spherical and cartesian works now --- src/equations/covariant_shallow_water.jl | 14 ++++++------ src/equations/equations.jl | 4 ++-- src/equations/reference_data.jl | 4 ++-- test/test_2d_shallow_water_covariant.jl | 27 ++++++++++++++++++++++-- 4 files changed, 36 insertions(+), 13 deletions(-) diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index 07dd9c2..8e4c44e 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -349,23 +349,23 @@ end # If the initial velocity field is defined in Cartesian coordinates and the chosen global # coordinate system is spherical, perform the appropriate conversion @inline function cartesian2global(u, x, - ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) + ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) h_vlon, h_vlat, h_vrad = cartesian2spherical(u[2], u[3], u[4], x) -return SVector(u[1], h_vlon, h_vlat, h_vrad) + return SVector(u[1], h_vlon, h_vlat, h_vrad) end # If the initial velocity field is defined in spherical coordinates and the chosen global # coordinate system is Cartesian, perform the appropriate conversion @inline function spherical2global(u, x, - ::CovariantShallowWaterEquations2D{GlobalCartesianCoordinates}) -h_vx, h_vy, h_vz = spherical2cartesian(u[2], u[3], u[4], x) -return SVector(u[1], h_vx, h_vy, h_vz) + ::CovariantShallowWaterEquations2D{GlobalCartesianCoordinates}) + h_vx, h_vy, h_vz = spherical2cartesian(u[2], u[3], u[4], x) + return SVector(u[1], h_vx, h_vy, h_vz) end # If the initial velocity field is defined in spherical coordinates and the chosen global # coordinate system is spherical, do not convert @inline function spherical2global(u, x, - ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) -return u + ::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) + return u end end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index dcfde48..a1ab5ff 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -94,7 +94,7 @@ function transform_initial_condition(initial_condition, ::AbstractCovariantEquat function initial_condition_transformed(x, t, aux_vars, equations) return Trixi.prim2cons(global2contravariant(initial_condition(x, t, equations), aux_vars, equations), aux_vars, - equations) + equations) end return initial_condition_transformed end @@ -327,7 +327,7 @@ end # return spherical components vlon = -sinlon * vx + coslon * vy vlat = -sinlat * coslon * vx - sinlat * sinlon * vy + coslat * vz - vrad = coslat * coslon * vx + coslat * sinlon * vy + sinlat * vz + vrad = coslat * coslon * vx + coslat * sinlon * vy + sinlat * vz return vlon, vlat, vrad end diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 4c531f6..5459a63 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -68,8 +68,8 @@ This problem is adapted from Case 1 of the test suite described in the following return cartesian2global(SVector(h, vx, vy, vz, 0.0f0), x, equations) end -function initial_condition_geostrophic_balance(x, t, - equations::CovariantShallowWaterEquations2D) +function initial_condition_geostrophic_balance(x, t, + equations::CovariantShallowWaterEquations2D) (; gravity, rotation_rate) = equations radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) diff --git a/test/test_2d_shallow_water_covariant.jl b/test/test_2d_shallow_water_covariant.jl index d5e54d6..48a24c5 100644 --- a/test/test_2d_shallow_water_covariant.jl +++ b/test/test_2d_shallow_water_covariant.jl @@ -7,7 +7,7 @@ include("test_trixiatmo.jl") EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") -@trixiatmo_testset "spherical shallow water, covariant weak form" begin +@trixiatmo_testset "Spherical shallow water, covariant weak form, global spherical coords" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_covariant_cubed_sphere.jl"), l2=[ @@ -19,7 +19,30 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") 1.4786544678290738, 0.0013754600033514114, 0.0007564014737144256 - ]) + ], global_coordinate_system=GlobalSphericalCoordinates()) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +@trixiatmo_testset "Spherical shallow water, covariant weak form, global Cartesian coords" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_covariant_cubed_sphere.jl"), + l2=[ + 0.3065314463985936, + 0.00019984467582352902, + 8.767819502807507e-5 + ], + linf=[ + 1.4786544678290738, + 0.0013754600033514114, + 0.0007564014737144256 + ], global_coordinate_system=GlobalCartesianCoordinates()) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 9c720dd276eb7d5cec774df953c4be5c2350d5d3 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Thu, 2 Jan 2025 00:29:15 -0500 Subject: [PATCH 12/21] add EC formulation --- ..._shallowwater_covariant_cubed_sphere_EC.jl | 75 ++++++++++++ src/TrixiAtmo.jl | 3 +- src/equations/covariant_shallow_water.jl | 108 +++++++++--------- src/equations/equations.jl | 7 +- .../dg_2d_manifold_in_3d_covariant.jl | 50 ++++++++ test/test_2d_shallow_water_covariant.jl | 19 +++ 6 files changed, 203 insertions(+), 59 deletions(-) create mode 100644 examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl diff --git a/examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl b/examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl new file mode 100644 index 0000000..9813eb0 --- /dev/null +++ b/examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl @@ -0,0 +1,75 @@ +############################################################################### +# DGSEM for the shallow water equations on the cubed sphere +############################################################################### + +using OrdinaryDiffEq, Trixi, TrixiAtmo + +############################################################################### +# Parameters + +initial_condition = initial_condition_geostrophic_balance +polydeg = 3 +cells_per_dimension = 5 +n_saves = 10 + +############################################################################### +# Spatial discretization + +tspan = (0.0, 1.0 * SECONDS_PER_DAY) + +mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg, + initial_refinement_level = 0, + element_local_mapping = true) + +equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, + EARTH_ROTATION_RATE, + global_coordinate_system = GlobalSphericalCoordinates()) + +# Create DG solver with polynomial degree = p +volume_flux = (flux_split_covariant, flux_nonconservative_split_covariant) +surface_flux = (flux_split_covariant, flux_nonconservative_split_covariant) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +initial_condition_transformed = transform_initial_condition(initial_condition, equations) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver, + source_terms = source_terms_split_covariant) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0 to T +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 50, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, + solution_variables = cons2cons) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.4) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 100.0, save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 7f3b433..7a4ea54 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -37,7 +37,8 @@ export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, export GlobalCartesianCoordinates, GlobalSphericalCoordinates export flux_chandrashekar, flux_LMARS, flux_split_covariant, flux_nonconservative_weak_form, - flux_nonconservative_split_covariant + flux_nonconservative_split_covariant, flux_split_covariant_lax_friedrichs, + source_terms_split_covariant export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal, lake_at_rest_error, source_terms_lagrange_multiplier, source_terms_weak_form, diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index 8e4c44e..0a7f2c0 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -1,10 +1,14 @@ @muladd begin #! format: noindent +""" + CovariantShallowWaterEquations2D{GlobalCoordinateSystem} <: + AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} +""" struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} <: AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} - gravity::RealT - rotation_rate::RealT + gravity::RealT # acceleration due to gravity + rotation_rate::RealT # rotation rate for Coriolis term global_coordinate_system::GlobalCoordinateSystem function CovariantShallowWaterEquations2D(gravity::RealT, rotation_rate::RealT; @@ -14,35 +18,37 @@ struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} < end end +# Our implementation of flux-differencing formulation uses nonconservative terms, but the +# standard weak form does not. To handle both options, we have defined a dummy kernel for +# the nonconservative terms that does nothing when VolumeIntegralWeakForm is used with a +# nonconservative system. +Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True() + # The conservative variables are the height and contravariant momentum components function Trixi.varnames(::typeof(cons2cons), ::CovariantShallowWaterEquations2D) return ("h", "h_vcon1", "h_vcon2") end -# Convenience function to extract the velocity -function velocity(u, ::CovariantShallowWaterEquations2D) - return SVector(u[2] / u[1], u[3] / u[1]) +# The primitive variables are the height and contravariant velocity components +function Trixi.varnames(::typeof(cons2prim), ::CovariantShallowWaterEquations2D) + return ("h", "vcon1", "vcon2") end -# Convenience function to extract the momentum -function momentum(u, ::CovariantShallowWaterEquations2D) - return SVector(u[2], u[3]) +# The change of variables contravariant2global converts the two local contravariant vector +# components u[2] and u[3] to the three global vector components specified by +# equations.global_coordinate_system (e.g. spherical or Cartesian). This transformation +# works for both primitive and conservative variables, although varnames refers +# specifically to transformations from conservative variables. +function Trixi.varnames(::typeof(contravariant2global), + ::CovariantShallowWaterEquations2D) + return ("h", "h_vglo1", "h_vglo2", "h_vglo3") end -# Our implementation of flux-differencing formulation uses nonconservative terms, but the -# standard weak form does not. To handle both options, we have defined a dummy kernel for -# the nonconservative terms that does nothing when VolumeIntegralWeakForm is used with a -# nonconservative system. -Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True() - -# Entropy function (total energy per unit volume) -@inline function Trixi.entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D) - h, h_vcon1, h_vcon2 = u - Gcov = metric_covariant(aux_vars, equations) - vcon = SVector(h_vcon1 / h, h_vcon2 / h) - vcov = Gcov * vcon - return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2) -end +# Convenience functions to extract variables from state vector +@inline waterheight(u, ::CovariantShallowWaterEquations2D) = u[1] +@inline velocity(u, ::CovariantShallowWaterEquations2D) = SVector(u[2] / u[1], + u[3] / u[1]) +@inline momentum(u, ::CovariantShallowWaterEquations2D) = SVector(u[2], u[3]) @inline function Trixi.cons2prim(u, aux_vars, ::CovariantShallowWaterEquations2D) h, h_vcon1, h_vcon2 = u @@ -54,21 +60,12 @@ end return SVector(h, h * vcon1, h * vcon2) end -# Entropy variables (partial derivatives of entropy with respect to conservative variables) @inline function Trixi.cons2entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D) - h, h_vcon1, h_vcon2 = u - Gcov = metric_covariant(aux_vars, equations) - vcon = SVector(h_vcon1 / h, h_vcon2 / h) - vcov = Gcov * vcon - w1 = equations.gravity * h - 0.5f0 * dot(vcov, vcon) - return SVector{3}(w1, vcov[1], vcov[2]) -end - -# Height and three global momentum components -function Trixi.varnames(::typeof(contravariant2global), - ::CovariantShallowWaterEquations2D) - return ("h", "h_vglo1", "h_vglo2", "h_vglo3") + h = waterheight(u, equations) + vcon = velocity(u, equations) + vcov = metric_covariant(aux_vars, equations) * vcon + return SVector{3}(equations.gravity * h - 0.5f0 * dot(vcov, vcon), vcov[1], vcov[2]) end # Convert contravariant momentum components to the global coordinate system @@ -84,6 +81,15 @@ end vcon1, vcon2 = basis_contravariant(aux_vars, equations) * SVector(u[2], u[3], u[4]) return SVector(u[1], vcon1, vcon2) end + +# Entropy function (total energy per unit volume) +@inline function Trixi.entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D) + h = waterheight(u, equations) + vcon = velocity(u, equations) + vcov = metric_covariant(aux_vars, equations) * vcon + return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2) +end + # The flux for the covariant form takes in the element container and node/element indices # in order to give the flux access to the geometric information @inline function Trixi.flux(u, aux_vars, orientation::Integer, @@ -171,16 +177,15 @@ end # Geometric and Coriolis sources for a rotating sphere with VolumeIntegralWeakForm @inline function source_terms_weak_form(u, x, t, aux_vars, equations::CovariantShallowWaterEquations2D) - h, h_vcon1, h_vcon2 = u - # Geometric variables Gcon = metric_contravariant(aux_vars, equations) (Gamma1, Gamma2) = christoffel_symbols(aux_vars, equations) J = area_element(aux_vars, equations) # Physical variables - h_vcon = SVector{2}(h_vcon1, h_vcon2) - v_con = SVector{2}(h_vcon1 / h, h_vcon2 / h) + h = waterheight(u, equations) + h_vcon = momentum(u, equations) + v_con = velocity(u, equations) # Doubly-contravariant flux tensor T = h_vcon * v_con' + 0.5f0 * equations.gravity * h^2 * Gcon @@ -192,8 +197,8 @@ end s_geo = SVector(sum(Gamma1 .* T), sum(Gamma2 .* T)) # Combined source terms - source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon1 - Gcon[1, 1] * h_vcon2) - source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon1 - Gcon[2, 1] * h_vcon2) + source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon[1] - Gcon[1, 1] * h_vcon[2]) + source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon[1] - Gcon[2, 1] * h_vcon[2]) # Do not scale by Jacobian since apply_jacobian! is called before this return SVector(zero(eltype(u)), -source_1, -source_2) @@ -202,7 +207,6 @@ end # Geometric and Coriolis sources for a rotating sphere with VolumeIntegralFluxDifferencing @inline function source_terms_split_covariant(u, x, t, aux_vars, equations::CovariantShallowWaterEquations2D) - h, h_vcon1, h_vcon2 = u # Geometric variables Gcov = metric_covariant(aux_vars, equations) @@ -211,8 +215,8 @@ end J = area_element(aux_vars, equations) # Physical variables - h_vcon = SVector{2}(h_vcon1, h_vcon2) - v_con = SVector{2}(h_vcon1 / h, h_vcon2 / h) + h_vcon = momentum(u, equations) + v_con = velocity(u, equations) # Doubly-contravariant and mixed inertial flux tensors h_vcon_vcon = h_vcon * v_con' @@ -226,8 +230,8 @@ end Gcon * (Gamma1 * h_vcov_vcon[1, :] + Gamma2 * h_vcov_vcon[2, :])) # Combined source terms - source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon1 - Gcon[1, 1] * h_vcon2) - source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon1 - Gcon[2, 1] * h_vcon2) + source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon[1] - Gcon[1, 1] * h_vcon[2]) + source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon[1] - Gcon[2, 1] * h_vcon[2]) # Do not scale by Jacobian since apply_jacobian! is called before this return SVector(zero(eltype(u)), -source_1, -source_2) @@ -246,13 +250,10 @@ end Gcon_ll = metric_contravariant(aux_vars_ll, equations) Gcon_rr = metric_contravariant(aux_vars_rr, equations) - phi_ll = max(h_ll * equations.gravity, 0) - phi_rr = max(h_rr * equations.gravity, 0) - return max(abs(h_vcon_ll[orientation] / h_ll) + - sqrt(Gcon_ll[orientation, orientation] * phi_ll), + sqrt(Gcon_ll[orientation, orientation] * h_ll * equations.gravity), abs(h_vcon_rr[orientation] / h_rr) + - sqrt(Gcon_rr[orientation, orientation] * phi_rr)) + sqrt(Gcon_rr[orientation, orientation] * h_rr * equations.gravity)) end # Maximum wave speeds with respect to the covariant basis @@ -260,9 +261,8 @@ end equations::CovariantShallowWaterEquations2D) h, h_vcon1, h_vcon2 = u Gcon = metric_contravariant(aux_vars, equations) - phi = max(h * equations.gravity, 0) - return abs(h_vcon1 / h) + sqrt(Gcon[1, 1] * phi), - abs(h_vcon2 / h) + sqrt(Gcon[2, 2] * phi) + return abs(h_vcon1 / h) + sqrt(Gcon[1, 1] * h * equations.gravity), + abs(h_vcon2 / h) + sqrt(Gcon[2, 2] * h * equations.gravity) end # Rossby-Haurwitz wave diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a1ab5ff..3b394da 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -184,7 +184,7 @@ end "basis_contravariant[1,2]", # e₁ ⋅ a² "basis_contravariant[2,2]", # e₂ ⋅ a² "basis_contravariant[3,2]", # e₃ ⋅ a² - "area_element", # √G = √(G₁₁G₂₂ - G₁₂G₂₁) = ||a₁ × a₂|| + "area_element", # J = √(G₁₁G₂₂ - G₁₂G₂₁) = ||a₁ × a₂|| "metric_covariant[1,1]", # G₁₁ "metric_covariant[1,2]", # G₁₂ = G₂₁ "metric_covariant[2,2]", # G₂₂ @@ -216,7 +216,7 @@ end aux_vars[11], aux_vars[12]) end -# Extract the area element √G = (det(AᵀA))^(1/2) from the auxiliary variables +# Extract the area element J = (det(AᵀA))^(1/2) from the auxiliary variables @inline function area_element(aux_vars, ::AbstractCovariantEquations{2}) return aux_vars[13] end @@ -274,10 +274,9 @@ end aux_vars_rr, orientation_or_normal_direction, equations::AbstractCovariantEquations) - sqrtG = area_element(aux_vars_ll, equations) λ = dissipation.max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr, orientation_or_normal_direction, equations) - return -0.5f0 * sqrtG * λ * (u_rr - u_ll) + return -0.5f0 * area_element(aux_vars_ll, equations) * λ * (u_rr - u_ll) end # Define the two-point nonconservative flux for weak form to be a no-op diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl index ccc1f77..81fa008 100644 --- a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl @@ -123,6 +123,56 @@ end Trixi.weak_form_kernel!(du, u, element, mesh, False(), equations, dg, cache, alpha) end +# Non-conservative flux differencing kernel which uses contravariant flux components, +# passing the geometric information contained in the auxiliary variables to the flux +# function +@inline function Trixi.flux_differencing_kernel!(du, u, element, mesh::P4estMesh{2}, + nonconservative_terms::True, + equations::AbstractCovariantEquations{2}, + volume_flux, dg::DGSEM, cache, + alpha = true) + (; derivative_split) = dg.basis + (; aux_node_vars) = cache.auxiliary_variables + symmetric_flux, nonconservative_flux = volume_flux + + # Apply the symmetric flux as usual + Trixi.flux_differencing_kernel!(du, u, element, mesh, False(), equations, + symmetric_flux, dg, cache, alpha) + + for j in eachnode(dg), i in eachnode(dg) + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + + # ξ¹ direction + integral_contribution = zero(u_node) + for ii in eachnode(dg) + u_node_ii = Trixi.get_node_vars(u, equations, dg, ii, j, element) + aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, + element) + flux1 = nonconservative_flux(u_node, u_node_ii, aux_node, aux_node_ii, 1, + equations) + integral_contribution = integral_contribution + + derivative_split[i, ii] * flux1 + end + + # ξ² direction + for jj in eachnode(dg) + u_node_jj = Trixi.get_node_vars(u, equations, dg, i, jj, element) + aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, + element) + flux2 = nonconservative_flux(u_node, u_node_jj, aux_node, aux_node_jj, 2, + equations) + integral_contribution = integral_contribution + + derivative_split[j, jj] * flux2 + end + + Trixi.multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution, + equations, dg, i, j, element) + end + + return nothing +end + # Pointwise interface flux in local coordinates for problems without nonconservative terms @inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, nonconservative_terms::False, diff --git a/test/test_2d_shallow_water_covariant.jl b/test/test_2d_shallow_water_covariant.jl index 48a24c5..0e0cb68 100644 --- a/test/test_2d_shallow_water_covariant.jl +++ b/test/test_2d_shallow_water_covariant.jl @@ -52,4 +52,23 @@ end @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 end end + +@trixiatmo_testset "Spherical shallow water, entropy-conservative covariant form" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_covariant_cubed_sphere_EC.jl"), + l2=[0.9995330728180628, + 0.000111592713364556, + 9.12269016730292e-5], + linf=[3.659066044233896, + 0.0012784722821582717, + 0.0012784722821565994]) + # 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 From 7317ed34a7409614009133169f48a77e2cd57bca Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Thu, 2 Jan 2025 17:08:19 -0500 Subject: [PATCH 13/21] add tests and improve elixirs/docs/cases --- ...lowwater_covariant_geostrophic_balance.jl} | 3 +- ...llowwater_covariant_rossby_haurwitz_EC.jl} | 10 +- src/TrixiAtmo.jl | 7 +- src/callbacks_step/analysis_covariant.jl | 19 ++++ src/equations/covariant_shallow_water.jl | 81 ---------------- src/equations/reference_data.jl | 96 +++++++++++++++---- test/test_2d_shallow_water_covariant.jl | 42 +++----- 7 files changed, 118 insertions(+), 140 deletions(-) rename examples/{elixir_shallowwater_covariant_cubed_sphere.jl => elixir_shallowwater_covariant_geostrophic_balance.jl} (99%) rename examples/{elixir_shallowwater_covariant_cubed_sphere_EC.jl => elixir_shallowwater_covariant_rossby_haurwitz_EC.jl} (94%) diff --git a/examples/elixir_shallowwater_covariant_cubed_sphere.jl b/examples/elixir_shallowwater_covariant_geostrophic_balance.jl similarity index 99% rename from examples/elixir_shallowwater_covariant_cubed_sphere.jl rename to examples/elixir_shallowwater_covariant_geostrophic_balance.jl index 825b70c..0eadcb6 100644 --- a/examples/elixir_shallowwater_covariant_cubed_sphere.jl +++ b/examples/elixir_shallowwater_covariant_geostrophic_balance.jl @@ -11,12 +11,11 @@ initial_condition = initial_condition_geostrophic_balance polydeg = 3 cells_per_dimension = 5 n_saves = 10 +tspan = (0.0, 1.0 * SECONDS_PER_DAY) ############################################################################### # Spatial discretization -tspan = (0.0, 1.0 * SECONDS_PER_DAY) - mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg, initial_refinement_level = 0, element_local_mapping = true) diff --git a/examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl similarity index 94% rename from examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl rename to examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl index 9813eb0..cf34d19 100644 --- a/examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl +++ b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl @@ -7,23 +7,22 @@ using OrdinaryDiffEq, Trixi, TrixiAtmo ############################################################################### # Parameters -initial_condition = initial_condition_geostrophic_balance +initial_condition = initial_condition_rossby_haurwitz polydeg = 3 cells_per_dimension = 5 n_saves = 10 +tspan = (0.0, 1.0 * SECONDS_PER_DAY) ############################################################################### # Spatial discretization -tspan = (0.0, 1.0 * SECONDS_PER_DAY) - mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg, initial_refinement_level = 0, element_local_mapping = true) equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, - global_coordinate_system = GlobalSphericalCoordinates()) + global_coordinate_system = GlobalCartesianCoordinates()) # Create DG solver with polynomial degree = p volume_flux = (flux_split_covariant, flux_nonconservative_split_covariant) @@ -51,7 +50,8 @@ summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results analysis_callback = AnalysisCallback(semi, interval = 50, save_analysis = true, - extra_analysis_errors = (:conservation_error,)) + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (entropy,)) # The SaveSolutionCallback allows to save the solution to a file in regular intervals save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 7a4ea54..a4364b8 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -34,6 +34,7 @@ include("callbacks_step/callbacks_step.jl") export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D, CovariantLinearAdvectionEquation2D, CovariantShallowWaterEquations2D + export GlobalCartesianCoordinates, GlobalSphericalCoordinates export flux_chandrashekar, flux_LMARS, flux_split_covariant, flux_nonconservative_weak_form, @@ -52,11 +53,11 @@ export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProdu export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, SECONDS_PER_DAY -export global2contravariant, contravariant2global, spherical2cartesian, +export global2contravariant, contravariant2global, spherical2cartesian, cartesian2spherical, transform_initial_condition -export initial_condition_gaussian, initial_condition_gaussian_cartesian, - initial_condition_geostrophic_balance +export initial_condition_gaussian, initial_condition_geostrophic_balance, + initial_condition_rossby_haurwitz export examples_dir end # module TrixiAtmo diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index 1c5849e..935bf69 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -1,6 +1,25 @@ @muladd begin #! format: noindent +# When the equations are of type AbstractCovariantEquations, the functions which we would +# like to integrate depend on the solution as well as the auxiliary variables +function Trixi.integrate(func::Func, u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + equations::AbstractCovariantEquations{2}, dg::DG, + cache; normalize = true) where {Func} + (; aux_node_vars) = cache.auxiliary_variables + + Trixi.integrate_via_indices(u, mesh, equations, dg, cache; + normalize = normalize) do u, i, j, element, equations, + dg + u_local = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_local = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + return func(u_local, aux_local, equations) + end +end + # For the covariant form, we want to integrate using the exact area element # √G = (det(AᵀA))^(1/2), which is stored in cache.auxiliary_variables, not the approximate # area element used in the Cartesian formulation, which stored in cache.elements diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index 0a7f2c0..c6af48b 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -265,87 +265,6 @@ end abs(h_vcon2 / h) + sqrt(Gcon[2, 2] * h * equations.gravity) end -# Rossby-Haurwitz wave -function initial_condition_rossby_haurwitz(x, t, - equations::CovariantShallowWaterEquations2D{<:GlobalSphericalCoordinates}) - (; gravity, rotation_rate) = equations - - radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - lon = atan(x[2], x[1]) - lat = asin(x[3] / radius) - - h_0 = 8.0f3 - K = 7.848f-6 - R = 4.0f0 - - A = 0.5f0 * K * (2 * rotation_rate + K) * (cos(lat))^2 + - 0.25f0 * K^2 * (cos(lat))^(2 * R) * - ((R + 1) * (cos(lat))^2 + - (2 * R^2 - R - 2) - 2 * R^2 / ((cos(lat))^2)) - B = 2 * (rotation_rate + K) * K / ((R + 1) * (R + 2)) * (cos(lat))^R * - ((R^2 + 2R + 2) - (R + 1)^2 * (cos(lat))^2) - C = 0.25f0 * K^2 * (cos(lat))^(2 * R) * ((R + 1) * (cos(lat))^2 - (R + 2)) - - h = h_0 + - (1 / gravity) * - (radius^2 * A + radius^2 * B * cos(R * lon) + radius^2 * C * cos(2 * R * lon)) - - vlon = radius * K * cos(lat) + - radius * K * (cos(lat))^(R - 1) * (R * (sin(lat))^2 - (cos(lat))^2) * - cos(R * lon) - vlat = -radius * K * R * (cos(lat))^(R - 1) * sin(lat) * sin(R * lon) - - # convert to conservative variables - return SVector(h, h * vlon, h * vlat) -end - -@inline function galewsky_velocity(θ, u_0, θ_0, θ_1) - if (θ_0 < θ) && (θ < θ_1) - u = u_0 / exp(-4 / (θ_1 - θ_0)^2) * exp(1 / (θ - θ_0) * 1 / (θ - θ_1)) - else - u = zero(θ) - end - return u -end - -@inline function galewsky_integrand(θ, u_0, θ_0, θ_1, a, - equations::CovariantShallowWaterEquations2D) - (; rotation_rate) = equations - u = galewsky_velocity(θ, u_0, θ_0, θ_1) - return u * (2 * rotation_rate * sin(θ) + u * tan(θ) / a) -end - -function initial_condition_barotropic_instability(x, t, - equations::CovariantShallowWaterEquations2D{GlobalSphericalCoordinates}) - (; gravity) = equations - realT = eltype(x) - radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - lon = atan(x[2], x[1]) - lat = asin(x[3] / radius) - - # compute zonal and meridional velocity components - u_0 = 80.0f0 - lat_0 = convert(realT, π / 7) - lat_1 = convert(realT, π / 2) - lat_0 - vlon = galewsky_velocity(lat, u_0, lat_0, lat_1) - vlat = zero(eltype(x)) - - # numerically integrate (here we use the QuadGK package) to get height - galewsky_integral, _ = quadgk(latp -> galewsky_integrand(latp, u_0, lat_0, lat_1, - radius, - equations), -π / 2, lat) - h = 10158.0f0 - radius / gravity * galewsky_integral - - # add perturbation to initiate instability - α, β = convert(realT, 1 / 3), convert(realT, 1 / 15) - lat_2 = convert(realT, π / 4) - if (-π < lon) && (lon < π) - h = h + 120.0f0 * cos(lat) * exp(-((lon / α)^2)) * exp(-((lat_2 - lat) / β)^2) - end - # convert to conservative variables - return SVector(h, h * vlon, h * vlat) -end - # If the initial velocity field is defined in Cartesian coordinates and the chosen global # coordinate system is spherical, perform the appropriate conversion @inline function cartesian2global(u, x, diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 5459a63..38f3297 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -1,8 +1,8 @@ @muladd begin #! format: noindent -# Physical constants -const EARTH_RADIUS = 6.37122e6 # m +# Physical constants in SI units, with values taken from the Williamson et al. test suite +const EARTH_RADIUS = 6.37122e6 # m const EARTH_GRAVITATIONAL_ACCELERATION = 9.80616 # m/s² const EARTH_ROTATION_RATE = 7.292e-5 # rad/s const SECONDS_PER_DAY = 8.64e4 @@ -11,7 +11,7 @@ const SECONDS_PER_DAY = 8.64e4 initial_condition_gaussian(x, t, equations) This Gaussian bell case is a smooth initial condition suitable for testing the convergence -of discretizations of the linear advection equation on a spherical domain of radius $6. +of discretizations of the linear advection equation on a spherical domain of radius $a = 6. 37122 \times 10^3\ \mathrm{m}$, representing the surface of the Earth. Denoting the Euclidean norm as $\lVert \cdot \rVert$, the initial height field is given by ```math @@ -33,9 +33,8 @@ This problem is adapted from Case 1 of the test suite described in the following spherical geometry. Journal of Computational Physics, 102(1):211-224. [DOI: 10.1016/S0021-9991(05)80016-6](https://doi.org/10.1016/S0021-9991(05)80016-6) """ -@inline function initial_condition_gaussian(x, t, equations::AbstractEquations) +@inline function initial_condition_gaussian(x, t, equations) RealT = eltype(x) - a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) # radius of the sphere omega = convert(RealT, 2π) / (12 * SECONDS_PER_DAY) # angular velocity alpha = convert(RealT, π / 4) # angle of rotation @@ -63,27 +62,88 @@ This problem is adapted from Case 1 of the test suite described in the following # get Cartesian velocity components vx, vy, vz = omega * Trixi.cross(axis, x) - # Prescribe the rotated bell shape and Cartesian velocity components. - # The last variable is the bottom topography, which we set to zero + # Convert primitive variables from Cartesian coordinates to the chosen global + # coordinate system, which depends on the equation type return cartesian2global(SVector(h, vx, vy, vz, 0.0f0), x, equations) end -function initial_condition_geostrophic_balance(x, t, - equations::CovariantShallowWaterEquations2D) - (; gravity, rotation_rate) = equations +@doc raw""" + initial_condition_geostrophic_balance(x, t, equations) - radius = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - lat = asin(x[3] / radius) +Steady geostrophic balance for the spherical shallow water equations, corresponding to a +purely zonal velocity field given as a function of the latitude $\theta$ by +$v_\lambda(\theta) = v_0 \cos\theta$, where we define +$v_0 = 2\pi a \cos(\theta) / 12 \ \mathrm{days}$ in terms of the Earth's radius +$a = 6.37122 \times 10^3\ \mathrm{m}$. The height field then varies with the latitude as +```math +h(\theta) = 1/g \Big(gh_0 - \Big(a \omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), +``` +where $gh_0 = 2.94 \times 10^4 \ \mathrm{m}^2/\mathrm{s}^2$, +$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and $\omega = 7.292e-5 \mathrm{s}^{-1}$. This +problem corresponds to Case 2 of the test suite described in the following paper: +- D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A + standard test set for numerical approximations to the shallow water equations in + spherical geometry. Journal of Computational Physics, 102(1):211-224. + [DOI: 10.1016/S0021-9991(05)80016-6](https://doi.org/10.1016/S0021-9991(05)80016-6) +""" +@inline function initial_condition_geostrophic_balance(x, t, equations) + RealT = eltype(x) + a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) # radius of the sphere + lat = asin(x[3] / a) # compute zonal and meridional components of the velocity - V = convert(eltype(x), 2π) * radius / (12 * SECONDS_PER_DAY) - vlon, vlat = V * cos(lat), zero(eltype(x)) + v_0 = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY) + vlon, vlat = v_0 * cos(lat), zero(eltype(x)) # compute geopotential height - h = 1 / gravity * - (2.94f4 - (radius * rotation_rate * V + 0.5f0 * V^2) * (sin(lat))^2) + h = 1 / EARTH_GRAVITATIONAL_ACCELERATION * + (2.94f4 - (a * EARTH_ROTATION_RATE * v_0 + 0.5f0 * v_0^2) * (sin(lat))^2) + + # Convert primitive variables from spherical coordinates to the chosen global + # coordinate system, which depends on the equation type + return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations) +end + +@doc raw""" + initial_condition_rossby_haurwitz(x, t, equations) + +Rossby-Haurwitz wave with wave number 4, corresponding to Case 5 of the test suite +described in the following paper: +- D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A + standard test set for numerical approximations to the shallow water equations in + spherical geometry. Journal of Computational Physics, 102(1):211-224. + [DOI: 10.1016/S0021-9991(05)80016-6](https://doi.org/10.1016/S0021-9991(05)80016-6) +""" +@inline function initial_condition_rossby_haurwitz(x, t, equations) + RealT = eltype(x) + a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) + lon = atan(x[2], x[1]) + lat = asin(x[3] / a) + + h_0 = 8.0f3 + K = 7.848f-6 + R = 4.0f0 + + A = 0.5f0 * K * (2 * EARTH_ROTATION_RATE + K) * (cos(lat))^2 + + 0.25f0 * K^2 * (cos(lat))^(2 * R) * + ((R + 1) * (cos(lat))^2 + (2 * R^2 - R - 2) - 2 * R^2 / ((cos(lat))^2)) + B = 2 * (EARTH_ROTATION_RATE + K) * K / ((R + 1) * (R + 2)) * (cos(lat))^R * + ((R^2 + 2R + 2) - (R + 1)^2 * (cos(lat))^2) + C = 0.25f0 * K^2 * (cos(lat))^(2 * R) * ((R + 1) * (cos(lat))^2 - (R + 2)) + + # compute geopotential height + h = h_0 + + (1 / EARTH_GRAVITATIONAL_ACCELERATION) * + (a^2 * A + a^2 * B * cos(R * lon) + a^2 * C * cos(2 * R * lon)) + + # compute zonal and meridional components of the velocity + vlon = a * K * cos(lat) + + a * K * (cos(lat))^(R - 1) * (R * (sin(lat))^2 - + (cos(lat))^2) * cos(R * lon) + vlat = -a * K * R * (cos(lat))^(R - 1) * sin(lat) * sin(R * lon) - # convert primitive variables from spherical coordinates to the chosen coordinate system - return spherical2global(SVector(h, vlon, vlat, zero(eltype(x))), x, equations) + # Convert primitive variables from spherical coordinates to the chosen global + # coordinate system, which depends on the equation type + return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations) end end # muladd diff --git a/test/test_2d_shallow_water_covariant.jl b/test/test_2d_shallow_water_covariant.jl index 0e0cb68..f9c4e6d 100644 --- a/test/test_2d_shallow_water_covariant.jl +++ b/test/test_2d_shallow_water_covariant.jl @@ -7,9 +7,9 @@ include("test_trixiatmo.jl") EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") -@trixiatmo_testset "Spherical shallow water, covariant weak form, global spherical coords" begin +@trixiatmo_testset "elixir_shallowwater_covariant_geostrophic_balance" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_covariant_cubed_sphere.jl"), + "elixir_shallowwater_covariant_geostrophic_balance.jl"), l2=[ 0.3065314463985936, 0.00019984467582352902, @@ -19,7 +19,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") 1.4786544678290738, 0.0013754600033514114, 0.0007564014737144256 - ], global_coordinate_system=GlobalSphericalCoordinates()) + ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -30,38 +30,18 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") end end -@trixiatmo_testset "Spherical shallow water, covariant weak form, global Cartesian coords" begin +@trixiatmo_testset "elixir_shallowwater_covariant_rossby_haurwitz_EC" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_covariant_cubed_sphere.jl"), + "elixir_shallowwater_covariant_rossby_haurwitz_EC.jl"), l2=[ - 0.3065314463985936, - 0.00019984467582352902, - 8.767819502807507e-5 + 265.98581316017135, + 0.17690835250329687, + 0.25383326240368487 ], linf=[ - 1.4786544678290738, - 0.0013754600033514114, - 0.0007564014737144256 - ], global_coordinate_system=GlobalCartesianCoordinates()) - # 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 shallow water, entropy-conservative covariant form" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_covariant_cubed_sphere_EC.jl"), - l2=[0.9995330728180628, - 0.000111592713364556, - 9.12269016730292e-5], - linf=[3.659066044233896, - 0.0012784722821582717, - 0.0012784722821565994]) + 579.0744773822007, + 0.5054767269383867, + 0.5628603103981209]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 4272119be61e6585b2762fd65f01b3fa2d7a58d6 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Thu, 2 Jan 2025 19:54:19 -0500 Subject: [PATCH 14/21] more detailed docs --- src/equations/reference_data.jl | 42 ++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 38f3297..0e3f404 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -72,15 +72,15 @@ end Steady geostrophic balance for the spherical shallow water equations, corresponding to a purely zonal velocity field given as a function of the latitude $\theta$ by -$v_\lambda(\theta) = v_0 \cos\theta$, where we define -$v_0 = 2\pi a \cos(\theta) / 12 \ \mathrm{days}$ in terms of the Earth's radius -$a = 6.37122 \times 10^3\ \mathrm{m}$. The height field then varies with the latitude as +$v_\lambda(\theta) = v_0 \cos\theta$, where we define $v_0 = 2\pi a / (12 \ \mathrm{days})$ +in terms of the Earth's radius $a = 6.37122 \times 10^3\ \mathrm{m}$. The height field +then varies with the latitude as ```math -h(\theta) = 1/g \Big(gh_0 - \Big(a \omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), +h(\theta) = \frac{1}{g} \Big(gh_0 - \Big(a \omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), ``` where $gh_0 = 2.94 \times 10^4 \ \mathrm{m}^2/\mathrm{s}^2$, -$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and $\omega = 7.292e-5 \mathrm{s}^{-1}$. This -problem corresponds to Case 2 of the test suite described in the following paper: +$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and +$\omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the test suite described in the following paper: - D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. @@ -107,8 +107,34 @@ end @doc raw""" initial_condition_rossby_haurwitz(x, t, equations) -Rossby-Haurwitz wave with wave number 4, corresponding to Case 5 of the test suite -described in the following paper: +Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and meridional velocity components are given by +```math +\begin{aligned} +v_\lambda(\lambda,\theta) &= a \omega \cos \theta+a K \cos ^{R-1} \theta +\left(R \sin ^2 \theta-\cos ^2 \theta\right) \cos (R \lambda)\\ +v_\theta(\lambda,\theta) &= -a K R \cos ^{R-1} \theta \sin \theta \sin (R \lambda) +\end{aligned} +``` +where $\omega = K = 7.848 \times 10^{-6} \ \mathrm{s}^{-1}$ and $R = 4$ are given +constants, and $a = 6.37122 \times 10^3\ \mathrm{m}$ is the Earth's radius. Taking +$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$, +and $h_0 = 8000 \ \mathrm{m}$ and defining the functions +```math +\begin{aligned} +A(\theta) &= \frac{\omega}{2}(2 \Omega+\omega) \cos ^2 \theta + +\frac{1}{4} K^2 \cos ^{2 R} \theta\Big((R+1) \cos ^2 \theta +\left(2 R^2-R-2\right) - +\big(2 R^2 / \cos^2 \theta\big) \Big), \\ +B(\theta) &= \frac{2(\Omega+\omega) K}{(R+1)(R+2)} \cos ^R \theta\big((R^2+2 R+2)- +(R+1)^2 \cos ^2 \theta\big), \\ +C(\theta) &= \frac{1}{4} K^2 \cos^{2 R} \theta\big((R+1) \cos ^2 \theta-(R+2)\big), +\end{aligned} +``` +the initial height field is given by +```math +h(\lambda,\theta) = h_0 + \frac{a^2}{g}\Big(A(\theta) + B(\theta)\cos(R\lambda) + +C(\theta)\cos(2R\lambda) \Big). +``` +This problem corresponds to Case 5 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. From 2a5d2a79581f19a25cb590ae62cab9783b987981 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Thu, 2 Jan 2025 19:58:56 -0500 Subject: [PATCH 15/21] fix typo in docs --- src/equations/reference_data.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 0e3f404..913bdab 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -76,11 +76,11 @@ $v_\lambda(\theta) = v_0 \cos\theta$, where we define $v_0 = 2\pi a / (12 \ \mat in terms of the Earth's radius $a = 6.37122 \times 10^3\ \mathrm{m}$. The height field then varies with the latitude as ```math -h(\theta) = \frac{1}{g} \Big(gh_0 - \Big(a \omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), +h(\theta) = \frac{1}{g} \Big(gh_0 - \Big(a \Omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), ``` where $gh_0 = 2.94 \times 10^4 \ \mathrm{m}^2/\mathrm{s}^2$, $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and -$\omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the test suite described in the following paper: +$\Omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the test suite described in the following paper: - D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. @@ -117,7 +117,7 @@ v_\theta(\lambda,\theta) &= -a K R \cos ^{R-1} \theta \sin \theta \sin (R \lambd ``` where $\omega = K = 7.848 \times 10^{-6} \ \mathrm{s}^{-1}$ and $R = 4$ are given constants, and $a = 6.37122 \times 10^3\ \mathrm{m}$ is the Earth's radius. Taking -$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$, +$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\Omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$, and $h_0 = 8000 \ \mathrm{m}$ and defining the functions ```math \begin{aligned} From 9021899757315063c4cb21d56706ec3bbd3b9400 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Thu, 2 Jan 2025 20:00:42 -0500 Subject: [PATCH 16/21] improve docs more --- src/equations/reference_data.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 913bdab..6fd1341 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -107,17 +107,17 @@ end @doc raw""" initial_condition_rossby_haurwitz(x, t, equations) -Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and meridional velocity components are given by +Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and meridional velocity components are given, respectively, by ```math \begin{aligned} v_\lambda(\lambda,\theta) &= a \omega \cos \theta+a K \cos ^{R-1} \theta -\left(R \sin ^2 \theta-\cos ^2 \theta\right) \cos (R \lambda)\\ -v_\theta(\lambda,\theta) &= -a K R \cos ^{R-1} \theta \sin \theta \sin (R \lambda) +\left(R \sin ^2 \theta-\cos ^2 \theta\right) \cos (R \lambda),\\ +v_\theta(\lambda,\theta) &= -a K R \cos ^{R-1} \theta \sin \theta \sin (R \lambda), \end{aligned} ``` where $\omega = K = 7.848 \times 10^{-6} \ \mathrm{s}^{-1}$ and $R = 4$ are given constants, and $a = 6.37122 \times 10^3\ \mathrm{m}$ is the Earth's radius. Taking -$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\Omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$, +$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\Omega = 7.292 \times 10^{-5} \ \mathrm{s}^{-1}$, and $h_0 = 8000 \ \mathrm{m}$ and defining the functions ```math \begin{aligned} From 0356b55240e9c9d7ea95a3ed862babfa11ff5da1 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Fri, 3 Jan 2025 16:05:22 -0500 Subject: [PATCH 17/21] add docstring for covariant SWE --- src/equations/covariant_advection.jl | 13 +++--- src/equations/covariant_shallow_water.jl | 53 +++++++++++++++++++++++- src/equations/reference_data.jl | 2 +- 3 files changed, 59 insertions(+), 9 deletions(-) diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index 5e08fea..da7a1e2 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -24,21 +24,20 @@ are spatially varying but assumed to be constant in time, so we do not apply any dissipation to such variables. The resulting system is then given on the reference element as ```math -\sqrt{G} \frac{\partial}{\partial t} +J \frac{\partial}{\partial t} \left[\begin{array}{c} h \\ v^1 \\ v^2 \end{array}\right] + \frac{\partial}{\partial \xi^1} -\left[\begin{array}{c} \sqrt{G} h v^1 \\ 0 \\ 0 \end{array}\right] +\left[\begin{array}{c} J h v^1 \\ 0 \\ 0 \end{array}\right] + \frac{\partial}{\partial \xi^2} -\left[\begin{array}{c} \sqrt{G} h v^2 \\ 0 \\ 0 \end{array}\right] +\left[\begin{array}{c} J h v^2 \\ 0 \\ 0 \end{array}\right] = \left[\begin{array}{c} 0 \\ 0 \\ 0 \end{array}\right], ``` -where $G$ is the determinant of the covariant metric tensor expressed as a 2 by 2 matrix -with entries $G_{ij} = \vec{a}_i \cdot \vec{a}_j$. Note that the variable advection -velocity components could alternatively be stored as auxiliary variables, similarly to the -geometric information. +where $J = \lVert\vec{a}^1 \times \vec{a}^2 \rVert$ is the area element. Note that the +variable advection velocity components could alternatively be stored as auxiliary +variables, similarly to the geometric information. """ struct CovariantLinearAdvectionEquation2D{GlobalCoordinateSystem} <: AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index c6af48b..29de7cc 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -1,9 +1,60 @@ @muladd begin #! format: noindent -""" +@doc raw""" CovariantShallowWaterEquations2D{GlobalCoordinateSystem} <: AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} + +Denoting the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) by +$\nabla_b$ and summing over repeated indices, the shallow water equations can be expressed +on a two-dimensional surface in three-dimensional ambient space as +```math +\begin{aligned} +\partial_t h + \nabla_b (hv^b) &= 0,\\ +\partial_t (hv^a) + \nabla_b (hv^av^b) + gh G^{ab}\partial_b(h + b) +&= -fJ G^{ab}\varepsilon_{bc} hv^c, +\end{aligned} +``` +where $h$ is the fluid height, $v^a$ and $G^{ab}$ are the contravariant velocity and metric +tensor components, $g$ is the gravitational constant, $f$ is the Coriolis parameter, and +$J$ is the area element. Combining the inertial and gravitational terms in order to define +the momentum flux components +```math +\tau^{ab} = hv^a v^b + \frac{1}{2}G^{ab}gh^2, +``` +the covariant shallow water equations can be expressed as a system of conservation laws +with a source term: +```math +J \frac{\partial}{\partial t} +\left[\begin{array}{c} h \\ hv^1 \\ hv^2 \end{array}\right] ++ +\frac{\partial}{\partial \xi^1} +\left[\begin{array}{c} J h v^1 \\ J \tau^{11} \\ J \tau^{12} \end{array}\right] ++ +\frac{\partial}{\partial \xi^2} +\left[\begin{array}{c} J h v^2 \\ J \tau^{21} \\ J \tau^{22} \end{array}\right] += J \left[\begin{array}{c} 0 \\ +-\Gamma^1_{ab}\tau^{ab} - f J \big(G^{12}hv^1 - G^{11}hv^2\big) \\ +-\Gamma^2_{ab}\tau^{ab} - f J \big(G^{22}hv^1 - G^{21}hv^2\big) + \end{array}\right]. +``` +TrixiAtmo.jl implements standard weak-form DG methods for both the above system, as well as +entropy-stable split forms based on a novel flux-differencing discretization of the +covariant derivative. + +## References +- M. Baldauf (2020). Discontinuous Galerkin solver for the shallow-water equations in + covariant form on the sphere and the ellipsoid. Journal of Computational Physics + 410:109384. [DOI: 10.1016/j.jcp.2020.109384](https://doi.org/10.1016/j.jcp.2020.109384) +- L. Bao, R. D. Nair, and H. M. Tufo (2014). A mass and momentum flux-form high-order + discontinuous Galerkin shallow water model on the cubed-sphere. A mass and momentum + flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. + Journal of Computational Physics 271:224-243. + [DOI: 10.1016/j.jcp.2013.11.033](https://doi.org/10.1016/j.jcp.2013.11.033) + +!!! warning "Experimental implementation" + The use of entropy-stable split-form/flux-differencing formulations for covariant + equations is an experimental feature and may change in future releases. """ struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} <: AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3} diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 6fd1341..1e3fd7a 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -1,7 +1,7 @@ @muladd begin #! format: noindent -# Physical constants in SI units, with values taken from the Williamson et al. test suite +# Physical constants in SI units (reference values from the Williamson et al. test suite) const EARTH_RADIUS = 6.37122e6 # m const EARTH_GRAVITATIONAL_ACCELERATION = 9.80616 # m/s² const EARTH_ROTATION_RATE = 7.292e-5 # rad/s From 192c84ba33ac537d3a59e01c52141bcb73761f3e Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Fri, 3 Jan 2025 16:08:11 -0500 Subject: [PATCH 18/21] sqrtG -> J --- src/callbacks_step/analysis_covariant.jl | 14 +++++++------- src/equations/covariant_advection.jl | 8 ++++---- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index 935bf69..c823594 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -21,7 +21,7 @@ function Trixi.integrate(func::Func, u, end # For the covariant form, we want to integrate using the exact area element -# √G = (det(AᵀA))^(1/2), which is stored in cache.auxiliary_variables, not the approximate +# J = (det(AᵀA))^(1/2), which is stored in cache.auxiliary_variables, not the approximate # area element used in the Cartesian formulation, which stored in cache.elements function Trixi.integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, @@ -42,10 +42,10 @@ function Trixi.integrate_via_indices(func::Func, u, for element in eachelement(dg, cache) for j in eachnode(dg), i in eachnode(dg) aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) - sqrtG = area_element(aux_node, equations) - integral += weights[i] * weights[j] * sqrtG * + J = area_element(aux_node, equations) + integral += weights[i] * weights[j] * J * func(u, i, j, element, equations, dg, args...) - total_volume += weights[i] * weights[j] * sqrtG + total_volume += weights[i] * weights[j] * J end end @@ -108,14 +108,14 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, diff = func(u_exact, equations) - func(u_numerical, equations) # For the L2 error, integrate with respect to area element stored in aux vars - sqrtG = area_element(aux_node, equations) - l2_error += diff .^ 2 * (weights[i] * weights[j] * sqrtG) + J = area_element(aux_node, equations) + l2_error += diff .^ 2 * (weights[i] * weights[j] * J) # Compute Linf error as usual linf_error = @. max(linf_error, abs(diff)) # Increment total volume according to the volume element stored in aux vars - total_volume += weights[i] * weights[j] * sqrtG + total_volume += weights[i] * weights[j] * J end end diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index da7a1e2..0347be2 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -93,9 +93,9 @@ end @inline function Trixi.flux(u, aux_vars, orientation::Integer, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u)) - sqrtG = area_element(aux_vars, equations) + J = area_element(aux_vars, equations) vcon = velocity(u, equations) - return SVector(sqrtG * u[1] * vcon[orientation], z, z) + return SVector(J * u[1] * vcon[orientation], z, z) end # Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity @@ -105,10 +105,10 @@ end orientation_or_normal_direction, equations::CovariantLinearAdvectionEquation2D) z = zero(eltype(u_ll)) - sqrtG = area_element(aux_vars_ll, equations) + J = area_element(aux_vars_ll, equations) λ = dissipation.max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr, orientation_or_normal_direction, equations) - return -0.5f0 * sqrtG * λ * SVector(u_rr[1] - u_ll[1], z, z) + return -0.5f0 * J * λ * SVector(u_rr[1] - u_ll[1], z, z) end # Maximum contravariant wave speed with respect to specific basis vector From bacf745bbd4134bc9183bdfcb16bf1b068ec21bb Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Fri, 3 Jan 2025 17:17:02 -0500 Subject: [PATCH 19/21] fix typo in quadrilaterals --- src/meshes/p4est_icosahedron_quads.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_icosahedron_quads.jl b/src/meshes/p4est_icosahedron_quads.jl index ed6e6f2..8cb9611 100644 --- a/src/meshes/p4est_icosahedron_quads.jl +++ b/src/meshes/p4est_icosahedron_quads.jl @@ -112,7 +112,7 @@ end # / -------->ξ | └------->ξ \ # /__________________|___________________\ # -# Each of those quadrilaterlas is subdivided into trees_per_face_dimension^2 trees. +# Each of those quadrilaterals is subdivided into trees_per_face_dimension^2 trees. # # We use the following numbering for the 12 vertices of the icosahedron # Fig 3: From 4622df45e91c7f3cc04504ad4c175d4e5da82aad Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Sat, 4 Jan 2025 23:38:45 -0500 Subject: [PATCH 20/21] add vorticity calculation for covariant form --- Project.toml | 1 - ...allowwater_covariant_rossby_haurwitz_EC.jl | 2 +- src/TrixiAtmo.jl | 2 - ...ve_solution_2d_manifold_in_3d_cartesian.jl | 13 ++-- src/callbacks_step/save_solution_covariant.jl | 63 +++++++++++++++---- src/equations/reference_data.jl | 30 +++++---- 6 files changed, 77 insertions(+), 34 deletions(-) diff --git a/Project.toml b/Project.toml index fb2eaf3..c99dfe2 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ version = "0.1.0-DEV" [deps] HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" -Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" diff --git a/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl index cf34d19..3081e39 100644 --- a/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl +++ b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl @@ -55,7 +55,7 @@ analysis_callback = AnalysisCallback(semi, interval = 50, # The SaveSolutionCallback allows to save the solution to a file in regular intervals save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, - solution_variables = cons2cons) + solution_variables = cons2prim_and_vorticity) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl = 0.4) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index a4364b8..76ff950 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -20,8 +20,6 @@ using Reexport: @reexport using LoopVectorization: @turbo using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace -using Infiltrator - @reexport using StaticArrays: SVector, SMatrix include("auxiliary/auxiliary.jl") diff --git a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl index 9807925..e782d82 100644 --- a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl +++ b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl @@ -3,7 +3,8 @@ # Convert to another set of variables using a solution_variables function function convert_variables(u, solution_variables, mesh::P4estMesh{2}, - equations::AbstractEquations{3}, dg, cache) + equations::AbstractEquations{3}, + dg, cache) (; contravariant_vectors) = cache.elements # Extract the number of solution variables to be output # (may be different than the number of conservative variables) @@ -129,7 +130,7 @@ end u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) - # Return th solution variables + # Return the solution variables return SVector(cons2prim(u_node, equations)..., relative_vorticity) end @@ -177,7 +178,9 @@ end end # Variable names for cons2prim_and_vorticity -Trixi.varnames(::typeof(cons2prim_and_vorticity), equations::ShallowWaterEquations3D) = (varnames(cons2prim, - equations)..., - "vorticity") +function Trixi.varnames(::typeof(cons2prim_and_vorticity), + equations::Union{ShallowWaterEquations3D, + AbstractCovariantEquations{2}}) + return (varnames(cons2prim, equations)..., "vorticity") end +end # muladd diff --git a/src/callbacks_step/save_solution_covariant.jl b/src/callbacks_step/save_solution_covariant.jl index b39da3a..2d9753b 100644 --- a/src/callbacks_step/save_solution_covariant.jl +++ b/src/callbacks_step/save_solution_covariant.jl @@ -1,32 +1,71 @@ @muladd begin #! format: noindent -# Convert to another set of variables using a solution_variables function that depends on -# auxiliary variables +# Convert to another set of variables using a solution_variables function function convert_variables(u, solution_variables, mesh::P4estMesh{2}, equations::AbstractCovariantEquations{2}, dg, cache) (; aux_node_vars) = cache.auxiliary_variables - # Extract the number of solution variables to be output # (may be different than the number of conservative variables) - n_vars = length(solution_variables(Trixi.get_node_vars(u, equations, dg, 1, 1, 1), - get_node_aux_vars(aux_node_vars, equations, dg, - 1, 1, - 1), equations)) + n_vars = length(Trixi.varnames(solution_variables, equations)) + # Allocate storage for output data = Array{eltype(u)}(undef, n_vars, nnodes(dg), nnodes(dg), nelements(dg, cache)) # Loop over all nodes and convert variables, passing in auxiliary variables Trixi.@threaded for element in eachelement(dg, cache) for j in eachnode(dg), i in eachnode(dg) - u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) - aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) - u_node_converted = solution_variables(u_node, aux_node, equations) + if applicable(solution_variables, u, equations, dg, cache, i, j, element) + # The solution_variables function depends on the solution on other nodes + data_node = solution_variables(u, equations, dg, cache, i, j, element) + else + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, + element) + data_node = solution_variables(u_node, aux_node, equations) + end + for v in 1:n_vars - data[v, i, j, element] = u_node_converted[v] + data[v, i, j, element] = data_node[v] end end end return data end -end # muladd + +# Calculate the primitive variables and the relative vorticity at a given node +@inline function cons2prim_and_vorticity(u, equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, i, j, element) + (; aux_node_vars) = cache.auxiliary_variables + u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + relative_vorticity = calc_vorticity_node(u, equations, dg, cache, i, j, element) + return SVector(cons2prim(u_node, aux_node, equations)..., relative_vorticity) +end + +# Calculate relative vorticity ζ = ∂₁v₂ - ∂₂v₁ for equations in covariant form +@inline function calc_vorticity_node(u, equations::AbstractCovariantEquations{2}, + dg::DGSEM, cache, i, j, element) + (; derivative_matrix) = dg.basis + (; aux_node_vars) = cache.auxiliary_variables + + dv2dxi1 = dv1dxi2 = zero(eltype(u)) + for ii in eachnode(dg) + u_node_ii = Trixi.get_node_vars(u, equations, dg, ii, j, element) + aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j, element) + vcov = metric_covariant(aux_node_ii, equations) * velocity(u_node_ii, equations) + dv2dxi1 = dv2dxi1 + derivative_matrix[i, ii] * vcov[2] + end + + for jj in eachnode(dg) + u_node_jj = Trixi.get_node_vars(u, equations, dg, i, jj, element) + aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj, element) + vcov = metric_covariant(aux_node_jj, equations) * velocity(u_node_jj, equations) + dv1dxi2 = dv1dxi2 + derivative_matrix[j, jj] * vcov[1] + end + + # compute the relative vorticity + aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element) + return (dv2dxi1 - dv1dxi2) / area_element(aux_node, equations) +end +end # @muladd diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index 1e3fd7a..ad90912 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -76,11 +76,13 @@ $v_\lambda(\theta) = v_0 \cos\theta$, where we define $v_0 = 2\pi a / (12 \ \mat in terms of the Earth's radius $a = 6.37122 \times 10^3\ \mathrm{m}$. The height field then varies with the latitude as ```math -h(\theta) = \frac{1}{g} \Big(gh_0 - \Big(a \Omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), +h(\theta) = \frac{1}{g} +\Big(gh_0 - \Big(a \Omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big), ``` where $gh_0 = 2.94 \times 10^4 \ \mathrm{m}^2/\mathrm{s}^2$, $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and -$\Omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the test suite described in the following paper: +$\Omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the +test suite described in the following paper: - D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. @@ -107,7 +109,8 @@ end @doc raw""" initial_condition_rossby_haurwitz(x, t, equations) -Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and meridional velocity components are given, respectively, by +Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and +meridional velocity components are given, respectively, by ```math \begin{aligned} v_\lambda(\lambda,\theta) &= a \omega \cos \theta+a K \cos ^{R-1} \theta @@ -121,18 +124,18 @@ $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\Omega = 7.292 \times 10^{-5} \ \mathr and $h_0 = 8000 \ \mathrm{m}$ and defining the functions ```math \begin{aligned} -A(\theta) &= \frac{\omega}{2}(2 \Omega+\omega) \cos ^2 \theta + -\frac{1}{4} K^2 \cos ^{2 R} \theta\Big((R+1) \cos ^2 \theta +\left(2 R^2-R-2\right) - +A(\theta) &= \frac{\omega}{2}(2 \Omega+\omega) \cos^2 \theta + +\frac{1}{4} K^2 \cos^{2 R} \theta\Big((R+1) \cos^2\theta +\left(2 R^2-R-2\right) - \big(2 R^2 / \cos^2 \theta\big) \Big), \\ -B(\theta) &= \frac{2(\Omega+\omega) K}{(R+1)(R+2)} \cos ^R \theta\big((R^2+2 R+2)- -(R+1)^2 \cos ^2 \theta\big), \\ -C(\theta) &= \frac{1}{4} K^2 \cos^{2 R} \theta\big((R+1) \cos ^2 \theta-(R+2)\big), +B(\theta) &= \frac{2(\Omega+\omega) K}{(R+1)(R+2)} \cos ^R \theta\big((R^2+2 R+2) - +(R+1)^2 \cos^2 \theta\big), \\ +C(\theta) &= \frac{1}{4} K^2 \cos^{2 R} \theta\big((R+1) \cos^2 \theta-(R+2)\big), \end{aligned} ``` the initial height field is given by ```math -h(\lambda,\theta) = h_0 + \frac{a^2}{g}\Big(A(\theta) + B(\theta)\cos(R\lambda) + -C(\theta)\cos(2R\lambda) \Big). +h(\lambda,\theta) = h_0 + +\frac{a^2}{g}\Big(A(\theta) + B(\theta)\cos(R\lambda) + C(\theta)\cos(2R\lambda) \Big). ``` This problem corresponds to Case 5 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 @@ -147,13 +150,14 @@ This problem corresponds to Case 5 of the test suite described in the following lat = asin(x[3] / a) h_0 = 8.0f3 + omega = 7.848f-6 K = 7.848f-6 R = 4.0f0 - A = 0.5f0 * K * (2 * EARTH_ROTATION_RATE + K) * (cos(lat))^2 + + A = 0.5f0 * K * (2 * EARTH_ROTATION_RATE + omega) * (cos(lat))^2 + 0.25f0 * K^2 * (cos(lat))^(2 * R) * ((R + 1) * (cos(lat))^2 + (2 * R^2 - R - 2) - 2 * R^2 / ((cos(lat))^2)) - B = 2 * (EARTH_ROTATION_RATE + K) * K / ((R + 1) * (R + 2)) * (cos(lat))^R * + B = 2 * (EARTH_ROTATION_RATE + omega) * K / ((R + 1) * (R + 2)) * (cos(lat))^R * ((R^2 + 2R + 2) - (R + 1)^2 * (cos(lat))^2) C = 0.25f0 * K^2 * (cos(lat))^(2 * R) * ((R + 1) * (cos(lat))^2 - (R + 2)) @@ -163,7 +167,7 @@ This problem corresponds to Case 5 of the test suite described in the following (a^2 * A + a^2 * B * cos(R * lon) + a^2 * C * cos(2 * R * lon)) # compute zonal and meridional components of the velocity - vlon = a * K * cos(lat) + + vlon = a * omega * cos(lat) + a * K * (cos(lat))^(R - 1) * (R * (sin(lat))^2 - (cos(lat))^2) * cos(R * lon) vlat = -a * K * R * (cos(lat))^(R - 1) * sin(lat) * sin(R * lon) From 18d9e91d6930861029a04647a09378eac1fa3d21 Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Mon, 6 Jan 2025 09:19:36 -0500 Subject: [PATCH 21/21] improve docs and comments --- ...llowwater_covariant_geostrophic_balance.jl | 31 ++++++++++++------- ...allowwater_covariant_rossby_haurwitz_EC.jl | 27 ++++++++++------ ...allowwater_cubed_sphere_shell_advection.jl | 1 + ...erical_advection_covariant_cubed_sphere.jl | 1 + ...al_advection_covariant_quad_icosahedron.jl | 20 +++++++----- src/callbacks_step/analysis_covariant.jl | 2 +- ...ve_solution_2d_manifold_in_3d_cartesian.jl | 2 +- src/callbacks_step/stepsize_dg2d.jl | 2 +- src/equations/reference_data.jl | 7 +++-- .../containers_2d_manifold_in_3d_covariant.jl | 2 +- .../dg_2d_manifold_in_3d_cartesian.jl | 2 +- test/test_2d_shallow_water_covariant.jl | 5 +-- 12 files changed, 64 insertions(+), 38 deletions(-) diff --git a/examples/elixir_shallowwater_covariant_geostrophic_balance.jl b/examples/elixir_shallowwater_covariant_geostrophic_balance.jl index 0eadcb6..a1040d4 100644 --- a/examples/elixir_shallowwater_covariant_geostrophic_balance.jl +++ b/examples/elixir_shallowwater_covariant_geostrophic_balance.jl @@ -1,5 +1,5 @@ ############################################################################### -# DGSEM for the shallow water equations on the cubed sphere +# DGSEM for the shallow water equations in covariant form on the cubed sphere ############################################################################### using OrdinaryDiffEq, Trixi, TrixiAtmo @@ -11,7 +11,7 @@ initial_condition = initial_condition_geostrophic_balance polydeg = 3 cells_per_dimension = 5 n_saves = 10 -tspan = (0.0, 1.0 * SECONDS_PER_DAY) +tspan = (0.0, 5.0 * SECONDS_PER_DAY) ############################################################################### # Spatial discretization @@ -24,11 +24,17 @@ equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, global_coordinate_system = GlobalSphericalCoordinates()) +# The covariant shallow water equations are treated as a nonconservative system in order to +# handle flux-differencing formulations of the covariant derivative. With +# VolumeIntegralWeakForm, there are actually no nonconservative terms, but we must still +# pass a no-op function "flux_nonconservative_weak_form" as the nonconservative surface flux +surface_flux = (flux_lax_friedrichs, flux_nonconservative_weak_form) + # Create DG solver with polynomial degree = p -solver = DGSEM(polydeg = polydeg, - surface_flux = (flux_lax_friedrichs, flux_nonconservative_weak_form), - volume_integral = VolumeIntegralWeakForm()) +solver = DGSEM(polydeg = polydeg, volume_integral = VolumeIntegralWeakForm(), + surface_flux = surface_flux) +# Transform the initial condition to the proper set of conservative variables initial_condition_transformed = transform_initial_condition(initial_condition, equations) # A semidiscretization collects data structures and functions for the spatial discretization @@ -41,12 +47,13 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transform # Create ODE problem with time span from 0 to T ode = semidiscretize(semi, tspan) -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers +# 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 = 50, +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the +# results +analysis_callback = AnalysisCallback(semi, interval = 200, save_analysis = true, extra_analysis_errors = (:conservation_error,)) @@ -57,14 +64,16 @@ save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl = 0.4) -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE +# solver callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) ############################################################################### # run the simulation -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed +# callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 100.0, save_everystep = false, callback = callbacks); diff --git a/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl index 3081e39..c9bf273 100644 --- a/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl +++ b/examples/elixir_shallowwater_covariant_rossby_haurwitz_EC.jl @@ -1,5 +1,6 @@ ############################################################################### -# DGSEM for the shallow water equations on the cubed sphere +# Entropy-conservative DGSEM for the shallow water equations in covariant form +# on the cubed sphere ############################################################################### using OrdinaryDiffEq, Trixi, TrixiAtmo @@ -11,7 +12,7 @@ initial_condition = initial_condition_rossby_haurwitz polydeg = 3 cells_per_dimension = 5 n_saves = 10 -tspan = (0.0, 1.0 * SECONDS_PER_DAY) +tspan = (0.0, 7.0 * SECONDS_PER_DAY) ############################################################################### # Spatial discretization @@ -24,13 +25,18 @@ equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, global_coordinate_system = GlobalCartesianCoordinates()) -# Create DG solver with polynomial degree = p +# Use entropy-conservative two-point fluxes for volume and surface terms volume_flux = (flux_split_covariant, flux_nonconservative_split_covariant) surface_flux = (flux_split_covariant, flux_nonconservative_split_covariant) +# The following flux should be used to add interface dissipation: +# surface_flux = (flux_split_covariant_lax_friedrichs, flux_nonconservative_split_covariant) + +# Create DG solver with polynomial degree = p solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) +# Transform the initial condition to the proper set of conservative variables initial_condition_transformed = transform_initial_condition(initial_condition, equations) # A semidiscretization collects data structures and functions for the spatial discretization @@ -43,12 +49,13 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transform # Create ODE problem with time span from 0 to T ode = semidiscretize(semi, tspan) -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers +# 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 = 50, +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the +# results. Note that entropy should be conserved at the semi-discrete level. +analysis_callback = AnalysisCallback(semi, interval = 200, save_analysis = true, extra_analysis_errors = (:conservation_error,), extra_analysis_integrals = (entropy,)) @@ -60,14 +67,16 @@ save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves, # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl = 0.4) -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE +# solver callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) ############################################################################### # run the simulation -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed +# callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 100.0, save_everystep = false, callback = callbacks); diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl index 667297f..727add3 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl @@ -46,6 +46,7 @@ mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = 3, #initial_refinement_level = 0, # Comment to use cells_per_dimension in the convergence test element_local_mapping = false) +# Transform the initial condition to the proper set of conservative variables initial_condition_transformed = transform_initial_condition(initial_condition, equations) # A semidiscretization collects data structures and functions for the spatial discretization diff --git a/examples/elixir_spherical_advection_covariant_cubed_sphere.jl b/examples/elixir_spherical_advection_covariant_cubed_sphere.jl index 8193213..1cc513e 100644 --- a/examples/elixir_spherical_advection_covariant_cubed_sphere.jl +++ b/examples/elixir_spherical_advection_covariant_cubed_sphere.jl @@ -25,6 +25,7 @@ mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = Trixi.polydeg(solver), element_local_mapping = true) +# Transform the initial condition to the proper set of conservative variables initial_condition_transformed = transform_initial_condition(initial_condition, equations) # A semidiscretization collects data structures and functions for the spatial discretization diff --git a/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl b/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl index 7699976..ddab0ab 100644 --- a/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl +++ b/examples/elixir_spherical_advection_covariant_quad_icosahedron.jl @@ -1,5 +1,5 @@ ############################################################################### -# DGSEM for the linear advection equation on the cubed sphere +# DGSEM for the linear advection equation on a quadrilateral icosahedral grid ############################################################################### # To run a convergence test, use # convergence_test("../examples/elixir_spherical_advection_covariant_quad_icosahedron.jl", 4, cells_per_dimension = (1,1)) @@ -18,12 +18,13 @@ equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = Global solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, volume_integral = VolumeIntegralWeakForm()) -# Create a 2D cubed sphere mesh the size of the Earth. For the covariant form to work -# properly, we currently need polydeg to equal that of the solver, and +# Create a 2D quadrilateral icosahedral mesh the size of the Earth. For the covariant form +# to work properly, we currently need polydeg to equal that of the solver, and # initial_refinement_level = 0 (default) mesh = P4estMeshQuadIcosahedron2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = Trixi.polydeg(solver)) +# Transform the initial condition to the proper set of conservative variables initial_condition_transformed = transform_initial_condition(initial_condition, equations) # A semidiscretization collects data structures and functions for the spatial discretization @@ -35,11 +36,12 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transform # Create ODE problem with time span from 0 to T ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY)) -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation +# setup and resets the timers summary_callback = SummaryCallback() -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the +# results analysis_callback = AnalysisCallback(semi, interval = 10, save_analysis = true, extra_analysis_errors = (:conservation_error,)) @@ -51,14 +53,16 @@ save_solution = SaveSolutionCallback(interval = 10, # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl = 0.7) -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE +# solver callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) ############################################################################### # run the simulation -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed +# callbacks sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0, save_everystep = false, callback = callbacks); diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index c823594..9b0c689 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/src/callbacks_step/analysis_covariant.jl @@ -124,4 +124,4 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2}, return l2_error, linf_error end -end # muladd +end # @muladd diff --git a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl index e782d82..fabf353 100644 --- a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl +++ b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl @@ -183,4 +183,4 @@ function Trixi.varnames(::typeof(cons2prim_and_vorticity), AbstractCovariantEquations{2}}) return (varnames(cons2prim, equations)..., "vorticity") end -end # muladd +end # @muladd diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 45530be..0d9c888 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -69,4 +69,4 @@ function Trixi.max_dt(u, t, mesh::P4estMesh{2}, constant_speed::False, end return 2 / (nnodes(dg) * max_scaled_speed) end -end # muladd +end # @muladd diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index ad90912..8c27a80 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -81,7 +81,7 @@ h(\theta) = \frac{1}{g} ``` where $gh_0 = 2.94 \times 10^4 \ \mathrm{m}^2/\mathrm{s}^2$, $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and -$\Omega = 7.292 \times 10^{-5} \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the +$\Omega = 7.292 \times 10^{-5}\ \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the test suite described in the following paper: - D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in @@ -110,7 +110,8 @@ end initial_condition_rossby_haurwitz(x, t, equations) Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and -meridional velocity components are given, respectively, by +meridional velocity components are given, respectively, as functions of the longitude +$\lambda$ and latitude $\theta$ by ```math \begin{aligned} v_\lambda(\lambda,\theta) &= a \omega \cos \theta+a K \cos ^{R-1} \theta @@ -176,4 +177,4 @@ This problem corresponds to Case 5 of the test suite described in the following # coordinate system, which depends on the equation type return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations) end -end # muladd +end # @muladd diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl index 9580f1a..59ce7d9 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl @@ -350,4 +350,4 @@ function calc_christoffel_symbols!(aux_node_vars, mesh::P4estMesh{2, 3}, Gcon[2, 2] * christoffel_firstkind_2[2, 2] end end -end # muladd +end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl index 1b494ef..c7028e5 100644 --- a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl @@ -152,4 +152,4 @@ function calc_sources_2d_manifold_in_3d!(du, u, t, source_terms, return nothing end -end # muladd +end # @muladd diff --git a/test/test_2d_shallow_water_covariant.jl b/test/test_2d_shallow_water_covariant.jl index f9c4e6d..7735de4 100644 --- a/test/test_2d_shallow_water_covariant.jl +++ b/test/test_2d_shallow_water_covariant.jl @@ -19,7 +19,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") 1.4786544678290738, 0.0013754600033514114, 0.0007564014737144256 - ]) + ], tspan=(0.0, 1.0 * SECONDS_PER_DAY)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -41,7 +41,8 @@ end linf=[ 579.0744773822007, 0.5054767269383867, - 0.5628603103981209]) + 0.5628603103981209 + ], tspan=(0.0, 1.0 * SECONDS_PER_DAY)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let