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