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