LCOV - code coverage report
Current view: top level - src/executioners - SecantSolve.C (source / functions) Hit Total Coverage
Test: idaholab/moose framework: 2bf808 Lines: 129 132 97.7 %
Date: 2025-07-17 01:28:37 Functions: 8 9 88.9 %
Legend: Lines: hit not hit

          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 :   allocateStorage(true);
      29             : 
      30         576 :   _transformed_pps_values.resize(_transformed_pps.size());
      31         676 :   for (size_t i = 0; i < _transformed_pps.size(); i++)
      32         100 :     _transformed_pps_values[i].resize(4);
      33         576 :   _secondary_transformed_pps_values.resize(_secondary_transformed_pps.size());
      34         644 :   for (size_t i = 0; i < _secondary_transformed_pps.size(); i++)
      35          68 :     _secondary_transformed_pps_values[i].resize(4);
      36         576 : }
      37             : 
      38             : void
      39         864 : 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         864 :   if (primary)
      46             :   {
      47         576 :     xn_m1_tagid = _problem.addVectorTag("xn_m1", Moose::VECTOR_TAG_SOLUTION);
      48         576 :     fxn_m1_tagid = _problem.addVectorTag("fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      49         576 :     xn_m2_tagid = _problem.addVectorTag("xn_m2", Moose::VECTOR_TAG_SOLUTION);
      50         576 :     fxn_m2_tagid = _problem.addVectorTag("fxn_m2", Moose::VECTOR_TAG_SOLUTION);
      51         576 :     _xn_m1_tagid = xn_m1_tagid;
      52         576 :     _fxn_m1_tagid = fxn_m1_tagid;
      53         576 :     _xn_m2_tagid = xn_m2_tagid;
      54         576 :     _fxn_m2_tagid = fxn_m2_tagid;
      55             :   }
      56             :   else
      57             :   {
      58         288 :     xn_m1_tagid = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION);
      59         288 :     fxn_m1_tagid = _problem.addVectorTag("secondary_fxn_m1", Moose::VECTOR_TAG_SOLUTION);
      60         288 :     xn_m2_tagid = _problem.addVectorTag("secondary_xn_m2", Moose::VECTOR_TAG_SOLUTION);
      61         288 :     fxn_m2_tagid = _problem.addVectorTag("secondary_fxn_m2", Moose::VECTOR_TAG_SOLUTION);
      62         288 :     _secondary_xn_m1_tagid = xn_m1_tagid;
      63         288 :     _secondary_fxn_m1_tagid = fxn_m1_tagid;
      64         288 :     _secondary_xn_m2_tagid = xn_m2_tagid;
      65         288 :     _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         864 :   _solver_sys.addVector(xn_m1_tagid, false, PARALLEL);
      72         864 :   _solver_sys.addVector(fxn_m1_tagid, false, PARALLEL);
      73         864 :   _solver_sys.addVector(xn_m2_tagid, false, PARALLEL);
      74         864 :   _solver_sys.addVector(fxn_m2_tagid, false, PARALLEL);
      75         864 : }
      76             : 
      77             : void
      78       42852 : 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       42852 :   if (primary)
      85             :   {
      86       30592 :     fxn_m1_tagid = _fxn_m1_tagid;
      87       30592 :     xn_m1_tagid = _xn_m1_tagid;
      88       30592 :     fxn_m2_tagid = _fxn_m2_tagid;
      89       30592 :     xn_m2_tagid = _xn_m2_tagid;
      90             :   }
      91             :   else
      92             :   {
      93       12260 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
      94       12260 :     xn_m1_tagid = _secondary_xn_m1_tagid;
      95       12260 :     fxn_m2_tagid = _secondary_fxn_m2_tagid;
      96       12260 :     xn_m2_tagid = _secondary_xn_m2_tagid;
      97             :   }
      98             : 
      99             :   // Save previous variable values
     100       42852 :   NumericVector<Number> & solution = _solver_sys.solution();
     101       42852 :   NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
     102       42852 :   NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
     103       42852 :   NumericVector<Number> & fxn_m2 = _solver_sys.getVector(fxn_m2_tagid);
     104       42852 :   NumericVector<Number> & xn_m2 = _solver_sys.getVector(xn_m2_tagid);
     105             : 
     106             :   // Advance one step
     107       42852 :   xn_m2 = xn_m1;
     108       42852 :   fxn_m2 = fxn_m1;
     109             : 
     110             :   // Before a solve, solution is a sequence term, after a solve, solution is the evaluated term
     111       42852 :   xn_m1 = solution;
     112       42852 : }
     113             : 
     114             : void
     115       42852 : SecantSolve::savePostprocessorValues(const bool primary)
     116             : {
     117             :   const std::vector<PostprocessorName> * transformed_pps;
     118             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     119       42852 :   if (primary)
     120             :   {
     121       30592 :     transformed_pps = &_transformed_pps;
     122       30592 :     transformed_pps_values = &_transformed_pps_values;
     123             :   }
     124             :   else
     125             :   {
     126       12260 :     transformed_pps = &_secondary_transformed_pps;
     127       12260 :     transformed_pps_values = &_secondary_transformed_pps_values;
     128             :   }
     129             : 
     130             :   // Save previous postprocessor values
     131       49652 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     132             :   {
     133             :     // Advance one step
     134        6800 :     (*transformed_pps_values)[i][3] = (*transformed_pps_values)[i][1];
     135        6800 :     (*transformed_pps_values)[i][2] = (*transformed_pps_values)[i][0];
     136             : 
     137             :     // Save current value
     138             :     // Primary: this is done before the timestep's solves and before timestep_begin transfers,
     139             :     // so the value is the result of the previous Secant update (xn_m1)
     140             :     // Secondary: this is done after the secondary solve, but before timestep_end postprocessors
     141             :     // are computed, or timestep_end transfers are received.
     142             :     // This value is the same as before the solve (xn_m1)
     143        6800 :     (*transformed_pps_values)[i][1] = getPostprocessorValueByName((*transformed_pps)[i]);
     144             :   }
     145       42852 : }
     146             : 
     147             : bool
     148       15831 : SecantSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary)
     149             : {
     150             :   // Need at least two evaluations to compute the Secant slope
     151       15831 :   if (primary)
     152        9730 :     return _fixed_point_it > 1;
     153             :   else
     154        6101 :     return _main_fixed_point_it > 1;
     155             : }
     156             : 
     157             : void
     158        5030 : SecantSolve::transformPostprocessors(const bool primary)
     159             : {
     160             :   Real relaxation_factor;
     161             :   const std::vector<PostprocessorName> * transformed_pps;
     162             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     163        5030 :   if (primary)
     164             :   {
     165        2940 :     relaxation_factor = _relax_factor;
     166        2940 :     transformed_pps = &_transformed_pps;
     167        2940 :     transformed_pps_values = &_transformed_pps_values;
     168             :   }
     169             :   else
     170             :   {
     171        2090 :     relaxation_factor = _secondary_relaxation_factor;
     172        2090 :     transformed_pps = &_secondary_transformed_pps;
     173        2090 :     transformed_pps_values = &_secondary_transformed_pps_values;
     174             :   }
     175             : 
     176             :   // Relax postprocessors for the main application
     177       10060 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     178             :   {
     179             :     // Get new postprocessor value
     180        5030 :     const Real fxn_m1 = getPostprocessorValueByName((*transformed_pps)[i]);
     181        5030 :     const Real xn_m1 = (*transformed_pps_values)[i][1];
     182        5030 :     const Real fxn_m2 = (*transformed_pps_values)[i][2];
     183        5030 :     const Real xn_m2 = (*transformed_pps_values)[i][3];
     184             : 
     185             :     // Save fxn_m1, received or computed before the solve
     186        5030 :     (*transformed_pps_values)[i][0] = fxn_m1;
     187             : 
     188             :     // Compute and set relaxed value
     189        5030 :     Real new_value = fxn_m1;
     190        5030 :     if (!MooseUtils::absoluteFuzzyEqual(fxn_m1 - xn_m1 - fxn_m2 + xn_m2, 0))
     191        5030 :       new_value = xn_m1 - (fxn_m1 - xn_m1) * (xn_m1 - xn_m2) / (fxn_m1 - xn_m1 - fxn_m2 + xn_m2);
     192             : 
     193             :     // Relax update if desired
     194        5030 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1;
     195             : 
     196        5030 :     _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
     197             :   }
     198        5030 : }
     199             : 
     200             : void
     201        5239 : SecantSolve::transformVariables(const std::set<dof_id_type> & target_dofs, const bool primary)
     202             : {
     203             :   Real relaxation_factor;
     204             :   TagID fxn_m1_tagid;
     205             :   TagID xn_m1_tagid;
     206             :   TagID fxn_m2_tagid;
     207             :   TagID xn_m2_tagid;
     208        5239 :   if (primary)
     209             :   {
     210        3688 :     relaxation_factor = _relax_factor;
     211        3688 :     fxn_m1_tagid = _fxn_m1_tagid;
     212        3688 :     xn_m1_tagid = _xn_m1_tagid;
     213        3688 :     fxn_m2_tagid = _fxn_m2_tagid;
     214        3688 :     xn_m2_tagid = _xn_m2_tagid;
     215             :   }
     216             :   else
     217             :   {
     218        1551 :     relaxation_factor = _secondary_relaxation_factor;
     219        1551 :     fxn_m1_tagid = _secondary_fxn_m1_tagid;
     220        1551 :     xn_m1_tagid = _secondary_xn_m1_tagid;
     221        1551 :     fxn_m2_tagid = _secondary_fxn_m2_tagid;
     222        1551 :     xn_m2_tagid = _secondary_xn_m2_tagid;
     223             :   }
     224             : 
     225        5239 :   NumericVector<Number> & solution = _solver_sys.solution();
     226        5239 :   NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
     227        5239 :   NumericVector<Number> & fxn_m2 = _solver_sys.getVector(fxn_m2_tagid);
     228        5239 :   NumericVector<Number> & xn_m2 = _solver_sys.getVector(xn_m2_tagid);
     229             : 
     230             :   // Save the most recent evaluation of the coupled problem
     231        5239 :   NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
     232        5239 :   fxn_m1 = solution;
     233             : 
     234      464555 :   for (const auto & dof : target_dofs)
     235             :   {
     236             :     // Avoid 0 denominator issue
     237      459316 :     Real new_value = fxn_m1(dof);
     238      459316 :     if (!MooseUtils::absoluteFuzzyEqual(solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof), 0))
     239      376662 :       new_value = xn_m1(dof) - (solution(dof) - xn_m1(dof)) * (xn_m1(dof) - xn_m2(dof)) /
     240      376662 :                                    (solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof));
     241             : 
     242             :     // Relax update
     243      459316 :     new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1(dof);
     244             : 
     245      459316 :     solution.set(dof, new_value);
     246             :   }
     247        5239 :   solution.close();
     248        5239 :   _solver_sys.update();
     249        5239 : }
     250             : 
     251             : void
     252       14911 : SecantSolve::printFixedPointConvergenceHistory(Real initial_norm,
     253             :                                                const std::vector<Real> & timestep_begin_norms,
     254             :                                                const std::vector<Real> & timestep_end_norms) const
     255             : {
     256       14911 :   _console << "\n 0 Secant initialization |R| = "
     257       14911 :            << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
     258             : 
     259       14911 :   Real max_norm_old = initial_norm;
     260       87037 :   for (unsigned int i = 0; i <= _fixed_point_it; ++i)
     261             :   {
     262       72126 :     Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
     263       72126 :     std::stringstream secant_prefix;
     264       72126 :     if (i < 2)
     265       27281 :       secant_prefix << " Secant initialization |R| = ";
     266             :     else
     267       44845 :       secant_prefix << " Secant step           |R| = ";
     268             : 
     269       72126 :     _console << std::setw(2) << i + 1 << secant_prefix.str()
     270       72126 :              << Console::outputNorm(max_norm_old, max_norm) << '\n';
     271       72126 :     max_norm_old = max_norm;
     272       72126 :   }
     273       14911 : }

Generated by: LCOV version 1.14