Skip to content

Commit

Permalink
[ci-skip] First commit of the branch to add a refined mesh test for t…
Browse files Browse the repository at this point in the history
…he Rosenbrock shallow water test case with the Williamson5 test configuration. This includes the Williamson 5 initial conditions and driver, as well as a change to the Rosenbrock solver to allow for bottom topography
  • Loading branch information
davelee2804 committed Nov 11, 2024
1 parent 6603998 commit 913a1bc
Show file tree
Hide file tree
Showing 3 changed files with 157 additions and 5 deletions.
29 changes: 24 additions & 5 deletions src/ShallowWaterRosenbrock.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,8 @@ function shallow_water_rosenbrock_time_step!(
model, dΩ, dω, Y, V, Q, R, S, f, g, y₁, y₀, # in args
RTMMchol, L2MMchol, Amat, Bchol, Blfchol, # more in args
dt, τ, leap_frog,
assem=SparseMatrixAssembler(SparseMatrixCSC{Float64,Int},Vector{Float64},R,S)) # ...yet more in args
assem=SparseMatrixAssembler(SparseMatrixCSC{Float64,Int},Vector{Float64},R,S),
topog=nothing)
# energetically balanced second order rosenbrock shallow water solver
# reference: eqns (24) and (39) of
# https://github.com/BOM-Monash-Collaborations/articles/blob/main/energetically_balanced_time_integration/EnergeticallyBalancedTimeIntegration_SW.tex
Expand All @@ -29,7 +30,11 @@ function shallow_water_rosenbrock_time_step!(
# 1.1: the mass flux
compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*h₁)
# 1.2: the bernoulli function
compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,u₁u₁,h₁,g)
if topog==nothing
compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,u₁u₁,h₁,g)
else
compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,u₁u₁,h₁+topog,g)
end
# 1.3: the potential vorticity
compute_potential_vorticity!(q₁,H1h,H1hchol,dΩ,R,S,h₁,u₁,f,n,assem)
# 1.4: assemble the momentum and continuity equation residuals
Expand All @@ -45,7 +50,11 @@ function shallow_water_rosenbrock_time_step!(
# 2.1: the mass flux
compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*(2.0*h₁ + h₂)/6.0+u₂*(h₁ + 2.0*h₂)/6.0)
# 2.2: the bernoulli function
compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,(u₁u₁ + u₁u₂ + u₂u₂)/3.0,0.5*(h₁ + h₂),g)
if topog==nothing
compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,(u₁u₁ + u₁u₂ + u₂u₂)/3.0,0.5*(h₁ + h₂),g)
else
compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,(u₁u₁ + u₁u₂ + u₂u₂)/3.0,0.5*(h₁ + h₂)+topog,g)
end
# 2.3: the potential vorticity
compute_potential_vorticity!(q₂,H1h,H1hchol,dΩ,R,S,h₂,u₂,f,n,assem)
# 2.4: assemble the momentum and continuity equation residuals
Expand Down Expand Up @@ -77,6 +86,7 @@ function shallow_water_rosenbrock_time_stepper(
λ, dt, τ, N;
mass_matrix_solver::Gridap.Algebra.LinearSolver=Gridap.Algebra.BackslashSolver(),
jacobian_matrix_solver::Gridap.Algebra.LinearSolver=Gridap.Algebra.BackslashSolver(),
tₒ=nothing,
leap_frog=false,
write_diagnostics=true,
write_diagnostics_freq=1,
Expand Down Expand Up @@ -134,6 +144,15 @@ function shallow_water_rosenbrock_time_stepper(
solve!(rhs2, Mchol, rhs2)
yn = FEFunction(Y, rhs2)

# project the bottom topography onto the L2 space
topog = nothing
if t₀ != nothing
b₄(q) = (q*t₀)dΩ
rhs3 = assemble_vector(b₄,Q)
topog = FEFunction(Q, copy(rhs3))
solve!(get_free_dof_values(topog),L2MMchol,get_free_dof_values(topog))
end

un, hn = yn

hnv, fv, ynv = get_free_dof_values(hn,f,yn)
Expand Down Expand Up @@ -171,7 +190,7 @@ function shallow_water_rosenbrock_time_stepper(
istep = 1
shallow_water_rosenbrock_time_step!(yn, ϕ, F, q1, q2, duh1, duh2, H1h, H1hchol, y_wrk,
model, dΩ, dω, Y, V, Q, R, S, f, g, ym1, ym2,
RTMMchol, L2MMchol, A, Bchol, Blfchol, dt, τ, false)
RTMMchol, L2MMchol, A, Bchol, Blfchol, dt, τ, false, topog)

if (write_diagnostics && write_diagnostics_freq>0 && mod(istep, write_diagnostics_freq) == 0)
compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(Ω))
Expand All @@ -189,7 +208,7 @@ function shallow_water_rosenbrock_time_stepper(
yn=aux
shallow_water_rosenbrock_time_step!(yn, ϕ, F, q1, q2, duh1, duh2, H1h, H1hchol, y_wrk,
model, dΩ, dω, Y, V, Q, R, S, f, g, ym1, ym2,
RTMMchol, L2MMchol, A, Bchol, Blfchol, dt, τ, leap_frog)
RTMMchol, L2MMchol, A, Bchol, Blfchol, dt, τ, leap_frog, topog)

# IMPORTANT NOTE: We need to extract un, hn out of yn at each iteration because
# the association of yn with its object instance changes at the beginning of
Expand Down
55 changes: 55 additions & 0 deletions test/mpi/Williamson5InitialConditions.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
# Initial conditions for the zonal flow over an isolated mountain ("Williamson 5")
# shallow water test case. Involves a balanced geostrophic flow with a Gaussian
# hill in the northern hemisphere.
# reference:
# D.L. Williamson, J.B. Drake, J.J. Hack, R. Jakob, P.N. Swarztrauber,
# J. Comput. Phys. 102 (1992) 211–224.

const Rₑ = 6371220.0
const gₑ = 9.80616
const U₀ = 20.0
const Ωₑ = 7.292e-5
const H₀ = 5960.0
const α₀ = 0.0

function (θϕr)
θ,ϕ,r = θϕr
U₀*(cos(ϕ)*cos(α₀) + cos(θ)*sin(ϕ)*sin(α₀))
end

function (θϕr)
θ,ϕ,r = θϕr
-U₀*sin(θ)*sin(α₀);
end

# Initial velocity
function u₀(xyz)
θϕr = xyz2θϕr(xyz)
u = (θϕr)
v = (θϕr)
spherical_to_cartesian_matrix(θϕr)VectorValue(u,v,0)
end

# Topography
function topography(xyz)
θϕr = xyz2θϕr(xyz)
θc = -π/2.0
ϕc = π/6.0
bo = 2000.0
rad = π/9.0
rsq =- ϕc)*- ϕc) +- θc)*- θc)
r = sqrt(rsq)
b = 0.0
if(r < rad)
b = bo*(1.0 - r/rad)
end
b
end

