Residual is nan

hansangseok · March 5, 2025, 12:32pm

import pyvista
from dolfinx import mesh, fem, plot, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import conus
import basix
from petsc4py import PETSc

# material define #
density = conus.Material.Material("MassDensity", mass_density = 24.5)
elastic = conus.Material.Material("Elastic", E = 200000, nu = 0.2)
plastic = conus.Material.Material("ElastoPlastic", fy = 400, Et = 100, beta = 0.3, status="Isotropic")
material = conus.Material.MaterialBehavior(density, elastic, plastic)

# modeling #
L = 1  # Length of box 
W = 0.2 # Height in cross section of box

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                         [5, 1, 1], cell_type=mesh.CellType.hexahedron)

V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))

deg_quad = 1  # quadrature degree for internal state variable representation
W0e = basix.ufl.quadrature_element(
    domain.basix_cell(), value_shape=(), scheme="default", degree=deg_quad
)
We = basix.ufl.quadrature_element(
    domain.basix_cell(), value_shape=(6,), scheme="default", degree=deg_quad
)

S_ = fem.functionspace(domain, We)
S = fem.functionspace(domain, W0e)

# boundary condition #
def clamped_boundary(x):
    return np.isclose(x[0], 0)

def free_end(x):
    return np.isclose(x[0], L)

fdim = domain.topology.dim - 1

boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)

facet_values = np.hstack([np.full(len(boundary_facets), 1, dtype=np.int32),
                          np.full(len(free_end_facets), 2, dtype=np.int32)])
facet_entities = np.hstack([boundary_facets, free_end_facets])

facet_tags = mesh.meshtags(domain, fdim, facet_entities, facet_values)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)

bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)

f = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Function definition
c2 = fem.Function(S)
n_elas = fem.Function(S_)
sig = fem.Function(S_)
sig_n = fem.Function(S_)
back_n = fem.Function(S_)
pl_strain = fem.Function(S)
dp = fem.Function(S)

u = fem.Function(V, name="Total_displacement")
u = fem.Function(V, name="Total_displacement")
du = fem.Function(V, name="Iteration_correction")
Du = fem.Function(V, name="Current_increment")

v = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)


# strain and stress definition
def epsilon(v): # 3 x 3 tensor form
    e =ufl.sym(ufl.grad(v))
    return ufl.as_tensor([[e[0,0], e[0,1], e[0,2]], [e[0,1], e[1,1], e[1,2]], [e[0,2], e[1,2], e[2,2]]]) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

def elastic_sigma(material, disp): # 3 x 3 tensor form
    return (
        material.properties["Elastic"]["lambda_"] * ufl.nabla_div(disp) * ufl.Identity(len(disp)) 
        + 2 * material.properties["Elastic"]["mu"] * epsilon(disp)
    ) # σ = λ tr(ε) I₃ + 2 μ ε


# operator of form transfer
def to_vector(N):
    return ufl.as_vector([N[0, 0], N[1, 1], N[2, 2], N[0, 1], N[1,2], N[0,2]])

def to_tensor(N):
    return ufl.as_tensor([[N[0],N[3],N[5]],
                          [N[3],N[1],N[4]],
                          [N[5],N[4],N[2]]])

# operator of selecting positive value
ppos = lambda x: ufl.max_value(x, 0)

# hardening algorithm
def return_mapping(material, disp, old_sig, old_p, old_back):
    if material.properties["ElastoPlastic"]["status"] == "Isotropic":
        return Isotropic_hardening(material, disp, old_sig, old_p)
    elif material.properties["ElastoPlastic"]["status"] == "Kinematic":
        return Kinematic_hardening(material, u, old_sig, old_p, old_back)
    elif material.properties["ElastoPlastic"]["status"] == "Combined":
        return Combined_hardening(material, u, old_sig, old_p, old_back)

