NaN values when assembling PSPG stabilization

General
1

Hi.

I’m using dolfinx version 0.9.0 on a docker container.

I’m trying to solve Navier-Stokes with a PETSc nonlinear solver, it diverges when u_sol = 0 due to the residual having norm = NaN. The NaN values are coming from the block of the jacobian containing the PSPG stabilization term: J = ufl.derivative(ufl.inner(tau_pspg * R, ufl.grad(q)) * ufl.dx, u_sol, du).
If I remove tau_pspg from this term, the NaN values disapear, it seems very weird to me, because projecting tau_psps shows that it’s value is aproximately 0.0024 in the entire domain and has no NaN nor Infinite values.

I would appreciate if someone could guide me towards discovering how are this NaN values arising or how could I get rid of them.

Here is a MWE:

from petsc4py import PETSc
import numpy as np
from basix.ufl import element
from dolfinx import fem
import ufl
from mpi4py import MPI
from dolfinx.mesh import create_unit_square
from dolfinx.fem.petsc import assemble_matrix

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
dt = fem.Constant(mesh, PETSc.ScalarType(1 / 200))
mu = fem.Constant(mesh, PETSc.ScalarType(1 / 100))

element_velocity = element("Lagrange", "triangle", 1, shape=(2,))
element_pressure = element("Lagrange", "triangle", 1)

V = fem.functionspace(mesh, element_velocity)
Q = fem.functionspace(mesh, element_pressure)

v = ufl.TestFunction(V)
q = ufl.TestFunction(Q)
du = ufl.TrialFunction(V)

u_sol = fem.Function(V)
p_sol = fem.Function(Q)
u_prev = fem.Function(V)

def sigma(u, p, mu):
    return 2 * mu * ufl.sym(ufl.nabla_grad(u)) - p * ufl.Identity(len(u))

R = (u_sol - u_prev) / dt + ufl.dot(u_sol, ufl.nabla_grad(u_sol))
R -= ufl.div(sigma(u_sol, p_sol, mu))

h = ufl.CellDiameter(mesh)
vnorm = ufl.sqrt(ufl.inner(u_sol, u_sol))

tau_supg1 = h / (
    (2.0 * vnorm) + fem.Constant(mesh, PETSc.ScalarType(1e-10))
)  # avoid division by zero
tau_supg2 = dt / 2.0
tau_supg3 = (h * h) / (4.0 * mu)
tau_supg = (
    1 / (tau_supg1**2) + 1 / (tau_supg2**2) + 1 / (tau_supg3**2)
) ** (-1 / 2)

tau_pspg = tau_supg


J = ufl.derivative(ufl.inner(tau_pspg * R, ufl.grad(q)) * ufl.dx, u_sol, du)
A = assemble_matrix(fem.form(J))
A.assemble()
print(f"pspg operator norm: {A.norm()}")

J = ufl.derivative(ufl.inner(R, ufl.grad(q)) * ufl.dx, u_sol, du)
A = assemble_matrix(fem.form(J))
A.assemble()
print(f"pspg operator norm (without tau): {A.norm()}")

tau_pspg_f = fem.Function(Q)
tau_pspg_f.interpolate(fem.Expression(tau_pspg, Q.element.interpolation_points()))
print(f"nan values in tau: {np.where(np.isnan(tau_pspg_f.x.array))}")
print(f"inf values in tau: {np.where(np.isinf(tau_pspg_f.x.array))}")
print(f"tau min: {tau_pspg_f.x.array.min()}")
print(f"tau max: {tau_pspg_f.x.array.max()}")
1 Like
2

The NAN’s are coming from the differntiation of your supg1 stabilization parameter. To avoid proliferation of nonlinearities, I’d recommend thaking the previous solution for that stabilization parameter:

vnorm = ufl.sqrt(ufl.inner(u_prev, u_prev))
tau_supg1 = h / (
    (2.0 * vnorm) + fem.Constant(mesh, PETSc.ScalarType(1e-10))
)
1 Like