Hello,
In order to stabilize a newton solver for Naviers stokes in fenicsx I would like to use the TS solver from PETSc to do pseudo time stepping. This would allow me to improve the convergence of the Newton. To do this, I’ve tried to create a TSProblem class that calculates the Jacobian and the residual vector. If I’ve understood the petsc doc correctly, J and F must not contain a time term, which will be added implicitly.
However, although the calculation proceeds correctly and time increases, the vector obtained at the end contains zero velocity everywhere except at the inlet (dirichlet bc). This problem seems stupid to me, but I must be missing something.
Here’s some functional code:
from mpi4py import MPI
import dolfinx
import os
from pathlib import Path
import time
from ufl import split, TestFunctions, TestFunction, derivative, adjoint, action, Identity, grad, transpose, inner, dx, replace, div, Measure, TrialFunction, sqrt, dot
import ufl
import basix.ufl
from petsc4py import PETSc
from petsc4py.PETSc import Options
from dolfinx.fem.petsc import NonlinearProblem, apply_lifting, assemble_vector, set_bc, assemble_matrix, create_matrix, create_vector
from dolfinx.fem import Function, functionspace, Constant, dirichletbc, locate_dofs_topological, locate_dofs_geometrical, form, assemble_scalar
from dolfinx.mesh import create_box, create_rectangle, locate_entities_boundary, CellType
from dolfinx import mesh
import numpy as np
from dolfinx.io import XDMFFile, VTXWriter
from dolfinx.io.gmshio import model_to_mesh
import scipy
import gmsh
from dolfinx.la import norm
comm = MPI.COMM_WORLD
gmsh.initialize()
tags = {}
tags["inlet"] = 1
tags["outlet"] = 2
tags["wall"] = 3
model = gmsh.model
model_rank = 0
comm = MPI.COMM_WORLD
if comm.rank == model_rank:
factory = model.occ
s1 = factory.addRectangle(-1, 0, 0, 1, 1)
s2 = factory.addRectangle(0, -1, 0, 7.5, 2)
s3 = 3
factory.fuse([(2, s1)], [(2, s2)], tag=s3)
factory.synchronize()
ov = model.getBoundary([(2, s3)])
l1, l2, l3, l4, l5, l6 = (val[1] for val in ov) # counterclockwise (l6 = inlet)
# Tag interior
model.addPhysicalGroup(2, [s3], tag=0)
# Tag boundaries
model.addPhysicalGroup(1, [l6], tag=tags["inlet"])
model.setPhysicalName(1, tags["inlet"], "inlet")
model.addPhysicalGroup(1, [l4], tag=tags["outlet"])
model.setPhysicalName(1, tags["outlet"], "outlet")
model.addPhysicalGroup(1, [l1, l2, l3, l5], tag=tags["wall"])
model.setPhysicalName(1, tags["wall"], "wall")
# Set uniform mesh size at all points
mesh_size = 0.05
model.mesh.setSize(model.getEntities(0), mesh_size)
# Generate 2D mesh
model.mesh.generate(dim=2)
# Convert Gmsh mesh to DOLFINx
mesh, _, facet_tags = model_to_mesh(model, MPI.COMM_WORLD, rank=model_rank, gdim=2)
facet_tags.name = "facet_tags"
gmsh.finalize()
metadata = {"quadrature_degree": 4}
dx = Measure("dx", domain=mesh, metadata=metadata)
Re = 100
v_max_inlet = Constant(mesh, PETSc.ScalarType(1))
rho_ = Constant(mesh, PETSc.ScalarType(1))
mu_ = Constant(mesh, PETSc.ScalarType(rho_*v_max_inlet*2/Re))
nu_ = Constant(mesh, PETSc.ScalarType(mu_/rho_))
V_element = basix.ufl.element("CG", mesh.basix_cell(), 2, shape = (mesh.topology.dim,))
P_element = basix.ufl.element("CG", mesh.basix_cell(), 1)
U_element = basix.ufl.mixed_element([V_element,P_element])
U = functionspace(mesh, U_element); U0 = U.sub(0);P0 = U.sub(1)
w_v, w_p = TestFunctions(U)
uh = Function(U)
v, p = split(uh)
uh_prev = Function(U)
v_prev, p_prev = split(uh_prev)
V_v = U.sub(0).collapse()[0]
def v_inflow_eval(x):
values = np.zeros((2, x.shape[1]))
values[0] = 4.0 * x[1] * (1.0 - x[1])
values[1] = np.zeros(x.shape[1])
return values
v_wall = Function(V_v)
v_wall.x.array[:] = 0.0
v_inflow = Function(V_v)
v_inflow.interpolate(v_inflow_eval)
fdim = mesh.topology.dim - 1
wall_dofs_v = locate_dofs_topological((U0,V_v), fdim, facet_tags.find(tags["wall"]))
inlet_dofs_v = locate_dofs_topological((U0,V_v), fdim, facet_tags.find(tags["inlet"]))
bcs = [dirichletbc(v_inflow, inlet_dofs_v, U0), dirichletbc(v_wall, wall_dofs_v, U0)]
F_stdy = (-div(v) * w_p * dx + nu_ * inner(grad(v),grad(w_v)) * dx + inner(grad(v) * v, w_v)*dx - p * div(w_v) * dx)
F_PV = F_stdy
class TSProblem:
def __init__(self, F, u, bc):
V = u.function_space
du = TrialFunction(V)
self.L = form(F)
self.a = form((derivative(F, u, du)))
self.bc = bc
self._F, self._J = None, None
self.u = u
def F(self, ts, t, x, F):
"""Assemble residual vector."""
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.x.petsc_vec)
self.u.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
with F.localForm() as f_local:
f_local.set(0.0)
assemble_vector(F, self.L)
apply_lifting(F, [self.a], bcs=[self.bc], x0=[x], alpha=-1.0)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(F, self.bc, x, -1.0)
def J(self, ts, t, x, J, P):
"""Assemble Jacobian matrix."""
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.x.petsc_vec)
self.u.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
J.zeroEntries()
assemble_matrix(J, self.a, bcs=self.bc)
J.assemble()
problem = TSProblem(F_PV, uh, bcs)
J = create_matrix(problem.a)
b = create_vector(problem.L)
# Create TS Newton solver
ts = PETSc.TS().create()
ts.setRHSFunction(problem.F, b)
ts.setRHSJacobian(problem.J, J)
ts.setType('pseudo')
problem_prefix = "ns_"
ts.setOptionsPrefix(problem_prefix)
opts = PETSc.Options()
opts.prefixPush(problem_prefix)
opts["ts_pseudo_increment"] = 1.2
opts["ts_monitor"] = None
opts['snes_monitor'] = None
opts["snes_type"] = "newtonls"
opts['snes_linesearch_type'] = 'bt'
opts['snes_linesearch_order']= 2
opts["ksp_converged_reason"] = None
opts["ksp_rtol"] = 1e-10
#opts["ksp_monitor"] = None
opts.prefixPop() # ns_
ts.setFromOptions()
snes = ts.getSNES()
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
x = uh.x.petsc_vec.copy()
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
ts.setTime(0.0)
ts.setTimeStep(0.05)
ts.setMaxSteps(20)
ts.setMaxTime(1e10)
ts.solve(x)
x.copy(uh.x.petsc_vec)
uh.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
uh.x.scatter_forward()
Thanks a lot !
