Skip to content

Commit

Permalink
fix initial condition
Browse files Browse the repository at this point in the history
  • Loading branch information
@tristanmontoya
tristanmontoya committed Dec 1, 2024
1 parent b86ed28 commit 68f458b
Show file tree
Hide file tree
Showing 4 changed files with 54 additions and 56 deletions.
2 changes: 1 addition & 1 deletion src/equations/covariant_advection.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ as
```
where $G$ is the determinant of the covariant metric tensor expressed as a matrix with
entries $G_{ij} = \vec{a}_i \cdot \vec{a}_j$. Note that the variable advection velocity
components could also be stored as auxiliary variables, similarly to the geometric
components could alternatively be stored as auxiliary variables, similarly to the geometric
information. The third velocity component $v^3$, representing the velocity component in the
direction normal to the surface, is set to zero when solving equations on two-dimensional
surfaces, and will remain zero throughout the simulation.
Expand Down
43 changes: 14 additions & 29 deletions src/equations/equations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,10 +73,19 @@ dispatching on the return type.
transform_initial_condition(initial_condition, equations)
Takes in a function with the signature `initial_condition(x, t)` which returns an initial
condition given in terms of global velocity or momentum components, and returns another
function with the signature `initial_condition_transformed(x, t, aux_vars, equations)`
which returns the same initial condition with the velocity or momentum vector given in
terms of contravariant components.
condition given in terms of global Cartesian or zonal/meridional velocity components, and
returns another function `initial_condition_transformed(x, t, equations)` or
`initial_condition_transformed(x, t, aux_vars, equations)` which returns the same initial
data, but transformed to the appropriate prognostic variables used internally by the
solver. For the covariant form, this involves a transformation of the global velocity
components to contravariant components using `aux_vars` as well as a conversion from
primitive to conservative variables. For standard Cartesian formulations, this simply
involves a conversion from primitive to conservative variables.
!!! note
When using the covariant formulation, the initial velocity components should be defined
in the coordinate system specified by the `GlobalCoordinateSystem` type parameter in
[`AbstractCovariantEquations`](@ref).
!!!
"""
function transform_initial_condition(initial_condition, ::AbstractCovariantEquations)
function initial_condition_transformed(x, t, aux_vars, equations)
Expand All @@ -87,6 +96,7 @@ function transform_initial_condition(initial_condition, ::AbstractCovariantEquat
return initial_condition_transformed
end

# Version of transform
function transform_initial_condition(initial_condition, ::AbstractEquations)
function initial_condition_transformed(x, t, equations)
return Trixi.prim2cons(initial_condition(x, t, equations), equations)
Expand Down Expand Up @@ -147,31 +157,6 @@ end
return aux_vars[13]
end

# Transform zonal and meridional velocity/momentum components to Cartesian components
function spherical2cartesian(vlon, vlat, x)
# Co-latitude
colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2))

# Longitude
if sign(x[2]) == 0.0
signy = 1.0
else
signy = sign(x[2])
end
r_xy = sqrt(x[1]^2 + x[2]^2)
if r_xy == 0.0
lon = pi / 2
else
lon = signy * acos(x[1] / r_xy)
end

v1 = -cos(colat) * cos(lon) * vlat - sin(lon) * vlon
v2 = -cos(colat) * sin(lon) * vlat + cos(lon) * vlon
v3 = sin(colat) * vlat

return SVector(v1, v2, v3)
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
Expand Down
33 changes: 23 additions & 10 deletions src/equations/reference_data.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,27 +47,32 @@ the test suite described in the following paper:
RealT = eltype(x)

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
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 / (a^2) # bump width
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))

# 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)

# 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))

# get zonal and meridional components of the velocity
lon, lat = atan(x[2], x[1]), asin(x[3] / a)
vlon = V * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha))
vlat = -V * sin(lon) * sin(alpha)
vx, vy, vz = spherical2cartesian(vlon, vlat, x)
# get Cartesian velocity components
vx, vy, vz = omega * Trixi.cross(axis, x)

# Prescribe the Cartesian velocity components
# Prescribe the rotated bell shape and Cartesian velocity components
return SVector(h, vx, vy, vz, 0.0f0)
end

Expand All @@ -78,24 +83,32 @@ end
RealT = eltype(x)

a = EARTH_RADIUS # radius of the sphere in metres
V = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY) # speed of rotation
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 / (a^2) # bump width
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))

# 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)

# 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))

# get zonal and meridional components of the velocity
lon, lat = atan(x[2], x[1]), asin(x[3] / a)
vlon = V * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha))
vlat = -V * sin(lon) * sin(alpha)
vlon = omega * a * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha))
vlat = -omega * a * sin(lon) * sin(alpha)

# Prescribe the spherical velocity components
return SVector(h, vlon, vlat, 0.0f0)
Expand Down
32 changes: 16 additions & 16 deletions test/test_spherical_advection.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,17 +11,17 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples")
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_shallowwater_cubed_sphere_shell_advection.jl"),
l2=[
0.7963221028128544,
20.317489255709344,
8.811382597221147,
20.317512894705885,
0.7963216338636225,
20.317829852541713,
8.810001095824774,
20.31782985254294,
0.0
],
linf=[
10.872101730749478,
289.65159632622454,
95.16499199537975,
289.65159632621726,
10.872101731796079,
289.65159635384043,
95.12887120427786,
289.65159635384407,
0.0
])
# and small reference values
Expand All @@ -39,17 +39,17 @@ end
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_shallowwater_cubed_sphere_shell_advection.jl"),
l2=[
0.8933384470346987,
22.84334163690791,
9.776600357597692,
22.843351937152637,
0.8933429672952714,
22.84887991902509,
9.758850586757735,
22.84887991902542,
0.0
],
linf=[
14.289456304666146,
380.6958334078372,
120.59259301636666,
380.6958334078299,
14.289456304624764,
380.6958334067349,
120.59259301602742,
380.69583340674217,
0.0
], element_local_mapping=true)
# and small reference values
Expand Down

0 comments on commit 68f458b

Please sign in to comment.