Newton Solver did not Converge or Report errors Fenicsx

dolfinx
1

Hi,

I am trying to move fenics code to fenicsx. However I am getting a “Newton solver will not converge because maximum number of iterations reached” error.
I’m using wsl on windows 11. I’ve tried using both docker and intalling fenicsx in ubuntu through wsl.
In both cases I can’t see what is happening at each iteration, as in the errors at each step are not being printed, so I am not sure why the code is not converging.

I am trying to solve the stokes equations in cyindrical coordinates with a free boundary where I have performed a coordinate sytem transformation to fix the free boundary.
Below is the code. Any help would be appreciated

import numpy as np
from dolfinx.mesh import locate_entities, meshtags, locate_entities_boundary, create_unit_square, CellType
from dolfinx.fem.petsc import NonlinearProblem 
from dolfinx import nls, log
from dolfinx.fem import  ( Function, VectorFunctionSpace, FunctionSpace, dirichletbc, locate_dofs_geometrical, 
                          locate_dofs_topological )
from petsc4py import PETSc
from ufl import (FiniteElement, VectorElement, TestFunctions, dx, 
                 MixedElement, split, Cell, Mesh, VectorElement, Measure, SpatialCoordinate)
from mpi4py import MPI 

ep = 0.1 # aspect ratio

msh = create_unit_square(MPI.COMM_WORLD, 50,50,  CellType.triangle)

# defining Function space, Trial & Test functions
v_cg2 = VectorElement("Lagrange", msh.ufl_cell(), 2) # velocities
q_cg1 = FiniteElement("Lagrange", msh.ufl_cell(), 1) # pressure
s_cg2 = FiniteElement("Lagrange", msh.ufl_cell(), 2) # free boundary

W = FunctionSpace(msh, MixedElement([v_cg2, q_cg1, s_cg2]))

V = FunctionSpace(msh, v_cg2)
Q = FunctionSpace(msh, q_cg1)
S = FunctionSpace(msh, s_cg2)

w = Function(W)
(u,p,r) = w.split()
(v,q,s) = TestFunctions(W)
# Boundary Conditions
x = SpatialCoordinate(msh)
fdim = msh.topology.dim -1

D = 30 # draw ratio 
# inlet 
def inlet(x):
    return np.isclose(x[1],0)

R0, _ = W.sub(2).collapse()
inlet_radius = Function(R0)

def inlet_radius_expression(x):
    return np.stack((np.ones(x.shape[1], dtype=PETSc.ScalarType)))

inlet_radius.interpolate(inlet_radius_expression)
inlet_dof_R = locate_dofs_geometrical((W.sub(2),S), inlet)
bcr_inflow = dirichletbc(inlet_radius, inlet_dof_R, W.sub(2))

# velocities at inlet
W0, _ = W.sub(0).collapse()

def inlet_velocity_expression(x): # setting constants as PETSc.Scalar type is meant to help compilation
    return np.stack(( np.zeros(x.shape[1]), np.ones(x.shape[1])) ).astype(PETSc.ScalarType)

inlet_velocity = Function(W0)
inlet_velocity.interpolate(inlet_velocity_expression)
inlet_facets = locate_entities_boundary(msh, fdim, inlet)
inlet_dof_W = locate_dofs_topological((W.sub(0),V), fdim, inlet_facets)
bcu_inflow = dirichletbc(inlet_velocity, inlet_dof_W, W.sub(0))

# outlet
def outlet(x):
    return np.isclose(x[1],1)

def outlet_velocity_expression(x):
    return np.stack(( np.zeros(x.shape[1]),  D*np.ones(x.shape[1]))).astype(PETSc.ScalarType)

outlet_velocity = Function(W0)
outlet_velocity.interpolate(outlet_velocity_expression)
outlet_facets = locate_entities_boundary(msh, fdim, outlet)
outlet_dof_W = locate_dofs_topological((W.sub(0),V), fdim, outlet_facets)
bcu_outflow = dirichletbc(outlet_velocity, outlet_dof_W, W.sub(0))    

