I have a nonlinear problem that solves for optimal displacement of a hyperelastic problem. The domain is an annulus with inner radius 1 and outer radius 2 (MWE given after the below output). I need to print the first absolute residual corresponding to the Newton Iteration 0 for each iteration of the loop (which are 1.779e+00 and 2.032e+00) of the output which is given as
output:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.078e+02 (tol = 1.000e-10) r (rel) = 6.061e+01 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 1.289e+00 (tol = 1.000e-10) r (rel) = 7.243e-01 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 3.506e-01 (tol = 1.000e-10) r (rel) = 1.971e-01 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 2.260e-01 (tol = 1.000e-10) r (rel) = 1.270e-01 (tol = 1.000e-09)
Newton iteration 5: r (abs) = 2.075e-04 (tol = 1.000e-10) r (rel) = 1.166e-04 (tol = 1.000e-09)
Newton iteration 6: r (abs) = 9.085e-01 (tol = 1.000e-10) r (rel) = 5.106e-01 (tol = 1.000e-09)
Newton iteration 7: r (abs) = 6.713e-05 (tol = 1.000e-10) r (rel) = 3.773e-05 (tol = 1.000e-09)
Newton iteration 8: r (abs) = 3.540e-03 (tol = 1.000e-10) r (rel) = 1.989e-03 (tol = 1.000e-09)
Newton iteration 9: r (abs) = 3.631e-07 (tol = 1.000e-10) r (rel) = 2.041e-07 (tol = 1.000e-09)
Newton iteration 10: r (abs) = 4.708e-07 (tol = 1.000e-10) r (rel) = 2.646e-07 (tol = 1.000e-09)
Newton iteration 11: r (abs) = 6.865e-11 (tol = 1.000e-10) r (rel) = 3.858e-11 (tol = 1.000e-09)
Newton solver finished in 11 iterations and 11 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.032e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.214e+02 (tol = 1.000e-10) r (rel) = 1.089e+02 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 1.290e+01 (tol = 1.000e-10) r (rel) = 6.349e+00 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 5.405e+00 (tol = 1.000e-10) r (rel) = 2.659e+00 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 3.252e+00 (tol = 1.000e-10) r (rel) = 1.600e+00 (tol = 1.000e-09)
Newton iteration 5: r (abs) = 1.077e-01 (tol = 1.000e-10) r (rel) = 5.300e-02 (tol = 1.000e-09)
Newton iteration 6: r (abs) = 2.892e-02 (tol = 1.000e-10) r (rel) = 1.423e-02 (tol = 1.000e-09)
Newton iteration 7: r (abs) = 1.086e-04 (tol = 1.000e-10) r (rel) = 5.343e-05 (tol = 1.000e-09)
Newton iteration 8: r (abs) = 1.161e-03 (tol = 1.000e-10) r (rel) = 5.710e-04 (tol = 1.000e-09)
Newton iteration 9: r (abs) = 1.548e-07 (tol = 1.000e-10) r (rel) = 7.618e-08 (tol = 1.000e-09)
Newton iteration 10: r (abs) = 3.063e-09 (tol = 1.000e-10) r (rel) = 1.507e-09 (tol = 1.000e-09)
Newton iteration 11: r (abs) = 1.206e-10 (tol = 1.000e-10) r (rel) = 5.933e-11 (tol = 1.000e-09)
Newton solver finished in 11 iterations and 11 linear solver iterations.
The MWE is as follows:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from ufl import *
from ufl import cofac
from dolfin import *
import numpy as np
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C = C2 - C1
mesh = generate_mesh(C, 100)
class inner(SubDomain):
def inside(self, x, on_boundary):
return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class outer(SubDomain):
def inside(self, x, on_boundary):
return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary
class annulus(SubDomain):
def inside(self, x, on_boundary):
return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
inner().mark(boundary_markers, 1)
outer().mark(boundary_markers, 2)
annulus().mark(surface_markers, 3)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(Constant(1)*ds(1)),assemble(Constant(1)*ds(2)))
print(assemble(Constant(1)*dx(3)))
V = VectorFunctionSpace(mesh, "Lagrange", 1)
n=FacetNormal(mesh)
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
T = Constant(5.0)
bcs=[]
# Kinematics
d = mesh.geometry().dim()
I = Identity(d)
F = I + grad(u)
C = F.T*F
M = cofac(F)
Finvt = (inv(F)).T
Ic = tr(C)
J = det(F)
E, nu = 10.0, 0.499
mu, lmbda = 27.9, Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) + ((nu)*(mu)/(1-2*nu))*(J-1)**2 - mu*(ln(J))
Coords = Expression(('x[0]','x[1]'),degree=0)
coords = project(Coords, V)
phi = coords+u
Pi = psi*dx - (0.5)*dot(M*T*(-n), phi)*ds(1)
ddPI = derivative(Pi, u, v)
solve(ddPI == 0, u, bcs,
form_compiler_parameters=ffc_options)
I tried the suggestion of @nate in the link Return relative residual from Newton solver and customized the solver as follows:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from ufl import *
from ufl import cofac
from dolfin import *
import numpy as np
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)
def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)
class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())
def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)
PETScOptions.set("NewtonSolver", "gmres")
PETScOptions.set("NewtonSolver_monitor")
PETScOptions.set("pc_type", "ilu")
self.linear_solver().set_from_options()
def converged(self, r, problem, iteration):
if iteration == 0:
self.r0 = r.norm("l2")
print(f"Iteration {iteration}, residual {r.norm('l2'):.6e}")
return super().converged(r, problem, iteration)
C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C = C2 - C1
mesh = generate_mesh(C, 100)
class inner(SubDomain):
def inside(self, x, on_boundary):
return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary
class outer(SubDomain):
def inside(self, x, on_boundary):
return near(x[0]**2 + x[1]**2, 4, eps=1e-3) and on_boundary
class annulus(SubDomain):
def inside(self, x, on_boundary):
return x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
inner().mark(boundary_markers, 1)
outer().mark(boundary_markers, 2)
annulus().mark(surface_markers, 3)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(Constant(1)*ds(1)),assemble(Constant(1)*ds(2)))
print(assemble(Constant(1)*dx(3)))
V = VectorFunctionSpace(mesh, "Lagrange", 1)
n=FacetNormal(mesh)
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
T = Constant(5.0)
bcs=[]
# Kinematics
d = mesh.geometry().dim()
I = Identity(d)
F = I + grad(u)
C = F.T*F
M = cofac(F)
Finvt = (inv(F)).T
Ic = tr(C)
J = det(F)
E, nu = 10.0, 0.499
mu, lmbda = 27.9, Constant(E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu/2)*(Ic - 3) + ((nu)*(mu)/(1-2*nu))*(J-1)**2 - mu*(ln(J))
Coords = Expression(('x[0]','x[1]'),degree=0)
coords = project(Coords, V)
phi = coords+u
for j in range(1,3):
Pi = psi*dx - (0.5)*dot(M*(T*j)*(-n), phi)*ds(1)
ddPI = derivative(Pi, u, v)
J = derivative(ddPI, u)
problem = Problem(J, ddPI, bcs)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())
This gives me the output:
Newton iteration 0: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 5: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 6: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 7: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 8: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 9: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 10: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 11: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 12: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 13: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 14: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 15: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 16: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 17: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 18: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 19: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 20: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 21: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 22: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 23: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 24: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 25: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 26: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 27: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 28: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 29: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 30: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 31: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 32: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 33: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 34: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 35: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 36: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 37: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 38: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 39: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 40: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 41: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 42: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 43: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 44: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 45: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 46: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 47: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 48: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 49: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 50: r (abs) = 1.779e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Iteration 0, residual 1.779230e+00
It seems to be stuck at the first Newton iteration only. I am not sure what I am missing here.
Also, I would like to find the values of the test function v and also save the test function in a pvd file to be later viewed in Paraview.
Is there a way I can print the values of the test function and save it as a PVD file?