Crouzeix Raviart SWE's formulation

variational formulation
1

Hi, i get the following error in my formulation: Error: Unable to solve linear system using PETSc Krylov solver. Here’s some simplified code:

from dolfin import *
#Parameters:
f=-1.028e-4;
g=9.81;
c_d=Constant(0.012); 
α_h=Constant(1.0);
R = as_matrix([[0, -1],[1, 0]])
T = 44400;
dt=10;
num_steps=4440;
k  = Constant(dt)
#dummy mesh, with this it works
#mesh=UnitSquareMesh(10,10)

#the real mesh (i think the problem is here) its loaded as
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
facets = MeshFunction("size_t", mesh, 1) 

V = VectorElement('CR',triangle,1)
Q = FiniteElement('Lagrange',triangle,1)
element = Q*V
Z = FunctionSpace(mesh,element)
z = TrialFunction(Z)
(η,u) = split(z)
(ϕ,w) = TestFunctions(Z)
z_n = Function(Z)
z_out = Function(Z)
z1 = Function(Z)
η_n, u_n = split(z_n)
η_out, u_out = split(z_out)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
tol = 1.E-4

class b_cerr(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]<8.5536e4-tol and x[1]>-3.0412e4+tol
b_cerr = b_cerr()
b_cerr.mark(boundary_markers,0)

class b_ab(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]>8.553e4-tol or x[1]<-3.0412e4+tol
b_ab = b_ab()
b_ab.mark(boundary_markers,1)

bcη_cerr = DirichletBC(Z.sub(0), Constant(0.0), b_cerr)
bcu_cerr = DirichletBC(Z.sub(1), Constant((0.0,0.0)), b_cerr)

#this is a dummy boundary condition, the actual i use its interpolated from data and updated every time step
bcη_ab = DirichletBC(Z.sub(0), Constant(1.0), b_ab)
bcu_ab = DirichletBC(Z.sub(1), Constant((2.0,3.0)), b_ab)

bcs=[bcη_cerr, bcu_cerr, bcη_ab, bcu_ab]

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
ν = FacetNormal(mesh)

from ufl import sign
def norma2(u):
    return sqrt(dot(u,u) + DOLFIN_EPS)
def avg_gamma(u, gamma):
     return (0.5+gamma("+"))*u("+") + (0.5-gamma("-"))*u("-")
lbda = 0.5*sign(dot(u_n,ν))
θ=0.5
H=η+h
H_n=η_n+h


B = inner((η-η_out)/(2*k),ϕ)*dx + inner((u-u_n)/k,w)*dx \
    -θ*inner(h*u,grad(ϕ))*dx \
    -(1-θ)*inner(h*u_out,grad(ϕ))*dx \
    +inner(η_n*u_n,grad(ϕ))*dx \
    -dot(u_n,div(outer(w,u_n)))*dx \
    +θ*(inner(f*R*u+1/H_n*(c_d*norma2(u_n)*u),w)*dx + g*inner(grad(η),w)*dx) \
    +(1-θ)*(inner(f*R*u_n+1/H_n*(c_d*norma2(u_n)*u_n),w)*dx + g*inner(grad(η_n),w)*dx) \
    +inner(avg_gamma(dot(u_n,ν)*u_n,lbda),jump(w))*dS

a1 = lhs(B)
L1 = rhs(B)

prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
parameters['krylov_solver']['nonzero_initial_guess'] = True

η_i=project(Constant(0.0),Z.sub(0).collapse())
u_i=project(Constant((0.0,0.0)),Z.sub(1).collapse())
assign(z_n.sub(0), η_i)
assign(z_n.sub(1), u_i)
assign(z_out.sub(0), η_i)
assign(z_out.sub(1), u_i)

t = dt
for m in range(num_steps):      
    A1 = assemble(a1)    
    [bc.apply(A1) for bc in bcs] 
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcs]
    solve(A1, z1.vector(), b1, "bicgstab", "default")    
    η1, u1 = z1.split()    
    # Move to next time step
    assign(z_n.sub(0), η1)
    assign(z_n.sub(1), u1)    
    η_n,u_n = z_n.split()    
    assign(z_out.sub(0), η_n)
    assign(z_out.sub(1), u_n)      
    t += dt

I don’t know how to do this properly so i upload the mesh and the h function
mesh_xdmf
mesh_h5
h
I can guess there is more than 1 error in this, can you spot some? Thanks!

2

Hey,

I feel like this could be anything, from insufficient boundary conditions to problems with the meshing…

Maybe try a different solver, like ‘gmres’ with ‘ilu’ preconditioning. This combination usually works great for me. In my opinion defining a dedicated solver also gives more control, something like

solver = KrylovSolver('gmres', 'ilu')
prm    = solver.parameters
prm['absolute_tolerance']  = 1E-10
prm['relative_tolerance']  = 1E-6
prm['maximum_iterations']  = 10000
prm['monitor_convergence'] = True

....

solver.solve(A, z1.vector(), b)

Also in general for MWE, I think, a simpler mesh like

mesh = UnitSquareMesh(64,64)

is better. Also I feel like your MWE can be simplified (leaving out the function h, and stuff like this); I can only speak for myself, but I personally am more inclined to try something out, when I can just copy paste the code into a temporary file, without downloading many files.

Best Regards

Emanuel

2 Likes
3

Thanks for your reply Emanuel, i’ll try what you say. If i use the UnitSquareMesh and some easy h function, the code runs with no error message so i coudn’t show the error like that.

4

ok, but that sounds great, maybe try to isolate the problem and find which one causes the error (the mesh or the h-function). For meshing problems you could go all out with a super fine mesh, problems with h can maybe solved with a proper preconditioning method, a finer mesh, different function spaces, different solver… :slight_smile:
Let us know what works and what does not!

1 Like