Home-cooking a NR solver: how to update solution field using numpy array

General
1

I would like to implement a arc-length solver.
For didactic reasons, I thought I would start reproducing a basic Newton-Raphson solver.

Say I have this setting, a 2D hyperelastic problem

import numpy as np
import ufl
import dolfinx
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx.plot import create_vtk_mesh
import pyvista
pyvista.set_jupyter_backend("panel")
L = 5.; H = 1;
cell_size = 0.1;

nx = 10
ny = 5

domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([L, H])], 
                               [nx, ny], mesh.CellType.quadrilateral)
ndim = domain.geometry.dim
V = fem.VectorFunctionSpace(domain, ("CG", 2))
def left(x):
    return np.isclose(x[0], 0)

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

fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full(len(left_facets), 1, dtype=np.int32), np.full(len(right_facets), 2, dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.indices[facet_tag.values==1])
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]
B = fem.Constant(domain, PETSc.ScalarType((0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0)))
v = ufl.TestFunction(V)
u = fem.Function(V)
# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))
# Elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2) 
problem = fem.petsc.NonlinearProblem(F, u, bcs)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

but now I want to write my own basic NR solver, to work on a local mesh.
To show where I get stuck, one load step is sufficient.
So I thought I would go on like this

from scipy import linalg
from scipy.linalg import norm
MAX_ITER = 10
MAX_ERROR = 10**-5
## function to get the Jacobian as numpy array
def ufl_form_2_numpy(F,u):
    J = ufl.derivative(F, u)
    z = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(J),bcs=bcs)
    z.assemble()
    z.convert("dense")
    return z.getDenseArray()

iter = 0
disp = 0
err = 10**5
tval0 = -1
while (err > MAX_ERROR) and iter <= MAX_ITER:
    print(5*"*"+f"iteration number: {iter}" + 3*"*"+ f"the error norm is {err}")
# set the traction to the vaue at the end of the first load step
    T.value[1] = 1 * tval0
# compute RHS and LHS from forms directly
    
    my_RHS = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(F)).getArray()
    my_Jacobian = ufl_form_2_numpy(F,u)
# solve linearised system
    x = linalg.solve(my_Jacobian, -my_RHS)
# update displacement
    disp =+ x
    
   err = norm(my_RHS)
   iter +=1 

### For this to work, I need to update the solution field u, for the new forms to be computed correctly.

I was unable to see how I could update the femFunction u, so that the forms are computed correctly based on the updated displacement for new iterations.

Thank you for your help and patience on such menial matters. Any other advice is of course greatly welcome. The approach above might be stone-age like, but it would make implementation of an arc-length procedure trivial, as all the necessary quantities are comfily handed out ready to be used.

2

The data of a dolfinx.fem.Function is stored as a numpy array in the attribute x.array, so if u is a dolfinx.fem.Function you can access the data and rewrite it by doing e.g.

u.x.array[:] = new_data

where new_data is a numpy array of correct shape. This would rewrite data for all degrees of freedom since I sliced it with :, but you could also index specific dofs if you wanted to.

Hope this is of help:)

3 Likes
3

Hello, thank you very much most helpful, I should have really known that.

So I edited my go at a NR solver, adding the update for the displacement field at the end of the loop

from scipy import linalg
from scipy.linalg import norm
MAX_ITER = 10
MAX_ERROR = 10**-5
# little function to get the Jacobian as a numpy array
def ufl_form_2_numpy(F,u):
    J = ufl.derivative(F, u)
    z = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(J),bcs=bcs)
    z.assemble()
    z.convert("dense")
    return z.getDenseArray()
iter = 0
disp = 0
err = 10**5
tval0 = -1


