Hi everyone,
I’m working on a non linear singular problem, where I have to build and apply the Nullspace in every iteration and solve it with a PETscKrylovSolver so I have to programm the Newton Method per hand. First I’m trying to adapt the nonlinear poisson demo with a NewtonSolver to get a working algorithm and then use the PETscKrylovSolver but my assembled system doesn’t seem to update in each iteration. When using a residual criterion the first absolute error coincides with the dema, while using a incremental criterion doesn’t. Appreciate any help!
Cheers
Extra Question: How do I get the NonLinearVariationalSolver to report the iteration number, absolute and relative errors?
from dolfin import *
import numpy as np
import sys
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
du = TrialFunction(V)
u = Function(V)
du_ = Function(V)
f = Expression("x[0]*sin(x[1])", degree=2)
bc = DirichletBC(V, Constant(1.0), DirichletBoundary())
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx
J = derivative(F, u, du)
def ConvergenceFunction(iteration, v, v0, abs_tol, rel_tol, convergence):
absolute = norm(v, "l2")
if (iteration == 0):
absolute0 = absolute
else:
absolute0 = v0
relative = absolute/absolute0
if absolute < abs_tol or relative < rel_tol:
convergence = True
return convergence, absolute, relative
OwnNewton = False
criterion = "incremental"
if OwnNewton:
u = interpolate(Constant(0.0), V)
relaxation = 1.0
iteration = 0
iteration_max = 5
absolute = 1.0
absolute_tol = 1.0E-6
relative_tol = 1.0E-6
convergence = False
while iteration < iteration_max and convergence != True:
A, b = assemble_system(J, -F, bc)
solve(A, du_.vector(), b)
if (criterion == "residual"):
convergence, absolute, relative = ConvergenceFunction(iteration, b, absolute, absolute_tol, relative_tol, convergence)
elif (criterion == "incremental"):
convergence, absolute, relative = ConvergenceFunction(iteration, du_, absolute, absolute_tol, relative_tol, convergence)
else:
print("Convergence criterion not valid")
sys.exit()
#norm_u0 = norm(u, "l2")
#norm_du = norm(du_, "l2")
u.vector()[:] += relaxation * du_.vector()
#norm_u = norm(u, "l2")
print("Newton iteration %d: r (abs) = %.3e (tol = %.3e) r (rel) = %.3e (tol = %.3e)" % (
iteration, absolute, absolute_tol, relative, relative_tol))
#print("Newton iteration %d: norm_u0 = %.3e, norm_du = %.3e, norm_u = %.3e" % (iteration, norm_u0, norm_du, norm_u))
iteration += 1
else:
solve(F == 0, u, bc, J = J,
solver_parameters={"newton_solver":{"relative_tolerance":1e-6, "convergence_criterion":criterion}})
#problem = NonlinearVariationalProblem(F, u, bc, J)
#solver = NonlinearVariationalSolver(problem)
#solver.parameters["newton_solver"]["convergence_criterion"] = criterion
#solver.parameters["newton_solver"]["relative_tolerance"] = 1e-6
#solver.parameters["newton_solver"]["report"] = True
#solver.solve
file = File("Results/mod_nonlinear_poisson.pvd")
file << u