Line data    Source code

       1              : // Copyright 2024, UChicago Argonne, LLC
       2              : // All Rights Reserved
       3              : // Software Name: NEML2 -- the New Engineering material Model Library, version 2
       4              : // By: Argonne National Laboratory
       5              : // OPEN SOURCE LICENSE (MIT)
       6              : //
       7              : // Permission is hereby granted, free of charge, to any person obtaining a copy
       8              : // of this software and associated documentation files (the "Software"), to deal
       9              : // in the Software without restriction, including without limitation the rights
      10              : // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
      11              : // copies of the Software, and to permit persons to whom the Software is
      12              : // furnished to do so, subject to the following conditions:
      13              : //
      14              : // The above copyright notice and this permission notice shall be included in
      15              : // all copies or substantial portions of the Software.
      16              : //
      17              : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
      18              : // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      19              : // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
      20              : // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      21              : // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
      22              : // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
      23              : // THE SOFTWARE.
      24              : 
      25              : #include <iostream>
      26              : #include <iomanip>
      27              : 
      28              : #include "neml2/solvers/Newton.h"
      29              : #include "neml2/tensors/Scalar.h"
      30              : #include "neml2/tensors/functions/linalg/vector_norm.h"
      31              : #include "neml2/tensors/functions/linalg/solve.h"
      32              : 
      33              : namespace neml2
      34              : {
      35              : register_NEML2_object(Newton);
      36              : 
      37              : OptionSet
      38           18 : Newton::expected_options()
      39              : {
      40           18 :   OptionSet options = NonlinearSolver::expected_options();
      41           18 :   options.doc() = "The standard Newton-Raphson solver which always takes the 'full' Newton step.";
      42              : 
      43           18 :   return options;
      44            0 : }
      45              : 
      46           21 : Newton::Newton(const OptionSet & options)
      47           21 :   : NonlinearSolver(options)
      48              : {
      49           21 : }
      50              : 
      51              : Newton::Result
      52          216 : Newton::solve(NonlinearSystem & system, const NonlinearSystem::Sol<false> & x0)
      53              : {
      54              :   // Solve always takes place in the "scaled" space
      55          216 :   auto x = system.scale(x0);
      56              : 
      57              :   // The initial residual for relative convergence check
      58          216 :   auto R = system.residual(x);
      59          216 :   auto nR = linalg::vector_norm(R);
      60          216 :   auto nR0 = nR.clone();
      61              : 
      62              :   // Check for initial convergence
      63          216 :   if (converged(0, nR0, nR0))
      64              :   {
      65              :     // The final update is only necessary if we use AD
      66            0 :     if (R.requires_grad())
      67            0 :       final_update(system, x, R, system.Jacobian<true>());
      68              : 
      69            0 :     return {RetCode::SUCCESS, system.unscale(x), 0};
      70              :   }
      71              : 
      72              :   // Prepare any solver internal data before the iterative update
      73          216 :   prepare(system, x);
      74              : 
      75              :   // Continuing iterating until one of:
      76              :   // 1. nR < atol (success)
      77              :   // 2. nR / nR0 < rtol (success)
      78              :   // 3. i > miters (failure)
      79         1322 :   for (size_t i = 1; i < miters; i++)
      80              :   {
      81         1321 :     auto J = system.Jacobian<true>();
      82         1321 :     update(system, x, R, J);
      83         1321 :     R = system.residual(x);
      84         1321 :     nR = linalg::vector_norm(R);
      85              : 
      86              :     // Check for convergence
      87         1321 :     if (converged(i, nR, nR0))
      88              :     {
      89              :       // The final update is only necessary if we use AD
      90          215 :       if (R.requires_grad())
      91           26 :         final_update(system, x, R, system.Jacobian<true>());
      92              : 
      93          215 :       return {RetCode::SUCCESS, system.unscale(x), i};
      94              :     }
      95         1321 :   }
      96              : 
      97            1 :   return {RetCode::MAXITER, system.unscale(x), miters};
      98          431 : }
      99              : 
     100              : bool
     101         1537 : Newton::converged(size_t itr, const ATensor & nR, const ATensor & nR0) const
     102              : {
     103              :   // LCOV_EXCL_START
     104              :   if (verbose)
     105              :     std::cout << "ITERATION " << std::setw(3) << itr << ", |R| = " << std::scientific
     106              :               << at::max(nR).item<double>() << ", |R0| = " << std::scientific
     107              :               << at::max(nR0).item<double>() << std::endl;
     108              :   // LCOV_EXCL_STOP
     109              : 
     110         1537 :   return at::all(at::logical_or(nR < atol, nR / nR0 < rtol)).item<bool>();
     111              : }
     112              : 
     113              : void
     114         1224 : Newton::update(NonlinearSystem & /*system*/,
     115              :                NonlinearSystem::Sol<true> & x,
     116              :                const NonlinearSystem::Res<true> & r,
     117              :                const NonlinearSystem::Jac<true> & J)
     118              : {
     119         1224 :   x = NonlinearSystem::Sol<true>(x.variable_data() + Tensor(solve_direction(r, J)));
     120         1224 : }
     121              : 
     122              : void
     123           26 : Newton::final_update(NonlinearSystem & /*system*/,
     124              :                      NonlinearSystem::Sol<true> & x,
     125              :                      const NonlinearSystem::Res<true> & r,
     126              :                      const NonlinearSystem::Jac<true> & J)
     127              : {
     128           26 :   x = NonlinearSystem::Sol<true>(Tensor(x) + Tensor(solve_direction(r, J)));
     129           26 : }
     130              : 
     131              : NonlinearSystem::Sol<true>
     132         1347 : Newton::solve_direction(const NonlinearSystem::Res<true> & r, const NonlinearSystem::Jac<true> & J)
     133              : {
     134         1347 :   if (r.base_dim() == 0 && J.base_dim() == 0)
     135           51 :     return NonlinearSystem::Sol<true>(-Tensor(r) / Tensor(J));
     136              : 
     137         1296 :   return NonlinearSystem::Sol<true>(-linalg::solve(J, r));
     138              : }
     139              : 
     140              : } // namespace neml2