Hey all,
I am attempting to solve a coupled Navier-Stokes problem, however once I make the meshh a bit more complex, from a straight tube to a bifurcatingg problem, the solver diverges after 0 iterations and I aam not sure why. I get the error:
Solving nonlinear variational problem.
0 SNES Function norm 9.750182999874e-09
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_LINEAR_SOLVE.
rho , eta , nu, dt, vdim = 1000 0.003 3e-06 0.001 2540388
t= 0.001
Traceback (most recent call last):
File “Patient_specific.py”, line 135, in
solver.solve()
RuntimeError:*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with PETScSNESSolver.
*** Reason: Solver did not converge.
*** Where: This error was encountered inside PETScSNESSolver.cpp.
*** Process: 0
*** DOLFIN version: 2018.1.0
*** Git changeset: 9e1777c0df7c762c7230189c87c13fd5ad33a4f0
*** -------------------------------------------------------------------------
I’m not sure what is causing this, I have tried a a few varientions with lower pressures input but not sure whats causing this. Here is MWE:
from dolfin import *
PETScOptions.set("snes_linesearch_monitor", "")
PETScOptions.set("snes_linesearch_type", "bt")
mesh_name = "P1_F1_smooth_600k.xdmf" #patient3_17mil.xdmf
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(), mesh_name) # Specify that the mesh is XDMF
xdmf.read(mesh) # Read the mesh from XDMF file
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0) # Read the boundary
indicators
xdmf.read(boundaries, 'boundaries')
inlet_vein = 11
outlet_vein = 10
walls = 13
inlet_graft = 12
################################
eta = 3e-3
rho = 1000
nu = eta/rho
dt = 0.01
T = .5
V_element = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q_element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W_element = MixedElement(V_element, Q_element) # Taylor-Hood
W = FunctionSpace(mesh, W_element)
vq = TestFunction(W) # Test function in the mixed space
delta_up = TrialFunction(W) # Trial function in the mixed space (XXX Note: for the increment!)
(delta_u, delta_p) = split(delta_up) # Function in each subspace to write the functional (XXX Note: for the increment!)
(v, q) = split( vq) # Test function in each subspace to write the functional
up = Function(W)
(u, p) = split(up)
up_prev = Function(W)
(u_prev, _) = split(up_prev)
n = FacetNormal(mesh)
P_out = Constant(0.)
P_in_v = Constant(0.05)
P_in_g = Constant(0.2)
ds = Measure("ds")(subdomain_data=boundaries)
F = ( inner(u, v)/Constant(dt)*dx # Implit Euler discretization
- inner(u_prev, v)/Constant(dt)*dx # Implit Euler discretization
+ nu*inner(grad(u), grad(v))*dx
+ inner(grad(u)*u, v)*dx
- div(v)*p*dx
+ div(u)*q*dx
+ (P_in_v * inner(v,n) * ds(inlet_vein))
+ (P_out * inner(v,n) * ds(outlet_vein))
+ (P_in_g * inner(v,n) * ds(inlet_graft))
)
J = derivative(F, up, delta_up)
walls_bc = DirichletBC(W.sub(0), Constant((0., 0., 0.)), boundaries, walls )
bc = [walls_bc]
snes_solver_parameters = {"nonlinear_solver": "snes",
"snes_solver": {"linear_solver": "mumps", #gmres
"maximum_iterations": 20,
"report": True,
"error_on_nonconvergence": True}}
problem = NonlinearVariationalProblem(F, up, bc, J)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
xdmffile_u = XDMFFile('Coupled/velocity.xdmf')
xdmffile_p = XDMFFile('Coupled/pressure.xdmf')
(u, p) = up.split()
K = int(T/dt)
for i in range(1, K):
# Compute the current time
t = i*dt
print("t= ",t)
# Solve the nonlinear problem
solver.solve()
# Store the solution in up_prev
assign(up_prev, up)
# Plot
(u, p) = up.split()
# Save solution to file (XDMF/HDF5)
xdmffile_u.write(u, t)
xdmffile_p.write(p, t)
Thanks in advance,
Jack