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 "SecantSolve.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 : SecantSolve::validParams()
      20             : {
      21           0 :   InputParameters params = FixedPointSolve::validParams();
      22             : 
      23           0 :   return params;
      24             : }
      25             : 
      26         576 : SecantSolve::SecantSolve(Executioner & ex) : FixedPointSolve(ex)
      27             : {
      28         576 :   _transformed_pps_values.resize(_transformed_pps.size());
      29         676 :   for (size_t i = 0; i < _transformed_pps.size(); i++)
      30         100 :     _transformed_pps_values[i].resize(4);
      31         576 :   _secondary_transformed_pps_values.resize(_secondary_transformed_pps.size());
      32         644 :   for (size_t i = 0; i < _secondary_transformed_pps.size(); i++)
      33          68 :     _secondary_transformed_pps_values[i].resize(4);
      34         576 : }
      35             : 
      36             : void
      37         864 : SecantSolve::allocateStorage(const bool primary)
      38             : {
      39         864 :   findTransformedSystem(primary);
      40         864 :   if (!_transformed_sys)
      41         656 :     return;
      42             : 
      43             :   // TODO: We would only need to store the solution for the degrees of freedom that
      44             :   // will be transformed, not the entire solution.
      45             :   // Store solution vectors for the two previous points and their evaluation
      46         208 :   if (primary)
      47             :   {
      48         150 :     _xn_m1_tagid =
      49         150 :         _problem.addVectorTag(Moose::PREVIOUS_FP_SOLUTION_TAG, Moose::VECTOR_TAG_SOLUTION);
      50         150 :     _fxn_m1_tagid = _problem.addVectorTag("fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      51         150 :     _xn_m2_tagid = _problem.addVectorTag("xn_m2", Moose::VECTOR_TAG_SOLUTION);
      52         150 :     _fxn_m2_tagid = _problem.addVectorTag("fxn_m2", Moose::VECTOR_TAG_SOLUTION);
      53             : 
      54         150 :     _transformed_sys->needSolutionState(1, Moose::SolutionIterationType::FixedPoint, PARALLEL);
      55         150 :     _transformed_sys->addVector(_fxn_m1_tagid, false, PARALLEL);
      56         150 :     _transformed_sys->addVector(_xn_m2_tagid, false, PARALLEL);
      57         150 :     _transformed_sys->addVector(_fxn_m2_tagid, false, PARALLEL);
      58             :   }
      59             :   else
      60             :   {
      61          58 :     _secondary_xn_m1_tagid = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION);
      62          58 :     _secondary_fxn_m1_tagid = _problem.addVectorTag("secondary_fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      63          58 :     _secondary_xn_m2_tagid = _problem.addVectorTag("secondary_xn_m2", Moose::VECTOR_TAG_SOLUTION);
      64          58 :     _secondary_fxn_m2_tagid = _problem.addVectorTag("secondary_fxn_m2", Moose::VECTOR_TAG_SOLUTION);
      65             : 
      66          58 :     _transformed_sys->addVector(_secondary_xn_m1_tagid, false, PARALLEL);
      67          58 :     _transformed_sys->addVector(_secondary_fxn_m1_tagid, false, PARALLEL);
      68          58 :     _transformed_sys->addVector(_secondary_xn_m2_tagid, false, PARALLEL);
      69          58 :     _transformed_sys->addVector(_secondary_fxn_m2_tagid, false, PARALLEL);
      70             :   }
      71             : }
      72             : 
      73             : void
      74        9545 : SecantSolve::saveVariableValues(const bool primary)
      75             : {
      76             :   TagID fxn_m1_tagid;
      77             :   TagID xn_m1_tagid;
      78             :   TagID fxn_m2_tagid;
      79             :   TagID xn_m2_tagid;
      80        9545 :   if (primary)
      81             :   {
      82        7401 :     fxn_m1_tagid = _fxn_m1_tagid;
      83        7401 :     xn_m1_tagid = _xn_m1_tagid;
      84        7401 :     fxn_m2_tagid = _fxn_m2_tagid;
      85        7401 :     xn_m2_tagid = _xn_m2_tagid;
      86             :   }
      87             :   else
      88             :   {
      89        2144 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
      90        2144 :     xn_m1_tagid = _secondary_xn_m1_tagid;
      91        2144 :     fxn_m2_tagid = _secondary_fxn_m2_tagid;
      92        2144 :     xn_m2_tagid = _secondary_xn_m2_tagid;
      93             :   }
      94             : 
      95             :   // Check to make sure allocateStorage has been called
      96             :   mooseAssert(fxn_m1_tagid != Moose::INVALID_TAG_ID,
      97             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
      98             :   mooseAssert(xn_m1_tagid != Moose::INVALID_TAG_ID,
      99             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
     100             :   mooseAssert(fxn_m2_tagid != Moose::INVALID_TAG_ID,
     101             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
     102             :   mooseAssert(xn_m2_tagid != Moose::INVALID_TAG_ID,
     103             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
     104             : 
     105             :   // Save previous variable values
     106        9545 :   NumericVector<Number> & solution = _transformed_sys->solution();
     107        9545 :   NumericVector<Number> & fxn_m1 = _transformed_sys->getVector(fxn_m1_tagid);
     108        9545 :   NumericVector<Number> & xn_m1 = _transformed_sys->getVector(xn_m1_tagid);
     109        9545 :   NumericVector<Number> & fxn_m2 = _transformed_sys->getVector(fxn_m2_tagid);
     110        9545 :   NumericVector<Number> & xn_m2 = _transformed_sys->getVector(xn_m2_tagid);
     111             : 
     112             :   // Advance one step
     113        9545 :   xn_m2 = xn_m1;
     114             : 
     115             :   // Before a solve, solution is a sequence term, after a solve, solution is the evaluated term
     116             :   // Primary is copied back by _transformed_sys->copyPreviousFixedPointSolutions()
     117        9545 :   if (!primary)
     118        2144 :     xn_m1 = solution;
     119             : 
     120             :   // Since we did not update on the 0th iteration, the solution is also the previous evaluated term
     121        9545 :   const unsigned int it = primary ? _fixed_point_it : _main_fixed_point_it;
     122        9545 :   if (it == 1)
     123        1021 :     fxn_m2 = solution;
     124             :   // Otherwise we just advance
     125             :   else
     126        8524 :     fxn_m2 = fxn_m1;
     127        9545 : }
     128             : 
     129             : void
     130       42273 : SecantSolve::savePostprocessorValues(const bool primary)
     131             : {
     132             :   const std::vector<PostprocessorName> * transformed_pps;
     133             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     134       42273 :   if (primary)
     135             :   {
     136       30162 :     transformed_pps = &_transformed_pps;
     137       30162 :     transformed_pps_values = &_transformed_pps_values;
     138             :   }
     139             :   else
     140             :   {
     141       12111 :     transformed_pps = &_secondary_transformed_pps;
     142       12111 :     transformed_pps_values = &_secondary_transformed_pps_values;
     143             :   }
     144       42273 :   const unsigned int it = primary ? _fixed_point_it : _main_fixed_point_it;
     145             : 
     146             :   // Save previous postprocessor values
     147       49493 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     148             :   {
     149             :     // Advance one step
     150        7220 :     (*transformed_pps_values)[i][3] = (*transformed_pps_values)[i][1];
     151             : 
     152             :     // Save current value
     153             :     // Primary: this is done before the timestep's solves and before timestep_begin transfers,
     154             :     // so the value is the result of the previous Secant update (xn_m1)
     155             :     // Secondary: this is done after the secondary solve, but before timestep_end postprocessors
     156             :     // are computed, or timestep_end transfers are received.
     157             :     // This value is the same as before the solve (xn_m1)
     158        7220 :     (*transformed_pps_values)[i][1] = getPostprocessorValueByName((*transformed_pps)[i]);
     159             : 
     160             :     // Since we did not update on the 1st iteration, the pp is also the previous evaluated term
     161        7220 :     if (it == 2)
     162         989 :       (*transformed_pps_values)[i][2] = (*transformed_pps_values)[i][1];
     163             :     // Otherwise we just advance
     164             :     else
     165        6231 :       (*transformed_pps_values)[i][2] = (*transformed_pps_values)[i][0];
     166             :   }
     167       42273 : }
     168             : 
     169             : bool
     170       15671 : SecantSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary)
     171             : {
     172             :   // Need at least two evaluations to compute the Secant slope
     173       15671 :   if (primary)
     174        9070 :     return _fixed_point_it > 0;
     175             :   else
     176        6601 :     return _main_fixed_point_it > 0;
     177             : }
     178             : 
     179             : void
     180        6439 : SecantSolve::transformPostprocessors(const bool primary)
     181             : {
     182        6439 :   if ((primary ? _fixed_point_it : _main_fixed_point_it) < 2)
     183        1009 :     return;
     184             : 
     185             :   Real relaxation_factor;
     186             :   const std::vector<PostprocessorName> * transformed_pps;
     187             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     188        5430 :   if (primary)
     189             :   {
     190        3214 :     relaxation_factor = _relax_factor;
     191        3214 :     transformed_pps = &_transformed_pps;
     192        3214 :     transformed_pps_values = &_transformed_pps_values;
     193             :   }
     194             :   else
     195             :   {
     196        2216 :     relaxation_factor = _secondary_relaxation_factor;
     197        2216 :     transformed_pps = &_secondary_transformed_pps;
     198        2216 :     transformed_pps_values = &_secondary_transformed_pps_values;
     199             :   }
     200             : 
     201             :   // Relax postprocessors for the main application
     202       10860 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     203             :   {
     204             :     // Get new postprocessor value
     205        5430 :     const Real fxn_m1 = getPostprocessorValueByName((*transformed_pps)[i]);
     206        5430 :     const Real xn_m1 = (*transformed_pps_values)[i][1];
     207        5430 :     const Real fxn_m2 = (*transformed_pps_values)[i][2];
     208        5430 :     const Real xn_m2 = (*transformed_pps_values)[i][3];
     209             : 
     210             :     // Save fxn_m1, received or computed before the solve
     211        5430 :     (*transformed_pps_values)[i][0] = fxn_m1;
     212             : 
     213             :     // Compute and set relaxed value
     214        5430 :     Real new_value = fxn_m1;
     215        5430 :     if (!MooseUtils::absoluteFuzzyEqual(fxn_m1 - xn_m1 - fxn_m2 + xn_m2, 0))
     216        5397 :       new_value = xn_m1 - (fxn_m1 - xn_m1) * (xn_m1 - xn_m2) / (fxn_m1 - xn_m1 - fxn_m2 + xn_m2);
     217             : 
     218             :     // Relax update if desired
     219        5430 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1;
     220             : 
     221        5430 :     _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
     222             :   }
     223             : }
     224             : 
     225             : void
     226        5758 : SecantSolve::transformVariables(const std::set<dof_id_type> & target_dofs, const bool primary)
     227             : {
     228             :   Real relaxation_factor;
     229             :   TagID fxn_m1_tagid;
     230             :   TagID xn_m1_tagid;
     231             :   TagID fxn_m2_tagid;
     232             :   TagID xn_m2_tagid;
     233        5758 :   if (primary)
     234             :   {
     235        3614 :     relaxation_factor = _relax_factor;
     236        3614 :     fxn_m1_tagid = _fxn_m1_tagid;
     237        3614 :     xn_m1_tagid = _xn_m1_tagid;
     238        3614 :     fxn_m2_tagid = _fxn_m2_tagid;
     239        3614 :     xn_m2_tagid = _xn_m2_tagid;
     240             :   }
     241             :   else
     242             :   {
     243        2144 :     relaxation_factor = _secondary_relaxation_factor;
     244        2144 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
     245        2144 :     xn_m1_tagid = _secondary_xn_m1_tagid;
     246        2144 :     fxn_m2_tagid = _secondary_fxn_m2_tagid;
     247        2144 :     xn_m2_tagid = _secondary_xn_m2_tagid;
     248             :   }
     249             : 
     250        5758 :   NumericVector<Number> & solution = _transformed_sys->solution();
     251        5758 :   NumericVector<Number> & xn_m1 = _transformed_sys->getVector(xn_m1_tagid);
     252        5758 :   NumericVector<Number> & fxn_m2 = _transformed_sys->getVector(fxn_m2_tagid);
     253        5758 :   NumericVector<Number> & xn_m2 = _transformed_sys->getVector(xn_m2_tagid);
     254             : 
     255             :   // Save the most recent evaluation of the coupled problem
     256        5758 :   NumericVector<Number> & fxn_m1 = _transformed_sys->getVector(fxn_m1_tagid);
     257        5758 :   fxn_m1 = solution;
     258             : 
     259      512627 :   for (const auto & dof : target_dofs)
     260             :   {
     261             :     // Avoid 0 denominator issue
     262      506869 :     Real new_value = fxn_m1(dof);
     263      506869 :     if (!MooseUtils::absoluteFuzzyEqual(solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof), 0))
     264      414953 :       new_value = xn_m1(dof) - (solution(dof) - xn_m1(dof)) * (xn_m1(dof) - xn_m2(dof)) /
     265      414953 :                                    (solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof));
     266             : 
     267             :     // Relax update
     268      506869 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1(dof);
     269             : 
     270      506869 :     solution.set(dof, new_value);
     271             :   }
     272        5758 :   solution.close();
     273        5758 :   _transformed_sys->update();
     274        5758 : }
     275             : 
     276             : void
     277       14651 : SecantSolve::printFixedPointConvergenceHistory(Real initial_norm,
     278             :                                                const std::vector<Real> & timestep_begin_norms,
     279             :                                                const std::vector<Real> & timestep_end_norms) const
     280             : {
     281       14651 :   _console << "\n 0 Secant initialization |R| = "
     282       14651 :            << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
     283             : 
     284       14651 :   Real max_norm_old = initial_norm;
     285       83952 :   for (unsigned int i = 0; i <= _fixed_point_it; ++i)
     286             :   {
     287       69301 :     Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
     288       69301 :     std::stringstream secant_prefix;
     289       69301 :     if (i < 1)
     290       14651 :       secant_prefix << " Secant initialization |R| = ";
     291             :     else
     292       54650 :       secant_prefix << " Secant step           |R| = ";
     293             : 
     294       69301 :     _console << std::setw(2) << i + 1 << secant_prefix.str()
     295       69301 :              << Console::outputNorm(max_norm_old, max_norm) << '\n';
     296       69301 :     max_norm_old = max_norm;
     297       69301 :   }
     298       14651 : }