def Isotropic_hardening(material, disp, old_sig, old_p):
    E = material.properties["Elastic"]["E"]
    Et = material.properties["ElastoPlastic"]["Et"]
    fy = material.properties["ElastoPlastic"]["fy"]
    beta = material.properties["ElastoPlastic"]["beta"]    
    mu = material.properties["Elastic"]["mu"]
    H = E * Et / (E - Et)
    
    sig_n = to_tensor(old_sig)
    sig_elas = sig_n + elastic_sigma(material, disp)
    s = ufl.dev(sig_elas)
    sig_eq = ufl.sqrt(ufl.inner(s, s))
    
    yield_func = sig_eq - ufl.sqrt(2/3)*(fy - (1-beta)*H*old_p)
    consistency_param = ppos(yield_func) / (2*mu + 2/3*H)
    n_elas = s / sig_eq * ppos(yield_func)/yield_func
    new_sig = sig_elas - 2*mu*consistency_param*n_elas
    dp = ufl.sqrt(2/3)*consistency_param
    c2 = consistency_param / sig_eq

    return to_vector(new_sig), to_vector(n_elas), dp, c2

def Kinematic_hardening(material, u, old_sig, old_p, old_back):
    E = material.properties["Elastic"]["E"]
    Et = material.properties["ElastoPlastic"]["Et"]
    fy = material.properties["ElastoPlastic"]["fy"]
    beta = material.properties["ElastoPlastic"]["beta"]    
    mu = material.properties["Elastic"]["mu"]
    lambda_ = material.properties["Elastic"]["lambda_"]
    H = E * Et / (E - Et)
    
    sig_n = old_sig
    sig_elas = sig_n + elastic_sigma(material, u)
    s = ufl.dev(sig_elas) - old_back
    sig_eq = ufl.sqrt(ufl.inner(s, s))   

    yield_func = sig_eq - ufl.sqrt(2/3)*(fy)
    consistency_param = ppos(yield_func) / (2*mu)
    n_elas = s / sig_eq * ppos(yield_func)/yield_func
    new_sig = sig_elas - 2*mu*consistency_param*n_elas
    new_back = old_back + 2/3*beta*H*consistency_param*n_elas
    dp = ufl.sqrt(2/3)*consistency_param     
    c2 = consistency_param / sig_eq
    
    return new_sig, n_elas, dp, new_back, c2

def Combined_hardening(material, u, old_sig, old_p, old_back):
    E = material.properties["Elastic"]["E"]
    Et = material.properties["ElastoPlastic"]["Et"]
    fy = material.properties["ElastoPlastic"]["fy"]
    beta = material.properties["ElastoPlastic"]["beta"]    
    mu = material.properties["Elastic"]["mu"]
    lambda_ = material.properties["Elastic"]["lambda_"]
    H = E * Et / (E - Et)
    
    sig_n = old_sig  
    sig_elas = sig_n + elastic_sigma(material, u)
    s = ufl.dev(sig_elas) - old_back
    sig_eq = ufl.sqrt(ufl.inner(s, s))
    
    yield_func = sig_eq - ufl.sqrt(2/3)*(fy - (1-beta)*H*old_p)
    consistency_param = ppos(yield_func) / (2*mu + 2/3*H)
    n_elas = s / sig_eq * ppos(yield_func)/yield_func
    new_sig = sig_elas - 2*mu*consistency_param*n_elas
    new_back = old_back + 2/3*beta*H*consistency_param*n_elas
    dp = ufl.sqrt(2/3)*consistency_param
    c2 = consistency_param / sig_eq
    
    return new_sig, n_elas, dp, new_back, c2


# consistent tangent stress
def tang_sig(material, disp):
    E = material.properties["Elastic"]["E"]
    Et = material.properties["ElastoPlastic"]["Et"]  
    mu = material.properties["Elastic"]["mu"]
    H = E * Et / (E - Et)
    
    k1 = 4*mu**2/(2*mu+2/3*H)
    k2 = 4*mu**2/c2
    N_elas = to_tensor(n_elas)
    tangent_sig = elastic_sigma(material, disp) - (k1-k2) * ufl.inner(N_elas, epsilon(disp)) * N_elas - k2 * mu * ufl.dev(epsilon(disp))
    return tangent_sig


print(epsilon(u_))
Residual = ufl.inner(epsilon(u_), to_tensor(sig)) * ufl.dx - ufl.dot(f, u_) * ufl.dx - ufl.dot(T, u_) * ds(2)
#print(Residual)
tangent_form = ufl.inner(epsilon(v), tang_sig(material, u_)) * ufl.dx



basix_celltype = getattr(basix.CellType, domain.topology.cell_type.name)
quadrature_points, weights = basix.make_quadrature(basix_celltype, 1)

