Skip to content

Commit

Permalink
Merge pull request #2167 from CliMA/dy/autodiff_options
Browse files Browse the repository at this point in the history 
Fix vertical metric terms, add SEM option for horizontal metric terms
  • Loading branch information
@dennisYatunin
dennisYatunin authored Feb 4, 2025
2 parents 83b6ccc + 68f07f4 commit 1d9f991
Show file tree
Hide file tree
Showing 8 changed files with 156 additions and 153 deletions.
2 changes: 1 addition & 1 deletion examples/hybrid/sphere/deformation_flow.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -405,7 +405,7 @@ compare_tracer_props(a, b; buffer = 1) = all(
@test compare_tracer_props(tracer_roughnesses(lim_third_upwind_sol), tracer_roughnesses(third_upwind_sol); buffer = 1.0)
@test compare_tracer_props(tracer_roughnesses(lim_fct_sol) , tracer_roughnesses(third_upwind_sol); buffer = 0.93)
@test compare_tracer_props(tracer_ranges(fct_sol) , tracer_ranges(third_upwind_sol); buffer = 1.0)
@test compare_tracer_props(tracer_ranges(lim_third_upwind_sol) , tracer_ranges(third_upwind_sol); buffer = 1.0)
@test compare_tracer_props(tracer_ranges(lim_third_upwind_sol) , tracer_ranges(third_upwind_sol); buffer = 1.2)
@test compare_tracer_props(tracer_ranges(lim_fct_sol) , tracer_ranges(third_upwind_sol); buffer = 1.0)
@test compare_tracer_props(tracer_ranges(lim_first_upwind_sol) , tracer_ranges(third_upwind_sol); buffer = 0.6)
@test compare_tracer_props(tracer_ranges(lim_centered_sol) , tracer_ranges(third_upwind_sol); buffer = 1.3)
Expand Down
94 changes: 37 additions & 57 deletions src/Grids/finitedifference.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,13 +56,21 @@ function _FiniteDifferenceGrid(topology::Topologies.IntervalTopology)
CT = Topologies.coordinate_type(mesh)
FT = Geometry.float_type(CT)
Nv_face = length(mesh.faces)
Nv_cent = Nv_face - 1
# construct on CPU, adapt to GPU
center_coordinates = DataLayouts.VF{CT, Nv_cent}(Array{FT}, Nv_cent)
face_coordinates = DataLayouts.VF{CT, Nv_face}(Array{FT}, Nv_face)
for v in 1:Nv_cent
center_coordinates[vindex(v)] = (mesh.faces[v] + mesh.faces[v + 1]) / 2
end
for v in 1:Nv_face
face_coordinates[vindex(v)] = mesh.faces[v]
end
center_local_geometry, face_local_geometry =
fd_geometry_data(face_coordinates, Val(Topologies.isperiodic(topology)))
center_local_geometry, face_local_geometry = fd_geometry_data(
center_coordinates,
face_coordinates,
Val(Topologies.isperiodic(topology)),
)

return FiniteDifferenceGrid(
topology,
Expand All @@ -74,84 +82,56 @@ end

# called by the FiniteDifferenceGrid constructor, and the ExtrudedFiniteDifferenceGrid constructor with Hypsography
function fd_geometry_data(
center_coordinates::DataLayouts.AbstractData{Geometry.ZPoint{FT}},
face_coordinates::DataLayouts.AbstractData{Geometry.ZPoint{FT}},
::Val{periodic},
) where {FT, periodic}
CT = Geometry.ZPoint{FT}
AIdx = (3,)
LG = Geometry.LocalGeometry{AIdx, CT, FT, SMatrix{1, 1, FT, 1}}
∂x∂ξ_axes = (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}())
(Ni, Nj, Nk, Nv, Nh) = size(face_coordinates)
Nv_face = Nv - periodic
Nv_cent = Nv - 1
center_local_geometry = similar(face_coordinates, LG, Val(Nv_cent))
center_local_geometry = similar(center_coordinates, LG, Val(Nv_cent))
face_local_geometry = similar(face_coordinates, LG, Val(Nv_face))
c1(args...) =
cent_coord(args...) =
Geometry.component(center_coordinates[CartesianIndex(args...)], 1)
face_coord(args...) =
Geometry.component(face_coordinates[CartesianIndex(args...)], 1)
for h in 1:Nh, k in 1:Nk, j in 1:Nj, i in 1:Ni
for v in 1:Nv_cent
# centers
coord⁻ = c1(i, j, k, v, h)
coord⁺ = c1(i, j, k, v + 1, h)
# use a "discrete Jacobian"
coord = (coord⁺ + coord⁻) / 2
Δcoord = coord⁺ - coord⁻
J = Δcoord
WJ = Δcoord
∂x∂ξ = SMatrix{1, 1}(J)
J = face_coord(i, j, k, v + 1, h) - face_coord(i, j, k, v, h)
WJ = J
x = CT(cent_coord(i, j, k, v, h))
∂x∂ξ = Geometry.AxisTensor(∂x∂ξ_axes, SMatrix{1, 1}(J))
center_local_geometry[CartesianIndex(i, j, k, v, h)] =
Geometry.LocalGeometry(
CT(coord),
J,
WJ,
Geometry.AxisTensor(
(
Geometry.LocalAxis{AIdx}(),
Geometry.CovariantAxis{AIdx}(),
),
∂x∂ξ,
),
)
Geometry.LocalGeometry(x, J, WJ, ∂x∂ξ)
end
for v in 1:Nv_face
coord = c1(i, j, k, v, h)
if v == 1
if periodic && v == 1
# periodic boundary face
J⁺ = face_coord(i, j, k, 2, h) - face_coord(i, j, k, 1, h)
J⁻ = face_coord(i, j, k, Nv, h) - face_coord(i, j, k, Nv - 1, h)
J = (J⁺ + J⁻) / 2
WJ = J
elseif !periodic && v == 1
# bottom face
if periodic
Δcoord⁺ = c1(i, j, k, 2, h) - c1(i, j, k, 1, h)
Δcoord⁻ = c1(i, j, k, Nv, h) - c1(i, j, k, Nv - 1, h)
J = (Δcoord⁺ + Δcoord⁻) / 2
WJ = J
else
coord⁺ = c1(i, j, k, 2, h)
J = coord⁺ - coord
WJ = J / 2
end
elseif v == Nv_cent + 1
@assert !periodic
J = face_coord(i, j, k, 2, h) - face_coord(i, j, k, 1, h)
WJ = J / 2
elseif v == Nv
# top face
coord⁻ = c1(i, j, k, Nv - 1, h)
J = coord - coord⁻
@assert !periodic
J = face_coord(i, j, k, Nv, h) - face_coord(i, j, k, Nv - 1, h)
WJ = J / 2
else
coord⁺ = c1(i, j, k, v + 1, h)
coord⁻ = c1(i, j, k, v - 1, h)
J = (coord⁺ - coord⁻) / 2
J = cent_coord(i, j, k, v, h) - cent_coord(i, j, k, v - 1, h)
WJ = J
end
∂x∂ξ = SMatrix{1, 1}(J)
x = CT(face_coord(i, j, k, v, h))
∂x∂ξ = Geometry.AxisTensor(∂x∂ξ_axes, SMatrix{1, 1}(J))
face_local_geometry[CartesianIndex(i, j, k, v, h)] =
Geometry.LocalGeometry(
CT(coord),
J,
WJ,
Geometry.AxisTensor(
(
Geometry.LocalAxis{AIdx}(),
Geometry.CovariantAxis{AIdx}(),
),
∂x∂ξ,
),
)
Geometry.LocalGeometry(x, J, WJ, ∂x∂ξ)
end
end
return (center_local_geometry, face_local_geometry)
Expand Down
Loading

0 comments on commit 1d9f991

Please sign in to comment.