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:
- Improve mesh quality; use finer mesh; use other geometry (imported via gmsh).
- Decrease solver relaxation parameter up to 0.1 (more damping).
- Use other KSP solvers (preonly, GAMG etc.)
- Change the initial guess for the NewtonSolver (via
MF.x.array[:]). - Change surface heat flux value (tried low ones such as 1 or 0.1 and high ones such as 10 or 100).
- 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
# Advecction
F2 = rho * cp * ufl.dot(u, ufl.grad(T)) * v_star * dx
# Diffusion
F2 += k * ufl.dot(ufl.grad(v), ufl.grad(T)) * dx
# 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