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