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 "CrankNicolson.h"
      11             : #include "NonlinearSystem.h"
      12             : #include "FEProblem.h"
      13             : 
      14             : registerMooseObject("MooseApp", CrankNicolson);
      15             : 
      16             : InputParameters
      17       14661 : CrankNicolson::validParams()
      18             : {
      19       14661 :   InputParameters params = TimeIntegrator::validParams();
      20       14661 :   params.addClassDescription("Crank-Nicolson time integrator.");
      21       14661 :   return params;
      22           0 : }
      23             : 
      24         264 : CrankNicolson::CrankNicolson(const InputParameters & parameters)
      25         264 :   : TimeIntegrator(parameters), _residual_old(addVector("residual_old", false, libMesh::GHOSTED))
      26             : {
      27         264 : }
      28             : 
      29             : void
      30       54960 : CrankNicolson::computeTimeDerivatives()
      31             : {
      32       54960 :   if (!_sys.solutionUDot())
      33           0 :     mooseError("CrankNicolson: Time derivative of solution (`u_dot`) is not stored. Please set "
      34             :                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
      35             : 
      36       54960 :   NumericVector<Number> & u_dot = *_sys.solutionUDot();
      37       54960 :   if (!_var_restriction)
      38             :   {
      39       38178 :     u_dot = *_solution;
      40       38178 :     computeTimeDerivativeHelper(u_dot, _solution_old);
      41             :   }
      42             :   else
      43             :   {
      44       16782 :     auto u_dot_sub = u_dot.get_subvector(_local_indices);
      45       16782 :     _solution->create_subvector(*_solution_sub, _local_indices, false);
      46       16782 :     _solution_old.create_subvector(*_solution_old_sub, _local_indices, false);
      47       16782 :     *u_dot_sub = *_solution_sub;
      48       16782 :     computeTimeDerivativeHelper(*u_dot_sub, *_solution_old_sub);
      49       16782 :     u_dot.restore_subvector(std::move(u_dot_sub), _local_indices);
      50             :     // Scatter info needed for ghosts
      51       16782 :     u_dot.close();
      52       16782 :   }
      53             : 
      54      128065 :   for (const auto i : index_range(_du_dot_du))
      55       73105 :     if (integratesVar(i))
      56       56323 :       _du_dot_du[i] = 2. / _dt;
      57       54960 : }
      58             : 
      59             : void
      60           0 : CrankNicolson::computeADTimeDerivatives(ADReal & ad_u_dot,
      61             :                                         const dof_id_type & dof,
      62             :                                         ADReal & /*ad_u_dotdot*/) const
      63             : {
      64           0 :   computeTimeDerivativeHelper(ad_u_dot, _solution_old(dof));
      65           0 : }
      66             : 
      67             : void
      68         124 : CrankNicolson::init()
      69             : {
      70         124 :   TimeIntegrator::init();
      71             : 
      72         124 :   if (!_sys.solutionUDot())
      73           0 :     mooseError("CrankNicolson: Time derivative of solution (`u_dot`) is not stored. Please set "
      74             :                "uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
      75             : 
      76             :   // time derivative is assumed to be zero on initial
      77         124 :   NumericVector<Number> & u_dot = *_sys.solutionUDot();
      78         124 :   u_dot.zero();
      79         124 :   computeDuDotDu();
      80             : 
      81         124 :   if (_nl)
      82             :   {
      83             :     // compute residual for the initial time step
      84             :     // Note: we can not directly pass _residual_old in computeResidualTag because
      85             :     //       the function will call postResidual, which will cause _residual_old
      86             :     //       to be added on top of itself prohibited by PETSc.
      87             :     //       Objects executed on initial have been executed by FEProblem,
      88             :     //       so we can and should directly call NonlinearSystem residual evaluation.
      89         124 :     _fe_problem.setCurrentResidualVectorTags({_nl->nonTimeVectorTag()});
      90         124 :     _nl->computeResidualTag(_nl->RHS(), _nl->nonTimeVectorTag());
      91         124 :     _fe_problem.clearCurrentResidualVectorTags();
      92             : 
      93         124 :     copyVector(_nl->RHS(), *_residual_old);
      94             :   }
      95         124 : }
      96             : 
      97             : void
      98       53050 : CrankNicolson::postResidual(NumericVector<Number> & residual)
      99             : {
     100             :   // PETSc 3.19 insists on having closed vectors when doing VecAXPY,
     101             :   // and that's probably a good idea with earlier versions too, but
     102             :   // we don't always get here with _Re_time closed.
     103             :   std::vector<unsigned char> inputs_closed = {
     104       53050 :       _Re_time->closed(), _Re_non_time->closed(), _residual_old->closed()};
     105             : 
     106             :   // We might have done work on one processor but not all processors,
     107             :   // so we have to sync our closed() checks.  Congrats to the BISON
     108             :   // folks for test coverage that caught that.
     109       53050 :   comm().min(inputs_closed);
     110             : 
     111       53050 :   if (!inputs_closed[0])
     112           0 :     _Re_time->close();
     113       53050 :   if (!inputs_closed[1])
     114           0 :     _Re_non_time->close();
     115       53050 :   if (!inputs_closed[2])
     116           0 :     _residual_old->close();
     117             : 
     118       53050 :   if (!_var_restriction)
     119             :   {
     120       36212 :     residual += *_Re_time;
     121       36212 :     residual += *_Re_non_time;
     122       36212 :     residual += *_residual_old;
     123             :   }
     124             :   else
     125             :   {
     126       16838 :     auto residual_sub = residual.get_subvector(_local_indices);
     127       16838 :     auto re_time_sub = _Re_time->get_subvector(_local_indices);
     128       16838 :     auto re_non_time_sub = _Re_non_time->get_subvector(_local_indices);
     129       16838 :     auto residual_old_sub = _residual_old->get_subvector(_local_indices);
     130       16838 :     *residual_sub += *re_time_sub;
     131       16838 :     *residual_sub += *re_non_time_sub;
     132       16838 :     *residual_sub += *residual_old_sub;
     133       16838 :     residual.restore_subvector(std::move(residual_sub), _local_indices);
     134       16838 :     _Re_time->restore_subvector(std::move(re_time_sub), _local_indices);
     135       16838 :     _Re_non_time->restore_subvector(std::move(re_non_time_sub), _local_indices);
     136       16838 :     _residual_old->restore_subvector(std::move(residual_old_sub), _local_indices);
     137       16838 :   }
     138       53050 : }
     139             : 
     140             : void
     141        9067 : CrankNicolson::postStep()
     142             : {
     143        9067 :   copyVector(*_Re_non_time, *_residual_old);
     144        9067 : }
     145             : 
     146             : Real
     147         124 : CrankNicolson::duDotDuCoeff() const
     148             : {
     149         124 :   return 2;
     150             : }