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 : // LibMesh includes
19 : #include "libmesh/equation_systems.h"
20 : #include "libmesh/continuation_system.h"
21 : #include "libmesh/linear_solver.h"
22 : #include "libmesh/time_solver.h"
23 : #include "libmesh/newton_solver.h"
24 : #include "libmesh/numeric_vector.h"
25 : #include "libmesh/sparse_matrix.h"
26 :
27 : namespace libMesh
28 : {
29 :
30 0 : ContinuationSystem::ContinuationSystem (EquationSystems & es,
31 : const std::string & name_in,
32 0 : const unsigned int number_in) :
33 : Parent(es, name_in, number_in),
34 0 : continuation_parameter(nullptr),
35 0 : quiet(true),
36 0 : continuation_parameter_tolerance(1.e-6),
37 0 : solution_tolerance(1.e-6),
38 0 : initial_newton_tolerance(0.01),
39 0 : old_continuation_parameter(0.),
40 0 : min_continuation_parameter(0.),
41 0 : max_continuation_parameter(0.),
42 0 : Theta(1.),
43 0 : Theta_LOCA(1.),
44 0 : n_backtrack_steps(5),
45 0 : n_arclength_reductions(5),
46 0 : ds_min(1.e-8),
47 0 : predictor(Euler),
48 0 : newton_stepgrowth_aggressiveness(1.),
49 0 : newton_progress_check(true),
50 0 : rhs_mode(Residual),
51 0 : tangent_initialized(false),
52 0 : newton_solver(nullptr),
53 0 : dlambda_ds(0.707),
54 0 : ds(0.1),
55 0 : ds_current(0.1),
56 0 : previous_dlambda_ds(0.),
57 0 : previous_ds(0.),
58 0 : newton_step(0)
59 : {
60 : // Warn about using untested code
61 : libmesh_experimental();
62 :
63 : // linear_solver is now in the ImplicitSystem base class, but we are
64 : // going to keep using it basically the way we did before it was
65 : // moved.
66 0 : linear_solver = LinearSolver<Number>::build(es.comm());
67 :
68 0 : if (this->prefix_with_name())
69 0 : linear_solver->init(this->prefix().c_str());
70 : else
71 0 : linear_solver->init();
72 :
73 0 : linear_solver->init_names(*this);
74 0 : }
75 :
76 :
77 :
78 :
79 0 : ContinuationSystem::~ContinuationSystem () = default;
80 :
81 :
82 :
83 0 : void ContinuationSystem::clear()
84 : {
85 : // Call the Parent's clear function
86 0 : Parent::clear();
87 0 : }
88 :
89 :
90 :
91 0 : void ContinuationSystem::init_data ()
92 : {
93 : // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
94 0 : du_ds = &(add_vector("du_ds"));
95 :
96 : // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
97 0 : previous_du_ds = &(add_vector("previous_du_ds"));
98 :
99 : // Add a vector to keep track of the previous nonlinear solution
100 : // at the old value of lambda.
101 0 : previous_u = &(add_vector("previous_u"));
102 :
103 : // Add a vector to keep track of the temporary solution "y" of Ay=G_{\lambda}.
104 0 : y = &(add_vector("y"));
105 :
106 : // Add a vector to keep track of the "old value" of "y" which is the solution of Ay=G_{\lambda}.
107 0 : y_old = &(add_vector("y_old"));
108 :
109 : // Add a vector to keep track of the temporary solution "z" of Az=-G.
110 0 : z = &(add_vector("z"));
111 :
112 : // Add a vector to keep track of the Newton update during the constrained PDE solves.
113 0 : delta_u = &(add_vector("delta_u"));
114 :
115 : // Call the Parent's initialization routine.
116 0 : Parent::init_data();
117 0 : }
118 :
119 :
120 :
121 :
122 0 : void ContinuationSystem::solve()
123 : {
124 : // Set the Residual RHS mode, and call the normal solve routine.
125 0 : rhs_mode = Residual;
126 0 : DifferentiableSystem::solve();
127 0 : }
128 :
129 :
130 :
131 :
132 0 : void ContinuationSystem::initialize_tangent()
133 : {
134 : // Be sure the tangent was not already initialized.
135 0 : libmesh_assert (!tangent_initialized);
136 :
137 : // Compute delta_s_zero, the initial arclength traveled during the
138 : // first step. Here we assume that previous_u and lambda_old store
139 : // the previous solution and control parameter. You may need to
140 : // read in an old solution (or solve the non-continuation system)
141 : // first and call save_current_solution() before getting here.
142 :
143 : // 1.) Compute delta_s_zero as ||u|| - ||u_old|| + ...
144 : // Compute norms of the current and previous solutions
145 : // Real norm_u = solution->l2_norm();
146 : // Real norm_previous_u = previous_u->l2_norm();
147 :
148 : // if (!quiet)
149 : // {
150 : // libMesh::out << "norm_u=" << norm_u << std::endl;
151 : // libMesh::out << "norm_previous_u=" << norm_previous_u << std::endl;
152 : // }
153 :
154 : // if (norm_u == norm_previous_u)
155 : // {
156 : // libMesh::err << "Warning, it appears u and previous_u are the "
157 : // << "same, are you sure this is correct?"
158 : // << "It's possible you forgot to set one or the other..."
159 : // << std::endl;
160 : // }
161 :
162 : // Real delta_s_zero = std::sqrt(
163 : // (norm_u - norm_previous_u)*(norm_u - norm_previous_u) +
164 : // (*continuation_parameter-old_continuation_parameter)*
165 : // (*continuation_parameter-old_continuation_parameter)
166 : // );
167 :
168 : // // 2.) Compute delta_s_zero as ||u -u_old|| + ...
169 : // *delta_u = *solution;
170 : // delta_u->add(-1., *previous_u);
171 : // delta_u->close();
172 : // Real norm_delta_u = delta_u->l2_norm();
173 : // Real norm_u = solution->l2_norm();
174 : // Real norm_previous_u = previous_u->l2_norm();
175 :
176 : // // Scale norm_delta_u by the bigger of either norm_u or norm_previous_u
177 : // norm_delta_u /= std::max(norm_u, norm_previous_u);
178 :
179 : // if (!quiet)
180 : // {
181 : // libMesh::out << "norm_u=" << norm_u << std::endl;
182 : // libMesh::out << "norm_previous_u=" << norm_previous_u << std::endl;
183 : // //libMesh::out << "norm_delta_u=" << norm_delta_u << std::endl;
184 : // libMesh::out << "norm_delta_u/max(|u|,|u_old|)=" << norm_delta_u << std::endl;
185 : // libMesh::out << "|norm_u-norm_previous_u|=" << std::abs(norm_u - norm_previous_u) << std::endl;
186 : // }
187 :
188 : // const Real dlambda = *continuation_parameter-old_continuation_parameter;
189 :
190 : // if (!quiet)
191 : // libMesh::out << "dlambda=" << dlambda << std::endl;
192 :
193 : // Real delta_s_zero = std::sqrt(
194 : // (norm_delta_u*norm_delta_u) +
195 : // (dlambda*dlambda)
196 : // );
197 :
198 : // if (!quiet)
199 : // libMesh::out << "delta_s_zero=" << delta_s_zero << std::endl;
200 :
201 : // 1.) + 2.)
202 : // // Now approximate the initial tangent d(lambda)/ds
203 : // this->dlambda_ds = (*continuation_parameter-old_continuation_parameter) / delta_s_zero;
204 :
205 :
206 : // // We can also approximate the deriv. wrt s by finite differences:
207 : // // du/ds = (u1 - u0) / delta_s_zero.
208 : // // FIXME: Use delta_u from above if we decide to keep that method.
209 : // *du_ds = *solution;
210 : // du_ds->add(-1., *previous_u);
211 : // du_ds->scale(1./delta_s_zero);
212 : // du_ds->close();
213 :
214 :
215 : // 3.) Treating (u-previous_u)/(lambda - lambda_old) as an approximation to du/d(lambda),
216 : // we follow the same technique as Carnes and Shadid.
217 : // const Real dlambda = *continuation_parameter-old_continuation_parameter;
218 : // libmesh_assert_greater (dlambda, 0.);
219 :
220 : // // Use delta_u for temporary calculation of du/d(lambda)
221 : // *delta_u = *solution;
222 : // delta_u->add(-1., *previous_u);
223 : // delta_u->scale(1. / dlambda);
224 : // delta_u->close();
225 :
226 : // // Determine initial normalization parameter
227 : // const Real solution_size = std::max(solution->l2_norm(), previous_u->l2_norm());
228 : // if (solution_size > 1.)
229 : // {
230 : // Theta = 1./solution_size;
231 :
232 : // if (!quiet)
233 : // libMesh::out << "Setting Normalization Parameter Theta=" << Theta << std::endl;
234 : // }
235 :
236 : // // Compute d(lambda)/ds
237 : // // The correct sign of d(lambda)/ds should be positive, since we assume that (lambda > lambda_old)
238 : // // but we could always double-check that as well.
239 : // Real norm_delta_u = delta_u->l2_norm();
240 : // this->dlambda_ds = 1. / std::sqrt(1. + Theta*Theta*norm_delta_u*norm_delta_u);
241 :
242 : // // Finally, compute du/ds = d(lambda)/ds * du/d(lambda)
243 : // *du_ds = *delta_u;
244 : // du_ds->scale(dlambda_ds);
245 : // du_ds->close();
246 :
247 :
248 : // 4.) Use normalized arclength formula to estimate delta_s_zero
249 : // // Determine initial normalization parameter
250 : // set_Theta();
251 :
252 : // // Compute (normalized) delta_s_zero
253 : // *delta_u = *solution;
254 : // delta_u->add(-1., *previous_u);
255 : // delta_u->close();
256 : // Real norm_delta_u = delta_u->l2_norm();
257 :
258 : // const Real dlambda = *continuation_parameter-old_continuation_parameter;
259 :
260 : // if (!quiet)
261 : // libMesh::out << "dlambda=" << dlambda << std::endl;
262 :
263 : // Real delta_s_zero = std::sqrt(
264 : // (Theta_LOCA*Theta_LOCA*Theta*norm_delta_u*norm_delta_u) +
265 : // (dlambda*dlambda)
266 : // );
267 : // *du_ds = *delta_u;
268 : // du_ds->scale(1./delta_s_zero);
269 : // dlambda_ds = dlambda / delta_s_zero;
270 :
271 : // if (!quiet)
272 : // {
273 : // libMesh::out << "delta_s_zero=" << delta_s_zero << std::endl;
274 : // libMesh::out << "initial d(lambda)/ds|_0 = " << dlambda_ds << std::endl;
275 : // libMesh::out << "initial ||du_ds||_0 = " << du_ds->l2_norm() << std::endl;
276 : // }
277 :
278 : // // FIXME: Also store the initial finite-differenced approximation to -du/dlambda as y.
279 : // // We stick to the convention of storing negative y, since that is what we typically
280 : // // solve for anyway.
281 : // *y = *delta_u;
282 : // y->scale(-1./dlambda);
283 : // y->close();
284 :
285 :
286 :
287 : // 5.) Assume dlambda/ds_0 ~ 1/sqrt(2) and determine the value of Theta_LOCA which
288 : // will satisfy this criterion
289 :
290 : // Initial change in parameter
291 0 : const Real dlambda = *continuation_parameter-old_continuation_parameter;
292 0 : libmesh_assert_not_equal_to (dlambda, 0.0);
293 :
294 : // Ideal initial value of dlambda_ds
295 0 : dlambda_ds = 1. / std::sqrt(2.);
296 0 : if (dlambda < 0.)
297 0 : dlambda_ds *= -1.;
298 :
299 : // This also implies the initial value of ds
300 0 : ds_current = dlambda / dlambda_ds;
301 :
302 0 : if (!quiet)
303 0 : libMesh::out << "Setting ds_current|_0=" << ds_current << std::endl;
304 :
305 : // Set y = -du/dlambda using finite difference approximation
306 0 : *y = *solution;
307 0 : y->add(-1., *previous_u);
308 0 : y->scale(-1./dlambda);
309 0 : y->close();
310 0 : const Real ynorm=y->l2_norm();
311 :
312 : // Finally, set the value of du_ds to be used in the upcoming
313 : // tangent calculation. du/ds = du/dlambda * dlambda/ds
314 0 : *du_ds = *y;
315 0 : du_ds->scale(-dlambda_ds);
316 0 : du_ds->close();
317 :
318 : // Determine additional solution normalization parameter
319 : // (Since we just set du/ds, it will be: ||du||*||du/ds||)
320 0 : set_Theta();
321 :
322 : // The value of Theta_LOCA which makes dlambda_ds = 1/sqrt(2),
323 : // assuming our Theta = ||du||^2.
324 : // Theta_LOCA = std::abs(dlambda);
325 :
326 : // Assuming general Theta
327 0 : Theta_LOCA = std::sqrt(1./Theta/ynorm/ynorm);
328 :
329 :
330 0 : if (!quiet)
331 : {
332 0 : libMesh::out << "Setting initial Theta_LOCA = " << Theta_LOCA << std::endl;
333 0 : libMesh::out << "Theta_LOCA^2*Theta = " << Theta_LOCA*Theta_LOCA*Theta << std::endl;
334 0 : libMesh::out << "initial d(lambda)/ds|_0 = " << dlambda_ds << std::endl;
335 0 : libMesh::out << "initial ||du_ds||_0 = " << du_ds->l2_norm() << std::endl;
336 : }
337 :
338 :
339 :
340 : // OK, we estimated the tangent at point u0.
341 : // Now, to estimate the tangent at point u1, we call the solve_tangent routine.
342 :
343 : // Set the flag which tells us the method has been initialized.
344 0 : tangent_initialized = true;
345 :
346 0 : solve_tangent();
347 :
348 : // Advance the solution and the parameter to the next value.
349 0 : update_solution();
350 0 : }
351 :
352 :
353 :
354 :
355 :
356 :
357 : // This is most of the "guts" of this class. This is where we implement
358 : // our custom Newton iterations and perform most of the solves.
359 0 : void ContinuationSystem::continuation_solve()
360 : {
361 : // Be sure the user has set the continuation parameter pointer
362 0 : libmesh_error_msg_if(!continuation_parameter,
363 : "You must set the continuation_parameter pointer "
364 : "to a member variable of the derived class, preferably in the "
365 : "Derived class's init_data function. This is how the ContinuationSystem "
366 : "updates the continuation parameter.");
367 :
368 : // Use extra precision for all the numbers printed in this function.
369 0 : std::streamsize old_precision = libMesh::out.precision();
370 0 : libMesh::out.precision(16);
371 0 : libMesh::out.setf(std::ios_base::scientific);
372 :
373 : // We can't start solving the augmented PDE system unless the tangent
374 : // vectors have been initialized. This only needs to occur once.
375 0 : if (!tangent_initialized)
376 0 : initialize_tangent();
377 :
378 : // Save the old value of -du/dlambda. This will be used after the Newton iterations
379 : // to compute the angle between previous tangent vectors. This cosine of this angle is
380 : //
381 : // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
382 : //
383 : // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
384 : // when we are approaching a turning point. Note that it can only shrink the step size.
385 0 : *y_old = *y;
386 :
387 : // Set pointer to underlying Newton solver
388 0 : if (!newton_solver)
389 0 : newton_solver = cast_ptr<NewtonSolver *> (this->time_solver->diff_solver().get());
390 :
391 : // A pair for catching return values from linear system solves.
392 0 : std::pair<unsigned int, Real> rval;
393 :
394 : // Convergence flag for the entire arcstep
395 0 : bool arcstep_converged = false;
396 :
397 : // Begin loop over arcstep reductions.
398 0 : for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
399 : {
400 0 : if (!quiet)
401 : {
402 0 : libMesh::out << "Current arclength stepsize, ds_current=" << ds_current << std::endl;
403 0 : libMesh::out << "Current parameter value, lambda=" << *continuation_parameter << std::endl;
404 : }
405 :
406 : // Upon exit from the nonlinear loop, the newton_converged flag
407 : // will tell us the convergence status of Newton's method.
408 0 : bool newton_converged = false;
409 :
410 : // The nonlinear residual before *any* nonlinear steps have been taken.
411 0 : Real nonlinear_residual_firststep = 0.;
412 :
413 : // The nonlinear residual from the current "k" Newton step, before the Newton step
414 0 : Real nonlinear_residual_beforestep = 0.;
415 :
416 : // The nonlinear residual from the current "k" Newton step, after the Newton step
417 0 : Real nonlinear_residual_afterstep = 0.;
418 :
419 : // The linear solver tolerance, can be updated dynamically at each Newton step.
420 0 : double current_linear_tolerance = 0.;
421 :
422 : // The nonlinear loop
423 0 : for (newton_step=0; newton_step<newton_solver->max_nonlinear_iterations; ++newton_step)
424 : {
425 0 : libMesh::out << "\n === Starting Newton step " << newton_step << " ===" << std::endl;
426 :
427 : // Set the linear system solver tolerance
428 : // // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
429 : // const Real residual_multiple = 1.e-4;
430 : // Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;
431 :
432 : // // But if the current residual isn't small, don't let the solver exit with zero iterations!
433 : // if (current_linear_tolerance > 1.)
434 : // current_linear_tolerance = residual_multiple;
435 :
436 : // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
437 0 : if (newton_step==0)
438 : {
439 : // At first step, only try reducing the residual by a small amount
440 0 : current_linear_tolerance = initial_newton_tolerance;//0.01;
441 : }
442 :
443 : else
444 : {
445 : // The new tolerance is based on the ratio of the most recent tolerances
446 0 : const Real alp=0.5*(1.+std::sqrt(5.));
447 0 : const Real gam=0.9;
448 :
449 0 : libmesh_assert_not_equal_to (nonlinear_residual_beforestep, 0.0);
450 0 : libmesh_assert_not_equal_to (nonlinear_residual_afterstep, 0.0);
451 :
452 0 : current_linear_tolerance =
453 0 : double(std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
454 0 : current_linear_tolerance*current_linear_tolerance));
455 :
456 : // Don't let it get ridiculously small!!
457 0 : if (current_linear_tolerance < 1.e-12)
458 0 : current_linear_tolerance = 1.e-12;
459 : }
460 :
461 0 : if (!quiet)
462 0 : libMesh::out << "Using current_linear_tolerance=" << current_linear_tolerance << std::endl;
463 :
464 :
465 : // Assemble the residual (and Jacobian).
466 0 : rhs_mode = Residual;
467 0 : assembly(true, // Residual
468 0 : true); // Jacobian
469 0 : rhs->close();
470 :
471 : // Save the current nonlinear residual. We don't need to recompute the residual unless
472 : // this is the first step, since it was already computed as part of the convergence check
473 : // at the end of the last loop iteration.
474 0 : if (newton_step==0)
475 : {
476 0 : nonlinear_residual_beforestep = rhs->l2_norm();
477 :
478 : // Store the residual before any steps have been taken. This will *not*
479 : // be updated at each step, and can be used to see if any progress has
480 : // been made from the initial residual at later steps.
481 0 : nonlinear_residual_firststep = nonlinear_residual_beforestep;
482 :
483 0 : const Real old_norm_u = solution->l2_norm();
484 0 : libMesh::out << " (before step) ||R||_{L2} = " << nonlinear_residual_beforestep << std::endl;
485 0 : libMesh::out << " (before step) ||R||_{L2}/||u|| = " << nonlinear_residual_beforestep / old_norm_u << std::endl;
486 :
487 : // In rare cases (very small arcsteps), it's possible that the residual is
488 : // already below our absolute linear tolerance.
489 0 : if (nonlinear_residual_beforestep < solution_tolerance)
490 : {
491 0 : if (!quiet)
492 0 : libMesh::out << "Initial guess satisfied linear tolerance, exiting with zero Newton iterations!" << std::endl;
493 :
494 : // Since we go straight from here to the solve of the next tangent, we
495 : // have to close the matrix before it can be assembled again.
496 0 : matrix->close();
497 0 : newton_converged=true;
498 0 : break; // out of Newton iterations, with newton_converged=true
499 : }
500 : }
501 :
502 : else
503 : {
504 0 : nonlinear_residual_beforestep = nonlinear_residual_afterstep;
505 : }
506 :
507 :
508 : // Solve the linear system G_u*z = G
509 : // Initial guess?
510 0 : z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
511 0 : z->close();
512 :
513 : // It's possible that we have selected the current_linear_tolerance so large that
514 : // a guess of z=zero yields a linear system residual |Az + R| small enough that the
515 : // linear solver exits in zero iterations. If this happens, we will reduce the
516 : // current_linear_tolerance until the linear solver does at least 1 iteration.
517 0 : do
518 : {
519 : rval =
520 0 : linear_solver->solve(*matrix,
521 0 : *z,
522 0 : *rhs,
523 : //1.e-12,
524 : current_linear_tolerance,
525 0 : newton_solver->max_linear_iterations); // max linear iterations
526 :
527 0 : if (rval.first==0)
528 : {
529 0 : if (newton_step==0)
530 : {
531 0 : libMesh::out << "Repeating initial solve with smaller linear tolerance!" << std::endl;
532 0 : current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
533 : }
534 : else
535 : // We shouldn't get here ... it means the linear solver did no work on a Newton
536 : // step other than the first one. If this happens, we need to think more about our
537 : // tolerance selection.
538 0 : libmesh_error_msg("Linear solver did no work!");
539 : }
540 :
541 0 : } while (rval.first==0);
542 :
543 :
544 0 : if (!quiet)
545 0 : libMesh::out << " G_u*z = G solver converged at step "
546 0 : << rval.first
547 0 : << " linear tolerance = "
548 0 : << rval.second
549 0 : << "."
550 0 : << std::endl;
551 :
552 : // Sometimes (I am not sure why) the linear solver exits after zero iterations.
553 : // Perhaps it is hitting PETSc's divergence tolerance dtol??? If this occurs,
554 : // we should break out of the Newton iteration loop because nothing further is
555 : // going to happen... Of course if the tolerance is already small enough after
556 : // zero iterations (how can this happen?!) we should not quit.
557 0 : if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
558 : {
559 0 : if (!quiet)
560 0 : libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
561 :
562 : // Try to find out the reason for convergence/divergence
563 0 : linear_solver->print_converged_reason();
564 :
565 0 : break; // out of Newton iterations
566 : }
567 :
568 : // Note: need to scale z by -1 since our code always solves Jx=R
569 : // instead of Jx=-R.
570 0 : z->scale(-1.);
571 0 : z->close();
572 :
573 :
574 :
575 :
576 :
577 :
578 : // Assemble the G_Lambda vector, skip residual.
579 0 : rhs_mode = G_Lambda;
580 :
581 : // Assemble both rhs and Jacobian
582 0 : assembly(true, // Residual
583 0 : false); // Jacobian
584 :
585 : // Not sure if this is really necessary
586 0 : rhs->close();
587 0 : const Real yrhsnorm=rhs->l2_norm();
588 0 : libmesh_error_msg_if(yrhsnorm == 0.0, "||G_Lambda|| = 0");
589 :
590 : // We select a tolerance for the y-system which is based on the inexact Newton
591 : // tolerance but scaled by an extra term proportional to the RHS (which is not -> 0 in this case)
592 0 : const double ysystemtol =
593 0 : double(current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm));
594 0 : if (!quiet)
595 0 : libMesh::out << "ysystemtol=" << ysystemtol << std::endl;
596 :
597 : // Solve G_u*y = G_{\lambda}
598 : // FIXME: Initial guess? This is really a solve for -du/dlambda so we could try
599 : // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda
600 : //*y = *solution;
601 : //y->add(-1., *previous_u);
602 : //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero...
603 : //y->close();
604 :
605 : // const unsigned int max_attempts=1;
606 : // unsigned int attempt=0;
607 : // do
608 : // {
609 : // if (!quiet)
610 : // libMesh::out << "Trying to solve tangent system, attempt " << attempt << std::endl;
611 :
612 : rval =
613 0 : linear_solver->solve(*matrix,
614 0 : *y,
615 0 : *rhs,
616 : //1.e-12,
617 : ysystemtol,
618 0 : newton_solver->max_linear_iterations); // max linear iterations
619 :
620 0 : if (!quiet)
621 0 : libMesh::out << " G_u*y = G_{lambda} solver converged at step "
622 0 : << rval.first
623 0 : << ", linear tolerance = "
624 0 : << rval.second
625 0 : << "."
626 0 : << std::endl;
627 :
628 : // Sometimes (I am not sure why) the linear solver exits after zero iterations.
629 : // Perhaps it is hitting PETSc's divergence tolerance dtol??? If this occurs,
630 : // we should break out of the Newton iteration loop because nothing further is
631 : // going to happen...
632 0 : if ((rval.first == 0) && (rval.second > ysystemtol))
633 : {
634 0 : if (!quiet)
635 0 : libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
636 :
637 0 : break; // out of Newton iterations
638 : }
639 :
640 : // ++attempt;
641 : // } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations));
642 :
643 :
644 :
645 :
646 :
647 : // Compute N, the residual of the arclength constraint eqn.
648 : // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds
649 : // We temporarily use the delta_u vector as a temporary vector for this calculation.
650 0 : *delta_u = *solution;
651 0 : delta_u->add(-1., *previous_u);
652 :
653 : // First part of the arclength constraint
654 0 : const Number N1 = Theta_LOCA*Theta_LOCA*Theta*delta_u->dot(*du_ds);
655 0 : const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds;
656 0 : const Number N3 = ds_current;
657 :
658 0 : if (!quiet)
659 : {
660 0 : libMesh::out << " N1=" << N1 << std::endl;
661 0 : libMesh::out << " N2=" << N2 << std::endl;
662 0 : libMesh::out << " N3=" << N3 << std::endl;
663 : }
664 :
665 : // The arclength constraint value
666 0 : const Number N = N1+N2-N3;
667 :
668 0 : if (!quiet)
669 0 : libMesh::out << " N=" << N << std::endl;
670 :
671 0 : const Number duds_dot_z = du_ds->dot(*z);
672 0 : const Number duds_dot_y = du_ds->dot(*y);
673 :
674 : //libMesh::out << "duds_dot_z=" << duds_dot_z << std::endl;
675 : //libMesh::out << "duds_dot_y=" << duds_dot_y << std::endl;
676 : //libMesh::out << "dlambda_ds=" << dlambda_ds << std::endl;
677 :
678 0 : const Number delta_lambda_numerator = -(N + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z);
679 0 : const Number delta_lambda_denominator = (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y);
680 :
681 0 : libmesh_assert_not_equal_to (delta_lambda_denominator, 0.0);
682 :
683 : // Now, we are ready to compute the step delta_lambda
684 0 : const Number delta_lambda_comp = delta_lambda_numerator /
685 : delta_lambda_denominator;
686 : // Lambda is real-valued
687 0 : const Real delta_lambda = libmesh_real(delta_lambda_comp);
688 :
689 : // Knowing delta_lambda, we are ready to update delta_u
690 : // delta_u = z - delta_lambda*y
691 0 : delta_u->zero();
692 0 : delta_u->add(1., *z);
693 0 : delta_u->add(-delta_lambda, *y);
694 0 : delta_u->close();
695 :
696 : // Update the system solution and the continuation parameter.
697 0 : solution->add(1., *delta_u);
698 0 : solution->close();
699 0 : *continuation_parameter += delta_lambda;
700 :
701 : // Did the Newton step actually reduce the residual?
702 0 : rhs_mode = Residual;
703 0 : assembly(true, // Residual
704 0 : false); // Jacobian
705 0 : rhs->close();
706 0 : nonlinear_residual_afterstep = rhs->l2_norm();
707 :
708 :
709 : // In a "normal" Newton step, ||du||/||R|| > 1 since the most recent
710 : // step is where you "just were" and the current residual is where
711 : // you are now. It can occur that ||du||/||R|| < 1, but these are
712 : // likely not good cases to attempt backtracking (?).
713 0 : const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep;
714 0 : if (!quiet)
715 0 : libMesh::out << " norm_du_norm_R=" << norm_du_norm_R << std::endl;
716 :
717 :
718 : // Factor to decrease the stepsize by for backtracking
719 0 : Real newton_stepfactor = 1.;
720 :
721 0 : const bool attempt_backtracking =
722 0 : (nonlinear_residual_afterstep > solution_tolerance)
723 0 : && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)
724 0 : && (n_backtrack_steps>0)
725 0 : && (norm_du_norm_R > 1.)
726 : ;
727 :
728 : // If residual is not reduced, do Newton back tracking.
729 0 : if (attempt_backtracking)
730 : {
731 0 : if (!quiet)
732 0 : libMesh::out << "Newton step did not reduce residual." << std::endl;
733 :
734 : // back off the previous step.
735 0 : solution->add(-1., *delta_u);
736 0 : solution->close();
737 0 : *continuation_parameter -= delta_lambda;
738 :
739 : // Backtracking: start cutting the Newton stepsize by halves until
740 : // the new residual is actually smaller...
741 0 : for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step)
742 : {
743 0 : newton_stepfactor *= 0.5;
744 :
745 0 : if (!quiet)
746 0 : libMesh::out << "Shrinking step size by " << newton_stepfactor << std::endl;
747 :
748 : // Take fractional step
749 0 : solution->add(newton_stepfactor, *delta_u);
750 0 : solution->close();
751 0 : *continuation_parameter += newton_stepfactor*delta_lambda;
752 :
753 0 : rhs_mode = Residual;
754 0 : assembly(true, // Residual
755 0 : false); // Jacobian
756 0 : rhs->close();
757 0 : nonlinear_residual_afterstep = rhs->l2_norm();
758 :
759 0 : if (!quiet)
760 0 : libMesh::out << "At shrink step "
761 0 : << backtrack_step
762 0 : << ", nonlinear_residual_afterstep="
763 0 : << nonlinear_residual_afterstep
764 0 : << std::endl;
765 :
766 0 : if (nonlinear_residual_afterstep < nonlinear_residual_beforestep)
767 : {
768 0 : if (!quiet)
769 0 : libMesh::out << "Backtracking succeeded!" << std::endl;
770 :
771 0 : break; // out of backtracking loop
772 : }
773 :
774 : else
775 : {
776 : // Back off that step
777 0 : solution->add(-newton_stepfactor, *delta_u);
778 0 : solution->close();
779 0 : *continuation_parameter -= newton_stepfactor*delta_lambda;
780 : }
781 :
782 : // Save a copy of the solution from before the Newton step.
783 : //std::unique_ptr<NumericVector<Number>> prior_iterate = solution->clone();
784 : }
785 : } // end if (attempte_backtracking)
786 :
787 :
788 : // If we tried backtracking but the residual is still not reduced, print message.
789 0 : if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep))
790 : {
791 : //libMesh::err << "Backtracking failed." << std::endl;
792 0 : libMesh::out << "Backtracking failed." << std::endl;
793 :
794 : // 1.) Quit, exit program.
795 : //libmesh_error_msg("Backtracking failed!");
796 :
797 : // 2.) Continue with last newton_stepfactor
798 0 : if (newton_step<3)
799 : {
800 0 : solution->add(newton_stepfactor, *delta_u);
801 0 : solution->close();
802 0 : *continuation_parameter += newton_stepfactor*delta_lambda;
803 0 : if (!quiet)
804 0 : libMesh::out << "Backtracking could not reduce residual ... continuing anyway!" << std::endl;
805 : }
806 :
807 : // 3.) Break out of Newton iteration loop with newton_converged = false,
808 : // reduce the arclength stepsize, and try again.
809 : else
810 : {
811 0 : break; // out of Newton iteration loop, with newton_converged=false
812 : }
813 : }
814 :
815 : // Another type of convergence check: suppose the residual has not been reduced
816 : // from its initial value after half of the allowed Newton steps have occurred.
817 : // In our experience, this typically means that it isn't going to converge and
818 : // we could probably save time by dropping out of the Newton iteration loop and
819 : // trying a smaller arcstep.
820 0 : if (this->newton_progress_check)
821 : {
822 0 : if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) &&
823 0 : (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations)))
824 : {
825 0 : libMesh::out << "Progress check failed: the current residual: "
826 0 : << nonlinear_residual_afterstep
827 0 : << ", is\n"
828 0 : << "larger than the initial residual, and half of the allowed\n"
829 0 : << "number of Newton iterations have elapsed.\n"
830 0 : << "Exiting Newton iterations with converged==false." << std::endl;
831 :
832 0 : break; // out of Newton iteration loop, newton_converged = false
833 : }
834 : }
835 :
836 : // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value
837 0 : if (*continuation_parameter < min_continuation_parameter)
838 : {
839 0 : libMesh::out << "Continuation parameter fell below min-allowable value." << std::endl;
840 0 : break; // out of Newton iteration loop, newton_converged = false
841 : }
842 :
843 : // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value
844 0 : if ( (max_continuation_parameter != 0.0) &&
845 0 : (*continuation_parameter > max_continuation_parameter) )
846 : {
847 0 : libMesh::out << "Current continuation parameter value: "
848 0 : << *continuation_parameter
849 0 : << " exceeded max-allowable value."
850 0 : << std::endl;
851 0 : break; // out of Newton iteration loop, newton_converged = false
852 : }
853 :
854 :
855 : // Check the convergence of the parameter and the solution. If they are small
856 : // enough, we can break out of the Newton iteration loop.
857 0 : const Real norm_delta_u = delta_u->l2_norm();
858 0 : const Real norm_u = solution->l2_norm();
859 0 : libMesh::out << " delta_lambda = " << delta_lambda << std::endl;
860 0 : libMesh::out << " newton_stepfactor*delta_lambda = " << newton_stepfactor*delta_lambda << std::endl;
861 0 : libMesh::out << " lambda_current = " << *continuation_parameter << std::endl;
862 0 : libMesh::out << " ||delta_u|| = " << norm_delta_u << std::endl;
863 0 : libMesh::out << " ||delta_u||/||u|| = " << norm_delta_u / norm_u << std::endl;
864 :
865 :
866 : // Evaluate the residual at the current Newton iterate. We don't want to detect
867 : // convergence due to a small Newton step when the residual is still not small.
868 0 : rhs_mode = Residual;
869 0 : assembly(true, // Residual
870 0 : false); // Jacobian
871 0 : rhs->close();
872 0 : const Real norm_residual = rhs->l2_norm();
873 0 : libMesh::out << " ||R||_{L2} = " << norm_residual << std::endl;
874 0 : libMesh::out << " ||R||_{L2}/||u|| = " << norm_residual / norm_u << std::endl;
875 :
876 :
877 : // FIXME: The norm_delta_u tolerance (at least) should be relative.
878 : // It doesn't make sense to converge a solution whose size is ~ 10^5 to
879 : // a tolerance of 1.e-6. Oh, and we should also probably check the
880 : // (relative) size of the residual as well, instead of just the step.
881 0 : if ((std::abs(delta_lambda) < continuation_parameter_tolerance) &&
882 : //(norm_delta_u < solution_tolerance) && // This is a *very* strict criterion we can probably skip
883 0 : (norm_residual < solution_tolerance))
884 : {
885 0 : if (!quiet)
886 0 : libMesh::out << "Newton iterations converged!" << std::endl;
887 :
888 0 : newton_converged = true;
889 0 : break; // out of Newton iterations
890 : }
891 : } // end nonlinear loop
892 :
893 0 : if (!newton_converged)
894 : {
895 0 : libMesh::out << "Newton iterations of augmented system did not converge!" << std::endl;
896 :
897 : // Reduce ds_current, recompute the solution and parameter, and continue to next
898 : // arcstep, if there is one.
899 0 : ds_current *= 0.5;
900 :
901 : // Go back to previous solution and parameter value.
902 0 : *solution = *previous_u;
903 0 : *continuation_parameter = old_continuation_parameter;
904 :
905 : // Compute new predictor with smaller ds
906 0 : apply_predictor();
907 : }
908 : else
909 : {
910 : // Set step convergence and break out
911 0 : arcstep_converged=true;
912 0 : break; // out of arclength reduction loop
913 : }
914 :
915 : } // end loop over arclength reductions
916 :
917 : // Check for convergence of the whole arcstep. If not converged at this
918 : // point, we have no choice but to quit.
919 0 : libmesh_error_msg_if(!arcstep_converged, "Arcstep failed to converge after max number of reductions! Exiting...");
920 :
921 : // Print converged solution control parameter and max value.
922 0 : libMesh::out << "lambda_current=" << *continuation_parameter << std::endl;
923 : //libMesh::out << "u_max=" << solution->max() << std::endl;
924 :
925 : // Reset old stream precision and flags.
926 0 : libMesh::out.precision(old_precision);
927 0 : libMesh::out.unsetf(std::ios_base::scientific);
928 :
929 : // Note: we don't want to go on to the next guess yet, since the user may
930 : // want to post-process this data. It's up to the user to call advance_arcstep()
931 : // when they are ready to go on.
932 0 : }
933 :
934 :
935 :
936 0 : void ContinuationSystem::advance_arcstep()
937 : {
938 : // Solve for the updated tangent du1/ds, d(lambda1)/ds
939 0 : solve_tangent();
940 :
941 : // Advance the solution and the parameter to the next value.
942 0 : update_solution();
943 0 : }
944 :
945 :
946 :
947 : // This function solves the tangent system:
948 : // [ G_u G_{lambda} ][(du/ds)_new ] = [ 0 ]
949 : // [ Theta*(du/ds)_old (dlambda/ds)_old ][(dlambda/ds)_new ] [-N_s]
950 : // The solution is given by:
951 : // .) Let G_u y = G_lambda, then
952 : // .) 2nd row yields:
953 : // (dlambda/ds)_new = 1.0 / ( (dlambda/ds)_old - Theta*(du/ds)_old*y )
954 : // .) 1st row yields
955 : // (du_ds)_new = -(dlambda/ds)_new * y
956 0 : void ContinuationSystem::solve_tangent()
957 : {
958 : // We shouldn't call this unless the current tangent already makes sense.
959 0 : libmesh_assert (tangent_initialized);
960 :
961 : // Set pointer to underlying Newton solver
962 0 : if (!newton_solver)
963 0 : newton_solver =
964 0 : cast_ptr<NewtonSolver *> (this->time_solver->diff_solver().get());
965 :
966 : // Assemble the system matrix AND rhs, with rhs = G_{\lambda}
967 0 : this->rhs_mode = G_Lambda;
968 :
969 : // Assemble Residual and Jacobian
970 0 : this->assembly(true, // Residual
971 0 : true); // Jacobian
972 :
973 : // Not sure if this is really necessary
974 0 : rhs->close();
975 :
976 : // Solve G_u*y = G_{\lambda}
977 : std::pair<unsigned int, Real> rval =
978 0 : linear_solver->solve(*matrix,
979 0 : *y,
980 0 : *rhs,
981 0 : 1.e-12, // relative linear tolerance
982 0 : 2*newton_solver->max_linear_iterations); // max linear iterations
983 :
984 : // FIXME: If this doesn't converge at all, the new tangent vector is
985 : // going to be really bad...
986 :
987 0 : if (!quiet)
988 0 : libMesh::out << "G_u*y = G_{lambda} solver converged at step "
989 0 : << rval.first
990 0 : << " linear tolerance = "
991 0 : << rval.second
992 0 : << "."
993 0 : << std::endl;
994 :
995 : // Save old solution and parameter tangents for possible use in higher-order
996 : // predictor schemes.
997 0 : previous_dlambda_ds = dlambda_ds;
998 0 : *previous_du_ds = *du_ds;
999 :
1000 :
1001 : // 1.) Previous, probably wrong, technique!
1002 : // // Solve for the updated d(lambda)/ds
1003 : // // denom = N_{lambda} - (du_ds)^t y
1004 : // // = d(lambda)/ds - (du_ds)^t y
1005 : // Real denom = dlambda_ds - du_ds->dot(*y);
1006 :
1007 : // //libMesh::out << "denom=" << denom << std::endl;
1008 : // libmesh_assert_not_equal_to (denom, 0.0);
1009 :
1010 : // dlambda_ds = 1.0 / denom;
1011 :
1012 :
1013 : // if (!quiet)
1014 : // libMesh::out << "dlambda_ds=" << dlambda_ds << std::endl;
1015 :
1016 : // // Compute the updated value of du/ds = -_dlambda_ds * y
1017 : // du_ds->zero();
1018 : // du_ds->add(-dlambda_ds, *y);
1019 : // du_ds->close();
1020 :
1021 :
1022 : // 2.) From Brian Carnes' paper...
1023 : // According to Carnes, y comes from solving G_u * y = -G_{\lambda}
1024 0 : y->scale(-1.);
1025 0 : const Real ynorm = y->l2_norm();
1026 0 : dlambda_ds = 1. / std::sqrt(1. + Theta_LOCA*Theta_LOCA*Theta*ynorm*ynorm);
1027 :
1028 : // Determine the correct sign for dlambda_ds.
1029 :
1030 : // We will use delta_u to temporarily compute this sign.
1031 0 : *delta_u = *solution;
1032 0 : delta_u->add(-1., *previous_u);
1033 0 : delta_u->close();
1034 :
1035 : const Real sgn_dlambda_ds =
1036 0 : libmesh_real(Theta_LOCA*Theta_LOCA*Theta*y->dot(*delta_u) +
1037 0 : (*continuation_parameter-old_continuation_parameter));
1038 :
1039 0 : if (sgn_dlambda_ds < 0.)
1040 : {
1041 0 : if (!quiet)
1042 0 : libMesh::out << "dlambda_ds is negative." << std::endl;
1043 :
1044 0 : dlambda_ds *= -1.;
1045 : }
1046 :
1047 : // Finally, set the new tangent vector, du/ds = dlambda/ds * y.
1048 0 : du_ds->zero();
1049 0 : du_ds->add(dlambda_ds, *y);
1050 0 : du_ds->close();
1051 :
1052 0 : if (!quiet)
1053 : {
1054 0 : libMesh::out << "d(lambda)/ds = " << dlambda_ds << std::endl;
1055 0 : libMesh::out << "||du_ds|| = " << du_ds->l2_norm() << std::endl;
1056 : }
1057 :
1058 : // Our next solve expects y ~ -du/dlambda, so scale it back by -1 again now.
1059 0 : y->scale(-1.);
1060 0 : y->close();
1061 0 : }
1062 :
1063 :
1064 :
1065 0 : void ContinuationSystem::set_Theta()
1066 : {
1067 : // // Use the norm of the latest solution, squared.
1068 : //const Real normu = solution->l2_norm();
1069 : //libmesh_assert_not_equal_to (normu, 0.0);
1070 : //Theta = 1./normu/normu;
1071 :
1072 : // // 1.) Use the norm of du, squared
1073 : // *delta_u = *solution;
1074 : // delta_u->add(-1, *previous_u);
1075 : // delta_u->close();
1076 : // const Real normdu = delta_u->l2_norm();
1077 :
1078 : // if (normdu < 1.) // don't divide by zero or make a huge scaling parameter.
1079 : // Theta = 1.;
1080 : // else
1081 : // Theta = 1./normdu/normdu;
1082 :
1083 : // 2.) Use 1.0, i.e. don't scale
1084 0 : Theta=1.;
1085 :
1086 : // 3.) Use a formula which attempts to make the "solution triangle" isosceles.
1087 : // libmesh_assert_less (std::abs(dlambda_ds), 1.);
1088 :
1089 : // *delta_u = *solution;
1090 : // delta_u->add(-1, *previous_u);
1091 : // delta_u->close();
1092 : // const Real normdu = delta_u->l2_norm();
1093 :
1094 : // Theta = std::sqrt(1. - dlambda_ds*dlambda_ds) / normdu * tau * ds;
1095 :
1096 :
1097 : // // 4.) Use the norm of du and the norm of du/ds
1098 : // *delta_u = *solution;
1099 : // delta_u->add(-1, *previous_u);
1100 : // delta_u->close();
1101 : // const Real normdu = delta_u->l2_norm();
1102 : // du_ds->close();
1103 : // const Real normduds = du_ds->l2_norm();
1104 :
1105 : // if (normduds < 1.e-12)
1106 : // {
1107 : // libMesh::out << "Setting initial Theta= 1./normdu/normdu" << std::endl;
1108 : // libMesh::out << "normdu=" << normdu << std::endl;
1109 :
1110 : // // Don't use this scaling if the solution delta is already O(1)
1111 : // if (normdu > 1.)
1112 : // Theta = 1./normdu/normdu;
1113 : // else
1114 : // Theta = 1.;
1115 : // }
1116 : // else
1117 : // {
1118 : // libMesh::out << "Setting Theta= 1./normdu/normduds" << std::endl;
1119 : // libMesh::out << "normdu=" << normdu << std::endl;
1120 : // libMesh::out << "normduds=" << normduds << std::endl;
1121 :
1122 : // // Don't use this scaling if the solution delta is already O(1)
1123 : // if ((normdu>1.) || (normduds>1.))
1124 : // Theta = 1./normdu/normduds;
1125 : // else
1126 : // Theta = 1.;
1127 : // }
1128 :
1129 0 : if (!quiet)
1130 0 : libMesh::out << "Setting Normalization Parameter Theta=" << Theta << std::endl;
1131 0 : }
1132 :
1133 :
1134 :
1135 0 : void ContinuationSystem::set_Theta_LOCA()
1136 : {
1137 : // We also recompute the LOCA normalization parameter based on the
1138 : // most recently computed value of dlambda_ds
1139 : // if (!quiet)
1140 : // libMesh::out << "(Theta_LOCA) dlambda_ds=" << dlambda_ds << std::endl;
1141 :
1142 : // Formula makes no sense if |dlambda_ds| > 1
1143 0 : libmesh_assert_less (std::abs(dlambda_ds), 1.);
1144 :
1145 : // 1.) Attempt to implement the method in LOCA paper
1146 : // const Real g = 1./std::sqrt(2.); // "desired" dlambda_ds
1147 :
1148 : // // According to the LOCA people, we only renormalize for
1149 : // // when |dlambda_ds| exceeds some pre-selected maximum (which they take to be zero, btw).
1150 : // if (std::abs(dlambda_ds) > .9)
1151 : // {
1152 : // // Note the *= ... This is updating the previous value of Theta_LOCA
1153 : // // Note: The LOCA people actually use Theta_LOCA^2 to normalize their arclength constraint.
1154 : // Theta_LOCA *= std::abs( (dlambda_ds/g)*std::sqrt( (1.-g*g) / (1.-dlambda_ds*dlambda_ds) ) );
1155 :
1156 : // // Suggested max-allowable value for Theta_LOCA
1157 : // if (Theta_LOCA > 1.e8)
1158 : // {
1159 : // Theta_LOCA = 1.e8;
1160 :
1161 : // if (!quiet)
1162 : // libMesh::out << "max Theta_LOCA=" << Theta_LOCA << " has been selected." << std::endl;
1163 : // }
1164 : // }
1165 : // else
1166 : // Theta_LOCA=1.0;
1167 :
1168 : // 2.) FIXME: Should we do *= or just =? This function is of dlambda_ds is
1169 : // < 1, |dlambda_ds| < 1/sqrt(2) ~~ .7071
1170 : // > 1, |dlambda_ds| > 1/sqrt(2) ~~ .7071
1171 0 : Theta_LOCA *= std::abs( dlambda_ds / std::sqrt( (1.-dlambda_ds*dlambda_ds) ) );
1172 :
1173 : // Suggested max-allowable value for Theta_LOCA. I've never come close
1174 : // to this value in my code.
1175 0 : if (Theta_LOCA > 1.e8)
1176 : {
1177 0 : Theta_LOCA = 1.e8;
1178 :
1179 0 : if (!quiet)
1180 0 : libMesh::out << "max Theta_LOCA=" << Theta_LOCA << " has been selected." << std::endl;
1181 : }
1182 :
1183 : // 3.) Use 1.0, i.e. don't scale
1184 : //Theta_LOCA=1.0;
1185 :
1186 0 : if (!quiet)
1187 0 : libMesh::out << "Setting Theta_LOCA=" << Theta_LOCA << std::endl;
1188 0 : }
1189 :
1190 :
1191 :
1192 0 : void ContinuationSystem::update_solution()
1193 : {
1194 : // Set some stream formatting flags
1195 0 : std::streamsize old_precision = libMesh::out.precision();
1196 0 : libMesh::out.precision(16);
1197 0 : libMesh::out.setf(std::ios_base::scientific);
1198 :
1199 : // We must have a tangent that makes sense before we can update the solution.
1200 0 : libmesh_assert (tangent_initialized);
1201 :
1202 : // Compute tau, the stepsize scaling parameter which attempts to
1203 : // reduce ds when the angle between the most recent two tangent
1204 : // vectors becomes large. tau is actually the (absolute value of
1205 : // the) cosine of the angle between these two vectors... so if tau ~
1206 : // 0 the angle is ~ 90 degrees, while if tau ~ 1 the angle is ~ 0
1207 : // degrees.
1208 0 : y_old->close();
1209 0 : y->close();
1210 0 : const Real yoldnorm = y_old->l2_norm();
1211 0 : const Real ynorm = y->l2_norm();
1212 0 : const Number yoldy = y_old->dot(*y);
1213 0 : const Real yold_over_y = yoldnorm/ynorm;
1214 :
1215 0 : if (!quiet)
1216 : {
1217 0 : libMesh::out << "yoldnorm=" << yoldnorm << std::endl;
1218 0 : libMesh::out << "ynorm=" << ynorm << std::endl;
1219 0 : libMesh::out << "yoldy=" << yoldy << std::endl;
1220 0 : libMesh::out << "yoldnorm/ynorm=" << yoldnorm/ynorm << std::endl;
1221 : }
1222 :
1223 : // Save the current value of ds before updating it
1224 0 : previous_ds = ds_current;
1225 :
1226 : // // 1.) Cosine method (for some reason this always predicts the angle is ~0)
1227 : // // Don't try dividing by zero
1228 : // if ((yoldnorm > 1.e-12) && (ynorm > 1.e-12))
1229 : // tau = std::abs(yoldy) / yoldnorm / ynorm;
1230 : // else
1231 : // tau = 1.;
1232 :
1233 : // // 2.) Relative size of old and new du/dlambda method with cutoff of 0.9
1234 : // if ((yold_over_y < 0.9) && (yold_over_y > 1.e-6))
1235 : // tau = yold_over_y;
1236 : // else
1237 : // tau = 1.;
1238 :
1239 : // 3.) Grow (or shrink) the arclength stepsize by the ratio of du/dlambda, but do not
1240 : // exceed the user-specified value of ds.
1241 0 : if (yold_over_y > 1.e-6)
1242 : {
1243 : // // 1.) Scale current ds by the ratio of successive tangents.
1244 : // ds_current *= yold_over_y;
1245 : // if (ds_current > ds)
1246 : // ds_current = ds;
1247 :
1248 : // 2.) Technique 1 tends to shrink the step fairly well (and even if it doesn't
1249 : // get very small, we still have step reduction) but it seems to grow the step
1250 : // very slowly. Another possible technique is step-doubling:
1251 : // if (yold_over_y > 1.)
1252 : // ds_current *= 2.;
1253 : // else
1254 : // ds_current *= yold_over_y;
1255 :
1256 : // 3.) Technique 2 may over-zealous when we are also using the Newton stepgrowth
1257 : // factor. For technique 3 we multiply by yold_over_y unless yold_over_y > 2
1258 : // in which case we use 2.
1259 : // if (yold_over_y > 2.)
1260 : // ds_current *= 2.;
1261 : // else
1262 : // ds_current *= yold_over_y;
1263 :
1264 : // 4.) Double-or-halve. We double the arc-step if the ratio of successive tangents
1265 : // is larger than 'double_threshold', halve it if it is less than 'halve_threshold'
1266 0 : const Real double_threshold = 0.5;
1267 0 : const Real halve_threshold = 0.5;
1268 0 : if (yold_over_y > double_threshold)
1269 0 : ds_current *= 2.;
1270 0 : else if (yold_over_y < halve_threshold)
1271 0 : ds_current *= 0.5;
1272 :
1273 :
1274 : // Also possibly use the number of Newton iterations required to compute the previous
1275 : // step (relative to the maximum-allowed number of Newton iterations) to grow the step.
1276 0 : if (newton_stepgrowth_aggressiveness > 0.)
1277 : {
1278 0 : libmesh_assert(newton_solver);
1279 0 : const unsigned int Nmax = newton_solver->max_nonlinear_iterations;
1280 :
1281 : // // The LOCA Newton step growth technique (note: only grows step length)
1282 : // const Real stepratio = static_cast<Real>(Nmax-(newton_step+1))/static_cast<Real>(Nmax-1.);
1283 : // const Real newtonstep_growthfactor = 1. + newton_stepgrowth_aggressiveness*stepratio*stepratio;
1284 :
1285 : // The "Nopt/N" method, may grow or shrink the step. Assume Nopt=Nmax/2.
1286 0 : const Real newtonstep_growthfactor =
1287 0 : newton_stepgrowth_aggressiveness * 0.5 *
1288 0 : static_cast<Real>(Nmax) / static_cast<Real>(newton_step+1);
1289 :
1290 0 : if (!quiet)
1291 0 : libMesh::out << "newtonstep_growthfactor=" << newtonstep_growthfactor << std::endl;
1292 :
1293 0 : ds_current *= newtonstep_growthfactor;
1294 : }
1295 : }
1296 :
1297 :
1298 : // Don't let the stepsize get above the user's maximum-allowed stepsize.
1299 0 : if (ds_current > ds)
1300 0 : ds_current = ds;
1301 :
1302 : // Check also for a minimum allowed stepsize.
1303 0 : if (ds_current < ds_min)
1304 : {
1305 0 : libMesh::out << "Enforcing minimum-allowed arclength stepsize of " << ds_min << std::endl;
1306 0 : ds_current = ds_min;
1307 : }
1308 :
1309 0 : if (!quiet)
1310 : {
1311 0 : libMesh::out << "Current step size: ds_current=" << ds_current << std::endl;
1312 : }
1313 :
1314 : // Recompute scaling factor Theta for
1315 : // the current solution before updating.
1316 0 : set_Theta();
1317 :
1318 : // Also, recompute the LOCA scaling factor, which attempts to
1319 : // maintain a reasonable value of dlambda/ds
1320 0 : set_Theta_LOCA();
1321 :
1322 0 : libMesh::out << "Theta*Theta_LOCA^2=" << Theta*Theta_LOCA*Theta_LOCA << std::endl;
1323 :
1324 : // Based on the asymptotic singular behavior of du/dlambda near simple turning points,
1325 : // we can compute a single parameter which may suggest that we are close to a singularity.
1326 0 : *delta_u = *solution;
1327 0 : delta_u->add(-1, *previous_u);
1328 0 : delta_u->close();
1329 0 : const Real normdu = delta_u->l2_norm();
1330 0 : const Real C = (std::log (Theta_LOCA*normdu) /
1331 0 : std::log (std::abs(*continuation_parameter-old_continuation_parameter))) - 1.0;
1332 0 : if (!quiet)
1333 0 : libMesh::out << "C=" << C << std::endl;
1334 :
1335 : // Save the current value of u and lambda before updating.
1336 0 : save_current_solution();
1337 :
1338 0 : if (!quiet)
1339 : {
1340 0 : libMesh::out << "Updating the solution with the tangent guess." << std::endl;
1341 0 : libMesh::out << "||u_old||=" << this->solution->l2_norm() << std::endl;
1342 0 : libMesh::out << "lambda_old=" << *continuation_parameter << std::endl;
1343 : }
1344 :
1345 : // Since we solved for the tangent vector, now we can compute an
1346 : // initial guess for the new solution, and an initial guess for the
1347 : // new value of lambda.
1348 0 : apply_predictor();
1349 :
1350 0 : if (!quiet)
1351 : {
1352 0 : libMesh::out << "||u_new||=" << this->solution->l2_norm() << std::endl;
1353 0 : libMesh::out << "lambda_new=" << *continuation_parameter << std::endl;
1354 : }
1355 :
1356 : // Unset previous stream flags
1357 0 : libMesh::out.precision(old_precision);
1358 0 : libMesh::out.unsetf(std::ios_base::scientific);
1359 0 : }
1360 :
1361 :
1362 :
1363 :
1364 0 : void ContinuationSystem::save_current_solution()
1365 : {
1366 : // Save the old solution vector
1367 0 : *previous_u = *solution;
1368 :
1369 : // Save the old value of lambda
1370 0 : old_continuation_parameter = *continuation_parameter;
1371 0 : }
1372 :
1373 :
1374 :
1375 0 : void ContinuationSystem::apply_predictor()
1376 : {
1377 0 : if (predictor == Euler)
1378 : {
1379 : // 1.) Euler Predictor
1380 : // Predict next the solution
1381 0 : solution->add(ds_current, *du_ds);
1382 0 : solution->close();
1383 :
1384 : // Predict next parameter value
1385 0 : *continuation_parameter += ds_current*dlambda_ds;
1386 : }
1387 :
1388 :
1389 0 : else if (predictor == AB2)
1390 : {
1391 : // 2.) 2nd-order explicit AB predictor
1392 0 : libmesh_assert_not_equal_to (previous_ds, 0.0);
1393 0 : const Real stepratio = ds_current/previous_ds;
1394 :
1395 : // Build up next solution value.
1396 0 : solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
1397 0 : solution->add(-0.5*ds_current*stepratio , *previous_du_ds);
1398 0 : solution->close();
1399 :
1400 : // Next parameter value
1401 0 : *continuation_parameter +=
1402 0 : 0.5*ds_current*((2.+stepratio)*dlambda_ds -
1403 0 : stepratio*previous_dlambda_ds);
1404 : }
1405 :
1406 : else
1407 0 : libmesh_error_msg("Unknown predictor!");
1408 0 : }
1409 :
1410 : } // namespace libMesh
|
