Skip to content

Commit c4394bb

Browse files
Add option to use ForwardDiff.jl to compute Christoffel symbols exactly and separate covariant form containers (#78)
* use ForwardDiff to compute Christoffel symbols exactly * update tests * formatting * restructure containers and remove unused fields from covariant solver * fix elixir and give file with metric terms descriptive name * ability to choose how to compute christoffel symbols * Apply suggestions from code review Co-authored-by: Andrés Rueda-Ramírez <[email protected]> * remove surface_integral from prolong2interfaces! call * fix surface integral issue, add verbose names for metric derivatives --------- Co-authored-by: Andrés Rueda-Ramírez <[email protected]>
1 parent 561a5d1 commit c4394bb

11 files changed

+452
-119
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 <[email protected]>", "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"
@@ -19,6 +20,7 @@ StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
1920
Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
2021

2122
[compat]
23+
ForwardDiff = "0.10.36"
2224
HDF5 = "0.16.10, 0.17"
2325
LinearAlgebra = "1"
2426
LoopVectorization = "0.12.118"

examples/elixir_shallowwater_covariant_rossby_haurwitz.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -37,8 +37,11 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
3737
# Transform the initial condition to the proper set of conservative variables
3838
initial_condition_transformed = transform_initial_condition(initial_condition, equations)
3939

40-
# A semidiscretization collects data structures and functions for the spatial discretization
40+
# A semidiscretization collects data structures and functions for the spatial discretization.
41+
# Even though `metric_terms = MetricTermsCovariantSphere()` is default, we pass it here
42+
# explicitly, such that `metric_terms` can be adjusted from the `trixi_include()` call in the tests
4143
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
44+
metric_terms = MetricTermsCovariantSphere(),
4245
source_terms = source_terms_geometric_coriolis)
4346

4447
###############################################################################

src/TrixiAtmo.jl

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ using LinearAlgebra: cross, norm, dot, det
1919
using Reexport: @reexport
2020
using LoopVectorization: @turbo
2121
using QuadGK: quadgk
22+
using ForwardDiff: derivative
2223
using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
2324

2425
@reexport using StaticArrays: SVector, SMatrix
@@ -29,7 +30,6 @@ include("equations/equations.jl")
2930
include("meshes/meshes.jl")
3031
include("semidiscretization/semidiscretization.jl")
3132
include("solvers/solvers.jl")
32-
include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
3333
include("callbacks_step/callbacks_step.jl")
3434

3535
export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D,
@@ -49,7 +49,8 @@ export source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier!
4949
export cons2prim_and_vorticity, contravariant2global
5050

5151
export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct,
52-
MetricTermsInvariantCurl
52+
MetricTermsInvariantCurl, MetricTermsCovariantSphere, ChristoffelSymbolsAutodiff,
53+
ChristoffelSymbolsCollocationDerivative
5354