# Initial fluid depth
function h₀(xyz)
θϕr = xyz2θϕr(xyz)
b = -cos(θ)*cos(ϕ)*sin(α₀) + sin(ϕ)*cos(α₀)
bt = topography(xyz)
H₀ - (Rₑ*Ωₑ*U₀ + 0.5*U₀*U₀)*b*b/gₑ - bt
end
78 changes: 78 additions & 0 deletions test/mpi/Williamson5ShallowWaterRosenbrock.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
module Williamson5ShallowWaterRosenbrock

using PartitionedArrays
using Test
using FillArrays
using Gridap
using GridapPETSc
using GridapGeosciences
using GridapDistributed
using GridapP4est

include("Williamson5InitialConditions.jl")

function petsc_gamg_options()
"""
-ksp_type gmres -ksp_rtol 1.0e-06 -ksp_atol 0.0
-ksp_monitor -pc_type gamg -pc_gamg_type agg
-mg_levels_esteig_ksp_type gmres -mg_coarse_sub_pc_type lu
-mg_coarse_sub_pc_factor_mat_ordering_type nd -pc_gamg_process_eq_limit 50
-pc_gamg_square_graph 9 pc_gamg_agg_nsmooths 1
"""
end
function petsc_mumps_options()
"""
-ksp_type preonly -ksp_error_if_not_converged true
-pc_type lu -pc_factor_mat_solver_type mumps
"""
end

order = 0
degree = 4

# magnitude of the descent direction of the implicit solve;
# neutrally stable for 0.5, L-stable for 1+sqrt(2)/2
λ = 1.0 + 0.5*sqrt(2.0)

dt = 480.0
nstep = Int(24*60^2*20/dt) # 20 days

function main(distribute,parts)
ranks = distribute(LinearIndices((prod(parts),)))

# Change directory to the location of the script, where the mesh data files are located
cd(@__DIR__)

coarse_model, cell_panels, coarse_cell_wise_vertex_coordinates=
parse_cubed_sphere_coarse_model("williamson-5-C12/connectivity-gridapgeo.txt","williamson-5-C12/geometry-gridapgeo.txt")

num_uniform_refinements=0

GridapPETSc.with(args=split(petsc_mumps_options())) do
model=CubedSphereDiscreteModel(ranks,
coarse_model,
coarse_cell_wise_vertex_coordinates,
cell_panels,
num_uniform_refinements;
radius=rₑ,
adaptive=true,
order=0)

hf, uf = shallow_water_rosenbrock_time_stepper(model, order, degree,
h₀, u₀, Ωₑ, gₑ, H₀,
λ, dt, 60.0, nstep;
t₀=topography,
leap_frog=true,
write_solution=true,
write_solution_freq=45,
write_diagnostics=true,
write_diagnostics_freq=1,
dump_diagnostics_on_screen=true)

end
end

with_mpi() do distribute
main(distribute,1)
end
end # module

0 comments on commit 913a1bc

Please sign in to comment.