Set Initial Condition to Nonlinear Problem

General
1

I am trying to solve lid driven stokes flow for low Reynolds numbers (Re <400), and to speed up computational time and to provide better convergence, I first solve stokes flow, and I am going to use the solution from the linear stokes as the initial guess/condition for the nonlinear incompressible Naiver-Stokes flow.

I do not know how to properly set the initial condition to the problem = NonlinearProblem(a, w, bcs=bcs, J=dF). It is my understanding that FEniCSx uses all zeros as the initial condition, and this causes the solver to diverge. I am trying to be similar to the “Cahn-Hilliard equation demo”, and I use w = U.interpolate like in the demo, where w is the function space, and U is the mixed solution from the Stokes flow.

What is the correct way to set an initial condition to the NonLinearProblem function in FEniCSx, using a previous solution as the initial condition?

I am using FEniCSx 0.9, running on WSL, on Ubuntu. Also, I have already validated the Stokes flow solver, and I am confident it works properly, and I am using Taylor-Hood elements to avoid stabilization issues.

Below is my code

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import log
from dolfinx.fem import Function, dirichletbc, functionspace, locate_dofs_topological, locate_dofs_geometrical
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner, dot, nabla_grad
from mpi4py import MPI

# Create mesh
msh = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
                       [16, 16], CellType.triangle)


# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)),
                         np.isclose(x[1], 0.0))


# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

# P2 = element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,)) # Velocity elements for P1-P1
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,)) # Velocity elements for Taylor-Hood (P2-P1)
P1 = element("Lagrange", msh.basix_cell(), 1) # Pressure elements
V, Q = functionspace(msh, P2), functionspace(msh, P1)
print(f"There are this Many Degrees of Freedom in the Pressure Nodes: {Q.dofmap.index_map.size_local}")

# Create the Taylor-Hood function space
TH = mixed_element([P2, P1])
W = functionspace(msh, TH)

# No slip boundary condition
W0, _ = W.sub(0).collapse()
noslip = Function(W0)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
bc0 = dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
dofs = locate_dofs_topological((W.sub(0), W0), 1, facets)
bc1 = dirichletbc(lid_velocity, dofs, W.sub(0))

zero = Function(Q)
dofs = locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = dirichletbc(zero, dofs, W.sub(1))

bcs = [bc0, bc1, bc2]

# Create stokes flow solution
W0 = W.sub(0)
Q, _ = W0.collapse()
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(Q)

# Stabilization parameters per Andre Massing
h = ufl.CellDiameter(msh)
Beta = 0.2
mu_T = Beta*h*h # Stabilization coefficient
a = inner(grad(u), grad(v)) * dx
a -= inner(p, div(v)) * dx
a += inner(div(u), q) * dx
a += mu_T*inner(grad(p), grad(q)) * dx # Stabilization term

L = inner(f, v) * dx
L -= mu_T * inner(f, grad(q)) * dx # Stabilization  term

from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, bcs = bcs, petsc_options={'ksp_type': 'bcgs', 'ksp_rtol':1e-10, 'ksp_atol':1e-10})
U = Function(W)
U = problem.solve() # Solve the problem
print('Solved Stokes Flow')

# ------ Create/Define weak form of Navier-Stokes Equations ------
W0 = W.sub(0)
Q, _ = W0.collapse()
w = Function(W)
(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)
f = Function(Q)

nu = 100 # Viscocity (Reynolds number equals 1/nu)
a = inner(dot(u, nabla_grad(u)),v) * dx # Advection
a += nu*inner(grad(u),grad(v)) * dx # Diffusion
a -= inner(p,div(v)) * dx # Pressure
a -= inner(q,div(u)) * dx # Incompressibility

dw = ufl.TrialFunction(W)
dF = ufl.derivative(a, w, dw)

w = U.interpolate

from dolfinx.fem.petsc import NonlinearProblem
problem = NonlinearProblem(a, w, bcs=bcs, J=dF)
from dolfinx.nls.petsc import NewtonSolver

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

