How to use solution of linear problem as initial guess of non linear problem

General
1

I followed the tutorial Stokes equations with Taylor-Hood elements in order to obtain solutions of the fields u and p.
I want now to introduce a non linear term in the Stokes equations, namely - (u \cdot \nabla) u , so I need a non linear solver. I tried to use dolfinx.nls.petsc.NewtonSolver in the fashion of Cahn-Hilliard equation , using the fields obtained previously as first guess.

I paste here the mwc, assuming u and p are the solution obtained with Stokes equations with Taylor-Hood elements

# get solution
u_initial=u
p_initial=p

up = Function(W)

u_size=len(u_initial.x.array)
p_size=len(p_initial.x.array)

# get previous solution and use it as initial condition
up.x.array[:u_size]=u_initial.x.array
up.x.array[u_size:(p_size+u_size)]=p_initial.x.array
up.x.scatter_forward()

uu, pp = split(up)

v, q = TestFunctions(W)

# define form
F= (1/Re)*inner(grad(uu), grad(v)) * dx \
    -inner(pp, div(v)) * dx -inner(div(uu), q) * dx \
    -inner(dot(uu, nabla_grad(uu)),v) * dx \
    -delta*inner(grad(pp),grad(q)) * dx 

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, up, bcs)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-3
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

r = solver.solve(up)

The function space W is the same I have used In previous part of the code to solve the Stokes equations

P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2*P1
W = FunctionSpace(mesh, TH)

I get RuntimeError: Newton solver did not converge because maximum number of iterations reached. I assume it should converge quickly since the first guess should not be too far away from the solution, as I’m trying to solve the same problem but with an additional term. What am I doing wrong?

2

Can you provide a minimal working example? E.g. it’s crucial to know what value you’re employing for Re and the initial approximation you’re computing for u. It may be worth considering a continuation method rather than simply using the solution of a different problem as an initial guess.

3

I can provide the rest of the code, but unfortunately not the original mesh. It is approximately a L-shaped 2D pipe with a chamber at the elbow, the coordinates range in [0,1]x[0,1]. The borders are smooth. The bc are parabolic inlet velocity, noslip on the sides and known pressure in the outlet. The solution I obtain with the following code has physical sense and it is robust for reynolds ranging from 10^2 to 10^4.

Here is the rest of the code, which comes before the one I have posted:

import dolfinx.fem.petsc
import numpy as np
import pyvista
from dolfinx import fem
import dolfinx.plot as dplt
from dolfinx.fem import (Constant,  VectorFunctionSpace, Function, FunctionSpace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological,extract_function_spaces)
from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, create_rectangle, CellType, Mesh
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc
from ufl import (VectorElement,FiniteElement,FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, TrialFunctions, TestFunctions,
                 div, dot, dx, grad, inner, lhs, rhs, MixedElement, nabla_grad, split)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
import matplotlib.pyplot as plt
from dolfinx import cpp as _cpp
from dolfinx.nls.petsc import NewtonSolver

###########################################
#options
reee=3000.0 #reynold value
c=1.0 #max velocity on inlet
###########################################

# coordinates of inlet
a=0.43902439024390244 # x1, y1=0
b=0.5365853658536586 # x2, y2=0
xm=(a+b)*0.5
# parabola parameters
def calc_parabola(x1, y1, x2, y2, x3, y3):
	denom = (x1-x2) * (x1-x3) * (x2-x3)
	A     = (x3 * (y2-y1) + x2 * (y1-y3) + x1 * (y3-y2)) / denom
	B     = (x3*x3 * (y1-y2) + x2*x2 * (y3-y1) + x1*x1 * (y2-y3)) / denom
	C     = (x2 * x3 * (x2-x3) * y1+x3 * x1 * (x3-x1) * y2+x1 * x2 * (x1-x2) * y3) / denom
	return A,B,C
    
alpha,beta,om=calc_parabola(a,0,b,0,xm,c)
# inlet velocity
def vel_in(x):
    val = alpha*x[0]*x[0] + beta*x[0] + om
    return np.stack((np.zeros(x.shape[1]), -val*np.ones(x.shape[1])))

