Skip to content

Commit

Permalink
make integrate_via_indices use exact Jacobian
Browse files Browse the repository at this point in the history
  • Loading branch information
tristanmontoya committed Dec 1, 2024
1 parent 3635746 commit b86ed28
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 3 deletions.
40 changes: 38 additions & 2 deletions src/callbacks_step/analysis_covariant.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,43 @@
@muladd begin
#! format: noindent

# For the covariant form, we want to integrate using the exact area element
# sqrtG = sqrt(det(G)), which is stored in cache.auxiliary_variables, not the approximate
# area element used in the Cartesian formulation, which stored in cache.elements.
function Trixi.integrate_via_indices(func::Func, u,
mesh::Union{StructuredMesh{2},
StructuredMeshView{2},
UnstructuredMesh2D, P4estMesh{2},
T8codeMesh{2}},
equations::AbstractCovariantEquations{2},
dg::DGSEM, cache, args...;
normalize = true) where {Func}
(; weights) = dg.basis
(; aux_node_vars) = cache.auxiliary_variables

# Initialize integral with zeros of the right shape
integral = zero(func(u, 1, 1, 1, equations, dg, args...))
total_volume = zero(real(mesh))

# Use quadrature to numerically integrate over entire domain
for element in eachelement(dg, cache)
for j in eachnode(dg), i in eachnode(dg)
aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
sqrtG = area_element(aux_node, equations)
integral += weights[i] * weights[j] * sqrtG *
func(u, i, j, element, equations, dg, args...)
total_volume += weights[i] * weights[j] * sqrtG
end
end

# Normalize with total volume
if normalize
integral = integral / total_volume
end

return integral
end

# Entropy time derivative which uses auxiliary variables
function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t,
mesh::P4estMesh{2},
Expand Down Expand Up @@ -28,7 +65,6 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
(; weights) = dg.basis
(; node_coordinates) = cache.elements
(; aux_node_vars) = cache.auxiliary_variables

# Set up data structures
l2_error = zero(func(Trixi.get_node_vars(u, equations, dg, 1, 1, 1), equations))
linf_error = copy(l2_error)
Expand All @@ -51,7 +87,7 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
u_numerical = Trixi.get_node_vars(u, equations, dg, i, j, element)
diff = func(u_exact, equations) - func(u_numerical, equations)

# For the L2 error, integrate with respect to volume element stored in aux vars
# For the L2 error, integrate with respect to area element stored in aux vars
sqrtG = area_element(aux_node, equations)
l2_error += diff .^ 2 * (weights[i] * weights[j] * sqrtG)

Expand Down
2 changes: 1 addition & 1 deletion src/equations/reference_data.jl
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ the test suite described in the following paper:
@inline function initial_condition_gaussian(x, t, ::AbstractEquations)
RealT = eltype(x)

a = EARTH_RADIUS # radius of the sphere in metres
a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) # radius of the sphere
V = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY) # speed of rotation
alpha = convert(RealT, π / 4) # angle of rotation
h_0 = 1000.0f0 # bump height in metres
Expand Down

0 comments on commit b86ed28

Please sign in to comment.