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 "PicardSolve.h"
      11             : 
      12             : #include "Executioner.h"
      13             : #include "FEProblemBase.h"
      14             : #include "NonlinearSystem.h"
      15             : #include "AllLocalDofIndicesThread.h"
      16             : #include "Console.h"
      17             : 
      18             : InputParameters
      19           0 : PicardSolve::validParams()
      20             : {
      21           0 :   InputParameters params = FixedPointSolve::validParams();
      22           0 :   return params;
      23             : }
      24             : 
      25       59923 : PicardSolve::PicardSolve(Executioner & ex) : FixedPointSolve(ex) {}
      26             : 
      27             : void
      28       67973 : PicardSolve::allocateStorage(const bool primary)
      29             : {
      30       67973 :   if (!performingRelaxation(primary))
      31       67770 :     return;
      32         203 :   findTransformedSystem(primary);
      33             : 
      34             :   const std::vector<PostprocessorName> * transformed_pps;
      35             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
      36         203 :   if (primary)
      37             :   {
      38         123 :     if (_transformed_sys)
      39             :     {
      40          24 :       _old_tag_id =
      41          24 :           _problem.addVectorTag(Moose::PREVIOUS_FP_SOLUTION_TAG, Moose::VECTOR_TAG_SOLUTION);
      42          24 :       _transformed_sys->needSolutionState(1, Moose::SolutionIterationType::FixedPoint, PARALLEL);
      43             :     }
      44         123 :     transformed_pps = &_transformed_pps;
      45         123 :     transformed_pps_values = &_transformed_pps_values;
      46             :   }
      47             :   else
      48             :   {
      49          80 :     if (_transformed_sys)
      50             :     {
      51          12 :       _secondary_old_tag_id = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION);
      52          12 :       _transformed_sys->addVector(_secondary_old_tag_id, false, PARALLEL);
      53             :     }
      54             : 
      55          80 :     transformed_pps = &_secondary_transformed_pps;
      56          80 :     transformed_pps_values = &_secondary_transformed_pps_values;
      57             :   }
      58             : 
      59             :   // Allocate storage for the previous postprocessor values
      60         203 :   (*transformed_pps_values).resize((*transformed_pps).size());
      61         370 :   for (const auto i : index_range(*transformed_pps))
      62         167 :     (*transformed_pps_values)[i].resize(1);
      63             : }
      64             : 
      65             : void
      66        1438 : PicardSolve::saveVariableValues(const bool primary)
      67             : {
      68             :   // Primary is copied back by _transformed_sys->copyPreviousFixedPointSolutions()
      69        1438 :   if (!performingRelaxation(primary) || primary)
      70        1146 :     return;
      71             : 
      72             :   // Check to make sure allocateStorage has been called
      73             :   mooseAssert(_secondary_old_tag_id != Moose::INVALID_TAG_ID,
      74             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
      75             : 
      76             :   // Save variable previous values
      77         292 :   NumericVector<Number> & solution = _transformed_sys->solution();
      78         292 :   NumericVector<Number> & transformed_old = _transformed_sys->getVector(_secondary_old_tag_id);
      79         292 :   transformed_old = solution;
      80             : }
      81             : 
      82             : void
      83      279243 : PicardSolve::savePostprocessorValues(const bool primary)
      84             : {
      85      279243 :   if (!performingRelaxation(primary))
      86      266355 :     return;
      87             : 
      88             :   const std::vector<PostprocessorName> * transformed_pps;
      89             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
      90       12888 :   if (primary)
      91             :   {
      92        9357 :     transformed_pps = &_transformed_pps;
      93        9357 :     transformed_pps_values = &_transformed_pps_values;
      94             :   }
      95             :   else
      96             :   {
      97        3531 :     transformed_pps = &_secondary_transformed_pps;
      98        3531 :     transformed_pps_values = &_secondary_transformed_pps_values;
      99             :   }
     100             : 
     101             :   // Save postprocessor previous values
     102       24674 :   for (const auto i : index_range(*transformed_pps))
     103       11786 :     (*transformed_pps_values)[i][0] = getPostprocessorValueByName((*transformed_pps)[i]);
     104             : }
     105             : 
     106             : bool
     107       14482 : PicardSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary)
     108             : {
     109             :   // unrelaxed Picard is the default update for fixed point iterations
     110             :   // old values are required for relaxation
     111       14482 :   const auto fixed_point_it = primary ? _fixed_point_it : _main_fixed_point_it;
     112       14482 :   return performingRelaxation(primary) && fixed_point_it > 0;
     113             : }
     114             : 
     115             : void
     116       10886 : PicardSolve::transformPostprocessors(const bool primary)
     117             : {
     118             :   Real relaxation_factor;
     119             :   const std::vector<PostprocessorName> * transformed_pps;
     120             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     121       10886 :   if (primary)
     122             :   {
     123        7652 :     relaxation_factor = _relax_factor;
     124        7652 :     transformed_pps = &_transformed_pps;
     125        7652 :     transformed_pps_values = &_transformed_pps_values;
     126             :   }
     127             :   else
     128             :   {
     129        3234 :     relaxation_factor = _secondary_relaxation_factor;
     130        3234 :     transformed_pps = &_secondary_transformed_pps;
     131        3234 :     transformed_pps_values = &_secondary_transformed_pps_values;
     132             :   }
     133             : 
     134             :   // Relax the postprocessors
     135       21772 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     136             :   {
     137             :     // Get new postprocessor value
     138       10886 :     const Real current_value = getPostprocessorValueByName((*transformed_pps)[i]);
     139       10886 :     const Real old_value = (*transformed_pps_values)[i][0];
     140             : 
     141             :     // Compute and set relaxed value
     142       10886 :     Real new_value = current_value;
     143       10886 :     new_value = relaxation_factor * current_value + (1 - relaxation_factor) * old_value;
     144       10886 :     _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
     145             :   }
     146       10886 : }
     147             : 
     148             : void
     149        1011 : PicardSolve::transformVariables(const std::set<dof_id_type> & transformed_dofs, const bool primary)
     150             : {
     151             :   Real relaxation_factor;
     152             :   TagID old_tag_id;
     153        1011 :   if (primary)
     154             :   {
     155         719 :     relaxation_factor = _relax_factor;
     156         719 :     old_tag_id = _old_tag_id;
     157             :   }
     158             :   else
     159             :   {
     160         292 :     relaxation_factor = _secondary_relaxation_factor;
     161         292 :     old_tag_id = _secondary_old_tag_id;
     162             :   }
     163             : 
     164        1011 :   NumericVector<Number> & solution = _transformed_sys->solution();
     165        1011 :   NumericVector<Number> & transformed_old = _transformed_sys->getVector(old_tag_id);
     166             : 
     167      129634 :   for (const auto & dof : transformed_dofs)
     168      257246 :     solution.set(dof,
     169      128623 :                  (transformed_old(dof) * (1.0 - relaxation_factor)) +
     170      128623 :                      (solution(dof) * relaxation_factor));
     171             : 
     172        1011 :   solution.close();
     173        1011 :   _transformed_sys->update();
     174        1011 : }
     175             : 
     176             : void
     177       30582 : PicardSolve::printFixedPointConvergenceHistory(Real initial_norm,
     178             :                                                const std::vector<Real> & timestep_begin_norms,
     179             :                                                const std::vector<Real> & timestep_end_norms) const
     180             : {
     181       30582 :   _console << "\n 0 Picard |R| = "
     182       30582 :            << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
     183             : 
     184       30582 :   Real max_norm_old = initial_norm;
     185      159345 :   for (unsigned int i = 0; i <= _fixed_point_it; ++i)
     186             :   {
     187      128763 :     Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
     188      128763 :     _console << std::setw(2) << i + 1
     189      128763 :              << " Picard |R| = " << Console::outputNorm(max_norm_old, max_norm) << '\n';
     190      128763 :     max_norm_old = max_norm;
     191             :   }
     192             : 
     193       30582 :   _console << std::endl;
     194       30582 : }