Discontinuous Galerkin 2-layer shallow water method does not converge.

[garrettdreyfus]

Hello everyone!

I have been trying to solve a simplified version of the steady-state 2-D two layer shallow water equations with a prescribed entertainment and rigid lid. The equation set is as follows:

Where h_i is the layer thickness, u_i is the layer velocity, w is the entrainment, pi is the rigid lid pressure, eta is the depth of the interface between the two layers, Cd is a friction coefficient and b is a constant density difference between the two layers.

I initially solved these equations very easily using Continuous Galerkin elements and found some success in a rectangular domain with four no flux boundaries, but the solution was not very stable (it would break with small changes in the parameters) and it showed weird alternating velocities at the boundaries. The image below shows y-velocity in the second layer along the eastern wall.

I have been trying to adapt my solution using Discrete Galerkin elements for pressure and interface height as well as BDM elements for the velocity. I have been heavily drawing from the FireDrake steady-state continuity tutorial and Compatible Finite Element Methods for Geophysical Flows: Automation and Implementation Using Firedrake by Gibson et al.

I have included the code below. It runs correctly if I use Continuous Galerkin elements for every field, but if I switch to BDM and DG it does not converge. My eventual goal is to have an open boundary as well.

Thank you for any and all help in advance, and let me know if there is anything I can clarify!

from firedrake import *  

# Model domain 200x300km at 1km resolution
mesh = RectangleMesh(100,150,200000,300000)
# Contains mesh coordinates
x = SpatialCoordinate(mesh)

# reduced gravity
gprime=Constant(0.003)

Cmelt=Constant(3.1709792e-8)
rhofrac = Constant(991/1025)

# Function spaces for all our unknowns
U1 = FunctionSpace(mesh, "BDM", 2)
U2 = FunctionSpace(mesh, "BDM", 2)
ETA_FRAC = FunctionSpace(mesh, "DG", 1)
PI = FunctionSpace(mesh, "DG", 1)

# Mixed Function space with all our unknowns in their
V = U1 * U2 * ETA_FRAC * PI

#Split up ouor function space 
Z = Function(V)
pi_coeff = pi
u1, u2, eta_frac,  pi = split(Z)
u1_t, u2_t, eta_frac_t,  pi_t = TestFunctions(V)

def perp(u):
    return as_vector((-u[1], u[0]))

f_cor = Constant(-10**-4)
Cd = 10**-4

H = Constant(-500.0)
z_b = -0.0005*(300000-x[1])

Vf = FunctionSpace(mesh, "CG", 1)
entfunc = -sin(2*pi_coeff*(x[1]/300000.0))*(10**-7)
entfunc = Function(Vf).interpolate(entfunc)
z_b = Function(Vf).interpolate(z_b)
H = Function(Vf).interpolate(H)
#entfunc = (x[1]*0)+Constant(10**(-7))
#entfunc = Constant(0)
#entfunc  = ((x[1]-150000)/300000)*(10**-7) +(10**-8)
h1 = (z_b - H)*eta_frac
h2 = (z_b - H)-h1
eta = H+h2

n = FacetNormal(mesh)
un1 = 0.5*(dot(u1, n) + abs(dot(u1, n)))
un2 = 0.5*(dot(u2, n) + abs(dot(u2, n)))



h_eq = (h1*inner(u1,grad(eta_frac_t)) - eta_frac_t*entfunc)*dx - dot(jump(eta_frac_t), un1('+')*h1('+') - un1('-')*h1('-'))*dS - dot(eta_frac_t, un1*h1)*ds

pi_eq =  (h2*inner(u2,grad(pi_t)) + pi_t*entfunc)*dx - dot(jump(pi_t), un2('+')*h2('+') - un2('-')*h2('-'))*dS - dot(pi_t, un2*h2)*ds

u1_eq = (f_cor*inner(u1_t,perp(u1)) + inner(div(u1_t), pi)+ inner(u1_t,Cd*u1/h1))*dx - (dot(u1_t,n*pi))*ds 

u2_eq = (f_cor*inner(u2_t,perp(u2)) + inner(div(u2_t), gprime*eta + pi) + inner(u2_t,Cd*u2/h2))*dx - (dot(u1_t,n*(gprime*eta + pi)))*ds  


