Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #ifndef LIBMESH_CONTINUATION_SYSTEM_H
21 : #define LIBMESH_CONTINUATION_SYSTEM_H
22 :
23 : // Local Includes
24 : #include "libmesh/fem_system.h"
25 :
26 : namespace libMesh
27 : {
28 :
29 : // Forward Declarations
30 : template <typename T> class LinearSolver;
31 : class NewtonSolver;
32 :
33 : /**
34 : * This class inherits from the FEMSystem. It can be
35 : * used to do arclength continuation. Most of the ideas and
36 : * the notation here come from HB Keller's 1977 paper:
37 : *
38 : * \verbatim
39 : * @InProceedings{Kell-1977,
40 : * author = {H.~B.~Keller},
41 : * title = {{Numerical solution of bifurcation and nonlinear eigenvalue problems}},
42 : * booktitle = {Applications of Bifurcation Theory, P.~H.~Rabinowitz (ed.)},
43 : * year = 1977,
44 : * publisher = {Academic Press},
45 : * pages = {359--389},
46 : * notes = {QA 3 U45 No.\ 38 (PMA)}
47 : * }
48 : * \endverbatim
49 : *
50 : * \author John W. Peterson
51 : * \date 2007
52 : */
53 0 : class ContinuationSystem : public FEMSystem
54 : {
55 : public:
56 : /**
57 : * Constructor. Optionally initializes required
58 : * data structures.
59 : */
60 : ContinuationSystem (EquationSystems & es,
61 : const std::string & name,
62 : const unsigned int number);
63 :
64 : /**
65 : * Special functions.
66 : * - This class has the same restrictions as its base class.
67 : * - The destructor is defaulted out-of-line.
68 : */
69 : ContinuationSystem (const ContinuationSystem &) = delete;
70 : ContinuationSystem & operator= (const ContinuationSystem &) = delete;
71 : ContinuationSystem (ContinuationSystem &&) = delete;
72 : ContinuationSystem & operator= (ContinuationSystem &&) = delete;
73 : virtual ~ContinuationSystem ();
74 :
75 : /**
76 : * The type of system.
77 : */
78 : typedef ContinuationSystem sys_type;
79 :
80 : /**
81 : * The type of the parent.
82 : */
83 : typedef FEMSystem Parent;
84 :
85 : /**
86 : * Clear all the data structures associated with
87 : * the system.
88 : */
89 : virtual void clear () override;
90 :
91 : /**
92 : * Perform a standard "solve" of the system, without doing continuation.
93 : */
94 : virtual void solve () override;
95 :
96 : /**
97 : * Perform a continuation solve of the system. In general, you can only
98 : * begin the continuation solves after either reading in or solving for two
99 : * previous values of the control parameter. The prior two solutions are
100 : * required for starting up the continuation method.
101 : */
102 : void continuation_solve();
103 :
104 : /**
105 : * Call this function after a continuation solve to compute the tangent and
106 : * get the next initial guess.
107 : */
108 : void advance_arcstep();
109 :
110 : /**
111 : * The continuation parameter must be a member variable of the
112 : * derived class, and the "continuation_parameter" pointer defined
113 : * here must be a pointer to that variable. This is how the
114 : * continuation system updates the derived class's continuation
115 : * parameter.
116 : *
117 : * Also sometimes referred to as "lambda" in the code comments.
118 : */
119 : Real * continuation_parameter;
120 :
121 : /**
122 : * If quiet==false, the System prints extra information about what
123 : * it is doing.
124 : */
125 : bool quiet;
126 :
127 : /**
128 : * Sets (initializes) the max-allowable ds value and the current ds value.
129 : * Call this before beginning arclength continuation. The default max stepsize
130 : * is 0.1
131 : */
132 : void set_max_arclength_stepsize(Real maxds) { ds=maxds; ds_current=maxds; }
133 :
134 : /**
135 : * How tightly should the Newton iterations attempt to converge delta_lambda.
136 : * Defaults to 1.e-6.
137 : */
138 : double continuation_parameter_tolerance;
139 :
140 : /**
141 : * How tightly should the Newton iterations attempt to converge ||delta_u||
142 : * Defaults to 1.e-6.
143 : */
144 : double solution_tolerance;
145 :
146 : /**
147 : * How much to try to reduce the residual by at the first (inexact) Newton step.
148 : * This is frequently something large like 1/2 in an inexact Newton method, to
149 : * prevent oversolving.
150 : */
151 : double initial_newton_tolerance;
152 :
153 : /**
154 : * Stores the current solution and continuation parameter
155 : * (as "previous_u" and "old_continuation_parameter") for later referral.
156 : * You may need to call this e.g. after the first regular solve, in order
157 : * to store the first solution, before computing a second solution and
158 : * beginning arclength continuation.
159 : */
160 : void save_current_solution();
161 :
162 : /**
163 : * The system also keeps track of the old value of the continuation parameter.
164 : */
165 : Real old_continuation_parameter;
166 :
167 : /**
168 : * The minimum-allowable value of the continuation parameter. The Newton iterations
169 : * will quit if the continuation parameter falls below this value.
170 : */
171 : Real min_continuation_parameter;
172 :
173 : /**
174 : * The maximum-allowable value of the continuation parameter. The Newton iterations
175 : * will quit if the continuation parameter goes above this value. If this value is zero,
176 : * there is no maximum value for the continuation parameter.
177 : */
178 : Real max_continuation_parameter;
179 :
180 : /**
181 : * Arclength normalization parameter. Defaults to 1.0 (no normalization). Used to
182 : * ensure that one term in the arclength constraint equation does not wash out all
183 : * the others.
184 : */
185 : Real Theta;
186 :
187 : /**
188 : * Another normalization parameter, which is described in the LOCA manual.
189 : * This one is designed to maintain a relatively "fixed" value of d(lambda)/ds.
190 : * It is initially 1.0 and is updated after each solve.
191 : */
192 : Real Theta_LOCA;
193 :
194 : /**
195 : * Another scaling parameter suggested by the LOCA people. This one attempts
196 : * to shrink the stepsize ds whenever the angle between the previous two
197 : * tangent vectors gets large.
198 : */
199 : //Real tau;
200 :
201 : /**
202 : * Number of (Newton) backtracking steps to attempt if a Newton step does not
203 : * reduce the residual. This is backtracking within a *single* Newton step,
204 : * if you want to try a smaller arcstep, set n_arclength_reductions > 0.
205 : */
206 : unsigned int n_backtrack_steps;
207 :
208 : /**
209 : * Number of arclength reductions to try when Newton fails to reduce
210 : * the residual. For each arclength reduction, the arcstep size is
211 : * cut in half.
212 : */
213 : unsigned int n_arclength_reductions;
214 :
215 : /**
216 : * The minimum-allowed steplength, defaults to 1.e-8.
217 : */
218 : Real ds_min;
219 :
220 : /**
221 : * The code provides the ability to select from different predictor
222 : * schemes for getting the initial guess for the solution at the next
223 : * point on the solution arc.
224 : */
225 : enum Predictor {
226 : /**
227 : * First-order Euler predictor
228 : */
229 : Euler,
230 :
231 : /**
232 : * Second-order explicit Adams-Bashforth predictor
233 : */
234 : AB2,
235 :
236 : /**
237 : * Invalid predictor
238 : */
239 : Invalid_Predictor
240 : };
241 :
242 : Predictor predictor;
243 :
244 : /**
245 : * A measure of how rapidly one should attempt to grow the arclength
246 : * stepsize based on the number of Newton iterations required to solve
247 : * the problem. Default value is 1.0, if set to zero, will not try to
248 : * grow or shrink the arclength stepsize based on the number of Newton
249 : * iterations required.
250 : */
251 : Real newton_stepgrowth_aggressiveness;
252 :
253 : /**
254 : * True by default, the Newton progress check allows the Newton loop
255 : * to exit if half the allowed iterations have elapsed without a reduction
256 : * in the *initial* residual. In our experience this usually means the
257 : * Newton steps are going to fail eventually and we could save some time
258 : * by quitting early.
259 : */
260 : bool newton_progress_check;
261 :
262 : protected:
263 : /**
264 : * Initializes the member data fields associated with
265 : * the system, so that, e.g., \p assemble() may be used.
266 : */
267 : virtual void init_data () override;
268 :
269 : /**
270 : * There are (so far) two different vectors which may be assembled using the assembly routine:
271 : * 1.) The Residual = the normal PDE weighted residual
272 : * 2.) G_Lambda = the derivative wrt the parameter lambda of the PDE weighted residual
273 : *
274 : * It is up to the derived class to handle writing separate assembly code for the different cases.
275 : * Usually something like:
276 : * switch (rhs_mode)
277 : * {
278 : * case Residual:
279 : * {
280 : * Fu(i) += ... // normal PDE residual
281 : * break;
282 : * }
283 : *
284 : * case G_Lambda:
285 : * {
286 : * Fu(i) += ... // derivative wrt control parameter
287 : * break;
288 : * }
289 : */
290 : enum RHS_Mode {Residual,
291 : G_Lambda};
292 :
293 : RHS_Mode rhs_mode;
294 :
295 :
296 :
297 : private:
298 : /**
299 : * Before starting arclength continuation, we need at least 2 prior solutions
300 : * (both solution and u_previous should be filled with meaningful values)
301 : * And we need to initialize the tangent vector. This only needs to be
302 : * called once.
303 : */
304 : void initialize_tangent();
305 :
306 : /**
307 : * Special solve algorithm for solving the tangent system.
308 : */
309 : void solve_tangent();
310 :
311 : /**
312 : * This function (which assumes the most recent tangent vector has been
313 : * computed) updates the solution and the control parameter with the initial
314 : * guess for the next point on the continuation path.
315 : */
316 : void update_solution();
317 :
318 : /**
319 : * A centralized function for setting the normalization parameter theta
320 : */
321 : void set_Theta();
322 :
323 : /**
324 : * A centralized function for setting the other normalization parameter, i.e.
325 : * the one suggested by the LOCA developers.
326 : */
327 : void set_Theta_LOCA();
328 :
329 : /**
330 : * Applies the predictor to the current solution to get a guess for the
331 : * next solution.
332 : */
333 : void apply_predictor();
334 :
335 : /**
336 : * Extra work vectors used by the continuation algorithm.
337 : * These are added to the system by the init_data() routine.
338 : *
339 : *
340 : * The "solution" tangent vector du/ds.
341 : */
342 : NumericVector<Number> * du_ds;
343 :
344 : /**
345 : * The value of du_ds from the previous solution
346 : */
347 : NumericVector<Number> * previous_du_ds;
348 :
349 : /**
350 : * The solution at the previous value of the continuation variable.
351 : */
352 : NumericVector<Number> * previous_u;
353 :
354 : /**
355 : * Temporary vector "y" ... stores -du/dlambda, the solution of \f$ Ay=G_{\lambda} \f$.
356 : */
357 : NumericVector<Number> * y;
358 :
359 : /**
360 : * Temporary vector "y_old" ... stores the previous value of -du/dlambda,
361 : * which is the solution of \f$ Ay=G_{\lambda} \f$.
362 : */
363 : NumericVector<Number> * y_old;
364 :
365 : /**
366 : * Temporary vector "z" ... the solution of \f$ Az = -G \f$
367 : */
368 : NumericVector<Number> * z;
369 :
370 : /**
371 : * Temporary vector "delta u" ... the Newton step update in our custom
372 : * augmented PDE solve.
373 : */
374 : NumericVector<Number> * delta_u;
375 :
376 : /**
377 : * False until initialize_tangent() is called
378 : */
379 : bool tangent_initialized;
380 :
381 : /**
382 : * A pointer to the underlying Newton solver used by the \p DiffSystem.
383 : * From this pointer, we can get access to all the parameters and options
384 : * which are available to the "normal" Newton solver.
385 : */
386 : NewtonSolver * newton_solver;
387 :
388 : /**
389 : * The most recent value of the derivative of the continuation parameter
390 : * with respect to s. We use "lambda" here for shortness of notation, lambda
391 : * always refers to the continuation parameter.
392 : */
393 : Real dlambda_ds;
394 :
395 : /**
396 : * The initial arclength stepsize, selected by the user. This is
397 : * the max-allowable arclength stepsize, but the algorithm may frequently
398 : * reduce ds near turning points.
399 : */
400 : Real ds;
401 :
402 : /**
403 : * Value of stepsize currently in use. Will not exceed user-provided maximum
404 : * arclength stepsize ds.
405 : */
406 : Real ds_current;
407 :
408 : /**
409 : * The old parameter tangent value.
410 : */
411 : Real previous_dlambda_ds;
412 :
413 : /**
414 : * The previous arcstep length used.
415 : */
416 : Real previous_ds;
417 :
418 : /**
419 : * Loop counter for nonlinear (Newton) iteration loop.
420 : */
421 : unsigned int newton_step;
422 : };
423 :
424 : } // namespace libMesh
425 :
426 : #endif // LIBMESH_CONTINUATION_SYSTEM_H
|