with XDMFFile(MPI.COMM_WORLD, "mesh_domain.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, "mesh_boundaries.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")

ds = Measure("ds", domain=mesh, subdomain_data=ft)
fdim = mesh.topology.dim - 1

###### BCs
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2*P1
W = FunctionSpace(mesh, TH)

V=W.sub(0).collapse()[0]
Q=W.sub(1).collapse()[0]

# parabolic inlet
u_in = Function(V)
u_in.interpolate(vel_in)
where=ft.indices[ft.values==3]
dofs = locate_dofs_topological(V, fdim, where)
bc0 = dirichletbc(u_in, dofs)

# noslip
noslip = Constant(mesh,(0.0,0.0))
where=ft.indices[ft.values==1]
dofs = locate_dofs_topological(V, fdim, where)
bc1 = dirichletbc(noslip, dofs, V)

# noslip
noslip = Constant(mesh,(0.0,0.0))
where=ft.indices[ft.values==2]
dofs = locate_dofs_topological(V, fdim, where)
bc2 = dirichletbc(noslip, dofs, V)

# pressure=0 out
pres=Constant(mesh,(0.0))
where=ft.indices[ft.values==4]
dofs = locate_dofs_topological(Q, fdim, where)
bc3 = dirichletbc(pres, dofs, Q)

bcs = [bc0,bc1,bc2,bc3]
########

######## Variational problem
(u, p) = TrialFunction(V), TrialFunction(Q)
(v, q) = TestFunction(V), TestFunction(Q)
f = Constant(mesh, (ScalarType(0), ScalarType(0)))
Re = Constant(mesh,ScalarType( reee ))
delta=1
a = form([[(1/Re)*inner(grad(u), grad(v)) * dx, -inner(p, div(v)) * dx],
          [-inner(div(u), q) * dx, -delta*inner(grad(p),grad(q)) * dx ]])
L = form([inner(f, v) * dx, delta*inner(f, grad(q)) * dx])
########

######## solver
a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
       [None, a_p11]]
A = fem.petsc.assemble_matrix_nest(a, bcs=bcs)
A.assemble()
P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()
b = fem.petsc.assemble_vector_nest(L)
fem.petsc.apply_lifting_nest(b, a, bcs=bcs)

for b_sub in b.getNestSubVecs():
    b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs)
fem.petsc.set_bc_nest(b, bcs0)
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(
    ("u", nested_IS[0][0]),
    ("p", nested_IS[0][1]))
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
ksp.setFromOptions()

u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([_cpp.la.petsc.create_vector_wrap(u.x), _cpp.la.petsc.create_vector_wrap(p.x)])
ksp.solve(b, x)

########

The code I’ve already posted goes after this. I hope this is enough.
ps: I will try to perform the same but with a simpler mesh that I can share if it’s necessary

4

I confirm that the newton solver do not converge even with a much simpler domain [0,1]x[0,5], with inlet={(x,y),y=5}, sides={(x,y),x=0 or x=1}, outlet={(x,y),y=0}

5

Considering we don’t have the mesh with corresponding facet tags, it’s very difficult to help.

My first avenues for debugging would be the following:

  • See if you can solve the problem with a very small Reynolds number \mathrm{Re} \leq 1. Ideally see if you can solve the problem with \mathrm{Re} = 0 but this is not compatible with your formulation.
  • If you can solve the problem with a small Reynolds number, use a continuation scheme to get up to your target of \mathrm{Re} = 3000 assuming that your domain’s geometry has appropriate length scale.
  • If you cannot solve the case where \mathrm{Re} = 0 try a direct solver rather than iterative scheme.
1 Like
6

Thank you very much. I will try that.
Where can I find a reference for the implementation of the continuation scheme?
Could you just confirm that the following actually initialize the first guess for the newton solver as I intend?

# get solution
u_initial=u
p_initial=p

up = Function(W)

u_size=len(u_initial.x.array)
p_size=len(p_initial.x.array)

