https://mooseframework.inl.gov
SecantSolve.C
Go to the documentation of this file.
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"
16 #include "Console.h"
17 
20 {
22 
23  return params;
24 }
25 
27 {
29  for (size_t i = 0; i < _transformed_pps.size(); i++)
30  _transformed_pps_values[i].resize(4);
32  for (size_t i = 0; i < _secondary_transformed_pps.size(); i++)
34 }
35 
36 void
37 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  if (primary)
43  {
44  _xn_m1_tagid =
49 
54  }
55  else
56  {
61 
66  }
67 }
68 
69 void
71 {
72  TagID fxn_m1_tagid;
73  TagID xn_m1_tagid;
74  TagID fxn_m2_tagid;
75  TagID xn_m2_tagid;
76  if (primary)
77  {
78  fxn_m1_tagid = _fxn_m1_tagid;
79  xn_m1_tagid = _xn_m1_tagid;
80  fxn_m2_tagid = _fxn_m2_tagid;
81  xn_m2_tagid = _xn_m2_tagid;
82  }
83  else
84  {
85  fxn_m1_tagid = _secondary_fxn_m1_tagid;
86  xn_m1_tagid = _secondary_xn_m1_tagid;
87  fxn_m2_tagid = _secondary_fxn_m2_tagid;
88  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
103  NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
104  NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
105  NumericVector<Number> & fxn_m2 = _solver_sys.getVector(fxn_m2_tagid);
106  NumericVector<Number> & xn_m2 = _solver_sys.getVector(xn_m2_tagid);
107 
108  // Advance one step
109  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  if (!primary)
114  xn_m1 = solution;
115 
116  // Since we did not update on the 0th iteration, the solution is also the previous evaluated term
117  const unsigned int it = primary ? _fixed_point_it : _main_fixed_point_it;
118  if (it == 1)
119  fxn_m2 = solution;
120  // Otherwise we just advance
121  else
122  fxn_m2 = fxn_m1;
123 }
124 
125 void
127 {
128  const std::vector<PostprocessorName> * transformed_pps;
129  std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
130  if (primary)
131  {
132  transformed_pps = &_transformed_pps;
133  transformed_pps_values = &_transformed_pps_values;
134  }
135  else
136  {
137  transformed_pps = &_secondary_transformed_pps;
138  transformed_pps_values = &_secondary_transformed_pps_values;
139  }
140  const unsigned int it = primary ? _fixed_point_it : _main_fixed_point_it;
141 
142  // Save previous postprocessor values
143  for (size_t i = 0; i < (*transformed_pps).size(); i++)
144  {
145  // Advance one step
146  (*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  (*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  if (it == 2)
158  (*transformed_pps_values)[i][2] = (*transformed_pps_values)[i][1];
159  // Otherwise we just advance
160  else
161  (*transformed_pps_values)[i][2] = (*transformed_pps_values)[i][0];
162  }
163 }
164 
165 bool
167 {
168  // Need at least two evaluations to compute the Secant slope
169  if (primary)
170  return _fixed_point_it > 0;
171  else
172  return _main_fixed_point_it > 0;
173 }
174 
175 void
177 {
178  if ((primary ? _fixed_point_it : _main_fixed_point_it) < 2)
179  return;
180 
181  Real relaxation_factor;
182  const std::vector<PostprocessorName> * transformed_pps;
183  std::vector<std::vector<PostprocessorValue>> * transformed_pps_values;
184  if (primary)
185  {
186  relaxation_factor = _relax_factor;
187  transformed_pps = &_transformed_pps;
188  transformed_pps_values = &_transformed_pps_values;
189  }
190  else
191  {
192  relaxation_factor = _secondary_relaxation_factor;
193  transformed_pps = &_secondary_transformed_pps;
194  transformed_pps_values = &_secondary_transformed_pps_values;
195  }
196 
197  // Relax postprocessors for the main application
198  for (size_t i = 0; i < (*transformed_pps).size(); i++)
199  {
200  // Get new postprocessor value
201  const Real fxn_m1 = getPostprocessorValueByName((*transformed_pps)[i]);
202  const Real xn_m1 = (*transformed_pps_values)[i][1];
203  const Real fxn_m2 = (*transformed_pps_values)[i][2];
204  const Real xn_m2 = (*transformed_pps_values)[i][3];
205 
206  // Save fxn_m1, received or computed before the solve
207  (*transformed_pps_values)[i][0] = fxn_m1;
208 
209  // Compute and set relaxed value
210  Real new_value = fxn_m1;
211  if (!MooseUtils::absoluteFuzzyEqual(fxn_m1 - xn_m1 - fxn_m2 + xn_m2, 0))
212  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  new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1;
216 
217  _problem.setPostprocessorValueByName((*transformed_pps)[i], new_value);
218  }
219 }
220 
221 void
222 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  if (primary)
230  {
231  relaxation_factor = _relax_factor;
232  fxn_m1_tagid = _fxn_m1_tagid;
233  xn_m1_tagid = _xn_m1_tagid;
234  fxn_m2_tagid = _fxn_m2_tagid;
235  xn_m2_tagid = _xn_m2_tagid;
236  }
237  else
238  {
239  relaxation_factor = _secondary_relaxation_factor;
240  fxn_m1_tagid = _secondary_fxn_m1_tagid;
241  xn_m1_tagid = _secondary_xn_m1_tagid;
242  fxn_m2_tagid = _secondary_fxn_m2_tagid;
243  xn_m2_tagid = _secondary_xn_m2_tagid;
244  }
245 
247  NumericVector<Number> & xn_m1 = _solver_sys.getVector(xn_m1_tagid);
248  NumericVector<Number> & fxn_m2 = _solver_sys.getVector(fxn_m2_tagid);
249  NumericVector<Number> & xn_m2 = _solver_sys.getVector(xn_m2_tagid);
250 
251  // Save the most recent evaluation of the coupled problem
252  NumericVector<Number> & fxn_m1 = _solver_sys.getVector(fxn_m1_tagid);
253  fxn_m1 = solution;
254 
255  for (const auto & dof : target_dofs)
256  {
257  // Avoid 0 denominator issue
258  Real new_value = fxn_m1(dof);
259  if (!MooseUtils::absoluteFuzzyEqual(solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof), 0))
260  new_value = xn_m1(dof) - (solution(dof) - xn_m1(dof)) * (xn_m1(dof) - xn_m2(dof)) /
261  (solution(dof) - xn_m1(dof) - fxn_m2(dof) + xn_m2(dof));
262 
263  // Relax update
264  new_value = relaxation_factor * new_value + (1 - relaxation_factor) * xn_m1(dof);
265 
266  solution.set(dof, new_value);
267  }
268  solution.close();
270 }
271 
272 void
274  const std::vector<Real> & timestep_begin_norms,
275  const std::vector<Real> & timestep_end_norms) const
276 {
277  _console << "\n 0 Secant initialization |R| = "
278  << Console::outputNorm(std::numeric_limits<Real>::max(), initial_norm) << '\n';
279 
280  Real max_norm_old = initial_norm;
281  for (unsigned int i = 0; i <= _fixed_point_it; ++i)
282  {
283  Real max_norm = std::max(timestep_begin_norms[i], timestep_end_norms[i]);
284  std::stringstream secant_prefix;
285  if (i < 1)
286  secant_prefix << " Secant initialization |R| = ";
287  else
288  secant_prefix << " Secant step |R| = ";
289 
290  _console << std::setw(2) << i + 1 << secant_prefix.str()
291  << Console::outputNorm(max_norm_old, max_norm) << '\n';
292  max_norm_old = max_norm;
293  }
294 }
std::vector< std::vector< PostprocessorValue > > _transformed_pps_values
Previous values of the relaxed postprocessors.
std::vector< std::vector< PostprocessorValue > > _secondary_transformed_pps_values
Previous values of the postprocessors relaxed outside of the fixed point iteration (used as a subapp)...
FEProblemBase & _problem
Reference to FEProblem.
Definition: SolveObject.h:47
TagID _secondary_xn_m1_tagid
Vector tag id for the solution variable before the latest solve, as a sub app.
Definition: SecantSolve.h:100
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
Definition: MooseUtils.h:380
unsigned int TagID
Definition: MooseTypes.h:210
NumericVector< Number > & solution()
Definition: SystemBase.h:196
TagID _secondary_fxn_m2_tagid
Vector tag id for the result of the last but one secondary solve, as a sub app.
Definition: SecantSolve.h:103
PARALLEL
static InputParameters validParams()
Definition: SecantSolve.C:19
void setPostprocessorValueByName(const PostprocessorName &name, const PostprocessorValue &value, std::size_t t_index=0)
Set the value of a PostprocessorValue.
TagID _secondary_xn_m2_tagid
Vector tag id for the solution variable two primary solves ago, as a sub app.
Definition: SecantSolve.h:106
virtual void saveVariableValues(const bool primary) override final
Saves the current values of the variables, and update the old(er) vectors.
Definition: SecantSolve.C:70
virtual TagID addVectorTag(const TagName &tag_name, const Moose::VectorTagType type=Moose::VECTOR_TAG_RESIDUAL)
Create a Tag.
Definition: SubProblem.C:92
TagID _xn_m2_tagid
Vector tag id for the solution variable two solves ago, as a main app.
Definition: SecantSolve.h:94
TagID _fxn_m1_tagid
Vector tag id for the most recent solution variable, pre-Secant transform, as a main app...
Definition: SecantSolve.h:85
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
const TagName PREVIOUS_FP_SOLUTION_TAG
Definition: MooseTypes.C:29
std::vector< PostprocessorName > _secondary_transformed_pps
Postprocessors to be relaxed outside of fixed point iteration (used as a subapp)
virtual void transformPostprocessors(const bool primary) override final
Use the fixed point algorithm to transform the postprocessors.
Definition: SecantSolve.C:176
auto max(const L &left, const R &right)
const std::vector< PostprocessorName > _transformed_pps
The postprocessors (transferred or not) that are going to be relaxed.
virtual void printFixedPointConvergenceHistory(Real initial_norm, const std::vector< Real > &timestep_begin_norms, const std::vector< Real > &timestep_end_norms) const override final
Print the convergence history of the coupling, at every fixed point iteration.
Definition: SecantSolve.C:273
NumericVector< Number > & addVector(const std::string &vector_name, const bool project, const libMesh::ParallelType type)
Adds a solution length vector to the system.
void update()
Update the system (doing libMesh magic)
Definition: SystemBase.C:1235
TagID _secondary_fxn_m1_tagid
Vector tag id for the most recent solution variable, pre-Secant transform, as a sub app...
Definition: SecantSolve.h:97
const Real _relax_factor
Relaxation factor for fixed point Iteration.
virtual void transformVariables(const std::set< dof_id_type > &transformed_dofs, const bool primary) override final
Use the fixed point algorithm to transform the variables.
Definition: SecantSolve.C:222
Executioners are objects that do the actual work of solving your problem.
Definition: Executioner.h:30
Real _secondary_relaxation_factor
Relaxation factor outside of fixed point iteration (used as a subapp)
virtual void savePostprocessorValues(const bool primary) override final
Saves the current values of the postprocessors, and update the old(er) vectors.
Definition: SecantSolve.C:126
std::string stringify(const T &t)
conversion to string
Definition: Conversion.h:64
virtual void needSolutionState(const unsigned int state, Moose::SolutionIterationType iteration_type=Moose::SolutionIterationType::Time, libMesh::ParallelType parallel_type=GHOSTED)
Registers that the solution state state is needed.
Definition: SystemBase.C:1442
virtual void close()=0
TagID _fxn_m2_tagid
Vector tag id for the result of the last but one solve, as a main app.
Definition: SecantSolve.h:91
virtual const PostprocessorValue & getPostprocessorValueByName(const PostprocessorName &name) const
Retrieve the value of the Postprocessor.
unsigned int _fixed_point_it
TagID _xn_m1_tagid
Vector tag id for the solution variable before the latest solve, as a main app.
Definition: SecantSolve.h:88
SystemBase & _solver_sys
Reference to a system for creating vectors as needed for the solve, etc.
Definition: SolveObject.h:55
static InputParameters validParams()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static std::string outputNorm(const Real &old_norm, const Real &norm, const unsigned int precision=6)
A helper function for outputting norms in color.
Definition: Console.C:623
const TagID INVALID_TAG_ID
Definition: MooseTypes.C:23
virtual void allocateStorage(const bool primary) override final
Allocate storage for the fixed point algorithm.
Definition: SecantSolve.C:37
virtual bool useFixedPointAlgorithmUpdateInsteadOfPicard(const bool primary) override final
Use the fixed point algorithm transform instead of simply using the Picard update.
Definition: SecantSolve.C:166
virtual void set(const numeric_index_type i, const Number value)=0
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
virtual NumericVector< Number > & getVector(const std::string &name)
Get a raw NumericVector by name.
Definition: SystemBase.C:925
SecantSolve(Executioner &ex)
Definition: SecantSolve.C:26
unsigned int _main_fixed_point_it
fixed point iteration counter for the main app