5455
export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
5556
EARTH_ROTATION_RATE, SECONDS_PER_DAY
Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
"""
2+
MetricTermsCrossProduct()
3+
4+
Struct used for multiple dispatch on functions that compute the metric terms.
5+
When the argument `metric_terms` is of type `MetricTermsCrossProduct`, the
6+
contravariant vectors are computed using the cross-product form.
7+
"""
8+
struct MetricTermsCrossProduct end
9+
10+
"""
11+
MetricTermsInvariantCurl()
12+
13+
Struct used for multiple dispatch on functions that compute the metric terms.
14+
When the argument `metric_terms` is of type `MetricTermsInvariantCurl`, the
15+
contravariant vectors are computed using the invariant curl form.
16+
## References
17+
- Kopriva, D. A. (2006). Metric identities and the discontinuous spectral element method on
18+
curvilinear meshes. Journal of Scientific Computing 26, 301-327.
19+
[DOI: 10.1007/s10915-005-9070-8](https://doi.org/10.1007/s10915-005-9070-8)
20+
- Vinokur, M. and Yee, H. C. (2001). Extension of efficient low dissipation high order
21+
schemes for 3-D curvilinear moving grids. In Caughey, D. A., and Hafez, M. M. (eds.),
22+
Frontiers of Computational Fluid Dynamics 2002, World Scientific, Singapore, pp. 129–164.
23+
[DOI: 10.1142/9789812810793_0008](https://doi.org/10.1142/9789812810793_0008)
24+
"""
25+
struct MetricTermsInvariantCurl end
26+
27+
"""
28+
MetricTermsCovariantSphere(christoffel_symbols = ChristoffelSymbolsAutodiff())
29+
30+
Struct specifying options for computing geometric information for discretizations in
31+
covariant form based on an exact representation of the spherical geometry. Currently, the
32+
only field is `christoffel_symbols`, specifying the approach used to compute the
33+
Christoffel symbols, for which the options are [`ChristoffelSymbolsAutodiff`](@ref) or
34+
[`ChristoffelSymbolsCollocationDerivative`](@ref).
35+
"""
36+
struct MetricTermsCovariantSphere{ChristoffelSymbols}
37+
christoffel_symbols::ChristoffelSymbols
38+
function MetricTermsCovariantSphere(;
39+
christoffel_symbols::ChristoffelSymbols = ChristoffelSymbolsAutodiff()) where {ChristoffelSymbols}
40+
return new{ChristoffelSymbols}(christoffel_symbols)
41+
end
42+
end
43+
44+
@doc raw"""
45+
ChristoffelSymbolsAutodiff()
46+
47+
Struct used for multiple dispatch on functions that compute the Christoffel symbols. This
48+
option uses [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) to compute
49+
```math
50+
\Gamma_{jk}^i =
51+
\frac{1}{2}G^{il}\big(\partial_j G_{kl} + \partial_k G_{jl} - \partial_l G_{jk}\big)
52+
```
53+
using forward-mode automatic differentiation.
54+
!!! warning
55+
Using this option with [`GlobalSphericalCoordinates`](@ref) is prone to `NaN` values as
56+
a result of the polar singularity. This is remedied through the use of
57+
[`GlobalCartesianCoordinates`](@ref).
58+
"""
59+
struct ChristoffelSymbolsAutodiff end
60+
61+
@doc raw"""
62+
ChristoffelSymbolsCollocationDerivative()
63+
64+
Struct used for multiple dispatch on functions that compute the Christoffel symbols.
65+
Letting $I^N$ denote the degree $N$ polynomial interpolation operator collocated with the
66+
scheme's quadrature nodes, this option computes the Christoffel symbols at each quadrature
67+
node using the approximation
68+
```math
69+
\Gamma_{jk}^i \approx
70+
\frac{1}{2}G^{il}\big(\partial_j I^N G_{kl} + \partial_k I^N G_{jl} - \partial_l I^N G_{jk}\big).
71+
```
72+
"""
73+
struct ChristoffelSymbolsCollocationDerivative end
Lines changed: 2 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,2 @@
1-
"""
2-
MetricTermsCrossProduct()
3-
4-
Struct used for multiple dispatch on functions that compute the metric terms.
5-
When the argument `metric_terms` is of type `MetricTermsCrossProduct`, the
6-
contravariant vectors are computed using the cross-product form.
7-
"""
8-
struct MetricTermsCrossProduct end
9-
10-
"""
11-
MetricTermsInvariantCurl()
12-
13-
Struct used for multiple dispatch on functions that compute the metric terms.
14-
When the argument `metric_terms` is of type `MetricTermsInvariantCurl`, the
15-
contravariant vectors are computed using the invariant curl form.
16-
## References
17-
- Kopriva, D. A. (2006). Metric identities and the discontinuous spectral element method on
18-
curvilinear meshes. Journal of Scientific Computing 26, 301-327.
19-
[DOI: 10.1007/s10915-005-9070-8](https://doi.org/10.1007/s10915-005-9070-8)
20-
- Vinokur, M. and Yee, H. C. (2001). Extension of efficient low dissipation high order
21-
schemes for 3-D curvilinear moving grids. In Caughey, D. A., and Hafez, M. M. (eds.),
22-
Frontiers of Computational Fluid Dynamics 2002, World Scientific, Singapore, pp. 129–164.
23-
[DOI: 10.1142/9789812810793_0008](https://doi.org/10.1142/9789812810793_0008)
24-
"""
25-
struct MetricTermsInvariantCurl end
1+
include("metric_terms.jl")
2+
include("semidiscretization_hyperbolic_2d_manifold_in_3d.jl")

