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             : // MOOSE includes
      11             : #include "ExplicitTimeIntegrator.h"
      12             : #include "NonlinearSystem.h"
      13             : #include "FEProblem.h"
      14             : #include "PetscSupport.h"
      15             : #include "KernelBase.h"
      16             : #include "DGKernelBase.h"
      17             : #include "ScalarKernelBase.h"
      18             : #include "FVElementalKernel.h"
      19             : #include "FVFluxKernel.h"
      20             : #include "NodalKernelBase.h"
      21             : 
      22             : // libMesh includes
      23             : #include "libmesh/enum_convergence_flags.h"
      24             : 
      25             : using namespace libMesh;
      26             : 
      27             : InputParameters
      28       10179 : ExplicitTimeIntegrator::validParams()
      29             : {
      30       10179 :   InputParameters params = TimeIntegrator::validParams();
      31             : 
      32       40716 :   MooseEnum solve_type("consistent lumped lump_preconditioned", "consistent");
      33             : 
      34       30537 :   params.addParam<MooseEnum>(
      35             :       "solve_type",
      36             :       solve_type,
      37             :       "The way to solve the system.  A 'consistent' solve uses the full mass matrix and actually "
      38             :       "needs to use a linear solver to solve the problem.  'lumped' uses a lumped mass matrix with "
      39             :       "a simple inversion - incredibly fast but may be less accurate.  'lump_preconditioned' uses "
      40             :       "the lumped mass matrix as a preconditioner for the 'consistent' solve");
      41             : 
      42       20358 :   return params;
      43       10179 : }
      44             : 
      45         664 : ExplicitTimeIntegrator::ExplicitTimeIntegrator(const InputParameters & parameters)
      46             :   : TimeIntegrator(parameters),
      47             :     MeshChangedInterface(parameters),
      48         664 :     _solve_type(getParam<MooseEnum>("solve_type")),
      49        1328 :     _explicit_residual(addVector("explicit_residual", false, PARALLEL)),
      50        1328 :     _solution_update(addVector("solution_update", true, PARALLEL)),
      51        1328 :     _mass_matrix_diag_inverted(addVector("mass_matrix_diag_inverted", true, GHOSTED)),
      52        1328 :     _ones(nullptr)
      53             : {
      54         664 :   _Ke_time_tag = _fe_problem.getMatrixTagID("TIME");
      55             : 
      56             :   // This effectively changes the default solve_type to LINEAR instead of PJFNK,
      57             :   // so that it is valid to not supply solve_type in the Executioner block:
      58         664 :   if (_nl)
      59         332 :     _fe_problem.solverParams(_nl->number())._type = Moose::ST_LINEAR;
      60             : 
      61         664 :   if (_solve_type == LUMPED || _solve_type == LUMP_PRECONDITIONED)
      62         516 :     _ones = addVector("ones", true, PARALLEL);
      63             :   // don't set any of the common SNES and KSP-related petsc options to prevent unused option
      64             :   // warnings
      65         664 :   Moose::PetscSupport::dontAddCommonSNESOptions(_fe_problem);
      66         664 :   Moose::PetscSupport::dontAddCommonKSPOptions(_fe_problem);
      67         664 : }
      68             : 
      69             : void
      70         326 : ExplicitTimeIntegrator::initialSetup()
      71             : {
      72         326 :   setupSolver();
      73             : 
      74         326 :   if (_nl)
      75             :   {
      76         326 :     std::unordered_set<unsigned int> vars_to_check;
      77         326 :     if (!_vars.empty())
      78           0 :       vars_to_check = _vars;
      79             :     else
      80         664 :       for (const auto i : make_range(_nl->nVariables()))
      81         338 :         vars_to_check.insert(i);
      82             : 
      83         326 :     const auto mass_matrix_tag_id = massMatrixTagID();
      84         652 :     std::set<TagID> matrix_tags = {mass_matrix_tag_id};
      85         326 :     auto fv_object_starting_query = _fe_problem.theWarehouse()
      86         326 :                                         .query()
      87         326 :                                         .template condition<AttribMatrixTags>(matrix_tags)
      88         326 :                                         .template condition<AttribSysNum>(_nl->number());
      89             : 
      90         661 :     for (const auto var_id : vars_to_check)
      91             :     {
      92         338 :       const auto & var_name = _nl->system().variable_name(var_id);
      93         338 :       const bool field_var = _nl->hasVariable(var_name);
      94         338 :       if (!field_var)
      95             :         mooseAssert(_nl->hasScalarVariable(var_name),
      96             :                     var_name << " should be either a field or scalar variable");
      97         338 :       if (field_var)
      98             :       {
      99         287 :         const auto & field_var = _nl->getVariable(/*tid=*/0, var_name);
     100         287 :         if (field_var.isFV())
     101             :         {
     102           0 :           std::vector<FVElementalKernel *> fv_elemental_kernels;
     103           0 :           auto var_query = fv_object_starting_query.clone().template condition<AttribVar>(var_id);
     104           0 :           auto var_query_clone = var_query.clone();
     105           0 :           var_query.template condition<AttribSystem>("FVElementalKernel")
     106           0 :               .queryInto(fv_elemental_kernels);
     107           0 :           if (fv_elemental_kernels.size())
     108           0 :             continue;
     109             : 
     110           0 :           std::vector<FVFluxKernel *> fv_flux_kernels;
     111           0 :           var_query_clone.template condition<AttribSystem>("FVFluxKernel")
     112           0 :               .queryInto(fv_flux_kernels);
     113           0 :           if (fv_flux_kernels.size())
     114           0 :             continue;
     115           0 :         }
     116             :         else
     117             :         {
     118             :           // We are a finite element variable
     119         287 :           if (_nl->getKernelWarehouse()
     120         287 :                   .getMatrixTagObjectWarehouse(mass_matrix_tag_id, 0)
     121         287 :                   .hasVariableObjects(var_id))
     122         267 :             continue;
     123          20 :           if (_nl->getDGKernelWarehouse()
     124          20 :                   .getMatrixTagObjectWarehouse(mass_matrix_tag_id, 0)
     125          20 :                   .hasVariableObjects(var_id))
     126           0 :             continue;
     127          20 :           if (_nl->getNodalKernelWarehouse()
     128          20 :                   .getMatrixTagObjectWarehouse(mass_matrix_tag_id, 0)
     129          20 :                   .hasVariableObjects(var_id))
     130           0 :             continue;
     131             : #ifdef MOOSE_KOKKOS_ENABLED
     132          19 :           if (_nl->getKokkosKernelWarehouse()
     133          19 :                   .getMatrixTagObjectWarehouse(mass_matrix_tag_id, 0)
     134          19 :                   .hasVariableObjects(var_id))
     135          17 :             continue;
     136           2 :           if (_nl->getKokkosNodalKernelWarehouse()
     137           2 :                   .getMatrixTagObjectWarehouse(mass_matrix_tag_id, 0)
     138           2 :                   .hasVariableObjects(var_id))
     139           0 :             continue;
     140             : #endif
     141             :         }
     142             :       }
     143          51 :       else if (_nl->getScalarKernelWarehouse()
     144          51 :                    .getMatrixTagObjectWarehouse(mass_matrix_tag_id, 0)
     145          51 :                    .hasVariableObjects(var_id))
     146          51 :         continue;
     147             : 
     148           3 :       mooseError("No objects contributing to the mass matrix were found for variable '",
     149             :                  var_name,
     150             :                  "'. Did you, e.g., forget a time derivative term?");
     151             :     }
     152         323 :   }
     153         323 : }
     154             : 
     155             : void
     156         299 : ExplicitTimeIntegrator::init()
     157             : {
     158         299 :   if (_nl && _fe_problem.solverParams(_nl->number())._type != Moose::ST_LINEAR)
     159           3 :     mooseError(
     160             :         "The chosen time integrator requires 'solve_type = LINEAR' in the Executioner block.");
     161         296 : }
     162             : 
     163             : void
     164        4335 : ExplicitTimeIntegrator::preSolve()
     165             : {
     166        4335 : }
     167             : 
     168             : void
     169          44 : ExplicitTimeIntegrator::meshChanged()
     170             : {
     171          44 :   setupSolver();
     172          44 : }
     173             : 
     174             : void
     175         370 : ExplicitTimeIntegrator::setupSolver()
     176             : {
     177             :   // Can only be done after the system is initialized
     178         370 :   if (_ones)
     179             :   {
     180          86 :     if (_solve_type == LUMPED || _solve_type == LUMP_PRECONDITIONED)
     181          86 :       *_ones = 1.;
     182             :   }
     183             : 
     184         370 :   if (_solve_type == CONSISTENT || _solve_type == LUMP_PRECONDITIONED)
     185         308 :     _linear_solver = LinearSolver<Number>::build(comm());
     186             : 
     187         370 :   if (_solve_type == LUMP_PRECONDITIONED)
     188             :   {
     189          24 :     _preconditioner = std::make_unique<LumpedPreconditioner>(*_mass_matrix_diag_inverted);
     190          24 :     _linear_solver->attach_preconditioner(_preconditioner.get());
     191          24 :     _linear_solver->init();
     192             :   }
     193             : 
     194         370 :   if (_solve_type == CONSISTENT || _solve_type == LUMP_PRECONDITIONED)
     195         308 :     Moose::PetscSupport::setLinearSolverDefaults(_fe_problem, *_linear_solver);
     196         370 : }
     197             : 
     198             : bool
     199        4467 : ExplicitTimeIntegrator::performExplicitSolve(SparseMatrix<Number> & mass_matrix)
     200             : {
     201        4467 :   bool converged = false;
     202             : 
     203        4467 :   switch (_solve_type)
     204             :   {
     205        1755 :     case CONSISTENT:
     206             :     {
     207        1755 :       converged = solveLinearSystem(mass_matrix);
     208             : 
     209        1755 :       break;
     210             :     }
     211        2465 :     case LUMPED:
     212             :     {
     213             :       // Computes the sum of each row (lumping)
     214             :       // Note: This is actually how PETSc does it
     215             :       // It's not "perfectly optimal" - but it will be fast (and universal)
     216        2465 :       mass_matrix.vector_mult(*_mass_matrix_diag_inverted, *_ones);
     217             : 
     218             :       // "Invert" the diagonal mass matrix
     219        2465 :       _mass_matrix_diag_inverted->reciprocal();
     220             : 
     221             :       // Multiply the inversion by the RHS
     222        2465 :       _solution_update->pointwise_mult(*_mass_matrix_diag_inverted, *_explicit_residual);
     223             : 
     224             :       // Check for convergence by seeing if there is a nan or inf
     225        2465 :       auto sum = _solution_update->sum();
     226        2465 :       converged = std::isfinite(sum);
     227             : 
     228             :       // The linear iteration count remains zero
     229        2465 :       _n_linear_iterations = 0;
     230             : 
     231        2465 :       break;
     232             :     }
     233         247 :     case LUMP_PRECONDITIONED:
     234             :     {
     235         247 :       mass_matrix.vector_mult(*_mass_matrix_diag_inverted, *_ones);
     236         247 :       _mass_matrix_diag_inverted->reciprocal();
     237             : 
     238         247 :       converged = solveLinearSystem(mass_matrix);
     239             : 
     240         247 :       break;
     241             :     }
     242           0 :     default:
     243           0 :       mooseError("Unknown solve_type in ExplicitTimeIntegrator.");
     244             :   }
     245             : 
     246        4467 :   return converged;
     247             : }
     248             : 
     249             : bool
     250        2002 : ExplicitTimeIntegrator::solveLinearSystem(SparseMatrix<Number> & mass_matrix)
     251             : {
     252        2002 :   auto & es = _fe_problem.es();
     253             : 
     254             :   const auto num_its_and_final_tol =
     255        4906 :       _linear_solver->solve(mass_matrix,
     256        2002 :                             *_solution_update,
     257        2002 :                             *_explicit_residual,
     258        1452 :                             es.parameters.get<Real>("linear solver tolerance"),
     259        1452 :                             es.parameters.get<unsigned int>("linear solver maximum iterations"));
     260             : 
     261        2002 :   _n_linear_iterations += num_its_and_final_tol.first;
     262             : 
     263        2002 :   const bool converged = checkLinearConvergence();
     264             : 
     265        2002 :   return converged;
     266             : }
     267             : 
     268             : bool
     269        2002 : ExplicitTimeIntegrator::checkLinearConvergence()
     270             : {
     271        2002 :   auto reason = _linear_solver->get_converged_reason();
     272             : 
     273        2002 :   switch (reason)
     274             :   {
     275        1894 :     case CONVERGED_RTOL_NORMAL:
     276             :     case CONVERGED_ATOL_NORMAL:
     277             :     case CONVERGED_RTOL:
     278             :     case CONVERGED_ATOL:
     279             :     case CONVERGED_ITS:
     280             :     case CONVERGED_CG_NEG_CURVE:
     281             :     case CONVERGED_CG_CONSTRAINED:
     282             :     case CONVERGED_STEP_LENGTH:
     283             :     case CONVERGED_HAPPY_BREAKDOWN:
     284        1894 :       return true;
     285         108 :     case DIVERGED_NULL:
     286             :     case DIVERGED_ITS:
     287             :     case DIVERGED_DTOL:
     288             :     case DIVERGED_BREAKDOWN:
     289             :     case DIVERGED_BREAKDOWN_BICG:
     290             :     case DIVERGED_NONSYMMETRIC:
     291             :     case DIVERGED_INDEFINITE_PC:
     292             :     case DIVERGED_NAN:
     293             :     case DIVERGED_INDEFINITE_MAT:
     294             :     case CONVERGED_ITERATING:
     295             :     case DIVERGED_PCSETUP_FAILED:
     296         108 :       return false;
     297           0 :     default:
     298           0 :       mooseError("Unknown convergence reason in ExplicitTimeIntegrator.");
     299             :   }
     300             : }