T.value[1] = 1 * tval0    
while (err > MAX_ERROR) and iter <= MAX_ITER:
    print(5*"*"+f"iteration number: {iter}" + 3*"*"+ f"the error norm is {err}")
    
    
    my_RHS = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(F)).getArray()
    #my_RHS = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(F)).array
    my_Jacobian = ufl_form_2_numpy(F,u)
    # solve for the displacement increments
    x = linalg.solve(my_Jacobian, -my_RHS)
    err = norm(my_RHS)
    # update displacement
    disp =+ x
    # update fem.Function u
    u.x.array[:]=disp
    
    
    iter +=1

which regretfully crashes after two iterations, with the error

ValueError: array must not contain infs or NaNs

the problem being that the residual vector my_RHScontains NaNs.
The vector is computed as

my_RHS = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(F)).getArray()

I cannot get it why it would contain NaNs.

I am wondering if the whole update procedure is missing its goal. If I plot the deformed shape during the iterations I get the following, from left to right for the initial configuration and the two completed iterations

image
image1994×458 18.9 KB

so the frst iteration is very reasonable but then it breaks down completely, it does not even respect the boundary conditions (this is basically a beam clamped on its left edge, with a vertical traction on the right).
Is there something that goes on with the dof ordering? Is updating u sufficient to update the form F? I guess there is a way to avoid recompiling the forms at each iteration.

It seems that dolfinx’s solve does reassamble the updated systems without recompiling, Form compilation — FEniCS Tutorial @ Sorbonne, is there a way to do anything similar in my case? Would like to learn

Thank you very much

4

Hi again,

Yes, it seems there is a problem with the setting of the boundary conditions (BCs). When looking at your code, there is nowhere inside the while loop where the BCs are set. You need to apply BCs to my_RHS, see e.g. Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial. In the solution for-loop in the code in that tutorial, the necessary steps of applying lifting (and performing a ghost update if you are running in parallel) and setting the BCs of the right-hand side (RHS) vector are demonstrated. This is something which is handled by dolfinx.fem.petsc.NonlinearProblem, so when assembling the system manually you will need to do these actions when assembling the RHS vector.

Regarding recompilation, it looks like you should be able to assemble the Jacobian only on the first timestep, and then reuse it on every timestep within the solution loop, since the system matrix does not change. The only thing that changes is your right-hand side. If you had time-dependent material parameters and BCs, it should be possible to only update the values of these (as long as the parameters are defined as dolfinx.fem.Constant’s), without having to reassemble the matrix. @dokken could correct me if I’m wrong here.

Cheers,
Halvor

1 Like
5

Hi @hherlyng ,

many thanks for this, I will look at your suggestions in detail. Indeed bcs are accounted for in the LHS, but nowhere in the RHS, looking into it.
I was actually also looking at the code for the NewtonSolver, I am actually confused about the difference between dolfinx.fem.petsc.create_matrix(problem.a) and dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(J),bcs=bcs)

# Create matrix and vector to be used for assembly # of the non-linear problem self._A = create_matrix(problem.a) self.setJ(problem.J, self._A) self._b = create_vector(problem.L) self.setF(problem.F, self._b) self.set_form(problem.form)

Interestingly, if I define a NonlinearProblem as

problem = fem.petsc.NonlinearProblem(F, u, bcs)

and follow the aproach in the source code

def ufl_form_2_numpy(F,u):
    J = ufl.derivative(F, u)
    #z = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(J),bcs=bcs)
    z = dolfinx.fem.petsc.create_matrix(problem.a)
    z.assemble()
    z.convert("dense")
    return z.getDenseArray()

the RHS is returned as singular!

Thanks a lot for your support.

6

No problem, happy to help!

I assume you mean that the LHS is singular? This is because dolfinx.fem.petsc.create_matrix() does just that, creates a matrix (without filling any values, hence the singularity), which is dimensionally compatible with the bilinear form passed as an argument. The function dolfinx.fem.petsc.assemble_matrix(), however, actually assembles the system matrix based on the bilinear form and the BCs passed as arguments.

2 Likes
7

I see, thanks, so in the source code presumably the assembly happens during the solve step, ok thanks a lot, will look into the bcs for the RHS, many thanks