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