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         624 : SteffensenSolve::SteffensenSolve(Executioner & ex) : FixedPointSolve(ex) { allocateStorage(true); }
      28             : 
      29             : void
      30         936 : 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         936 :   if (primary)
      37             :   {
      38         624 :     xn_m1_tagid = _problem.addVectorTag("xn_m1", Moose::VECTOR_TAG_SOLUTION);
      39         624 :     fxn_m1_tagid = _problem.addVectorTag("fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      40         624 :     transformed_pps = &_transformed_pps;
      41         624 :     transformed_pps_values = &_transformed_pps_values;
      42         624 :     _xn_m1_tagid = xn_m1_tagid;
      43         624 :     _fxn_m1_tagid = fxn_m1_tagid;
      44             :   }
      45             :   else
      46             :   {
      47         312 :     xn_m1_tagid = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION);
      48         312 :     fxn_m1_tagid = _problem.addVectorTag("secondary_fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      49         312 :     transformed_pps = &_secondary_transformed_pps;
      50         312 :     transformed_pps_values = &_secondary_transformed_pps_values;
      51         312 :     _secondary_xn_m1_tagid = xn_m1_tagid;
      52         312 :     _secondary_fxn_m1_tagid = fxn_m1_tagid;
      53             :   }
      54             : 
      55             :   // Store a copy of the previous solution here
      56         936 :   _solver_sys.addVector(xn_m1_tagid, false, PARALLEL);
      57         936 :   _solver_sys.addVector(fxn_m1_tagid, false, PARALLEL);
      58             : 
      59             :   // Allocate storage for the previous postprocessor values
      60         936 :   (*transformed_pps_values).resize((*transformed_pps).size());
      61        1118 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
      62         182 :     (*transformed_pps_values)[i].resize(2);
      63         936 : }
      64             : 
      65             : void
      66         624 : SteffensenSolve::initialSetup()
      67             : {
      68         624 :   FixedPointSolve::initialSetup();
      69             : 
      70         624 :   auto & convergence = _problem.getConvergence(_problem.getMultiAppFixedPointConvergenceName());
      71         624 :   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         624 : }
      76             : 
      77             : void
      78       34038 : SteffensenSolve::saveVariableValues(const bool primary)
      79             : {
      80             :   unsigned int iteration;
      81             :   TagID fxn_m1_tagid;
      82             :   TagID xn_m1_tagid;
      83       34038 :   if (primary)
      84             :   {
      85       24718 :     iteration = _fixed_point_it;
      86       24718 :     fxn_m1_tagid = _fxn_m1_tagid;
      87       24718 :     xn_m1_tagid = _xn_m1_tagid;
      88             :   }
      89             :   else
      90             :   {
      91        9320 :     iteration = _main_fixed_point_it;
      92        9320 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
      93        9320 :     xn_m1_tagid = _secondary_xn_m1_tagid;
      94             :   }
      95             : 
      96             :   // Check to make sure allocateStorage has been called
      97             :   mooseAssert(fxn_m1_tagid != Moose::INVALID_TAG_ID,
      98             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
      99             :   mooseAssert(xn_m1_tagid != Moose::INVALID_TAG_ID,
     100             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
     101             : 
     102             :   // Save previous variable values
     103       34038 :   NumericVector<Number> & solution = _solver_sys.solution();
     104       34038 :   NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
     105       34038 :   NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
     106             : 
     107             :   // What 'solution' is with regards to the Steffensen solve depends on the step
     108       34038 :   if (iteration % 2 == 1)
     109       10684 :     xn_m1 = solution;
     110             :   else
     111       23354 :     fxn_m1 = solution;
     112       34038 : }
     113             : 
     114             : void
     115       34080 : SteffensenSolve::savePostprocessorValues(const bool primary)
     116             : {
     117             :   unsigned int iteration;
     118             :   const std::vector<PostprocessorName> * transformed_pps;
     119             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     120       34080 :   if (primary)
     121             :   {
     122       24718 :     iteration = _fixed_point_it;
     123       24718 :     transformed_pps = &_transformed_pps;
     124       24718 :     transformed_pps_values = &_transformed_pps_values;
     125             :   }
     126             :   else
     127             :   {
     128        9362 :     iteration = _main_fixed_point_it;
     129        9362 :     transformed_pps = &_secondary_transformed_pps;
     130        9362 :     transformed_pps_values = &_secondary_transformed_pps_values;
     131             :   }
     132             : 
     133             :   // Save previous postprocessor values
     134       39663 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     135             :   {
     136        5583 :     if (iteration % 2 == 0)
     137        2848 :       (*transformed_pps_values)[i][1] = getPostprocessorValueByName((*transformed_pps)[i]);
     138             :     else
     139        2735 :       (*transformed_pps_values)[i][0] = getPostprocessorValueByName((*transformed_pps)[i]);
     140             :   }
     141       34080 : }
     142             : 
     143             : bool
     144       12788 : SteffensenSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary)
     145             : {
     146             :   // Need at least two values to compute the Steffensen update, and the update is only performed
     147             :   // every other iteration as two evaluations of the coupled problem are necessary
     148       12788 :   if (primary)
     149        7494 :     return _fixed_point_it > 1 && (_fixed_point_it % 2 == 0);
     150             :   else
     151        5294 :     return _main_fixed_point_it > 1 && (_main_fixed_point_it % 2 == 0);
     152             : }
     153             : 
     154             : void
     155        1996 : SteffensenSolve::transformPostprocessors(const bool primary)
     156             : {
     157             :   Real relaxation_factor;
     158             :   const std::vector<PostprocessorName> * transformed_pps;
     159             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     160        1996 :   if (primary)
     161             :   {
     162        1160 :     relaxation_factor = _relax_factor;
     163        1160 :     transformed_pps = &_transformed_pps;
     164        1160 :     transformed_pps_values = &_transformed_pps_values;
     165             :   }
     166             :   else
     167             :   {
     168         836 :     relaxation_factor = _secondary_relaxation_factor;
     169         836 :     transformed_pps = &_secondary_transformed_pps;
     170         836 :     transformed_pps_values = &_secondary_transformed_pps_values;
     171             :   }
     172             : 
     173             :   // Relax postprocessors for the main application
     174        3992 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     175             :   {
     176             :     // Get new postprocessor value
     177        1996 :     const Real current_value = getPostprocessorValueByName((*transformed_pps)[i]);
     178        1996 :     const Real fxn_m1 = (*transformed_pps_values)[i][0];
     179        1996 :     const Real xn_m1 = (*transformed_pps_values)[i][1];
     180             : 
     181             :     // Compute and set relaxed value
     182        1996 :     Real new_value = current_value;
     183        1996 :     if (!MooseUtils::absoluteFuzzyEqual(current_value + xn_m1 - 2 * fxn_m1, 0))
     184        1996 :       new_value =
     185        1996 :           xn_m1 - (fxn_m1 - xn_m1) * (fxn_m1 - xn_m1) / (current_value + xn_m1 - 2 * fxn_m1);
     186             : 
     187             :     // Relax update
     188        1996 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * fxn_m1;
     189             : 
     190        1996 :     _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
     191             :   }
     192        1996 : }
     193             : 
     194             : void
     195        2037 : SteffensenSolve::transformVariables(const std::set<dof_id_type> & transformed_dofs,
     196             :                                     const bool primary)
     197             : {
     198             :   Real relaxation_factor;
     199             :   TagID fxn_m1_tagid;
     200             :   TagID xn_m1_tagid;
     201        2037 :   if (primary)
     202             :   {
     203        1335 :     relaxation_factor = _relax_factor;
     204        1335 :     fxn_m1_tagid = _fxn_m1_tagid;
     205        1335 :     xn_m1_tagid = _xn_m1_tagid;
     206             :   }
     207             :   else
     208             :   {
     209         702 :     relaxation_factor = _secondary_relaxation_factor;
     210         702 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
     211         702 :     xn_m1_tagid = _secondary_xn_m1_tagid;
     212             :   }
     213             : 
     214        2037 :   NumericVector<Number> & solution = _solver_sys.solution();
     215        2037 :   NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
     216        2037 :   NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
     217             : 
     218      186804 :   for (const auto & dof : transformed_dofs)
     219             :   {
     220             :     // Avoid 0 denominator issue
     221      184767 :     Real new_value = solution(dof);
     222      184767 :     if (!MooseUtils::absoluteFuzzyEqual(solution(dof) + xn_m1(dof) - 2 * fxn_m1(dof), 0))
     223      151173 :       new_value = xn_m1(dof) - (fxn_m1(dof) - xn_m1(dof)) * (fxn_m1(dof) - xn_m1(dof)) /
     224      151173 :                                    (solution(dof) + xn_m1(dof) - 2 * fxn_m1(dof));
     225             : 
     226             :     // Relax update
     227      184767 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * fxn_m1(dof);
     228             : 
     229      184767 :     solution.set(dof, new_value);
     230             :   }
     231        2037 :   solution.close();
     232        2037 :   _solver_sys.update();
     233        2037 : }
     234             : 
     235             : void
     236       12051 : SteffensenSolve::printFixedPointConvergenceHistory(
     237             :     const Real initial_norm,
     238             :     const std::vector<Real> & timestep_begin_norms,
     239             :     const std::vector<Real> & timestep_end_norms) const
     240             : {
     241       12051 :   _console << "\n 0 Steffensen initialization |R| = "
     242       12051 :            << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
     243             : 
     244       12051 :   Real max_norm_old = initial_norm;
     245       55074 :   for (unsigned int i = 0; i <= _fixed_point_it; ++i)
     246             :   {
     247       43023 :     Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
     248       43023 :     std::stringstream steffensen_prefix;
     249       43023 :     if (i == 0)
     250       12051 :       steffensen_prefix << " Steffensen initialization |R| = ";
     251       30972 :     else if (i % 2 == 0)
     252       12804 :       steffensen_prefix << " Steffensen half-step      |R| = ";
     253             :     else
     254       18168 :       steffensen_prefix << " Steffensen step           |R| = ";
     255             : 
     256       43023 :     _console << std::setw(2) << i + 1 << steffensen_prefix.str()
     257       43023 :              << Console::outputNorm(max_norm_old, max_norm) << '\n';
     258             : 
     259       43023 :     max_norm_old = max_norm;
     260       43023 :   }
     261       12051 : }