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!