map_c = domain.topology.index_map(domain.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)


def interpolate_quadrature(ufl_expr, function):
    expr_expr = fem.Expression(ufl_expr, quadrature_points)
    expr_eval = expr_expr.eval(domain, cells)
    function.x.array[:] = expr_eval.flatten()[:]

class CustomLinearProblem(fem.petsc.LinearProblem):
    def assemble_rhs(self, u=None):
        """Assemble right-hand side and lift Dirichlet bcs.
        Parameters
        ----------
        u : dolfinx.fem.Function, optional
            For non-zero Dirichlet bcs u_D, use this function to assemble rhs with the value u_D - u_{bc}
            where u_{bc} is the value of the given u at the corresponding. Typically used for custom Newton methods
            with non-zero Dirichlet bcs.
        """
        # Assemble rhs
        with self._b.localForm() as b_loc:
            b_loc.set(0)
        fem.petsc.assemble_vector(self._b, self._L)

        # Apply boundary conditions to the rhs
        x0 = [] if u is None else [u.vector]
        fem.petsc.apply_lifting(self._b, [self._a], bcs=[self.bcs], x0=x0)
        self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        x0 = None if u is None else u.vector
        fem.petsc.set_bc(self._b, self.bcs, x0)

    def assemble_lhs(self):
        self._A.zeroEntries()
        fem.petsc.assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
        self._A.assemble()

    def solve_system(self):
        # Solve linear system and update ghost values in the solution
        self._solver.solve(self._b, self._x)
        self.u.x.scatter_forward()

tangent_problem = CustomLinearProblem(
    tangent_form,
    -Residual,
    u=du,
    bcs=[bc],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)

Nitermax, tol = 200, 1e-6  # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr + 1)[1:] ** 0.5
results = np.zeros((Nincr + 1, 3))

sig.x.petsc_vec.set(0.0)
sig_n.x.petsc_vec.set(0.0)
pl_strain.x.petsc_vec.set(0.0)
u.x.petsc_vec.set(0.0)
n_elas.x.petsc_vec.set(0.0)
c2.x.petsc_vec.set(0.0)
back_n.x.petsc_vec.set(0.0)

sig_, n_elas_, dp_, c2_ = return_mapping(material, Du, sig_n, pl_strain, None)


for i, t in enumerate(load_steps):
    T.value = np.array([1000, 0, 0]) * t
    
    # compute the residual norm at the beginning of the load step
    print(f"========== Load step : {i} / {Nincr} ==========")

    tangent_problem.assemble_rhs()
    print(f"After assembling rhs, _b: {tangent_problem._b.getArray()}")
    #print(f"After assembling lhs, _A: {tangent_problem._A.getValuesCSR()}")
    #print(f"Solution vector _x: {tangent_problem._x.getArray()}")
    nRes0 = tangent_problem._b.norm()
    print(f"Initial residual norm: {nRes0}")
    if np.isnan(nRes0):
        print(f"Warning: Residual norm is NaN at load step {i}")
        break

    nRes = nRes0
    Du.x.array[:] = 0
    print(nRes)
    niter = 0

when i print residual print(f"After assembling rhs, _b:{tangent_problem._b.getArray()}")

the value is nan
why this happen?

nate · March 5, 2025, 4:36pm

I don’t have the time to go through your code, but you could try and hunt down singularities in your formulation. For example, make sure you don’t have operators like \mathrm{det}( \cdot ), (\nabla \cdot)^{-1}, \sqrt{\cdot}, etc, with a zero or constant initial guess of your solution varaible.

See also, Default absolute tolerance and relative tolerance - #4 by nate.

hansangseok · March 6, 2025, 3:55am

Thanks ! It seems that some parts of my code had division by zero, which caused the nan values to appear.

Topic		Replies	Views
Assemble_scalar of hyperelasticity PDE yields NaN dolfinx	1	375	May 11, 2023
Porting assemble_system and assemble to dolfinx0.7 dolfinx boundary-conditions , assembly	1	197	January 2, 2024
Getting residual "nan" from the start variational formulation	5	643	October 4, 2023
Hyperelasticity example issues General	9	246	May 3, 2024
ValueError: too many values to unpack (expected 1), while solving 2D elasticity problem Errors	4	519	September 27, 2022

Residual is nan

Related topics