log.set_log_level(log.LogLevel.INFO)

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# Compute the solution
solver.solve(w)

# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()

The solver does not converge even when for Reynolds numbers close to zero (Re = 0.01), which makes me think I am not properly applying the initial condition. Below is the output which quickly diverges.

[2025-01-06 21:01:56.489] [info] PETSc Krylov solver starting to solve system.
[2025-01-06 21:01:56.515] [info] PETSc Krylov solver starting to solve system.
[2025-01-06 21:01:56.525] [info] Newton iteration 2: r (abs) = 20006.926307338606 (tol = 1e-10), r (rel) = 1.800340334930685 (tol = 1e-09)
[2025-01-06 21:01:56.529] [info] PETSc Krylov solver starting to solve system.
[2025-01-06 21:01:56.539] [info] Newton iteration 3: r (abs) = 91450.1616653592 (tol = 1e-10), r (rel) = 8.22922082847314 (tol = 1e-09)
[2025-01-06 21:01:56.542] [info] PETSc Krylov solver starting to solve system.
[2025-01-06 21:01:56.552] [info] Newton iteration 4: r (abs) = 119053.0616316309 (tol = 1e-10), r (rel) = 10.71309133446424 (tol = 1e-09)
[2025-01-06 21:01:56.556] [info] PETSc Krylov solver starting to solve system.
[2025-01-06 21:01:56.568] [info] Newton iteration 5: r (abs) = 850273.3197059985 (tol = 1e-10), r (rel) = 76.51257018028937 (tol = 1e-09)
[2025-01-06 21:01:56.572] [info] PETSc Krylov solver starting to solve system.
[2025-01-06 21:01:56.583] [info] Newton iteration 6: r (abs) = 917642.0480645928 (tol = 1e-10), r (rel) = 82.57480268486343 (tol = 1e-09)
[2025-01-06 21:01:56.587] [info] PETSc Krylov solver starting to solve system.
[2025-01-06 21:01:56.598] [info] Newton iteration 7: r (abs) = 1490986.1047592827 (tol = 1e-10), r (rel) = 134.1676568396576 (tol = 1e-09)
[2025-01-06 21:01:56.601] [info] PETSc Krylov solver starting to solve system.
2

Change

opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"

to

opts[f"{option_prefix}ksp_type"] = "preonly"

Conjugate gradient solvers are for symmetric problems.
The basic rule of thumb is to always start with a direct solver, and only move to iterative solvers if necessary for your problem size and if you have a MWE that works with a direct solver.

There are quite a few inconsistencies in your code. I’m not sure to what degree this is the code you’re actually working with, but things you may want to revisit are:

  • w = U.interpolate → w.interpolate(U)
  • You’re repeatedly reusing / overwriting variable names, and sometimes repeating operations. This is confusing and error-prone.
  • You’re using Q as the pressure space and later as the velocity space.
  • Your stokes problem does not use the same viscosity as the NS problem, so it is a bad initial guess in any case
  • Are you sure your stabilized formulation of your stokes problem is required, given that you’re using Taylor-Hood elements?
  • You’re missing advection stabilization in your NS problem, so for higher Reynolds numbers this will break.
1 Like
3

Thank you so much for your help, I have changed some variable names and added a del line to before the Navier-Stokes formulation to clear variable such as u, v, p, q, a that are used multiple times to prevent any errors. Your solution of using w.interpolate(U) worked, and I am now using the same viscosity for Stokes flow and Navier-Stokes flow. Becuase I am using Taylor-Hood elements I do not need stabilization, however I want to be able to use linear pressure and velocity elements to speed up computation (P1-P1), so I am using Galerkin least squares stabilization to allow P1-P1 elements to be used for Stokes flow. Eventually I want to be able to use P1-P1 elements for Navier-Stokes flow of Reynolds number up to 100, which will require SUPG and PSPG stabilization.