Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PetscContactLineSearch.h"
      11             : #include "FEProblem.h"
      12             : #include "NonlinearSystem.h"
      13             : #include "libmesh/petsc_nonlinear_solver.h"
      14             : #include "libmesh/petsc_solver_exception.h"
      15             : #include <petscdm.h>
      16             : 
      17             : registerMooseObject("ContactApp", PetscContactLineSearch);
      18             : 
      19             : InputParameters
      20         164 : PetscContactLineSearch::validParams()
      21             : {
      22         164 :   return ContactLineSearchBase::validParams();
      23             : }
      24             : 
      25          82 : PetscContactLineSearch::PetscContactLineSearch(const InputParameters & parameters)
      26          82 :   : ContactLineSearchBase(parameters)
      27             : {
      28          82 :   _solver = dynamic_cast<PetscNonlinearSolver<Real> *>(
      29          82 :       _fe_problem.getNonlinearSystem(/*nl_sys_num=*/0).nonlinearSolver());
      30          82 :   if (!_solver)
      31           0 :     mooseError(
      32             :         "This line search operates only with Petsc, so Petsc must be your nonlinear solver.");
      33          82 : }
      34             : 
      35             : void
      36        2139 : PetscContactLineSearch::lineSearch()
      37             : {
      38        2139 :   PetscBool changed_y = PETSC_FALSE, changed_w = PETSC_FALSE;
      39             :   Vec X, F, Y, W, G, W1;
      40             :   SNESLineSearch line_search;
      41             :   PetscReal fnorm, xnorm, ynorm, gnorm;
      42             :   PetscBool domainerror;
      43             :   PetscReal ksp_rtol, ksp_abstol, ksp_dtol;
      44             :   PetscInt ksp_maxits;
      45             :   KSP ksp;
      46        2139 :   SNES snes = _solver->snes();
      47             : 
      48        2139 :   LibmeshPetscCall(SNESGetLineSearch(snes, &line_search));
      49        2139 :   LibmeshPetscCall(SNESLineSearchGetVecs(line_search, &X, &F, &Y, &W, &G));
      50        2139 :   LibmeshPetscCall(SNESLineSearchGetNorms(line_search, &xnorm, &fnorm, &ynorm));
      51        2139 :   LibmeshPetscCall(SNESLineSearchGetSNES(line_search, &snes));
      52        2139 :   LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_SUCCEEDED));
      53        2139 :   LibmeshPetscCall(SNESGetKSP(snes, &ksp));
      54        2139 :   LibmeshPetscCall(KSPGetTolerances(ksp, &ksp_rtol, &ksp_abstol, &ksp_dtol, &ksp_maxits));
      55        2139 :   LibmeshPetscCall(VecDuplicate(W, &W1));
      56             : 
      57        2139 :   if (!_user_ksp_rtol_set)
      58             :   {
      59          82 :     _user_ksp_rtol = ksp_rtol;
      60          82 :     _user_ksp_rtol_set = true;
      61             :   }
      62             : 
      63        2139 :   ++_nl_its;
      64             : 
      65             :   /* precheck */
      66        2139 :   LibmeshPetscCall(SNESLineSearchPreCheck(line_search, X, Y, &changed_y));
      67             : 
      68             :   /* temporary update */
      69        2139 :   _contact_lambda = 1.;
      70        2139 :   LibmeshPetscCall(VecWAXPY(W, -_contact_lambda, Y, X));
      71             : 
      72             :   /* compute residual to determine whether contact state has changed since the last non-linear
      73             :    * residual evaluation */
      74             :   _current_contact_state.clear();
      75        2139 :   LibmeshPetscCall(SNESComputeFunction(snes, W, F));
      76        2139 :   LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
      77        2139 :   if (domainerror)
      78           0 :     LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
      79             : 
      80        2139 :   LibmeshPetscCall(VecNorm(F, NORM_2, &fnorm));
      81             :   std::set<dof_id_type> contact_state_stored = _current_contact_state;
      82             :   _current_contact_state.clear();
      83        2139 :   printContactInfo(contact_state_stored);
      84             : 
      85        2139 :   if (_affect_ltol)
      86             :   {
      87          56 :     if (contact_state_stored != _old_contact_state)
      88             :     {
      89          18 :       LibmeshPetscCall(KSPSetTolerances(ksp, _contact_ltol, ksp_abstol, ksp_dtol, ksp_maxits));
      90          18 :       _console << "Contact set changed since previous non-linear iteration!" << std::endl;
      91             :     }
      92             :     else
      93          38 :       LibmeshPetscCall(KSPSetTolerances(ksp, _user_ksp_rtol, ksp_abstol, ksp_dtol, ksp_maxits));
      94             :   }
      95             : 
      96             :   size_t ls_its = 0;
      97        3298 :   while (ls_its < _allowed_lambda_cuts)
      98             :   {
      99        2779 :     _contact_lambda *= 0.5;
     100             :     /* update */
     101        2779 :     LibmeshPetscCall(VecWAXPY(W1, -_contact_lambda, Y, X));
     102             : 
     103        2779 :     LibmeshPetscCall(SNESComputeFunction(snes, W1, G));
     104        2779 :     LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
     105        2779 :     if (domainerror)
     106           0 :       LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
     107             : 
     108        2779 :     LibmeshPetscCall(VecNorm(G, NORM_2, &gnorm));
     109        2779 :     if (gnorm < fnorm)
     110             :     {
     111        1159 :       LibmeshPetscCall(VecCopy(G, F));
     112        1159 :       LibmeshPetscCall(VecCopy(W1, W));
     113        1159 :       fnorm = gnorm;
     114             :       contact_state_stored.swap(_current_contact_state);
     115             :       _current_contact_state.clear();
     116        1159 :       printContactInfo(contact_state_stored);
     117        1159 :       ++ls_its;
     118             :     }
     119             :     else
     120             :       break;
     121             :   }
     122             : 
     123        2139 :   LibmeshPetscCall(VecScale(Y, _contact_lambda));
     124             :   /* postcheck */
     125        2139 :   LibmeshPetscCall(SNESLineSearchPostCheck(line_search, X, Y, W, &changed_y, &changed_w));
     126             : 
     127        2139 :   if (changed_y)
     128           0 :     LibmeshPetscCall(VecWAXPY(W, -1., Y, X));
     129             : 
     130        2139 :   if (changed_w || changed_y)
     131             :   {
     132           0 :     LibmeshPetscCall(SNESComputeFunction(snes, W, F));
     133           0 :     LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
     134           0 :     if (domainerror)
     135           0 :       LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
     136             : 
     137             :     contact_state_stored.swap(_current_contact_state);
     138             :     _current_contact_state.clear();
     139           0 :     printContactInfo(contact_state_stored);
     140             :   }
     141             : 
     142             :   /* copy the solution over */
     143        2139 :   LibmeshPetscCall(VecCopy(W, X));
     144             : 
     145        2139 :   LibmeshPetscCall(SNESLineSearchComputeNorms(line_search));
     146             : 
     147        2139 :   LibmeshPetscCall(VecDestroy(&W1));
     148             : 
     149             :   _old_contact_state = std::move(contact_state_stored);
     150        2139 : }