# NewtonSolver convergence problem for Reduced-Order-Model

Hello,

I am trying to solve a non-linear potential flow problem with a heat transfer (under the Boussinesq approximation). Briefly, the problem includes solving the system of two tightly coupled equations (see [1], Section 2 and Appendix B for more info).

Strong formulation
u_i = - \dfrac{\kappa}{\mu} \left( \dfrac{\partial P}{\partial x_i} + \rho_0 \beta g_i (T-T_0) \right)
\dfrac{\partial u_i}{\partial x_i} = 0
\rho_0 c_p u_i \dfrac{\partial T}{\partial x_i} - \dfrac{\partial}{\partial x_i} \left(k \dfrac{\partial T}{\partial x_i} \right) - Q = 0
T = T_0 on \Gamma_T
k\dfrac{\partial T}{\partial x_i} n_i = q_h on \Gamma_h
P = P_0 on \Gamma_p
u_i n_i = q_u on \Gamma_u

Finite element formulation (v and w - weighting functions)
\displaystyle \int_\Omega \dfrac{\partial v}{\partial x_i} \dfrac{\kappa}{\mu} \left( \dfrac{\partial P}{\partial x_i} + \rho_0 \beta g_i (T-T_0) \right) dV + \int_{S_u} v q_u dS = 0
\displaystyle \int_\Omega w^* \left[ -\rho_0 c_p \dfrac{\kappa}{\mu} \left( \dfrac{\partial P}{\partial x_i} + \rho_0 \beta g_i (T-T_0) \right) \dfrac{\partial T}{\partial x_i} - Q \right] dV + \int_\Omega \dfrac{\partial w}{\partial x_i} k \dfrac{\partial T}{\partial x_i} dV - \int_{S_h} w q_h dS = 0
w^* = w + \tau u_i \dfrac{\partial w}{\partial x_i}
where \tau is a cell-wise stabilization parameter for SUPG scheme.

As the author of [1;2;3] suggests I use a damped-step NewtonSolver with GMRES linear solver and Jacobi preconditioner. However, the residual blows up to infinity and the system cannot be solved. Also, the suggested [1;2] SUPG stabilization does not help a lot (it slows down solver divergence a little). Other tweaks do not seem to help too.
What I have tried:

1. Improve mesh quality; use finer mesh; use other geometry (imported via gmsh).
2. Decrease solver relaxation parameter up to 0.1 (more damping).
3. Use other KSP solvers (preonly, GAMG etc.)
4. Change the initial guess for the NewtonSolver (via MF.x.array[:]).
5. Change surface heat flux value (tried low ones such as 1 or 0.1 and high ones such as 10 or 100).
6. Change value of a kappa parameter (decrease up to 1e-8).

Note on the last statement:
The problem seems to be very sensitive to the value of kappa which is defined in [1;2] as a fictious permeability. As the intuition suggests NewtonSolver diverges slower for lower kappa values. However, I was unable to find one for which it converges (or converges in finite number of steps).

The test case for this problem is similar to a well known “hot room” problem. The geometry is a box with a floor heated up and top and side walls being kept at constant temperature. The equations are solved for P that is defined as p-\rho_0 g h and for T defined as T-T_0 since these quantities are defined up to a constant in this problem. Also, the pressure point constraint is introduced as suggested in [1] to have at least one Dirichlet BC.
See the following example code:

import numpy as np

from mpi4py import MPI
from petsc4py import PETSc
import ufl

from dolfinx.fem import (Constant, Function, functionspace, dirichletbc, form,
locate_dofs_geometrical, Expression)
from dolfinx.io import VTXWriter, XDMFFile
from dolfinx.mesh import create_box, locate_entities, meshtags
from dolfinx import cpp, default_scalar_type, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

comm = MPI.COMM_WORLD

x1, x2 = 0, 0.1
y1, y2 = 0, 0.1
z1, z2 = 0, 0.12
nx, ny, nz = 20, 20, 25
domain = create_box(comm,
[np.array([x1, y1, z1]), np.array([x2, y2, z2])],
[nx, ny, nz], cpp.mesh.CellType.tetrahedron)
tdim = domain.topology.dim

def locate_sides_and_top(x):
return np.isclose(x[0], x1) & np.isclose(x[0], x2) \
& np.isclose(x[1], y1) * np.isclose(x[1], y2) \
& np.isclose(x[2], z2)

def locate_heater(x):
return np.isclose(x[2], z1)

### Mesh Tags
fdim = tdim - 1
heater_tag = 2

marked_facets = locate_entities(domain, fdim, locate_heater)
entities_tags = np.full(marked_facets.shape[0], heater_tag, dtype=marked_facets.dtype)
face_tags = meshtags(domain, fdim, marked_facets, entities_tags)
face_tags.name = 'Facet markers'
domain.topology.create_connectivity(fdim, tdim)

#################### Function spaces ####################

# For cell-wise velocity
v_dg0 = ufl.VectorElement("DG", domain.ufl_cell(), 0)
V0 = functionspace(domain, v_dg0)

# For saving velocity to a file
v_cg1 = ufl.VectorElement("CG", domain.ufl_cell(), 1)
V1 = functionspace(domain, v_cg1)

# For cell-wise stabilization constants
s_dg0 = ufl.FiniteElement("DG", domain.ufl_cell(), 0)
Q0 = functionspace(domain, s_dg0)

# Unknown fields
s_cg1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
me = ufl.MixedElement([s_cg1, s_cg1])
MS = functionspace(domain, me)
Q, Q_dofs = MS.sub(0).collapse()
S, S_dofs = MS.sub(1).collapse()

