Hi all, I am trying to compute the 2-norm of the residual when the code is running in parallel. The reason to do this is I would like to design the line search algorithm based on this residual.
Here is the MWE and the residual can be computed correctly when the code is running in serial.
import numpy as np
import dolfinx
from dolfinx.cpp.mesh import CellType
from mpi4py import MPI
from petsc4py import PETSc
import ufl
mesh = dolfinx.BoxMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([1, 1, 1])], \
[2, 2, 2], CellType.hexahedron)
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
E = 1e4 # elastic modulus
nu = 0.3 # Poisson's ratio
lambda_ = dolfinx.Constant(mesh, E*nu / (1+nu) / (1-2*nu))
mu = dolfinx.Constant(mesh, E / (2*(1+nu)))
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], 1)
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
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 = dolfinx.MeshTags(mesh, mesh.topology.dim-1, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = dolfinx.Function(V)
with u_bc.vector.localForm() as loc:
loc.set(0)
left_dofs = dolfinx.fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.indices[facet_tag.values==1])
bcs = [dolfinx.DirichletBC(u_bc, left_dofs)]
B = dolfinx.Constant(mesh, (0, 0, 0)) # body force
T = dolfinx.Constant(mesh, (0, 0, 0)) # traction
u = dolfinx.Function(V) # test function
v = ufl.TestFunction(V) # trial function
d = len(u) # spatial dimension
I = ufl.variable(ufl.Identity(d)) # identity tensor
F = ufl.variable(I + ufl.grad(u)) # deformation gradient
C = ufl.variable(F.T * F) # right Cauchy-Green tensor
Ic = ufl.variable(ufl.tr(C)) # first invariant
J = ufl.variable(ufl.det(F)) # third invariant
# strain energy density of compressible neo-Hookean model
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lambda_ / 2) * (ufl.ln(J))**2
P = ufl.diff(psi, F) # first Piola-Kirchhoff stress
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", metadata=metadata)
R = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2)
problem = dolfinx.fem.NonlinearProblem(R, u, bcs)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
T.value[2] = -5
num_its, converged = solver.solve(u)
u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
rhs = dolfinx.fem.assemble_vector(R)
dolfinx.fem.set_bc(rhs, bcs)
rhs.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
print(rhs.norm())
The output indicates the computed residual 5.35121e-13 is the same as the one given by the backend.
2021-11-28 18:43:03.983 ( 0.197s) [main ]topologycomputation.cpp:628 INFO| Computing mesh entities of dimension 1
2021-11-28 18:43:03.983 ( 0.197s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:43:03.983 ( 0.197s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 0: r (abs) = 1.875 (tol = 1e-10) r (rel) = inf(tol = 1e-09)
2021-11-28 18:43:03.983 ( 0.197s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:43:03.985 ( 0.199s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:43:03.985 ( 0.199s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 1: r (abs) = 0.0438781 (tol = 1e-10) r (rel) = 4.53887(tol = 1e-09)
2021-11-28 18:43:03.986 ( 0.200s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:43:03.986 ( 0.200s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:43:03.986 ( 0.200s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 2: r (abs) = 1.55256e-07 (tol = 1e-10) r (rel) = 1.60601e-05(tol = 1e-09)
2021-11-28 18:43:03.987 ( 0.201s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:43:03.987 ( 0.201s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:43:03.988 ( 0.201s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 3: r (abs) = 5.35121e-13 (tol = 1e-10) r (rel) = 5.53543e-11(tol = 1e-09)
2021-11-28 18:43:03.988 ( 0.201s) [main ] NewtonSolver.cpp:252 INFO| Newton solver finished in 3 iterations and 3 linear solver iterations.
2021-11-28 18:43:03.990 ( 0.204s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:43:03.990 ( 0.204s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
5.351209992524201e-13
2021-11-28 18:43:04.023 ( 0.237s) [main ] loguru.cpp:526 INFO| atexit
However, when running this code in parallel with mpirun -n 3 python3 ResidualTest.py, we can see the residuals are different:
2021-11-28 18:46:38.936 ( 0.233s) [main ]topologycomputation.cpp:628 INFO| Computing mesh entities of dimension 1
2021-11-28 18:46:38.936 ( 0.238s) [main ]topologycomputation.cpp:628 INFO| Computing mesh entities of dimension 1
2021-11-28 18:46:38.936 ( 0.233s) [main ]topologycomputation.cpp:628 INFO| Computing mesh entities of dimension 1
2021-11-28 18:46:38.936 ( 0.233s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.936 ( 0.239s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.936 ( 0.233s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.936 ( 0.233s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 0: r (abs) = 1.875 (tol = 1e-10) r (rel) = inf(tol = 1e-09)
2021-11-28 18:46:38.937 ( 0.234s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.937 ( 0.239s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.937 ( 0.234s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.939 ( 0.236s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.939 ( 0.242s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.939 ( 0.237s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.940 ( 0.236s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 1: r (abs) = 0.0438781 (tol = 1e-10) r (rel) = 4.53887(tol = 1e-09)
2021-11-28 18:46:38.940 ( 0.237s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.940 ( 0.242s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.940 ( 0.237s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.941 ( 0.238s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.941 ( 0.243s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.941 ( 0.238s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.941 ( 0.238s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 2: r (abs) = 1.55256e-07 (tol = 1e-10) r (rel) = 1.60601e-05(tol = 1e-09)
2021-11-28 18:46:38.942 ( 0.238s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.942 ( 0.244s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.942 ( 0.239s) [main ] PETScKrylovSolver.cpp:98 INFO| PETSc Krylov solver starting to solve system.
2021-11-28 18:46:38.942 ( 0.239s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.942 ( 0.245s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.942 ( 0.240s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.943 ( 0.240s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 3: r (abs) = 5.39295e-13 (tol = 1e-10) r (rel) = 5.57861e-11(tol = 1e-09)
2021-11-28 18:46:38.943 ( 0.240s) [main ] NewtonSolver.cpp:252 INFO| Newton solver finished in 3 iterations and 3 linear solver iterations.
2021-11-28 18:46:38.945 ( 0.242s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.947 ( 0.244s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.947 ( 0.250s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.947 ( 0.244s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.947 ( 0.250s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
2021-11-28 18:46:38.947 ( 0.244s) [main ]topologycomputation.cpp:669 INFO| Requesting connectivity 2 - 3
1.9425172391655907
1.9425172391655907
1.9425172391655907
2021-11-28 18:46:39.003 ( 0.300s) [main ] loguru.cpp:526 INFO| atexit
2021-11-28 18:46:39.003 ( 0.305s) [main ] loguru.cpp:526 INFO| atexit
2021-11-28 18:46:39.003 ( 0.300s) [main ] loguru.cpp:526 INFO| atexit
After many tries, I still cannot solve this problem. So could anyone please help me with this problem?
Thanks in advance!