# get previous solution and use it as initial condition
up.x.array[:u_size]=u_initial.x.array
up.x.array[u_size:(p_size+u_size)]=p_initial.x.array
up.x.scatter_forward()

uu, pp = split(up)

v, q = TestFunctions(W)

# define form
F= (1/Re)*inner(grad(uu), grad(v)) * dx \
    -inner(pp, div(v)) * dx -inner(div(uu), q) * dx \
    -inner(dot(uu, nabla_grad(uu)),v) * dx \
    -delta*inner(grad(pp),grad(q)) * dx 

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, up, bcs)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-3
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

r = solver.solve(up)
7

A simple example:

...
Re = dolfinx.fem.Constant(mesh, 0.0)
...
Re_vals = [0.0, 10.0, 100.0, ..., target_val]
for Re_val in Re_vals:
    Re.value = Re_val
    ...
    # solve problem
    ...

Since I can’t run your code I can’t test it. But it looks close enough. Consult the Stokes demo for examples of block formulations.

8

I have written a minimal working example with a simpler domain, which should produce a solution, but the newton solver do not converge. If you have time, please take a look:

import dolfinx.fem.petsc

import numpy as np
import pyvista

from dolfinx import fem

import dolfinx.plot as dplt

from dolfinx.fem import (Constant,  VectorFunctionSpace, Function, FunctionSpace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological,extract_function_spaces)

from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, create_rectangle, CellType, Mesh
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc

from ufl import (VectorElement,FiniteElement,FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, TrialFunctions, TestFunctions,
                 div, dot, dx, grad, inner, lhs, rhs, MixedElement, nabla_grad, split)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh

import matplotlib.pyplot as plt

from dolfinx import cpp as _cpp

from dolfinx.nls.petsc import NewtonSolver

################# options #####################
save_stokes_results=True
reee=100.0 #reynolds
c=1.0 # inlet max velocity
###############################################


a=0.0 
b=1.0 
xm=(a+b)*0.5

def calc_parabola(x1, y1, x2, y2, x3, y3):
	denom = (x1-x2) * (x1-x3) * (x2-x3)
	A     = (x3 * (y2-y1) + x2 * (y1-y3) + x1 * (y3-y2)) / denom
	B     = (x3*x3 * (y1-y2) + x2*x2 * (y3-y1) + x1*x1 * (y2-y3)) / denom
	C     = (x2 * x3 * (x2-x3) * y1+x3 * x1 * (x3-x1) * y2+x1 * x2 * (x1-x2) * y3) / denom
	return A,B,C
    
alpha,beta,om=calc_parabola(a,0,b,0,xm,c)

def vel_in(x):
    val = alpha*x[0]*x[0] + beta*x[0] + om
    return np.stack((np.zeros(x.shape[1]), -val*np.ones(x.shape[1])))

mesh = create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (1.0, 5.0)), n=(50, 200),
                            cell_type=CellType.triangle)

ds = Measure("ds", domain=mesh)

fdim = mesh.topology.dim - 1

P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2*P1
W = FunctionSpace(mesh, TH)

V=W.sub(0).collapse()[0]
Q=W.sub(1).collapse()[0]

# parabolic inlet
u_in = Function(V)
u_in.interpolate(vel_in)
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[1], 5.0))
dofs = locate_dofs_topological(V, fdim, where)
bc0 = dirichletbc(u_in, dofs)

# noslip
noslip = Constant(mesh,(0.0,0.0))
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[0], 0.0))
dofs = locate_dofs_topological(V, fdim, where)
bc1 = dirichletbc(noslip, dofs, V)

# noslip
noslip = Constant(mesh,(0.0,0.0))
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[0], 1.0))
dofs = locate_dofs_topological(V, fdim, where)
bc2 = dirichletbc(noslip, dofs, V)

# pressure=0 out
pres=Constant(mesh,(0.0))
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[1], 0.0))
dofs = locate_dofs_topological(Q, fdim, where)
bc3 = dirichletbc(pres, dofs, Q)

bcs = [bc0,bc1,bc2,bc3]

