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 "SteffensenSolve.h"
      11             : 
      12             : #include "Executioner.h"
      13             : #include "FEProblemBase.h"
      14             : #include "NonlinearSystem.h"
      15             : #include "AllLocalDofIndicesThread.h"
      16             : #include "Console.h"
      17             : #include "DefaultMultiAppFixedPointConvergence.h"
      18             : 
      19             : InputParameters
      20           0 : SteffensenSolve::validParams()
      21             : {
      22           0 :   InputParameters params = FixedPointSolve::validParams();
      23             : 
      24           0 :   return params;
      25             : }
      26             : 
      27         576 : SteffensenSolve::SteffensenSolve(Executioner & ex) : FixedPointSolve(ex) { allocateStorage(true); }
      28             : 
      29             : void
      30         864 : SteffensenSolve::allocateStorage(const bool primary)
      31             : {
      32             :   TagID fxn_m1_tagid;
      33             :   TagID xn_m1_tagid;
      34             :   const std::vector<PostprocessorName> * transformed_pps;
      35             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
      36         864 :   if (primary)
      37             :   {
      38         576 :     xn_m1_tagid = _problem.addVectorTag("xn_m1", Moose::VECTOR_TAG_SOLUTION);
      39         576 :     fxn_m1_tagid = _problem.addVectorTag("fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      40         576 :     transformed_pps = &_transformed_pps;
      41         576 :     transformed_pps_values = &_transformed_pps_values;
      42         576 :     _xn_m1_tagid = xn_m1_tagid;
      43         576 :     _fxn_m1_tagid = fxn_m1_tagid;
      44             :   }
      45             :   else
      46             :   {
      47         288 :     xn_m1_tagid = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION);
      48         288 :     fxn_m1_tagid = _problem.addVectorTag("secondary_fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      49         288 :     transformed_pps = &_secondary_transformed_pps;
      50         288 :     transformed_pps_values = &_secondary_transformed_pps_values;
      51         288 :     _secondary_xn_m1_tagid = xn_m1_tagid;
      52         288 :     _secondary_fxn_m1_tagid = fxn_m1_tagid;
      53             :   }
      54             : 
      55             :   // Store a copy of the previous solution here
      56         864 :   _solver_sys.addVector(xn_m1_tagid, false, PARALLEL);
      57         864 :   _solver_sys.addVector(fxn_m1_tagid, false, PARALLEL);
      58             : 
      59             :   // Allocate storage for the previous postprocessor values
      60         864 :   (*transformed_pps_values).resize((*transformed_pps).size());
      61        1032 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
      62         168 :     (*transformed_pps_values)[i].resize(2);
      63         864 : }
      64             : 
      65             : void
      66           0 : SteffensenSolve::initialSetup()
      67             : {
      68           0 :   FixedPointSolve::initialSetup();
      69             : 
      70           0 :   auto & convergence = _problem.getConvergence(_problem.getMultiAppFixedPointConvergenceName());
      71           0 :   if (!dynamic_cast<DefaultMultiAppFixedPointConvergence *>(&convergence))
      72           0 :     mooseError(
      73             :         "Only DefaultMultiAppFixedPointConvergence objects may be used for "
      74             :         "'multiapp_fixed_point_convergence' when using the Steffensen fixed point algorithm.");
      75           0 : }
      76             : 
      77             : void
      78       31001 : SteffensenSolve::saveVariableValues(const bool primary)
      79             : {
      80             :   unsigned int iteration;
      81             :   TagID fxn_m1_tagid;
      82             :   TagID xn_m1_tagid;
      83       31001 :   if (primary)
      84             :   {
      85       22514 :     iteration = _fixed_point_it;
      86       22514 :     fxn_m1_tagid = _fxn_m1_tagid;
      87       22514 :     xn_m1_tagid = _xn_m1_tagid;
      88             :   }
      89             :   else
      90             :   {
      91        8487 :     iteration = _main_fixed_point_it;
      92        8487 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
      93        8487 :     xn_m1_tagid = _secondary_xn_m1_tagid;
      94             :   }
      95             : 
      96             :   // Save previous variable values
      97       31001 :   NumericVector<Number> & solution = _solver_sys.solution();
      98       31001 :   NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
      99       31001 :   NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
     100             : 
     101             :   // What 'solution' is with regards to the Steffensen solve depends on the step
     102       31001 :   if (iteration % 2 == 1)
     103        9736 :     xn_m1 = solution;
     104             :   else
     105       21265 :     fxn_m1 = solution;
     106       31001 : }
     107             : 
     108             : void
     109       31001 : SteffensenSolve::savePostprocessorValues(const bool primary)
     110             : {
     111             :   unsigned int iteration;
     112             :   const std::vector<PostprocessorName> * transformed_pps;
     113             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     114       31001 :   if (primary)
     115             :   {
     116       22514 :     iteration = _fixed_point_it;
     117       22514 :     transformed_pps = &_transformed_pps;
     118       22514 :     transformed_pps_values = &_transformed_pps_values;
     119             :   }
     120             :   else
     121             :   {
     122        8487 :     iteration = _main_fixed_point_it;
     123        8487 :     transformed_pps = &_secondary_transformed_pps;
     124        8487 :     transformed_pps_values = &_secondary_transformed_pps_values;
     125             :   }
     126             : 
     127             :   // Save previous postprocessor values
     128       36101 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     129             :   {
     130        5100 :     if (iteration % 2 == 0)
     131        2600 :       (*transformed_pps_values)[i][1] = getPostprocessorValueByName((*transformed_pps)[i]);
     132             :     else
     133        2500 :       (*transformed_pps_values)[i][0] = getPostprocessorValueByName((*transformed_pps)[i]);
     134             :   }
     135       31001 : }
     136             : 
     137             : bool
     138       11647 : SteffensenSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary)
     139             : {
     140             :   // Need at least two values to compute the Steffensen update, and the update is only performed
     141             :   // every other iteration as two evaluations of the coupled problem are necessary
     142       11647 :   if (primary)
     143        6837 :     return _fixed_point_it > 1 && (_fixed_point_it % 2 == 0);
     144             :   else
     145        4810 :     return _main_fixed_point_it > 1 && (_main_fixed_point_it % 2 == 0);
     146             : }
     147             : 
     148             : void
     149        1820 : SteffensenSolve::transformPostprocessors(const bool primary)
     150             : {
     151             :   Real relaxation_factor;
     152             :   const std::vector<PostprocessorName> * transformed_pps;
     153             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     154        1820 :   if (primary)
     155             :   {
     156        1060 :     relaxation_factor = _relax_factor;
     157        1060 :     transformed_pps = &_transformed_pps;
     158        1060 :     transformed_pps_values = &_transformed_pps_values;
     159             :   }
     160             :   else
     161             :   {
     162         760 :     relaxation_factor = _secondary_relaxation_factor;
     163         760 :     transformed_pps = &_secondary_transformed_pps;
     164         760 :     transformed_pps_values = &_secondary_transformed_pps_values;
     165             :   }
     166             : 
     167             :   // Relax postprocessors for the main application
     168        3640 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     169             :   {
     170             :     // Get new postprocessor value
     171        1820 :     const Real current_value = getPostprocessorValueByName((*transformed_pps)[i]);
     172        1820 :     const Real fxn_m1 = (*transformed_pps_values)[i][0];
     173        1820 :     const Real xn_m1 = (*transformed_pps_values)[i][1];
     174             : 
     175             :     // Compute and set relaxed value
     176        1820 :     Real new_value = current_value;
     177        1820 :     if (!MooseUtils::absoluteFuzzyEqual(current_value + xn_m1 - 2 * fxn_m1, 0))
     178        1820 :       new_value =
     179        1820 :           xn_m1 - (fxn_m1 - xn_m1) * (fxn_m1 - xn_m1) / (current_value + xn_m1 - 2 * fxn_m1);
     180             : 
     181             :     // Relax update
     182        1820 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * fxn_m1;
     183             : 
     184        1820 :     _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
     185             :   }
     186        1820 : }
     187             : 
     188             : void
     189        1849 : SteffensenSolve::transformVariables(const std::set<dof_id_type> & transformed_dofs,
     190             :                                     const bool primary)
     191             : {
     192             :   Real relaxation_factor;
     193             :   TagID fxn_m1_tagid;
     194             :   TagID xn_m1_tagid;
     195        1849 :   if (primary)
     196             :   {
     197        1211 :     relaxation_factor = _relax_factor;
     198        1211 :     fxn_m1_tagid = _fxn_m1_tagid;
     199        1211 :     xn_m1_tagid = _xn_m1_tagid;
     200             :   }
     201             :   else
     202             :   {
     203         638 :     relaxation_factor = _secondary_relaxation_factor;
     204         638 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
     205         638 :     xn_m1_tagid = _secondary_xn_m1_tagid;
     206             :   }
     207             : 
     208        1849 :   NumericVector<Number> & solution = _solver_sys.solution();
     209        1849 :   NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
     210        1849 :   NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
     211             : 
     212      163868 :   for (const auto & dof : transformed_dofs)
     213             :   {
     214             :     // Avoid 0 denominator issue
     215      162019 :     Real new_value = solution(dof);
     216      162019 :     if (!MooseUtils::absoluteFuzzyEqual(solution(dof) + xn_m1(dof) - 2 * fxn_m1(dof), 0))
     217      132561 :       new_value = xn_m1(dof) - (fxn_m1(dof) - xn_m1(dof)) * (fxn_m1(dof) - xn_m1(dof)) /
     218      132561 :                                    (solution(dof) + xn_m1(dof) - 2 * fxn_m1(dof));
     219             : 
     220             :     // Relax update
     221      162019 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * fxn_m1(dof);
     222             : 
     223      162019 :     solution.set(dof, new_value);
     224             :   }
     225        1849 :   solution.close();
     226        1849 :   _solver_sys.update();
     227        1849 : }
     228             : 
     229             : void
     230       10977 : SteffensenSolve::printFixedPointConvergenceHistory(
     231             :     const Real initial_norm,
     232             :     const std::vector<Real> & timestep_begin_norms,
     233             :     const std::vector<Real> & timestep_end_norms) const
     234             : {
     235       10977 :   _console << "\n 0 Steffensen initialization |R| = "
     236       10977 :            << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
     237             : 
     238       10977 :   Real max_norm_old = initial_norm;
     239       50106 :   for (unsigned int i = 0; i <= _fixed_point_it; ++i)
     240             :   {
     241       39129 :     Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
     242       39129 :     std::stringstream steffensen_prefix;
     243       39129 :     if (i == 0)
     244       10977 :       steffensen_prefix << " Steffensen initialization |R| = ";
     245       28152 :     else if (i % 2 == 0)
     246       11632 :       steffensen_prefix << " Steffensen half-step      |R| = ";
     247             :     else
     248       16520 :       steffensen_prefix << " Steffensen step           |R| = ";
     249             : 
     250       39129 :     _console << std::setw(2) << i + 1 << steffensen_prefix.str()
     251       39129 :              << Console::outputNorm(max_norm_old, max_norm) << '\n';
     252             : 
     253       39129 :     max_norm_old = max_norm;
     254       39129 :   }
     255       10977 : }