# Changing relaxation parameter while running Newton solver

Hi, I am trying to solve a Mixed Nonlinear Variational Problem.

Using the NonlinearVariationalSolver I don’t get convergence. The idea is to decrease the relaxation paramter in each newton iteration, if the residual is not decreasing.
My first question is, is there any implementation with varying relaxation parameter?

Otherwise, I would implement it myself. The problem here is the use of python and cpp code.

So far I could successfully modify update_solution(self, x, dx, relaxation, problem, k) as part of the NewtonSolver.
The next step would be the modification of the Newton method which I suppose gets called inside the NonlinearVariationalSolver method.
How can I modify the Newton method here? I had a look at the cpp code here: https://bitbucket.org/fenics-project/dolfin/src/master/dolfin/fem/NonlinearVariationalSolver.cpp

For the custom newton solver this and
https://fenicsproject.org/qa/13805/possible-certain-values-displacement-newton-solver-running/
was helpful.

Thank You

I don’t really like using `NonlinearVariationalSolver` as it doesn’t allow for much low-level control.

I’m also not exactly sure what you’re asking for. But you can see an example of setting a custom relaxation parameter here, which is inspired by Amrein & Wihler 2014.

Thank you for your response.

I am working on a time-dependent mixed formulation and used the solution here to formulate and solve the problem in the linear case.

Adding a nonlinearity to my formulation results in the Newton method not converging. In the code below I have labelled this as Method1.

My goal is to change the relaxation_parameter in each newton iteration.

Method2 follows the example you have given. Unfortunately I get

``````Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
``````

in each time step and I can’t explain why. It seem like the method update_solution doesn’t get called, when I run custom_solver.solve().

Another approach is Method3. Here I get the error

``````*** Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
``````

This post sasys the problem is resulting from the wrong use of split(). I tried both vuh.split() and split(vuh) but it doesnt make a difference

``````    P1 = VectorElement('Lagrange', cell=triangle, degree=degree)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

tol = 1E-14

def clamped_boundary(x, on_boundary):  # beam is only fixed on the left end
return on_boundary and near(x[0], 0.0, tol) #x[0] < tol

u_D1 = Constant((0, 0))

bc1 = DirichletBC(V.sub(0), u_D1, clamped_boundary)
bc2 = DirichletBC(V.sub(1), u_D1, clamped_boundary)

def force_on_rhs_boundary(x, on_boundary):
return on_boundary and near(x[0], 60.0, tol)

class topbot(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and 0.0 + tol < x[0] < 60.0 - tol

kappa = 0.5

tt = np.linspace(0, 10, 100)
S = 1.0
SW = 0

fct = str(S) + '-(' + str(S) + '-' + str(SW) + ')*exp(-kappa*t)'

u_D22 = Expression((fct, '0'), degree=1, kappa=kappa, t=0)
u_D21 = Constant((0, 0))

bc3 = DirichletBC(V.sub(0), u_D21, force_on_rhs_boundary)
bc4 = DirichletBC(V.sub(1), u_D22, force_on_rhs_boundary)

# Mark boundaries as "ds"
boundaries = MeshFunction("size_t", mesh, dim=1)

boundaries.set_all(0)

#Mark boundaries where u*n=0
#tlt.mark(boundaries, 1)
topbot = topbot()

topbot.mark(boundaries, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

#Combine Dirichlet bcs
bcs = [bc3, bc4]

# Define strain and stress
def epsilon(z):  # strain
return 0.5*(nabla_grad(z) + nabla_grad(z).T)
# return sym(nabla_grad(z))

def sigma(z):  # stress
return lambda_*nabla_div(z) * Identity(d) + 2 * mu * epsilon(z)

# Define test functions
vu = Function(V)
v, u = split(vu)

vnun = Function(V)
v_n, u_n = split(vnun)

(phi1, phi2) = TestFunctions(V)

d = v.geometric_dimension()

# Define source terms
f_1 = Constant((0, 0))
f_2 = Constant((0, 0))

# Define variational problem
gamma = Constant(gamma_nitsche)
h = CellDiameter(mesh)
n = FacetNormal(mesh)

def ppos(x):
return (x+abs(x))/2

# https://github.com/MiroK/fenics-nitsche/blob/master/poisson/poisson_circle_dirichlet.py
# https://pastebin.com/3zPiJd9m # Bsp von 2018

# Weak formulation
F = (dot(v-v_n, phi2) * dx
- dt * inner(sigma(u), epsilon(phi2)) * dx
- dt * dot(v, phi1) * dx
+ dot(u-u_n, phi1) * dx)

# Add penalty term
F += gamma/h**2 * dot(ppos(dot(u, n)), dot(phi2, n)) * ds(1)

# Add Nitsche term:
F -= dot(dot(n, dot(sigma(u), n)), dot(phi2, n)) * ds(1)

class Problem(NonlinearProblem):
def __init__(self, a, L):
self.a = a
self.L = L
NonlinearProblem.__init__(self)

def F(self, b, x):
assemble(self.L, tensor=b)

def J(self, A, x):
assemble(self.a, tensor=A)

vuh = Function(V)

du = TrialFunction(V)

# Write custom Newton Solver:

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("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps")

self.linear_solver().set_from_options()

def update_solution(self, x, dx, relaxation_parameter, nonlinear_problem,
iteration):
tau = 1.0
theta = min(sqrt(2.0*tau/norm(dx, norm_type="l2", mesh=V.mesh())), 1.0)
print(theta)
#info(f"Newton damping parameter: {theta:.3e}")
x.axpy(-theta, dx)

# Time stepping
for n in range(num_steps):

# Update current time
t += dt
print('t= ', t)
# Update bc
u_D22.t = t
J = derivative(F, vuh, du)

# Compute derivative of F for Newton method
R = action(F, vuh)

DR = derivative(R, vuh)

# Method1:
#problem = NonlinearVariationalProblem(R, vuh, bcs, DR)
#solver = NonlinearVariationalSolver(problem)

# Method2:
problem = Problem(DR, F)
custom_solver = CustomSolver()
custom_solver.solve(problem, vuh.vector())

# Method3:
"""
a_tol, r_tol = 1e-7, 1e-10

bcs_du = []
for bc in bcs:
bcs_du.append(bc)

print(bcs_du)
U_inc = Function(V)
nIter = 0
eps = 1

while eps > 1e-10 and nIter < 10:              # Newton iterations
nIter += 1
A, b = assemble_system(J, -F, bcs_du)
solve(A, U_inc.vector(), b)     # Determine step direction
eps = np.linalg.norm(U_inc.vector().array(), ord=2)

a = assemble(F)
for bc in bcs_du:
bc.apply(a)
#print 'b.norm("l2") = ', b.norm('l2'), 'np.linalg.norm(a.array(), ord = 2) = ', np.linalg.norm(a.array(), ord = 2)
#fnorm = b.norm('l2')
lmbda = 1.0     # step size, initially 1

vuh.vector()[:] += lmbda*U_inc.vector()    # New u vector

print(      {0:2d}  {1:3.2E}  {2:5e}'.format(nIter, eps, fnorm))

"""

# Split solutions
vh, uh = vuh.split(True)

# Update previous solution
v_n = vh
u_n = uh
``````