# Define variational problem
(u, p) = TrialFunction(V), TrialFunction(Q)
(v, q) = TestFunction(V), TestFunction(Q)
f = Constant(mesh, (ScalarType(0), ScalarType(0)))

Re = Constant(mesh,ScalarType( reee ))
delta=1

a = form([[(1/Re)*inner(grad(u), grad(v)) * dx, -inner(p, div(v)) * dx],
          [-inner(div(u), q) * dx, -delta*inner(grad(p),grad(q)) * dx ]])


L = form([inner(f, v) * dx, delta*inner(f, grad(q)) * dx])

a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
       [None, a_p11]]

# nested

A = fem.petsc.assemble_matrix_nest(a, bcs=bcs)
A.assemble()

P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()

b = fem.petsc.assemble_vector_nest(L)

# Modify ('lift') the RHS for Dirichlet boundary conditions
fem.petsc.apply_lifting_nest(b, a, bcs=bcs)

# Sum contributions from ghost entries on the owner
for b_sub in b.getNestSubVecs():
    b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

# Set Dirichlet boundary condition values in the RHS
bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs)
fem.petsc.set_bc_nest(b, bcs0)

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

# Define the matrix blocks in the preconditioner with the velocity and
# pressure matrix index sets
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(
    ("u", nested_IS[0][0]),
    ("p", nested_IS[0][1]))

# Set the preconditioners for each block
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")

# Monitor the convergence of the KSP
ksp.setFromOptions()

u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([_cpp.la.petsc.create_vector_wrap(u.x), _cpp.la.petsc.create_vector_wrap(p.x)])
ksp.solve(b, x)

norm_u_0 = u.x.norm()
norm_p_0 = p.x.norm()
if MPI.COMM_WORLD.rank == 0:
    print("(A) Norm of velocity coefficient vector (nested, iterative): {}".format(norm_u_0))
    print("(A) Norm of pressure coefficient vector (nested, iterative): {}".format(norm_p_0))
    
## Write the solution to file

if save_stokes_results:
    with XDMFFile(MPI.COMM_WORLD, "out_stokes/velocity.xdmf", "w") as ufile_xdmf:
        u.x.scatter_forward()
        ufile_xdmf.write_mesh(mesh)
        ufile_xdmf.write_function(u)
    
    with XDMFFile(MPI.COMM_WORLD, "out_stokes/pressure.xdmf", "w") as pfile_xdmf:
        p.x.scatter_forward()
        pfile_xdmf.write_mesh(mesh)
        pfile_xdmf.write_function(p)

########### non linear problem ##########

# get previous solution
u_initial=u
p_initial=p

up = Function(W)

u_size=len(u_initial.x.array)
p_size=len(p_initial.x.array)

# use previous as initial condition
up.x.array[:u_size]=u_initial.x.array
up.x.array[u_size:(p_size+u_size)]=p_initial.x.array

up.x.scatter_forward()

uu, pp = split(up)

v, q = TestFunctions(W)

# define the new problem
F= (1/Re)*inner(grad(uu), grad(v)) * dx \
    -inner(pp, div(v)) * dx -inner(div(uu), q) * dx \
    -inner(dot(uu, nabla_grad(uu)),v) * dx \
    -delta*inner(grad(pp),grad(q)) * dx 

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, up, bcs)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5
solver.report = True
solver.max_it = 1000

# helper function to print norm
def normal2(arr):
    val=0
    for a in arr:
        val+=a**2
    return np.sqrt(val)

r = solver.solve(up)

print(f"Num iterations: {r[0]}")

curr_u_array=up.x.array[:u_size]
curr_p_array=up.x.array[u_size:(u_size+p_size)]

norm_u_0 = normal2(curr_u_array)
norm_p_0 = normal2(curr_p_array)

# check if norm of solution is comparable to the one of the stokes problem
if MPI.COMM_WORLD.rank == 0:
    print("(A) Norm2 of velocity coefficient vector (nested, iterative): {}".format(norm_u_0))
    print("(A) Norm2 of pressure coefficient vector (nested, iterative): {}".format(norm_p_0))