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         574 : SteffensenSolve::SteffensenSolve(Executioner & ex) : FixedPointSolve(ex) {}
      28             : 
      29             : void
      30         861 : SteffensenSolve::allocateStorage(const bool primary)
      31             : {
      32         861 :   findTransformedSystem(primary);
      33             : 
      34             :   TagID fxn_m1_tagid;
      35             :   TagID xn_m1_tagid;
      36             :   const std::vector<PostprocessorName> * transformed_pps;
      37             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
      38         861 :   if (primary)
      39             :   {
      40         574 :     xn_m1_tagid = _problem.addVectorTag("xn_m1", Moose::VECTOR_TAG_SOLUTION);
      41         574 :     fxn_m1_tagid = _problem.addVectorTag("fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      42         574 :     transformed_pps = &_transformed_pps;
      43         574 :     transformed_pps_values = &_transformed_pps_values;
      44         574 :     _xn_m1_tagid = xn_m1_tagid;
      45         574 :     _fxn_m1_tagid = fxn_m1_tagid;
      46             :   }
      47             :   else
      48             :   {
      49         287 :     xn_m1_tagid = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION);
      50         287 :     fxn_m1_tagid = _problem.addVectorTag("secondary_fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      51         287 :     transformed_pps = &_secondary_transformed_pps;
      52         287 :     transformed_pps_values = &_secondary_transformed_pps_values;
      53         287 :     _secondary_xn_m1_tagid = xn_m1_tagid;
      54         287 :     _secondary_fxn_m1_tagid = fxn_m1_tagid;
      55             :   }
      56             : 
      57             :   // Store a copy of the previous solution here
      58             :   // If we don't have a transformed system, we are not accelerating variables
      59         861 :   if (_transformed_sys)
      60             :   {
      61         208 :     _transformed_sys->addVector(xn_m1_tagid, false, PARALLEL);
      62         208 :     _transformed_sys->addVector(fxn_m1_tagid, false, PARALLEL);
      63             :   }
      64             : 
      65             :   // Allocate storage for the previous postprocessor values
      66         861 :   (*transformed_pps_values).resize((*transformed_pps).size());
      67        1028 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
      68         167 :     (*transformed_pps_values)[i].resize(2);
      69         861 : }
      70             : 
      71             : void
      72         574 : SteffensenSolve::initialSetup()
      73             : {
      74         574 :   FixedPointSolve::initialSetup();
      75             : 
      76         574 :   auto & convergence = _problem.getConvergence(_problem.getMultiAppFixedPointConvergenceName());
      77         574 :   if (!dynamic_cast<DefaultMultiAppFixedPointConvergence *>(&convergence))
      78           0 :     mooseError(
      79             :         "Only DefaultMultiAppFixedPointConvergence objects may be used for "
      80             :         "'multiapp_fixed_point_convergence' when using the Steffensen fixed point algorithm.");
      81         574 : }
      82             : 
      83             : void
      84        7312 : SteffensenSolve::saveVariableValues(const bool primary)
      85             : {
      86             :   unsigned int iteration;
      87             :   TagID fxn_m1_tagid;
      88             :   TagID xn_m1_tagid;
      89        7312 :   if (primary)
      90             :   {
      91        5758 :     iteration = _fixed_point_it;
      92        5758 :     fxn_m1_tagid = _fxn_m1_tagid;
      93        5758 :     xn_m1_tagid = _xn_m1_tagid;
      94             :   }
      95             :   else
      96             :   {
      97        1554 :     iteration = _main_fixed_point_it;
      98        1554 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
      99        1554 :     xn_m1_tagid = _secondary_xn_m1_tagid;
     100             :   }
     101             : 
     102             :   // Check to make sure allocateStorage has been called
     103             :   mooseAssert(fxn_m1_tagid != Moose::INVALID_TAG_ID,
     104             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
     105             :   mooseAssert(xn_m1_tagid != Moose::INVALID_TAG_ID,
     106             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
     107             : 
     108             :   // Save previous variable values
     109        7312 :   NumericVector<Number> & solution = _transformed_sys->solution();
     110        7312 :   NumericVector<Number> & fxn_m1 = _transformed_sys->getVector(fxn_m1_tagid);
     111        7312 :   NumericVector<Number> & xn_m1 = _transformed_sys->getVector(xn_m1_tagid);
     112             : 
     113             :   // What 'solution' is with regards to the Steffensen solve depends on the step
     114        7312 :   if (iteration % 2 == 1)
     115        2420 :     xn_m1 = solution;
     116             :   else
     117        4892 :     fxn_m1 = solution;
     118        7312 : }
     119             : 
     120             : void
     121       31203 : SteffensenSolve::savePostprocessorValues(const bool primary)
     122             : {
     123             :   unsigned int iteration;
     124             :   const std::vector<PostprocessorName> * transformed_pps;
     125             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     126       31203 :   if (primary)
     127             :   {
     128       22628 :     iteration = _fixed_point_it;
     129       22628 :     transformed_pps = &_transformed_pps;
     130       22628 :     transformed_pps_values = &_transformed_pps_values;
     131             :   }
     132             :   else
     133             :   {
     134        8575 :     iteration = _main_fixed_point_it;
     135        8575 :     transformed_pps = &_secondary_transformed_pps;
     136        8575 :     transformed_pps_values = &_secondary_transformed_pps_values;
     137             :   }
     138             : 
     139             :   // Save previous postprocessor values
     140       36309 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     141             :   {
     142        5106 :     if (iteration % 2 == 0)
     143        2605 :       (*transformed_pps_values)[i][1] = getPostprocessorValueByName((*transformed_pps)[i]);
     144             :     else
     145        2501 :       (*transformed_pps_values)[i][0] = getPostprocessorValueByName((*transformed_pps)[i]);
     146             :   }
     147       31203 : }
     148             : 
     149             : bool
     150       11704 : SteffensenSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary)
     151             : {
     152             :   // Need at least two values to compute the Steffensen update, and the update is only performed
     153             :   // every other iteration as two evaluations of the coupled problem are necessary
     154       11704 :   if (primary)
     155        6857 :     return _fixed_point_it > 1 && (_fixed_point_it % 2 == 0);
     156             :   else
     157        4847 :     return _main_fixed_point_it > 1 && (_main_fixed_point_it % 2 == 0);
     158             : }
     159             : 
     160             : void
     161        1825 : SteffensenSolve::transformPostprocessors(const bool primary)
     162             : {
     163             :   Real relaxation_factor;
     164             :   const std::vector<PostprocessorName> * transformed_pps;
     165             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     166        1825 :   if (primary)
     167             :   {
     168        1059 :     relaxation_factor = _relax_factor;
     169        1059 :     transformed_pps = &_transformed_pps;
     170        1059 :     transformed_pps_values = &_transformed_pps_values;
     171             :   }
     172             :   else
     173             :   {
     174         766 :     relaxation_factor = _secondary_relaxation_factor;
     175         766 :     transformed_pps = &_secondary_transformed_pps;
     176         766 :     transformed_pps_values = &_secondary_transformed_pps_values;
     177             :   }
     178             : 
     179             :   // Relax postprocessors for the main application
     180        3650 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     181             :   {
     182             :     // Get new postprocessor value
     183        1825 :     const Real current_value = getPostprocessorValueByName((*transformed_pps)[i]);
     184        1825 :     const Real fxn_m1 = (*transformed_pps_values)[i][0];
     185        1825 :     const Real xn_m1 = (*transformed_pps_values)[i][1];
     186             : 
     187             :     // Compute and set relaxed value
     188        1825 :     Real new_value = current_value;
     189        1825 :     if (!MooseUtils::absoluteFuzzyEqual(current_value + xn_m1 - 2 * fxn_m1, 0))
     190        1825 :       new_value =
     191        1825 :           xn_m1 - (fxn_m1 - xn_m1) * (fxn_m1 - xn_m1) / (current_value + xn_m1 - 2 * fxn_m1);
     192             : 
     193             :     // Relax update
     194        1825 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * fxn_m1;
     195             : 
     196        1825 :     _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
     197             :   }
     198        1825 : }
     199             : 
     200             : void
     201        1867 : SteffensenSolve::transformVariables(const std::set<dof_id_type> & transformed_dofs,
     202             :                                     const bool primary)
     203             : {
     204             :   Real relaxation_factor;
     205             :   TagID fxn_m1_tagid;
     206             :   TagID xn_m1_tagid;
     207        1867 :   if (primary)
     208             :   {
     209        1223 :     relaxation_factor = _relax_factor;
     210        1223 :     fxn_m1_tagid = _fxn_m1_tagid;
     211        1223 :     xn_m1_tagid = _xn_m1_tagid;
     212             :   }
     213             :   else
     214             :   {
     215         644 :     relaxation_factor = _secondary_relaxation_factor;
     216         644 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
     217         644 :     xn_m1_tagid = _secondary_xn_m1_tagid;
     218             :   }
     219             : 
     220        1867 :   NumericVector<Number> & solution = _transformed_sys->solution();
     221        1867 :   NumericVector<Number> & fxn_m1 = _transformed_sys->getVector(fxn_m1_tagid);
     222        1867 :   NumericVector<Number> & xn_m1 = _transformed_sys->getVector(xn_m1_tagid);
     223             : 
     224      166064 :   for (const auto & dof : transformed_dofs)
     225             :   {
     226             :     // Avoid 0 denominator issue
     227      164197 :     Real new_value = solution(dof);
     228      164197 :     if (!MooseUtils::absoluteFuzzyEqual(solution(dof) + xn_m1(dof) - 2 * fxn_m1(dof), 0))
     229      134343 :       new_value = xn_m1(dof) - (fxn_m1(dof) - xn_m1(dof)) * (fxn_m1(dof) - xn_m1(dof)) /
     230      134343 :                                    (solution(dof) + xn_m1(dof) - 2 * fxn_m1(dof));
     231             : 
     232             :     // Relax update
     233      164197 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * fxn_m1(dof);
     234             : 
     235      164197 :     solution.set(dof, new_value);
     236             :   }
     237        1867 :   solution.close();
     238        1867 :   _transformed_sys->update();
     239        1867 : }
     240             : 
     241             : void
     242       11034 : SteffensenSolve::printFixedPointConvergenceHistory(
     243             :     const Real initial_norm,
     244             :     const std::vector<Real> & timestep_begin_norms,
     245             :     const std::vector<Real> & timestep_end_norms) const
     246             : {
     247       11034 :   _console << "\n 0 Steffensen initialization |R| = "
     248       11034 :            << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
     249             : 
     250       11034 :   Real max_norm_old = initial_norm;
     251       50415 :   for (unsigned int i = 0; i <= _fixed_point_it; ++i)
     252             :   {
     253       39381 :     Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
     254       39381 :     std::stringstream steffensen_prefix;
     255       39381 :     if (i == 0)
     256       11034 :       steffensen_prefix << " Steffensen initialization |R| = ";
     257       28347 :     else if (i % 2 == 0)
     258       11718 :       steffensen_prefix << " Steffensen half-step      |R| = ";
     259             :     else
     260       16629 :       steffensen_prefix << " Steffensen step           |R| = ";
     261             : 
     262       39381 :     _console << std::setw(2) << i + 1 << steffensen_prefix.str()
     263       39381 :              << Console::outputNorm(max_norm_old, max_norm) << '\n';
     264             : 
     265       39381 :     max_norm_old = max_norm;
     266       39381 :   }
     267       11034 : }