bcs = [bcu_outflow, bcu_inflow, bcr_inflow] 
# need to integrate across free boundary to apply no Flux
facet_indicies = locate_entities(msh,fdim, lambda x: np.isclose(x[0],1))
                                          # ds(1) denotes free boundary
facet_markers = np.full_like(facet_indicies,1)

facet_indicies = np.hstack(facet_indicies).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)

sorted_facets = np.argsort(facet_indicies)
facet_tag = meshtags(msh, fdim, facet_indicies[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=msh, subdomain_data=facet_tag)

# Giving initial guesses using asyptotics for w and r
r.interpolate(lambda x:  0.5*np.exp(0*x[1]) ) 
u.sub(1).interpolate((lambda x:  np.exp(np.log(D)*x[1]) ) )
# Weak Formulation
ep = PETSc.ScalarType(ep)

sig_r = - p + 2/r * u[0].dx(0)
sig_z = - p + 2*u[1].dx(1)  - 2*x[0]/r*r.dx(1)*u[1].dx(0)
sig_rz = ep*(u[0].dx(1) - x[0]*r.dx(1)/r * u[0].dx(0)) + 1/ep *(1/r * u[1].dx(0))

rm = (- v[0].dx(0)*r*x[0]*sig_r                   \
      + ep*sig_rz*(-x[0]*(r**2 * v[0]).dx(1)   \
      + r.dx(1)*r*(x[0]**2 * v[0]).dx(0) )        \
     - r*v[0]*(-p + 2*u[0]/(r*x[0])) )*dx                      

zm = (- v[1].dx(0)*r*x[0]*sig_rz                  \
      + ep*sig_z*( -x[0]*(r**2 * v[1]).dx(1)    \
      + r.dx(1)*r*(x[0]**2*v[1]).dx(0)) )*dx   

ce = ( q*(r*(u[0]*x[0]).dx(0) + r**2*x[0]*u[1].dx(1) - x[0]**2*r*r.dx(1)*u[1].dx(0) ) ) *dx 

fb = (s*(r.dx(0))*r**2*x[0] )*dx + ( s*(u[0] - u[1]*r.dx(1))*r**2*x[0]  )*ds(1)

F = rm + zm + ce + fb

problem = NonlinearProblem(F,w, bcs)

solver = nls.petsc.NewtonSolver(msh.comm, problem)
solver.convergence_criterion = "incremental"
solver.atol = 1e-8
solver.rtol = 1e-8
solver.report = True

solver.solve(w)
2

If you add

log.set_log_level(log.LogLevel.INFO)

prior to solving, you will get the information of the NewtonSolver:

2023-02-07 07:19:43.612 (   0.471s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.618 (   0.477s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.624 (   0.483s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.624 (   0.483s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.630 (   0.489s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 3: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.631 (   0.490s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.636 (   0.495s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 4: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.637 (   0.496s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.644 (   0.503s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 5: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.644 (   0.503s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.650 (   0.509s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 6: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.651 (   0.510s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.656 (   0.515s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 7: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.657 (   0.516s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.663 (   0.522s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 8: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.664 (   0.523s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.669 (   0.528s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 9: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.670 (   0.529s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.676 (   0.535s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 10: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.676 (   0.535s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.682 (   0.541s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 11: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.683 (   0.542s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.688 (   0.547s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 12: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.689 (   0.548s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.694 (   0.553s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 13: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.695 (   0.554s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.701 (   0.560s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 14: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.701 (   0.560s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.707 (   0.566s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 15: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.708 (   0.567s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.713 (   0.572s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 16: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.714 (   0.573s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.720 (   0.579s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 17: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.720 (   0.579s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.726 (   0.585s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 18: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.727 (   0.586s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.732 (   0.591s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 19: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.733 (   0.592s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.738 (   0.597s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 20: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.739 (   0.598s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.745 (   0.604s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 21: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.745 (   0.604s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.751 (   0.610s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 22: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.752 (   0.611s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.757 (   0.616s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 23: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.758 (   0.617s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.763 (   0.623s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 24: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.764 (   0.623s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.770 (   0.629s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 25: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.771 (   0.630s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.776 (   0.635s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 26: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.777 (   0.636s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.783 (   0.642s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 27: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.783 (   0.642s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.789 (   0.648s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 28: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.789 (   0.648s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.795 (   0.654s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 29: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.796 (   0.655s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.801 (   0.660s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 30: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.802 (   0.661s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.808 (   0.667s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 31: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.808 (   0.667s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.814 (   0.673s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 32: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.815 (   0.674s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.820 (   0.679s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 33: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.821 (   0.680s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.826 (   0.685s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 34: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.827 (   0.686s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.833 (   0.692s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 35: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.833 (   0.692s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.839 (   0.698s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 36: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.840 (   0.699s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.845 (   0.704s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 37: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.846 (   0.705s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.851 (   0.710s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 38: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.852 (   0.711s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.858 (   0.717s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 39: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.858 (   0.717s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.864 (   0.723s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 40: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.865 (   0.724s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.870 (   0.729s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 41: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.871 (   0.730s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.877 (   0.736s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 42: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.877 (   0.736s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.883 (   0.742s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 43: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.884 (   0.743s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.889 (   0.748s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 44: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.890 (   0.749s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.895 (   0.755s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 45: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.896 (   0.755s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.902 (   0.761s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 46: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.903 (   0.762s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.908 (   0.767s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 47: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.909 (   0.768s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.914 (   0.773s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 48: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.915 (   0.774s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.920 (   0.780s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 49: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)
2023-02-07 07:19:43.921 (   0.780s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-02-07 07:19:43.927 (   0.786s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 50: r (abs) = inf (tol = 1e-08) r (rel) = -nan(tol = 1e-08)

which indicates that something is wrong with your variational form.
In particular

should be
(u,p,r) = ufl.split(w)
Following is a code that fixes this, + some interpolation issues:

import ufl
import numpy as np
from dolfinx import log, nls
from dolfinx.fem import (Function, FunctionSpace, VectorFunctionSpace,
                         dirichletbc, locate_dofs_geometrical,
                         locate_dofs_topological)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.mesh import (CellType, create_unit_square, locate_entities,
                          locate_entities_boundary, meshtags)
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (Cell, FiniteElement, Measure, Mesh, MixedElement,
                 SpatialCoordinate, TestFunctions, VectorElement, dx, split)

ep = 0.1  # aspect ratio

msh = create_unit_square(MPI.COMM_WORLD, 50, 50,  CellType.triangle)

# defining Function space, Trial & Test functions
v_cg2 = VectorElement("Lagrange", msh.ufl_cell(), 2)  # velocities
q_cg1 = FiniteElement("Lagrange", msh.ufl_cell(), 1)  # pressure
s_cg2 = FiniteElement("Lagrange", msh.ufl_cell(), 2)  # free boundary

W = FunctionSpace(msh, MixedElement([v_cg2, q_cg1, s_cg2]))

V = FunctionSpace(msh, v_cg2)
Q = FunctionSpace(msh, q_cg1)
S = FunctionSpace(msh, s_cg2)
w = Function(W)
(u, p, r) = ufl.split(w)
(v, q, s) = TestFunctions(W)
# Boundary Conditions
x = SpatialCoordinate(msh)
fdim = msh.topology.dim - 1

D = 30  # draw ratio
# inlet


def inlet(x):
    return np.isclose(x[1], 0)


R0, _ = W.sub(2).collapse()
inlet_radius = Function(R0)


def inlet_radius_expression(x):
    return np.stack((np.ones(x.shape[1], dtype=PETSc.ScalarType)))


inlet_radius.interpolate(inlet_radius_expression)
inlet_dof_R = locate_dofs_geometrical((W.sub(2), S), inlet)
bcr_inflow = dirichletbc(inlet_radius, inlet_dof_R, W.sub(2))

# velocities at inlet
W0, _ = W.sub(0).collapse()


# setting constants as PETSc.Scalar type is meant to help compilation
def inlet_velocity_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.ones(x.shape[1]))).astype(PETSc.ScalarType)


inlet_velocity = Function(W0)
inlet_velocity.interpolate(inlet_velocity_expression)
inlet_facets = locate_entities_boundary(msh, fdim, inlet)
inlet_dof_W = locate_dofs_topological((W.sub(0), V), fdim, inlet_facets)
bcu_inflow = dirichletbc(inlet_velocity, inlet_dof_W, W.sub(0))

# outlet


def outlet(x):
    return np.isclose(x[1], 1)


def outlet_velocity_expression(x):
    return np.stack((np.zeros(x.shape[1]),  D*np.ones(x.shape[1]))).astype(PETSc.ScalarType)


outlet_velocity = Function(W0)
outlet_velocity.interpolate(outlet_velocity_expression)
outlet_facets = locate_entities_boundary(msh, fdim, outlet)
outlet_dof_W = locate_dofs_topological((W.sub(0), V), fdim, outlet_facets)
bcu_outflow = dirichletbc(outlet_velocity, outlet_dof_W, W.sub(0))

bcs = [bcu_outflow, bcu_inflow, bcr_inflow]
# need to integrate across free boundary to apply no Flux
facet_indicies = locate_entities(msh, fdim, lambda x: np.isclose(x[0], 1))
# ds(1) denotes free boundary
facet_markers = np.full_like(facet_indicies, 1)

facet_indicies = np.hstack(facet_indicies).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)

sorted_facets = np.argsort(facet_indicies)
facet_tag = meshtags(
    msh, fdim, facet_indicies[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=msh, subdomain_data=facet_tag)

# Giving initial guesses using asyptotics for w and r
w.sub(2).interpolate(lambda x:  0.5*np.exp(0*x[1]))
w.sub(0).sub(1).interpolate((lambda x:  np.exp(np.log(D)*x[1])))
# Weak Formulation
ep = PETSc.ScalarType(ep)

sig_r = - p + 2/r * u[0].dx(0)
sig_z = - p + 2*u[1].dx(1) - 2*x[0]/r*r.dx(1)*u[1].dx(0)
sig_rz = ep*(u[0].dx(1) - x[0]*r.dx(1)/r * u[0].dx(0)) + \
    1/ep * (1/r * u[1].dx(0))

rm = (- v[0].dx(0)*r*x[0]*sig_r
      + ep*sig_rz*(-x[0]*(r**2 * v[0]).dx(1)
      + r.dx(1)*r*(x[0]**2 * v[0]).dx(0))
      - r*v[0]*(-p + 2*u[0]/(r*x[0])))*dx

zm = (- v[1].dx(0)*r*x[0]*sig_rz
      + ep*sig_z*(-x[0]*(r**2 * v[1]).dx(1)
      + r.dx(1)*r*(x[0]**2*v[1]).dx(0)))*dx

ce = (q*(r*(u[0]*x[0]).dx(0) + r**2*x[0] *
      u[1].dx(1) - x[0]**2*r*r.dx(1)*u[1].dx(0))) * dx

fb = (s*(r.dx(0))*r**2*x[0])*dx + (s*(u[0] - u[1]*r.dx(1))*r**2*x[0])*ds(1)

F = rm + zm + ce + fb
problem = NonlinearProblem(F, w, bcs)
log.set_log_level(log.LogLevel.INFO)
solver = nls.petsc.NewtonSolver(msh.comm, problem)
solver.convergence_criterion = "residual"
solver.atol = 1e-8
solver.rtol = 1e-8
solver.report = True

solver.solve(w)

The code does not converge due to your struct tolerance, but it gets down to:

2023-02-07 07:22:44.548 (  16.464s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 50: r (abs) = 4.60795e-05 (tol = 1e-08) r (rel) = 7.40067e-08(tol = 1e-08)
1 Like
3

Thanks so much for you’re help @dokken.