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 "PetscProjectSolutionOntoBounds.h"
      11             : #include "FEProblem.h"
      12             : #include "DisplacedProblem.h"
      13             : #include "NonlinearSystem.h"
      14             : #include "GeometricSearchData.h"
      15             : #include "PenetrationLocator.h"
      16             : 
      17             : #include "libmesh/petsc_nonlinear_solver.h"
      18             : #include "libmesh/petsc_solver_exception.h"
      19             : #include "libmesh/petsc_vector.h"
      20             : #include <petscdm.h>
      21             : 
      22             : registerMooseObject("ContactApp", PetscProjectSolutionOntoBounds);
      23             : 
      24             : InputParameters
      25          24 : PetscProjectSolutionOntoBounds::validParams()
      26             : {
      27          24 :   return LineSearch::validParams();
      28             : }
      29             : 
      30          12 : PetscProjectSolutionOntoBounds::PetscProjectSolutionOntoBounds(const InputParameters & parameters)
      31          12 :   : LineSearch(parameters), _nl(_fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0))
      32             : {
      33          12 :   _solver = dynamic_cast<PetscNonlinearSolver<Real> *>(
      34          12 :       _fe_problem.getNonlinearSystem(/*nl_sys_num=*/0).nonlinearSolver());
      35          12 :   if (!_solver)
      36           0 :     mooseError(
      37             :         "This line search operates only with Petsc, so Petsc must be your nonlinear solver.");
      38          12 : }
      39             : 
      40             : void
      41          12 : PetscProjectSolutionOntoBounds::initialSetup()
      42             : {
      43          12 :   _displaced_problem = _fe_problem.getDisplacedProblem().get();
      44          12 :   if (!_displaced_problem)
      45           0 :     mooseError("This line search only makes sense in a displaced context");
      46          12 : }
      47             : 
      48             : void
      49          24 : PetscProjectSolutionOntoBounds::lineSearch()
      50             : {
      51          24 :   PetscBool changed_y = PETSC_FALSE;
      52             :   Vec X, F, Y, W, G;
      53             :   SNESLineSearch line_search;
      54             :   PetscReal fnorm, xnorm, ynorm;
      55          24 :   SNES snes = _solver->snes();
      56             : 
      57          24 :   LibmeshPetscCall(SNESGetLineSearch(snes, &line_search));
      58          24 :   LibmeshPetscCall(SNESLineSearchGetVecs(line_search, &X, &F, &Y, &W, &G));
      59          24 :   LibmeshPetscCall(SNESLineSearchGetNorms(line_search, &xnorm, &fnorm, &ynorm));
      60          24 :   LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_SUCCEEDED));
      61             : 
      62          24 :   LibmeshPetscCall(SNESLineSearchPreCheck(line_search, X, Y, &changed_y));
      63             : 
      64             :   // apply the full newton step
      65          24 :   LibmeshPetscCall(VecWAXPY(W, -1., Y, X));
      66             : 
      67             :   {
      68          24 :     PetscVector<Number> solution(W, this->comm());
      69             : 
      70          24 :     _nl.setSolution(solution);
      71             : 
      72             :     // Displace the mesh and update the displaced geometric search objects
      73          24 :     _displaced_problem->updateMesh();
      74          24 :   }
      75             : 
      76             :   // A reference to the displaced MeshBase object
      77          24 :   const auto & mesh = _displaced_problem->mesh().getMesh();
      78             : 
      79          24 :   const auto & pen_locs = _displaced_problem->geomSearchData().getPenetrationLocators();
      80             : 
      81             :   // Keep track of the secondary nodes that we push back onto the primary face. We'll eventually
      82             :   // check to make sure that we didn't have a corner node at the intersection of two secondary faces
      83             :   // that we tried to displace twice. As this stands now this won't cover the case wherethe
      84             :   // intersection happens only across processes
      85             :   std::set<dof_id_type> nodes_displaced;
      86             : 
      87             :   std::vector<unsigned int> disp_nums;
      88             : 
      89             :   // Generate the displaced variable numbers
      90          72 :   for (const auto & disp_name : _displaced_problem->getDisplacementVarNames())
      91          48 :     disp_nums.push_back(_nl.system().variable_number(disp_name));
      92             : 
      93          48 :   for (const auto & pen_loc : pen_locs)
      94          72 :     for (const auto & pinfo_pair : pen_loc.second->_penetration_info)
      95             :     {
      96          48 :       auto node_id = pinfo_pair.first;
      97          48 :       auto pen_info = pinfo_pair.second;
      98             : 
      99             :       // We have penetration
     100          48 :       if (pen_info && pen_info->_distance > 0)
     101             :       {
     102             : // Avoid warning in optimized modes about unused variables
     103             : #ifndef NDEBUG
     104             :         // Make sure we haven't done this node before
     105             :         auto pair = nodes_displaced.insert(node_id);
     106             :         mooseAssert(pair.second, "Node id " << node_id << " has already been displaced");
     107             : #endif
     108             : 
     109          36 :         const auto & node = mesh.node_ref(node_id);
     110             : 
     111             :         // If this is not a local node, we will let displacement happen on another process
     112          36 :         if (node.processor_id() != this->processor_id())
     113          12 :           continue;
     114             : 
     115             :         // The vector that we need to displace by
     116             :         auto required_solution_change = pen_info->_distance * pen_info->_normal;
     117             : 
     118             :         unsigned component = 0;
     119             :         std::vector<PetscInt> indices;
     120             :         std::vector<PetscScalar> values;
     121          72 :         for (auto disp_num : disp_nums)
     122             :         {
     123          48 :           auto dof_number = node.dof_number(/*sys=*/0, disp_num, /*component=*/0);
     124          48 :           indices.push_back(static_cast<PetscInt>(dof_number));
     125          48 :           values.push_back(static_cast<PetscScalar>(required_solution_change(component++)));
     126             :         }
     127          24 :         LibmeshPetscCall(VecSetValues(
     128             :             W, static_cast<PetscInt>(indices.size()), indices.data(), values.data(), ADD_VALUES));
     129          24 :       }
     130             :     }
     131             : 
     132          24 :   LibmeshPetscCall(VecAssemblyBegin(W));
     133          24 :   LibmeshPetscCall(VecAssemblyEnd(W));
     134             : 
     135          24 :   LibmeshPetscCall(SNESComputeFunction(snes, W, F));
     136             : #if PETSC_VERSION_LESS_THAN(3, 25, 0)
     137             :   PetscBool domainerror;
     138             :   LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
     139             :   if (domainerror)
     140             :     LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
     141             : #else
     142             :   SNESConvergedReason reason;
     143          24 :   LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
     144          24 :   if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
     145           0 :     LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
     146             : #endif
     147             : 
     148          24 :   LibmeshPetscCall(VecNorm(F, NORM_2, &fnorm));
     149             : 
     150          24 :   LibmeshPetscCall(VecCopy(W, X));
     151             : 
     152          24 :   LibmeshPetscCall(SNESLineSearchComputeNorms(line_search));
     153          48 : }