MF = Function(MS)
p, T = ufl.split(MF)
q, v = ufl.TestFunctions(MS)

#################### Environment properties ####################

k = Constant(domain, default_scalar_type(24e-3))     # Thermal conductivity
mu = Constant(domain, default_scalar_type(1.83e-5))  # Dynamic viscosity
rho = Constant(domain, default_scalar_type(1.2))     # Density
beta = Constant(domain, default_scalar_type(1/300))  # Thermal expansion coef.
cp = Constant(domain, default_scalar_type(1000))     # Specific heat capacity at const. pressure
kappa = Constant(domain, default_scalar_type(1e-7))

g = Constant(domain, (default_scalar_type(0.0),
default_scalar_type(0.0),
default_scalar_type(-9.81)))   # Gravity

#################### Boundary Conditions ####################

bcT0 = dirichletbc(default_scalar_type(0), locate_dofs_geometrical(S, locate_sides_and_top), S)

bcT = [bcT0]

# Pressure point constraint
bcP_point = dirichletbc(
default_scalar_type(0),
locate_dofs_geometrical(Q, lambda x: np.isclose(x[0], x2) & np.isclose(x[1], y2) & np.isclose(x[2], z2)),
Q
)

bcp = [bcP_point]

bcs = bcp + bcT

#################### Equations ####################

ds_T = ufl.Measure('ds', domain, heater_tag, None, face_tags)
dx = ufl.dx

# --------------        Momentum balance (ROM)     -------------- #

# Velocity
u = - kappa / mu * (ufl.grad(p) + rho * beta * T * g)

# --------------             Mass balance           --------------#

# Continuity
F1 = - ufl.dot(ufl.grad(q), u) * dx

# -------------- Energy balance (SUPG stabilization) --------------#

### Cell-wise velocity vector
u0 = Function(V0)
u0_expr = Expression(u, V0.element.interpolation_points())
# u0.interpolate(u0_expr)
# u0.x.scatter_forward()

### Characteristic element size (computed once)
# Hexahedron: Characteristic element size = Diameter of Volume-equivalent sphere
# h = (6 * ufl.CellVolume(mesh) / np.pi)**(1/3)
# Tetrahedron: Characteristic element size = Height of Volume-equivalent regular tetrahedron
h = 2 / 3**(1/6) * ufl.CellVolume(domain)**(1/3)

h0 = Function(Q0)
h0_expr = Expression(h, Q0.element.interpolation_points())
h0.interpolate(h0_expr)
h0.x.scatter_forward()

### Characteristic timescale
tau1_invsq = 4 * ufl.dot(u0, u0) / h0**2
tau3_invsq = 9 * 16 * (k / (rho * cp))**2 / h0**4
tau_SUT = 10 * ufl.sqrt( 1 / ( tau1_invsq + tau3_invsq) )

tau0 = Function(Q0)
tau0_expr = Expression(tau_SUT, Q0.element.interpolation_points())
# tau0.interpolate(tau0_expr)
# tau0.x.scatter_forward()

### SUPG modified weighting function
v_star = v + tau0 * ufl.dot(u0, ufl.grad(v))

### Weak form residual
F2 = rho * cp * ufl.dot(u, ufl.grad(T)) * v_star * dx
# Diffusion
# Volumetric flux
# F2 -= v_star * s_T * dx
# Surface flux
f_n = Constant(domain, default_scalar_type(5))
F2 -= v * f_n * ds_T

# Combine to one form since TestFunctions are independent
F = F1 + F2

def ufl_eval(ufl_expr, V=None, f=None, upd=True):
""" Evaluate ufl epression on a given FunctionSpace or into a given Function."""
if f is None and V is None:
return
elif f is None:
f = Function(V)
elif V is None:
V = f.function_space()

f_expr = Expression(ufl_expr, V.element.interpolation_points())
f.interpolate(f_expr)

if upd:
f.x.scatter_forward()

return f

#################### Solving ####################

# Non-zero initial guess
MF.x.array[:] = 1.0

class MyNonlinearProblem(NonlinearProblem):

def __init__(self, *args, **kwargs) -> None:
super().__init__(*args, **kwargs)

def form(self, x: PETSc.Vec) -> None:
"""Redefine function to compute stabilization parameters every solver's step."""
super().form(x)
# return

u0.interpolate(u0_expr)
u0.x.scatter_forward()

tau0.interpolate(tau0_expr)
tau0.x.scatter_forward()

# Print once in parallel mode
if comm.rank == 0:
# MPI aware max()
u0max = u0.vector.max()
tau0max = tau0.vector.max()
log.log(log.LogLevel.INFO,
f"\n{u0max=}, \n{tau0max=}")

problem = MyNonlinearProblem(F, MF, bcs=bcs)
solver = NewtonSolver(comm, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-3
solver.report = True
solver.relaxation_parameter = 0.3
solver.max_it = 500

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
# opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "jacobi"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
solver.solve(MF)
log.set_log_level(log.LogLevel.WARNING)


[1] Joe Alexandersen et al., A “poor man’s” approach for high-resolution three-dimensional topology design for natural convection problems.
[2] Joe Alexandersen et al., A “poor man’s” approach to topology optimization of natural convection problems.
[3] Joe Alexandersen et al., Large scale three-dimensional topology optimisation of heat sinks cooled by natural convection

Forgot to mention that I am using FEniCSx from conda-forge and dolfinx v0.7.3, PETSc 3.20.5.