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