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 "PicardSolve.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 : PicardSolve::validParams()
      20             : {
      21           0 :   InputParameters params = FixedPointSolve::validParams();
      22           0 :   return params;
      23             : }
      24             : 
      25       61636 : PicardSolve::PicardSolve(Executioner & ex) : FixedPointSolve(ex) {}
      26             : 
      27             : void
      28       69688 : PicardSolve::allocateStorage(const bool primary)
      29             : {
      30       69688 :   if (!performingRelaxation(primary))
      31       69467 :     return;
      32             : 
      33             :   const std::vector<PostprocessorName> * transformed_pps;
      34             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
      35         221 :   if (primary)
      36             :   {
      37         186 :     _old_tag_id =
      38         186 :         _problem.addVectorTag(Moose::PREVIOUS_FP_SOLUTION_TAG, Moose::VECTOR_TAG_SOLUTION);
      39         186 :     _solver_sys.needSolutionState(1, Moose::SolutionIterationType::FixedPoint, PARALLEL);
      40             : 
      41         186 :     transformed_pps = &_transformed_pps;
      42         186 :     transformed_pps_values = &_transformed_pps_values;
      43             :   }
      44             :   else
      45             :   {
      46          35 :     _secondary_old_tag_id = _problem.addVectorTag("secondary_xn_m1", Moose::VECTOR_TAG_SOLUTION);
      47          35 :     _solver_sys.addVector(_secondary_old_tag_id, false, PARALLEL);
      48             : 
      49          35 :     transformed_pps = &_secondary_transformed_pps;
      50          35 :     transformed_pps_values = &_secondary_transformed_pps_values;
      51             :   }
      52             : 
      53             :   // Allocate storage for the previous postprocessor values
      54         221 :   (*transformed_pps_values).resize((*transformed_pps).size());
      55         351 :   for (const auto i : index_range(*transformed_pps))
      56         130 :     (*transformed_pps_values)[i].resize(1);
      57             : }
      58             : 
      59             : void
      60      275549 : PicardSolve::saveVariableValues(const bool primary)
      61             : {
      62             :   // Primary is copied back by _solver_sys.copyPreviousFixedPointSolutions()
      63      275549 :   if (!performingRelaxation(primary) || primary)
      64      275220 :     return;
      65             : 
      66             :   // Check to make sure allocateStorage has been called
      67             :   mooseAssert(_secondary_old_tag_id != Moose::INVALID_TAG_ID,
      68             :               "allocateStorage has not been called with primary = " + Moose::stringify(primary));
      69             : 
      70             :   // Save variable previous values
      71         329 :   NumericVector<Number> & solution = _solver_sys.solution();
      72         329 :   NumericVector<Number> & transformed_old = _solver_sys.getVector(_secondary_old_tag_id);
      73         329 :   transformed_old = solution;
      74             : }
      75             : 
      76             : void
      77      299827 : PicardSolve::savePostprocessorValues(const bool primary)
      78             : {
      79      299827 :   if (!performingRelaxation(primary))
      80      286598 :     return;
      81             : 
      82             :   const std::vector<PostprocessorName> * transformed_pps;
      83             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
      84       13229 :   if (primary)
      85             :   {
      86       12899 :     transformed_pps = &_transformed_pps;
      87       12899 :     transformed_pps_values = &_transformed_pps_values;
      88             :   }
      89             :   else
      90             :   {
      91         330 :     transformed_pps = &_secondary_transformed_pps;
      92         330 :     transformed_pps_values = &_secondary_transformed_pps_values;
      93             :   }
      94             : 
      95             :   // Save postprocessor previous values
      96       22596 :   for (const auto i : index_range(*transformed_pps))
      97        9367 :     (*transformed_pps_values)[i][0] = getPostprocessorValueByName((*transformed_pps)[i]);
      98             : }
      99             : 
     100             : bool
     101       14498 : PicardSolve::useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary)
     102             : {
     103             :   // unrelaxed Picard is the default update for fixed point iterations
     104             :   // old values are required for relaxation
     105       14498 :   const auto fixed_point_it = primary ? _fixed_point_it : _main_fixed_point_it;
     106       14498 :   return performingRelaxation(primary) && fixed_point_it > 0;
     107             : }
     108             : 
     109             : void
     110        8386 : PicardSolve::transformPostprocessors(const bool primary)
     111             : {
     112             :   Real relaxation_factor;
     113             :   const std::vector<PostprocessorName> * transformed_pps;
     114             :   std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
     115        8386 :   if (primary)
     116             :   {
     117        8375 :     relaxation_factor = _relax_factor;
     118        8375 :     transformed_pps = &_transformed_pps;
     119        8375 :     transformed_pps_values = &_transformed_pps_values;
     120             :   }
     121             :   else
     122             :   {
     123          11 :     relaxation_factor = _secondary_relaxation_factor;
     124          11 :     transformed_pps = &_secondary_transformed_pps;
     125          11 :     transformed_pps_values = &_secondary_transformed_pps_values;
     126             :   }
     127             : 
     128             :   // Relax the postprocessors
     129       16772 :   for (size_t i = 0; i < (*transformed_pps).size(); i++)
     130             :   {
     131             :     // Get new postprocessor value
     132        8386 :     const Real current_value = getPostprocessorValueByName((*transformed_pps)[i]);
     133        8386 :     const Real old_value = (*transformed_pps_values)[i][0];
     134             : 
     135             :     // Compute and set relaxed value
     136        8386 :     Real new_value = current_value;
     137        8386 :     new_value = relaxation_factor * current_value + (1 - relaxation_factor) * old_value;
     138        8386 :     _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
     139             :   }
     140        8386 : }
     141             : 
     142             : void
     143        1101 : PicardSolve::transformVariables(const std::set<dof_id_type> & transformed_dofs, const bool primary)
     144             : {
     145             :   Real relaxation_factor;
     146             :   TagID old_tag_id;
     147        1101 :   if (primary)
     148             :   {
     149         783 :     relaxation_factor = _relax_factor;
     150         783 :     old_tag_id = _old_tag_id;
     151             :   }
     152             :   else
     153             :   {
     154         318 :     relaxation_factor = _secondary_relaxation_factor;
     155         318 :     old_tag_id = _secondary_old_tag_id;
     156             :   }
     157             : 
     158        1101 :   NumericVector<Number> & solution = _solver_sys.solution();
     159        1101 :   NumericVector<Number> & transformed_old = _solver_sys.getVector(old_tag_id);
     160             : 
     161      145333 :   for (const auto & dof : transformed_dofs)
     162      288464 :     solution.set(dof,
     163      144232 :                  (transformed_old(dof) * (1.0 - relaxation_factor)) +
     164      144232 :                      (solution(dof) * relaxation_factor));
     165             : 
     166        1101 :   solution.close();
     167        1101 :   _solver_sys.update();
     168        1101 : }
     169             : 
     170             : void
     171       32294 : PicardSolve::printFixedPointConvergenceHistory(Real initial_norm,
     172             :                                                const std::vector<Real> & timestep_begin_norms,
     173             :                                                const std::vector<Real> & timestep_end_norms) const
     174             : {
     175       32294 :   _console << "\n 0 Picard |R| = "
     176       32294 :            << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
     177             : 
     178       32294 :   Real max_norm_old = initial_norm;
     179      161801 :   for (unsigned int i = 0; i <= _fixed_point_it; ++i)
     180             :   {
     181      129507 :     Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
     182      129507 :     _console << std::setw(2) << i + 1
     183      129507 :              << " Picard |R| = " << Console::outputNorm(max_norm_old, max_norm) << '\n';
     184      129507 :     max_norm_old = max_norm;
     185             :   }
     186             : 
     187       32294 :   _console << std::endl;
     188       32294 : }