src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl

Lines changed: 36 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,14 @@
11
# The SemidiscretizationHyperbolic constructor has been modified to remove the assertion
22
# that checks the compatibility between the mesh dimensionality and the equations'
33
# dimensionality. Instead, we now directly dispatch using the specific mesh type
4-
# (P4estMesh{2}) for 2D meshes and AbstractEquations{3} for 3D equations or
5-
# AbstractCovariantEquations{2, 3} for 2D covariant equations in three-dimensional space.
4+
# (P4estMesh{2}) for 2D meshes and AbstractEquations{3} for 3D equations.
65
# This change is necessary to support the implementation of a 2D manifold embedded in a 3D
76
# space. We also pass in the keyword arguments "metric_terms" and "auxiliary_field", which
87
# specify the approach used to compute the metric terms as well as any auxiliary fields
98
# (e.g. bottom topography or a background state) that may be needed by the solver but are
109
# not stored in the solution variables or elsewhere in the cache.
1110
function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{2},
12-
equations::Union{AbstractEquations{3},
13-
AbstractCovariantEquations{2,
14-
3}},
11+
equations::AbstractEquations{3},
1512
initial_condition,
1613
solver;
1714
source_terms = nothing,
@@ -38,3 +35,37 @@ function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{2},
3835
source_terms, solver,
3936
cache)
4037
end
38+
39+
# Constructor for SemidiscretizationHyperbolic for the covariant form. Requires
40+
# compatibility between the mesh and equations (i.e. the same `NDIMS` and `NDIMS_AMBIENT`)
41+
# and sets the default metric terms to MetricTermsCovariantSphere.
42+
function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{NDIMS, NDIMS_AMBIENT},
43+
equations::AbstractCovariantEquations{NDIMS,
44+
NDIMS_AMBIENT},
45+
initial_condition,
46+
solver;
47+
source_terms = nothing,
48+
boundary_conditions = boundary_condition_periodic,
49+
# `RealT` is used as real type for node locations etc.
50+
# while `uEltype` is used as element type of solutions etc.
51+
RealT = real(solver), uEltype = RealT,
52+
initial_cache = NamedTuple(),
53+
metric_terms = MetricTermsCovariantSphere(),
54+
auxiliary_field = nothing) where {NDIMS,
55+
NDIMS_AMBIENT}
56+
cache = (;
57+
Trixi.create_cache(mesh, equations, solver, RealT, metric_terms,
58+
auxiliary_field, uEltype)..., initial_cache...)
59+
_boundary_conditions = Trixi.digest_boundary_conditions(boundary_conditions, mesh,
60+
solver,
61+
cache)
62+
63+
SemidiscretizationHyperbolic{typeof(mesh), typeof(equations),
64+
typeof(initial_condition),
65+
typeof(_boundary_conditions), typeof(source_terms),
66+
typeof(solver), typeof(cache)}(mesh, equations,
67+
initial_condition,
68+
_boundary_conditions,
69+
source_terms, solver,
70+
cache)
71+
end

src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -65,8 +65,7 @@ end
6565
# This function dispatches on the dimensions of the mesh and the equation (AbstractEquations{3})
6666
function Trixi.init_elements(mesh::Union{P4estMesh{2, 3, RealT},
6767
T8codeMesh{2}},
68-
equations::Union{AbstractEquations{3},
69-
AbstractCovariantEquations{2, 3}},
68+
equations::AbstractEquations{3},
7069
basis,
7170
metric_terms,
7271
::Type{uEltype}) where {RealT <: Real, uEltype <: Real}
@@ -140,7 +139,7 @@ function init_elements_2d_manifold_in_3d!(elements,
140139
contravariant_vectors, inverse_jacobian) = elements
141140

142141
# The standard calc_node_coordinates! can be used, since Trixi.jl now dispatches on
143-
# P4estMesh{NDIMS, NDIMS_AMBIENT}, so it can be used here.
142+
# P4estMesh{NDIMS, NDIMS_AMBIENT}.
144143
Trixi.calc_node_coordinates!(node_coordinates, mesh, basis)
145144

146145
for element in 1:Trixi.ncells(mesh)

0 commit comments

Comments
 (0)