Getting residual "nan" from the start

variational formulation
1

Hi,

I am trying to solve a simple 1d visco-plasticity problem which is already solved using Abaqus here (chapter 2):
Lele, S. P. (2008). On a class of strain gradient plasticity theories: formulation and numerical implementation (Doctoral dissertation, Massachusetts Institute of Technology).

I’m getting RuntimeError: Newton solver did not converge because maximum number of iterations reached

When I check the log file, I see my residual is “nan” from the very start, and it stays the same.

Could anyone please help? Let me know for any clarification. Here is the implementation:

## Set-up
from dolfinx import fem, log, mesh, nls
from mpi4py import MPI
import numpy as np
from petsc4py import PETSc
import ufl

log.set_output_file("log.txt")

## Parameters
n_elems = 1000
x_left = 0.0
x_right = 10.0

dt = 0.02
t_final = 0.5  # Final time
ud_target = 1.0  # Displacement target
loading_rate = ud_target / t_final

μ = 100e9  # Shear modulus
S0 = 100e6  # Flow resistance
d0 = 0.1  # Reference flow rate
m = 0.02  # Rate sensitivity parameter

## Mesh
msh = mesh.create_interval(MPI.COMM_WORLD, n_elems, (x_left, x_right))

## Function space
P2 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 2)
ME = fem.FunctionSpace(msh, P2 * P2)

# # Boundary conditions
left_facets = mesh.locate_entities_boundary(msh, dim=0,
                                            marker=lambda x: np.isclose(x[0], x_left))
right_facets = mesh.locate_entities_boundary(msh, dim=0,
                                             marker=lambda x: np.isclose(x[0], x_right))
end_facets = mesh.locate_entities_boundary(msh, dim=0,
                                           marker=lambda x: np.logical_or(np.isclose(x[0], x_left),
                                                                          np.isclose(x[0], x_right)))                                              
bc_ud_left = fem.dirichletbc(PETSc.ScalarType(0.0),
                            fem.locate_dofs_topological(V=ME.sub(0), entity_dim=0, entities=left_facets),
                            ME.sub(0))
uD = fem.Constant(msh, PETSc.ScalarType(0.0))  # Displacement loading on the right end
bc_ud_right = fem.dirichletbc(uD,
                              fem.locate_dofs_topological(V=ME.sub(0), entity_dim=0, entities=right_facets),
                              ME.sub(0))
bc_γp = fem.dirichletbc(PETSc.ScalarType(0.0),
                        fem.locate_dofs_topological(V=ME.sub(1), entity_dim=0, entities=end_facets),
                        ME.sub(1))

bcs = [bc_ud_left, bc_ud_right, bc_γp]

## Variational form
v, q = ufl.TestFunctions(ME)
u = fem.Function(ME)  # Current solution
u0 = fem.Function(ME)  # Solution from the previous converged step

ud, γp = ufl.split(u)  # ud = displacement
ud0, γp0 = ufl.split(u0)

# Initial conditions
u.x.array[:] = 0.0

τ = μ * (ufl.grad(ud)[0] - γp)

γp_dot = (γp - γp0) / dt
τp = S0 * (ufl.algebra.Abs(γp_dot) / d0)**m * ufl.sign(γp_dot)

F0 = τ * ufl.grad(v)[0] * ufl.dx
F1 = (τ - τp) *  q * ufl.dx
F = F0 + F1

# # Solver
problem = fem.petsc.NonlinearProblem(F, u, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)

# # Solution
t = 0.0
u0.x.array[:] = u.x.array
while (t < t_final):
    t += dt
    uD.value = t * loading_rate
    r = solver.solve(u)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array
2

Perhaps Default absolute tolerance and relative tolerance - #4 by nate can help.

3

Thanks Nate. It definitely helped.

I have an exponent m, which is taken as 0.02. In the gradient it becomes negative (m - 1), which leads to division by zero, thus to nan. Just to check I changed it to 1.02, and the nan issue was resolved.

But I have to solve for m=0.02. So, I gave different initial values to u using —

u.sub(1).interpolate(lambda x: np.ones(x.shape[1], dtype=PETSc.ScalarType)*1e-7)
u.x.scatter_forward()

Now, it converges for first step only. The nan issue comes back for the second time step. I also tried different relaxation_parameter, but to no avail. I don’t know what else to try.

I have one question. I’m using 1000 2nd degree line elements. Each node has 2 dofs. So, my mixed element vector u should have 4002 dofs, and my split functions ud and γp should have 2001 dofs each. How can I check that this is indeed the case?

Thanks.

4

Consider Numerical continuation - Wikipedia to get to your desired value of m.

To print the number of DoFs, if V is your function space:

print(V.dofmap.index_map.size_global * V.dofmap.index_map_bs)
5

Thanks Nate. I’ll have a look at “Numerical continuation”. But I wonder how Abaqus solves it without resorting to such techniques. I think my mistakes lie somewhere else.

You provided the functions for dofs in the function space; I was asking about my split functions ud and γp. I think problem lies here. The current solution u is over the mixed function space ME, which has 4002 dofs. When I split u into ud and γp, how do I know that each got corresponding 2001 dofs?

Could you please have a cursory look at the code, and see if there is any obvious mistake?

Thanks.

6

Inserting into your code at the appropriate point:

def count_dofs(V):
	print(V.dofmap.index_map.size_global * V.dofmap.index_map_bs)

count_dofs(ME)
count_dofs(ME.sub(0).collapse()[0])
count_dofs(ME.sub(1).collapse()[0])

Gives me the expected

4002
2001
2001
1 Like