Skip to content

Commit

Permalink
Remove Gusto-specific test
Browse files Browse the repository at this point in the history
Not needed now that we actually test Gusto properly.
  • Loading branch information
connorjward committed Feb 7, 2025
1 parent d2df60f commit 0d556ac
Showing 1 changed file with 0 additions and 135 deletions.
135 changes: 0 additions & 135 deletions tests/firedrake/slate/test_slate_hybridization.py
Original file line number Diff line number Diff line change
Expand Up @@ -464,141 +464,6 @@ def test_mixed_poisson_approximated_schur_jacobi_prec(setup_poisson):
assert u_err < 1e-8


def test_slate_hybridization_gusto():

# ---------------------------------------------------------------------------- #
# Set up core objects
# ---------------------------------------------------------------------------- #

degree = 1
radius = Constant(6371220.)
dt = Constant(3000.0)
H = Constant(3e4/9.80616)
g = Constant(9.80616)
Omega = Constant(7.292e-5)
alpha = Constant(0.5)
beta_u = alpha*dt
beta_d = alpha*dt

mesh = IcosahedralSphereMesh(float(radius), refinement_level=0, degree=2)
x = SpatialCoordinate(mesh)

# Function spaces
VDG = FunctionSpace(mesh, "DG", degree)
VHDiv = FunctionSpace(mesh, "BDM", degree+1)
mesh.init_cell_orientations(x)

VCG1 = FunctionSpace(mesh, 'CG', 1)
Vmixed = MixedFunctionSpace((VHDiv, VDG, VDG))
Vreduced = MixedFunctionSpace((VHDiv, VDG))

# Mixed RHS and LHS functions
x_rhs = Function(Vmixed)

# Components of various mixed functions
u_rhs, D_rhs, b_rhs = split(x_rhs)

uD_lhs = Function(Vreduced)
u_test, D_test = TestFunctions(Vreduced)
u_trial, D_trial = TrialFunctions(Vreduced)
u_lhsr, D_lhsr = uD_lhs.subfunctions

# Reference profiles
uDb_bar = Function(Vmixed)
D_bar = split(uDb_bar)[1]
b_bar = split(uDb_bar)[2]
D_bar_subfunc = uDb_bar.subfunctions[1]
b_bar_subfunc = uDb_bar.subfunctions[2]

# Set up perp function
sphere_degree = mesh.coordinates.function_space().ufl_element().degree()
VecDG = VectorFunctionSpace(mesh, 'DG', sphere_degree)
N = Function(VecDG).interpolate(CellNormal(mesh))
perp = lambda u: cross(N, u)

# Coriolis
f = Function(VCG1)

# ---------------------------------------------------------------------------- #
# Set up problem and solvers
# ---------------------------------------------------------------------------- #

# We have a 3-function mixed space for wind, depth and buoyancy
# Strategy is to approximately eliminate buoyancy, leaving a wind-depth problem
# The wind-depth problem is solved using the hydridization preconditioner

# Approximate elimination of b
b = -0.5 * dt * dot(u_trial, grad(b_bar)) + b_rhs

solver_parameters = {
'ksp_type': 'preonly',
'mat_type': 'matfree',
'pc_type': 'python',
'pc_python_type': 'firedrake.HybridizationPC',
'hybridization': {
'ksp_type': 'cg',
'pc_type': 'gamg',
'ksp_rtol': 1e-1,
'mg_levels': {
'ksp_type': 'chebyshev',
'ksp_max_it': 1,
'pc_type': 'bjacobi',
'sub_pc_type': 'ilu'
}
}
}

n = FacetNormal(mesh)

# Linear thermal shallow water problem
uD_eqn = (
inner(u_test, (u_trial - u_rhs)) * dx
- beta_u * (D_trial - D_bar) * div(u_test * b_bar) * dx
+ beta_u * jump(u_test * b_bar, n) * avg(D_trial - D_bar) * dS
- beta_u * 0.5 * D_bar * b_bar * div(u_test) * dx
- beta_u * 0.5 * D_bar * b * div(u_test) * dx
- beta_u * 0.5 * b_bar * div(u_test*(D_trial - D_bar)) * dx
+ beta_u * 0.5 * jump((D_trial - D_bar)*u_test, n) * avg(b_bar) * dS
+ inner(D_test, (D_trial - D_rhs)) * dx
+ beta_d * D_test * div(D_bar*u_trial) * dx
+ beta_u * f * inner(u_test, perp(u_trial)) * dx
)

# Boundary conditions
bcs = []

# Solver for u, D
uD_problem = LinearVariationalProblem(
lhs(uD_eqn), rhs(uD_eqn), uD_lhs, bcs=bcs, constant_jacobian=True
)

# Provide callback for the nullspace of the trace system
def trace_nullsp(T):
return VectorSpaceBasis(constant=True)

appctx = {"trace_nullspace": trace_nullsp}
uD_solver = LinearVariationalSolver(
uD_problem, solver_parameters=solver_parameters, appctx=appctx
)

# ---------------------------------------------------------------------------- #
# Set some initial conditions
# ---------------------------------------------------------------------------- #

f.interpolate(2*Omega*x[2]/radius)
D_bar_subfunc.assign(H)
b_bar_subfunc.assign(g)

# Simplest test is to give a right-hand side that is zero
x_rhs.assign(1.0)

# ---------------------------------------------------------------------------- #
# Run
# ---------------------------------------------------------------------------- #

uD_solver.solve()


@pytest.mark.parametrize('counts', [(10001, 10002), (10002, 10001)])
def test_slate_hybridization_count_safe(counts):
g_count, c_count = counts
Expand Down

0 comments on commit 0d556ac

Please sign in to comment.