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