L = h_eq + u1_eq + u2_eq + pi_eq# + h_c  # + u1_visc + u2_visc

V_p = V.sub(2)

mixed_lb = Function(V)
mixed_ub = Function(V)


with mixed_lb.sub(0).dat.vec as vec:
    vec.set(float("-inf"))  # No lower bound for velocity
with mixed_ub.sub(0).dat.vec as vec:
    vec.set(float("inf"))  # No lower bound for velocity
with mixed_lb.sub(1).dat.vec as vec:
    vec.set(float("-inf"))  # No lower bound for velocity
with mixed_ub.sub(1).dat.vec as vec:
    vec.set(float("inf"))  # No lower bound for velocity

# Create Functions on the pressure subspace
lb = Function(V_p)
ub = Function(V_p)
# Assign the bounds
lower_bound_expr = Constant(0.01)
upper_bound_expr = Constant(0.99)
lb.interpolate(lower_bound_expr)
ub.interpolate(upper_bound_expr)
# Interpolate to functions
mixed_lb.sub(2).assign(lb)
mixed_ub.sub(2).assign(ub)
with mixed_lb.sub(3).dat.vec as vec:
    vec.set(float("-inf"))  # No lower bound for velocity
with mixed_ub.sub(3).dat.vec as vec:
    vec.set(float("inf"))  # No lower bound for velocity

solver_parameters={"ksp_type": "preonly",
    #'snes_type': 'vinewtonrsls',
    "mat_type": "aij",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    'snes_monitor': None,
    'snes_view':None,
    }


#V = U1 * U2 * ETA_FRAC * PI
bcs_u1 = [DirichletBC(V.sub(0), Constant((0.0,0.0)), 1),DirichletBC(V.sub(0), Constant((0.0,0.0)), 2),DirichletBC(V.sub(0), Constant((0.0,0.0)), 3),DirichletBC(V.sub(0), Constant((0.0,-0.00)), 4)]
bcs_u2 = [DirichletBC(V.sub(1), Constant((0.0,0.0)), 1),DirichletBC(V.sub(1), Constant((0.0,0.0)), 2),DirichletBC(V.sub(1), Constant((0.0,0.0)), 3),DirichletBC(V.sub(1), Constant((0.0,-0.00)), 4)]

#bcs_u1 = []
#bcs_u2 = []

bcs_eta_frac = []
bcs_ent = []
bcs_pi = []
bcs = bcs_u1+bcs_u2+bcs_eta_frac+bcs_ent+bcs_pi

Z.sub(0).interpolate(Constant((0.000,-0.0001)))
Z.sub(1).interpolate(Constant((0.000,0.001)))
Z.sub(2).interpolate(Constant(frac_enforce))
Z.sub(3).interpolate(Constant(0))

problem = NonlinearVariationalProblem(L , Z, bcs)
solver = NonlinearVariationalSolver(problem,solver_parameters=solver_parameters)
print("Solving began")
#solver.solve(bounds=(mixed_lb, mixed_ub)) # Solve for updated fields
solver.solve() # Solve for updated fields

wn1 = Function(V)
# mixed func. for both fields (n+1)
u1_out = Function(U1, name="Velocity")
div1 = Function(ETA_FRAC, name="DIV1")
div2 = Function(ETA_FRAC, name="DIV2")
u2_out = Function(U2, name="Velocity 2")
eta_frac_out = Function(ETA_FRAC, name="Interface depth")
pi_out = Function(PI, name="Rigid Lid Pressure")
ent_out = Function(PI, name="Entrainment")

outfile = VTKFile("results/wholething/out.pvd")

h1 = (z_b - H)*Z.sub(2)
h2 = (z_b - H)-h1
eta = H+h2

u1_out.assign(Z.sub(0))
u2_out.assign(Z.sub(1))
eta_frac_out.interpolate(eta)
div1.interpolate(div(Z.sub(0)*h1))
div2.interpolate(div(Z.sub(1)*h2))
pi_out.assign(Z.sub(3))
ent_out.assign(entfunc)

#V = U1 * U2 * ETA_FRAC * ENT * PI

outfile.write(u1_out,u2_out,eta_frac_out,pi_out,div1,div2,ent_out)