Skip to content

Commit 8121b38

Browse files
tristanmontoyatristanmontoya
committed
use ForwardDiff to compute Christoffel symbols exactly
1 parent a496119 commit 8121b38

File tree

4 files changed

+62
-75
lines changed

4 files changed

+62
-75
lines changed

‎Project.toml

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ authors = ["Benedict Geihe <bgeihe@uni-koeln.de>", "Tristan Montoya <montoya.tri
44
version = "0.1.0-DEV"
55

66
[deps]
7+
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
78
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
89
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
910
LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
@@ -18,6 +19,7 @@ StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
1819
Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
1920

2021
[compat]
22+
ForwardDiff = "0.10.36"
2123
HDF5 = "0.16.10, 0.17"
2224
LinearAlgebra = "1"
2325
LoopVectorization = "0.12.118"

‎src/TrixiAtmo.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ using StaticArrayInterface: static_size
1818
using LinearAlgebra: cross, norm, dot, det
1919
using Reexport: @reexport
2020
using LoopVectorization: @turbo
21+
using ForwardDiff: derivative
2122
using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
2223

2324
@reexport using StaticArrays: SVector, SMatrix

‎src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_covariant.jl

Lines changed: 49 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -213,10 +213,16 @@ function init_auxiliary_node_variables!(auxiliary_variables, mesh::P4estMesh{2,
213213
aux_node_vars[17:19, i, j, element] = SVector(metric_contravariant[1, 1],
214214
metric_contravariant[1, 2],
215215
metric_contravariant[2, 2])
216-
end
217216

218-
# Christoffel symbols of the second kind (aux_node_vars[20:25, :, :, element])
219-
calc_christoffel_symbols!(aux_node_vars, mesh, equations, dg, element)
217+
# Christoffel symbols of the second kind
218+
aux_node_vars[20:25, i, j, element] = calc_christoffel_symbols(v1, v2, v3,
219+
v4,
220+
dg.basis.nodes[i],
221+
dg.basis.nodes[j],
222+
radius,
223+
metric_contravariant,
224+
equations)
225+
end
220226
end
221227

222228
return nothing
@@ -290,68 +296,46 @@ end
290296
return SMatrix{3, 2}(A[1, 1], A[2, 1], 0.0f0, A[1, 2], A[2, 2], 0.0f0)
291297
end
292298

293-
# Calculate Christoffel symbols approximately using the collocation derivative. Note that
294-
# they could alternatively be computed exactly without affecting the entropy stability
295-
# properties of the scheme.
296-
function calc_christoffel_symbols!(aux_node_vars, mesh::P4estMesh{2, 3},
297-
equations::AbstractCovariantEquations{2, 3}, dg,
298-
element)
299-
(; derivative_matrix) = dg.basis
300-
301-
for j in eachnode(dg), i in eachnode(dg)
302-
303-
# Numerically differentiate covariant metric components with respect to ξ¹
304-
dG11dxi1 = zero(eltype(aux_node_vars))
305-
dG12dxi1 = zero(eltype(aux_node_vars))
306-
dG22dxi1 = zero(eltype(aux_node_vars))
307-
for ii in eachnode(dg)
308-
aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j,
309-
element)
310-
Gcov_ii = metric_covariant(aux_node_ii, equations)
311-
dG11dxi1 = dG11dxi1 + derivative_matrix[i, ii] * Gcov_ii[1, 1]
312-
dG12dxi1 = dG12dxi1 + derivative_matrix[i, ii] * Gcov_ii[1, 2]
313-
dG22dxi1 = dG22dxi1 + derivative_matrix[i, ii] * Gcov_ii[2, 2]
314-
end
315-
316-
# Numerically differentiate covariant metric components with respect to ξ²
317-
dG11dxi2 = zero(eltype(aux_node_vars))
318-
dG12dxi2 = zero(eltype(aux_node_vars))
319-
dG22dxi2 = zero(eltype(aux_node_vars))
320-
for jj in eachnode(dg)
321-
aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj,
322-
element)
323-
Gcov_jj = metric_covariant(aux_node_jj, equations)
324-
dG11dxi2 = dG11dxi2 + derivative_matrix[j, jj] * Gcov_jj[1, 1]
325-
dG12dxi2 = dG12dxi2 + derivative_matrix[j, jj] * Gcov_jj[1, 2]
326-
dG22dxi2 = dG22dxi2 + derivative_matrix[j, jj] * Gcov_jj[2, 2]
327-
end
299+
# Calculate the covariant metric tensor components G₁₁, G₁₂ (= G₂₁), and G₂₂ and return in
300+
# that order as an SVector of length 3
301+
function calc_metric_covariant(v1, v2, v3, v4, xi1, xi2, radius, equations)
302+
A = calc_basis_covariant(v1, v2, v3, v4, xi1, xi2, radius,
303+
equations.global_coordinate_system)
304+
Gcov = A' * A
305+
return SVector(Gcov[1, 1], Gcov[1, 2], Gcov[2, 2])
306+
end
328307

329-
# Compute Christoffel symbols of the first kind
330-
christoffel_firstkind_1 = SMatrix{2, 2}(0.5f0 * dG11dxi1,
331-
0.5f0 * dG11dxi2,
332-
0.5f0 * dG11dxi2,
333-
dG12dxi2 - 0.5f0 * dG22dxi1)
334-
christoffel_firstkind_2 = SMatrix{2, 2}(dG12dxi1 - 0.5f0 * dG11dxi2,
335-
0.5f0 * dG22dxi1,
336-
0.5f0 * dG22dxi1,
337-
0.5f0 * dG22dxi2)
338-
339-
# Raise indices to get Christoffel symbols of the second kind
340-
aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
341-
Gcon = metric_contravariant(aux_node, equations)
342-
aux_node_vars[20, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 1] +
343-
Gcon[1, 2] * christoffel_firstkind_2[1, 1]
344-
aux_node_vars[21, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[1, 2] +
345-
Gcon[1, 2] * christoffel_firstkind_2[1, 2]
346-
aux_node_vars[22, i, j, element] = Gcon[1, 1] * christoffel_firstkind_1[2, 2] +
347-
Gcon[1, 2] * christoffel_firstkind_2[2, 2]
348-
349-
aux_node_vars[23, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[1, 1] +
350-
Gcon[2, 2] * christoffel_firstkind_2[1, 1]
351-
aux_node_vars[24, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[1, 2] +
352-
Gcon[2, 2] * christoffel_firstkind_2[1, 2]
353-
aux_node_vars[25, i, j, element] = Gcon[2, 1] * christoffel_firstkind_1[2, 2] +
354-
Gcon[2, 2] * christoffel_firstkind_2[2, 2]
355-
end
308+
# Calculate the Christoffel symbols of the second kind using ForwardDiff to automatically
309+
# differentiate the metric, and return the components Γ¹₁₁, Γ¹₁₂ (= Γ¹₂₁), Γ¹₂₂, Γ²₁₁,
310+
# Γ²₁₂ (= Γ²₂₁), and Γ²₂₂ as an SVector of length 6
311+
function calc_christoffel_symbols(v1, v2, v3, v4, xi1, xi2, radius, Gcon, equations)
312+
313+
# Use ForwardDiff to differentiate the covariant metric tensor components
314+
dGdxi1 = derivative(x -> calc_metric_covariant(v1, v2, v3, v4, x, xi2, radius,
315+
equations), xi1)
316+
dGdxi2 = derivative(x -> calc_metric_covariant(v1, v2, v3, v4, xi1, x, radius,
317+
equations), xi2)
318+
319+
# Extract SVector components from derivatives of vectors
320+
dG11dxi1, dG12dxi1, dG22dxi1 = dGdxi1
321+
dG11dxi2, dG12dxi2, dG22dxi2 = dGdxi2
322+
323+
# Compute Christoffel symbols of the first kind
324+
Gamma_1 = SMatrix{2, 2}(0.5f0 * dG11dxi1, # Γ₁₁₁
325+
0.5f0 * dG11dxi2, # Γ₁₂₁
326+
0.5f0 * dG11dxi2, # Γ₁₁₂
327+
dG12dxi2 - 0.5f0 * dG22dxi1) # Γ₁₂₂
328+
Gamma_2 = SMatrix{2, 2}(dG12dxi1 - 0.5f0 * dG11dxi2, # Γ₂₁₁
329+
0.5f0 * dG22dxi1, # Γ₂₂₁
330+
0.5f0 * dG22dxi1, # Γ₂₁₂
331+
0.5f0 * dG22dxi2) # Γ₂₂₂
332+
333+
# Raise indices to get Christoffel symbols of the second kind
334+
return SVector(Gcon[1, 1] * Gamma_1[1, 1] + Gcon[1, 2] * Gamma_2[1, 1], # Γ¹₁₁
335+
Gcon[1, 1] * Gamma_1[1, 2] + Gcon[1, 2] * Gamma_2[1, 2], # Γ¹₁₂ (= Γ¹₂₁)
336+
Gcon[1, 1] * Gamma_1[2, 2] + Gcon[1, 2] * Gamma_2[2, 2], # Γ¹₂₂
337+
Gcon[2, 1] * Gamma_1[1, 1] + Gcon[2, 2] * Gamma_2[1, 1], # Γ²₁₁
338+
Gcon[2, 1] * Gamma_1[1, 2] + Gcon[2, 2] * Gamma_2[1, 2], # Γ²₁₂ (= Γ²₂₁)
339+
Gcon[2, 1] * Gamma_1[2, 2] + Gcon[2, 2] * Gamma_2[2, 2]) # Γ²₂₂
356340
end
357341
end # @muladd

‎test/test_2d_shallow_water_covariant.jl

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -11,14 +11,14 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples")
1111
@test_trixi_include(joinpath(EXAMPLES_DIR,
1212
"elixir_shallowwater_covariant_geostrophic_balance.jl"),
1313
l2=[
14-
0.30653144639858293,
15-
0.00019984467582353295,
16-
8.767819502806826e-5
14+
0.27823189465180176,
15+
0.00021164119121947655,
16+
8.588414721470682e-5
1717
],
1818
linf=[
19-
1.4786544678290738,
20-
0.001375460003351342,
21-
0.0007564014737142868
19+
1.818810315677183,
20+
0.001140237382440748,
21+
0.0010284231183847609
2222
],
2323
tspan=(0.0, 1.0 * SECONDS_PER_DAY))
2424
# Ensure that we do not have excessive memory allocations
@@ -34,8 +34,8 @@ end
3434
@trixiatmo_testset "elixir_shallowwater_covariant_rossby_haurwitz_EC" begin
3535
@test_trixi_include(joinpath(EXAMPLES_DIR,
3636
"elixir_shallowwater_covariant_rossby_haurwitz_EC.jl"),
37-
l2=[265.9858131601718, 0.1769083525032968, 0.2538332624036849],
38-
linf=[579.0744773821989, 0.5054767269383715, 0.5628603103981238],
37+
l2=[265.982080624612, 0.17691380660353082, 0.25383558815062035],
38+
linf=[579.0665426686955, 0.5054650879522844, 0.5629724540607254],
3939
tspan=(0.0, 1.0 * SECONDS_PER_DAY))
4040
# Ensure that we do not have excessive memory allocations
4141
# (e.g., from type instabilities)
@@ -50,8 +50,8 @@ end
5050
@trixiatmo_testset "elixir_shallowwater_covariant_rossby_haurwitz_EC (LLF surface flux)" begin
5151
@test_trixi_include(joinpath(EXAMPLES_DIR,
5252
"elixir_shallowwater_covariant_rossby_haurwitz_EC.jl"),
53-
l2=[265.9818409531758, 0.17644362975403768, 0.2535621643485035],
54-
linf=[574.6722628324133, 0.5155374806433987, 0.5497046315783312],
53+
l2=[265.9809675383955, 0.17644521491170614, 0.2535629213762403],
54+
linf=[574.6705006776556, 0.5155467030841373, 0.549735873681727],
5555
tspan=(0.0, 1.0 * SECONDS_PER_DAY),
5656
surface_flux=(flux_lax_friedrichs, flux_nonconservative_ec))
5757
# Ensure that we do not have excessive memory allocations

0 commit comments

Comments
 (0)