libMesh
continuation_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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/sparse_matrix.h"
25 
26 namespace libMesh
27 {
28 
30  const std::string & name_in,
31  const unsigned int number_in) :
32  Parent(es, name_in, number_in),
33  continuation_parameter(libmesh_nullptr),
34  quiet(true),
35  continuation_parameter_tolerance(1.e-6),
36  solution_tolerance(1.e-6),
37  initial_newton_tolerance(0.01),
38  old_continuation_parameter(0.),
39  min_continuation_parameter(0.),
40  max_continuation_parameter(0.),
41  Theta(1.),
42  Theta_LOCA(1.),
43  n_backtrack_steps(5),
44  n_arclength_reductions(5),
45  ds_min(1.e-8),
46  predictor(Euler),
47  newton_stepgrowth_aggressiveness(1.),
48  newton_progress_check(true),
49  rhs_mode(Residual),
50  linear_solver(LinearSolver<Number>::build(es.comm())),
51  tangent_initialized(false),
52  newton_solver(libmesh_nullptr),
53  dlambda_ds(0.707),
54  ds(0.1),
55  ds_current(0.1),
56  previous_dlambda_ds(0.),
57  previous_ds(0.),
58  newton_step(0)
59 {
60  // Warn about using untested code
61  libmesh_experimental();
62 
63  if (libMesh::on_command_line("--solver_system_names"))
64  linear_solver->init((this->name()+"_").c_str());
65  else
66  linear_solver->init();
67 }
68 
69 
70 
71 
73 {
74  this->clear();
75 }
76 
77 
78 
79 
81 {
82  // FIXME: Do anything here, e.g. zero vectors, etc?
83 
84  // Call the Parent's clear function
85  Parent::clear();
86 }
87 
88 
89 
91 {
92  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
93  du_ds = &(add_vector("du_ds"));
94 
95  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
96  previous_du_ds = &(add_vector("previous_du_ds"));
97 
98  // Add a vector to keep track of the previous nonlinear solution
99  // at the old value of lambda.
100  previous_u = &(add_vector("previous_u"));
101 
102  // Add a vector to keep track of the temporary solution "y" of Ay=G_{\lambda}.
103  y = &(add_vector("y"));
104 
105  // Add a vector to keep track of the "old value" of "y" which is the solution of Ay=G_{\lambda}.
106  y_old = &(add_vector("y_old"));
107 
108  // Add a vector to keep track of the temporary solution "z" of Az=-G.
109  z = &(add_vector("z"));
110 
111  // Add a vector to keep track of the Newton update during the constrained PDE solves.
112  delta_u = &(add_vector("delta_u"));
113 
114  // Call the Parent's initialization routine.
116 }
117 
118 
119 
120 
122 {
123  // Set the Residual RHS mode, and call the normal solve routine.
124  rhs_mode = Residual;
126 }
127 
128 
129 
130 
132 {
133  // Be sure the tangent was not already initialized.
135 
136  // Compute delta_s_zero, the initial arclength traveled during the
137  // first step. Here we assume that previous_u and lambda_old store
138  // the previous solution and control parameter. You may need to
139  // read in an old solution (or solve the non-continuation system)
140  // first and call save_current_solution() before getting here.
141 
142  // 1.) Compute delta_s_zero as ||u|| - ||u_old|| + ...
143  // Compute norms of the current and previous solutions
144  // Real norm_u = solution->l2_norm();
145  // Real norm_previous_u = previous_u->l2_norm();
146 
147  // if (!quiet)
148  // {
149  // libMesh::out << "norm_u=" << norm_u << std::endl;
150  // libMesh::out << "norm_previous_u=" << norm_previous_u << std::endl;
151  // }
152 
153  // if (norm_u == norm_previous_u)
154  // {
155  // libMesh::err << "Warning, it appears u and previous_u are the "
156  // << "same, are you sure this is correct?"
157  // << "It's possible you forgot to set one or the other..."
158  // << std::endl;
159  // }
160 
161  // Real delta_s_zero = std::sqrt(
162  // (norm_u - norm_previous_u)*(norm_u - norm_previous_u) +
163  // (*continuation_parameter-old_continuation_parameter)*
164  // (*continuation_parameter-old_continuation_parameter)
165  // );
166 
167  // // 2.) Compute delta_s_zero as ||u -u_old|| + ...
168  // *delta_u = *solution;
169  // delta_u->add(-1., *previous_u);
170  // delta_u->close();
171  // Real norm_delta_u = delta_u->l2_norm();
172  // Real norm_u = solution->l2_norm();
173  // Real norm_previous_u = previous_u->l2_norm();
174 
175  // // Scale norm_delta_u by the bigger of either norm_u or norm_previous_u
176  // norm_delta_u /= std::max(norm_u, norm_previous_u);
177 
178  // if (!quiet)
179  // {
180  // libMesh::out << "norm_u=" << norm_u << std::endl;
181  // libMesh::out << "norm_previous_u=" << norm_previous_u << std::endl;
182  // //libMesh::out << "norm_delta_u=" << norm_delta_u << std::endl;
183  // libMesh::out << "norm_delta_u/max(|u|,|u_old|)=" << norm_delta_u << std::endl;
184  // libMesh::out << "|norm_u-norm_previous_u|=" << std::abs(norm_u - norm_previous_u) << std::endl;
185  // }
186 
187  // const Real dlambda = *continuation_parameter-old_continuation_parameter;
188 
189  // if (!quiet)
190  // libMesh::out << "dlambda=" << dlambda << std::endl;
191 
192  // Real delta_s_zero = std::sqrt(
193  // (norm_delta_u*norm_delta_u) +
194  // (dlambda*dlambda)
195  // );
196 
197  // if (!quiet)
198  // libMesh::out << "delta_s_zero=" << delta_s_zero << std::endl;
199 
200  // 1.) + 2.)
201  // // Now approximate the initial tangent d(lambda)/ds
202  // this->dlambda_ds = (*continuation_parameter-old_continuation_parameter) / delta_s_zero;
203 
204 
205  // // We can also approximate the deriv. wrt s by finite differences:
206  // // du/ds = (u1 - u0) / delta_s_zero.
207  // // FIXME: Use delta_u from above if we decide to keep that method.
208  // *du_ds = *solution;
209  // du_ds->add(-1., *previous_u);
210  // du_ds->scale(1./delta_s_zero);
211  // du_ds->close();
212 
213 
214  // 3.) Treating (u-previous_u)/(lambda - lambda_old) as an approximation to du/d(lambda),
215  // we follow the same technique as Carnes and Shadid.
216  // const Real dlambda = *continuation_parameter-old_continuation_parameter;
217  // libmesh_assert_greater (dlambda, 0.);
218 
219  // // Use delta_u for temporary calculation of du/d(lambda)
220  // *delta_u = *solution;
221  // delta_u->add(-1., *previous_u);
222  // delta_u->scale(1. / dlambda);
223  // delta_u->close();
224 
225  // // Determine initial normalization parameter
226  // const Real solution_size = std::max(solution->l2_norm(), previous_u->l2_norm());
227  // if (solution_size > 1.)
228  // {
229  // Theta = 1./solution_size;
230 
231  // if (!quiet)
232  // libMesh::out << "Setting Normalization Parameter Theta=" << Theta << std::endl;
233  // }
234 
235  // // Compute d(lambda)/ds
236  // // The correct sign of d(lambda)/ds should be positive, since we assume that (lambda > lambda_old)
237  // // but we could always double-check that as well.
238  // Real norm_delta_u = delta_u->l2_norm();
239  // this->dlambda_ds = 1. / std::sqrt(1. + Theta*Theta*norm_delta_u*norm_delta_u);
240 
241  // // Finally, compute du/ds = d(lambda)/ds * du/d(lambda)
242  // *du_ds = *delta_u;
243  // du_ds->scale(dlambda_ds);
244  // du_ds->close();
245 
246 
247  // 4.) Use normalized arclength formula to estimate delta_s_zero
248  // // Determine initial normalization parameter
249  // set_Theta();
250 
251  // // Compute (normalized) delta_s_zero
252  // *delta_u = *solution;
253  // delta_u->add(-1., *previous_u);
254  // delta_u->close();
255  // Real norm_delta_u = delta_u->l2_norm();
256 
257  // const Real dlambda = *continuation_parameter-old_continuation_parameter;
258 
259  // if (!quiet)
260  // libMesh::out << "dlambda=" << dlambda << std::endl;
261 
262  // Real delta_s_zero = std::sqrt(
263  // (Theta_LOCA*Theta_LOCA*Theta*norm_delta_u*norm_delta_u) +
264  // (dlambda*dlambda)
265  // );
266  // *du_ds = *delta_u;
267  // du_ds->scale(1./delta_s_zero);
268  // dlambda_ds = dlambda / delta_s_zero;
269 
270  // if (!quiet)
271  // {
272  // libMesh::out << "delta_s_zero=" << delta_s_zero << std::endl;
273  // libMesh::out << "initial d(lambda)/ds|_0 = " << dlambda_ds << std::endl;
274  // libMesh::out << "initial ||du_ds||_0 = " << du_ds->l2_norm() << std::endl;
275  // }
276 
277  // // FIXME: Also store the initial finite-differenced approximation to -du/dlambda as y.
278  // // We stick to the convention of storing negative y, since that is what we typically
279  // // solve for anyway.
280  // *y = *delta_u;
281  // y->scale(-1./dlambda);
282  // y->close();
283 
284 
285 
286  // 5.) Assume dlambda/ds_0 ~ 1/sqrt(2) and determine the value of Theta_LOCA which
287  // will satisfy this criterion
288 
289  // Initial change in parameter
291  libmesh_assert_not_equal_to (dlambda, 0.0);
292 
293  // Ideal initial value of dlambda_ds
294  dlambda_ds = 1. / std::sqrt(2.);
295  if (dlambda < 0.)
296  dlambda_ds *= -1.;
297 
298  // This also implies the initial value of ds
299  ds_current = dlambda / dlambda_ds;
300 
301  if (!quiet)
302  libMesh::out << "Setting ds_current|_0=" << ds_current << std::endl;
303 
304  // Set y = -du/dlambda using finite difference approximation
305  *y = *solution;
306  y->add(-1., *previous_u);
307  y->scale(-1./dlambda);
308  y->close();
309  const Real ynorm=y->l2_norm();
310 
311  // Finally, set the value of du_ds to be used in the upcoming
312  // tangent calculation. du/ds = du/dlambda * dlambda/ds
313  *du_ds = *y;
315  du_ds->close();
316 
317  // Determine additional solution normalization parameter
318  // (Since we just set du/ds, it will be: ||du||*||du/ds||)
319  set_Theta();
320 
321  // The value of Theta_LOCA which makes dlambda_ds = 1/sqrt(2),
322  // assuming our Theta = ||du||^2.
323  // Theta_LOCA = std::abs(dlambda);
324 
325  // Assuming general Theta
326  Theta_LOCA = std::sqrt(1./Theta/ynorm/ynorm);
327 
328 
329  if (!quiet)
330  {
331  libMesh::out << "Setting initial Theta_LOCA = " << Theta_LOCA << std::endl;
332  libMesh::out << "Theta_LOCA^2*Theta = " << Theta_LOCA*Theta_LOCA*Theta << std::endl;
333  libMesh::out << "initial d(lambda)/ds|_0 = " << dlambda_ds << std::endl;
334  libMesh::out << "initial ||du_ds||_0 = " << du_ds->l2_norm() << std::endl;
335  }
336 
337 
338 
339  // OK, we estimated the tangent at point u0.
340  // Now, to estimate the tangent at point u1, we call the solve_tangent routine.
341 
342  // Set the flag which tells us the method has been initialized.
343  tangent_initialized = true;
344 
345  solve_tangent();
346 
347  // Advance the solution and the parameter to the next value.
348  update_solution();
349 }
350 
351 
352 
353 
354 
355 
356 // This is most of the "guts" of this class. This is where we implement
357 // our custom Newton iterations and perform most of the solves.
359 {
360  // Be sure the user has set the continuation parameter pointer
362  libmesh_error_msg("You must set the continuation_parameter pointer " \
363  << "to a member variable of the derived class, preferably in the " \
364  << "Derived class's init_data function. This is how the ContinuationSystem " \
365  << "updates the continuation parameter.");
366 
367  // Use extra precision for all the numbers printed in this function.
368  std::streamsize old_precision = libMesh::out.precision();
370  libMesh::out.setf(std::ios_base::scientific);
371 
372  // We can't start solving the augmented PDE system unless the tangent
373  // vectors have been initialized. This only needs to occur once.
374  if (!tangent_initialized)
376 
377  // Save the old value of -du/dlambda. This will be used after the Newton iterations
378  // to compute the angle between previous tangent vectors. This cosine of this angle is
379  //
380  // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
381  //
382  // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
383  // when we are approaching a turning point. Note that it can only shrink the step size.
384  *y_old = *y;
385 
386  // Set pointer to underlying Newton solver
387  if (!newton_solver)
388  newton_solver = cast_ptr<NewtonSolver *> (this->time_solver->diff_solver().get());
389 
390  // A pair for catching return values from linear system solves.
391  std::pair<unsigned int, Real> rval;
392 
393  // Convergence flag for the entire arcstep
394  bool arcstep_converged = false;
395 
396  // Begin loop over arcstep reductions.
397  for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
398  {
399  if (!quiet)
400  {
401  libMesh::out << "Current arclength stepsize, ds_current=" << ds_current << std::endl;
402  libMesh::out << "Current parameter value, lambda=" << *continuation_parameter << std::endl;
403  }
404 
405  // Upon exit from the nonlinear loop, the newton_converged flag
406  // will tell us the convergence status of Newton's method.
407  bool newton_converged = false;
408 
409  // The nonlinear residual before *any* nonlinear steps have been taken.
410  Real nonlinear_residual_firststep = 0.;
411 
412  // The nonlinear residual from the current "k" Newton step, before the Newton step
413  Real nonlinear_residual_beforestep = 0.;
414 
415  // The nonlinear residual from the current "k" Newton step, after the Newton step
416  Real nonlinear_residual_afterstep = 0.;
417 
418  // The linear solver tolerance, can be updated dynamically at each Newton step.
419  Real current_linear_tolerance = 0.;
420 
421  // The nonlinear loop
423  {
424  libMesh::out << "\n === Starting Newton step " << newton_step << " ===" << std::endl;
425 
426  // Set the linear system solver tolerance
427  // // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
428  // const Real residual_multiple = 1.e-4;
429  // Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;
430 
431  // // But if the current residual isn't small, don't let the solver exit with zero iterations!
432  // if (current_linear_tolerance > 1.)
433  // current_linear_tolerance = residual_multiple;
434 
435  // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
436  if (newton_step==0)
437  {
438  // At first step, only try reducing the residual by a small amount
439  current_linear_tolerance = initial_newton_tolerance;//0.01;
440  }
441 
442  else
443  {
444  // The new tolerance is based on the ratio of the most recent tolerances
445  const Real alp=0.5*(1.+std::sqrt(5.));
446  const Real gam=0.9;
447 
448  libmesh_assert_not_equal_to (nonlinear_residual_beforestep, 0.0);
449  libmesh_assert_not_equal_to (nonlinear_residual_afterstep, 0.0);
450 
451  current_linear_tolerance = std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
452  current_linear_tolerance*current_linear_tolerance
453  );
454 
455  // Don't let it get ridiculously small!!
456  if (current_linear_tolerance < 1.e-12)
457  current_linear_tolerance = 1.e-12;
458  }
459 
460  if (!quiet)
461  libMesh::out << "Using current_linear_tolerance=" << current_linear_tolerance << std::endl;
462 
463 
464  // Assemble the residual (and Jacobian).
465  rhs_mode = Residual;
466  assembly(true, // Residual
467  true); // Jacobian
468  rhs->close();
469 
470  // Save the current nonlinear residual. We don't need to recompute the residual unless
471  // this is the first step, since it was already computed as part of the convergence check
472  // at the end of the last loop iteration.
473  if (newton_step==0)
474  {
475  nonlinear_residual_beforestep = rhs->l2_norm();
476 
477  // Store the residual before any steps have been taken. This will *not*
478  // be updated at each step, and can be used to see if any progress has
479  // been made from the initial residual at later steps.
480  nonlinear_residual_firststep = nonlinear_residual_beforestep;
481 
482  const Real old_norm_u = solution->l2_norm();
483  libMesh::out << " (before step) ||R||_{L2} = " << nonlinear_residual_beforestep << std::endl;
484  libMesh::out << " (before step) ||R||_{L2}/||u|| = " << nonlinear_residual_beforestep / old_norm_u << std::endl;
485 
486  // In rare cases (very small arcsteps), it's possible that the residual is
487  // already below our absolute linear tolerance.
488  if (nonlinear_residual_beforestep < solution_tolerance)
489  {
490  if (!quiet)
491  libMesh::out << "Initial guess satisfied linear tolerance, exiting with zero Newton iterations!" << std::endl;
492 
493  // Since we go straight from here to the solve of the next tangent, we
494  // have to close the matrix before it can be assembled again.
495  matrix->close();
496  newton_converged=true;
497  break; // out of Newton iterations, with newton_converged=true
498  }
499  }
500 
501  else
502  {
503  nonlinear_residual_beforestep = nonlinear_residual_afterstep;
504  }
505 
506 
507  // Solve the linear system G_u*z = G
508  // Initial guess?
509  z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
510  z->close();
511 
512  // It's possible that we have selected the current_linear_tolerance so large that
513  // a guess of z=zero yields a linear system residual |Az + R| small enough that the
514  // linear solver exits in zero iterations. If this happens, we will reduce the
515  // current_linear_tolerance until the linear solver does at least 1 iteration.
516  do
517  {
518  rval =
519  linear_solver->solve(*matrix,
520  *z,
521  *rhs,
522  //1.e-12,
523  current_linear_tolerance,
524  newton_solver->max_linear_iterations); // max linear iterations
525 
526  if (rval.first==0)
527  {
528  if (newton_step==0)
529  {
530  libMesh::out << "Repeating initial solve with smaller linear tolerance!" << std::endl;
531  current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
532  }
533  else
534  // We shouldn't get here ... it means the linear solver did no work on a Newton
535  // step other than the first one. If this happens, we need to think more about our
536  // tolerance selection.
537  libmesh_error_msg("Linear solver did no work!");
538  }
539 
540  } while (rval.first==0);
541 
542 
543  if (!quiet)
544  libMesh::out << " G_u*z = G solver converged at step "
545  << rval.first
546  << " linear tolerance = "
547  << rval.second
548  << "."
549  << std::endl;
550 
551  // Sometimes (I am not sure why) the linear solver exits after zero iterations.
552  // Perhaps it is hitting PETSc's divergence tolerance dtol??? If this occurs,
553  // we should break out of the Newton iteration loop because nothing further is
554  // going to happen... Of course if the tolerance is already small enough after
555  // zero iterations (how can this happen?!) we should not quit.
556  if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
557  {
558  if (!quiet)
559  libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
560 
561  // Try to find out the reason for convergence/divergence
562  linear_solver->print_converged_reason();
563 
564  break; // out of Newton iterations
565  }
566 
567  // Note: need to scale z by -1 since our code always solves Jx=R
568  // instead of Jx=-R.
569  z->scale(-1.);
570  z->close();
571 
572 
573 
574 
575 
576 
577  // Assemble the G_Lambda vector, skip residual.
578  rhs_mode = G_Lambda;
579 
580  // Assemble both rhs and Jacobian
581  assembly(true, // Residual
582  false); // Jacobian
583 
584  // Not sure if this is really necessary
585  rhs->close();
586  const Real yrhsnorm=rhs->l2_norm();
587  if (yrhsnorm == 0.0)
588  libmesh_error_msg("||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  const Real ysystemtol=current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm);
593  if (!quiet)
594  libMesh::out << "ysystemtol=" << ysystemtol << std::endl;
595 
596  // Solve G_u*y = G_{\lambda}
597  // FIXME: Initial guess? This is really a solve for -du/dlambda so we could try
598  // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda
599  //*y = *solution;
600  //y->add(-1., *previous_u);
601  //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero...
602  //y->close();
603 
604  // const unsigned int max_attempts=1;
605  // unsigned int attempt=0;
606  // do
607  // {
608  // if (!quiet)
609  // libMesh::out << "Trying to solve tangent system, attempt " << attempt << std::endl;
610 
611  rval =
612  linear_solver->solve(*matrix,
613  *y,
614  *rhs,
615  //1.e-12,
616  ysystemtol,
617  newton_solver->max_linear_iterations); // max linear iterations
618 
619  if (!quiet)
620  libMesh::out << " G_u*y = G_{lambda} solver converged at step "
621  << rval.first
622  << ", linear tolerance = "
623  << rval.second
624  << "."
625  << std::endl;
626 
627  // Sometimes (I am not sure why) the linear solver exits after zero iterations.
628  // Perhaps it is hitting PETSc's divergence tolerance dtol??? If this occurs,
629  // we should break out of the Newton iteration loop because nothing further is
630  // going to happen...
631  if ((rval.first == 0) && (rval.second > ysystemtol))
632  {
633  if (!quiet)
634  libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
635 
636  break; // out of Newton iterations
637  }
638 
639  // ++attempt;
640  // } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations));
641 
642 
643 
644 
645 
646  // Compute N, the residual of the arclength constraint eqn.
647  // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds
648  // We temporarily use the delta_u vector as a temporary vector for this calculation.
649  *delta_u = *solution;
650  delta_u->add(-1., *previous_u);
651 
652  // First part of the arclength constraint
654  const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds;
655  const Number N3 = ds_current;
656 
657  if (!quiet)
658  {
659  libMesh::out << " N1=" << N1 << std::endl;
660  libMesh::out << " N2=" << N2 << std::endl;
661  libMesh::out << " N3=" << N3 << std::endl;
662  }
663 
664  // The arclength constraint value
665  const Number N = N1+N2-N3;
666 
667  if (!quiet)
668  libMesh::out << " N=" << N << std::endl;
669 
670  const Number duds_dot_z = du_ds->dot(*z);
671  const Number duds_dot_y = du_ds->dot(*y);
672 
673  //libMesh::out << "duds_dot_z=" << duds_dot_z << std::endl;
674  //libMesh::out << "duds_dot_y=" << duds_dot_y << std::endl;
675  //libMesh::out << "dlambda_ds=" << dlambda_ds << std::endl;
676 
677  const Number delta_lambda_numerator = -(N + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z);
678  const Number delta_lambda_denominator = (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y);
679 
680  libmesh_assert_not_equal_to (delta_lambda_denominator, 0.0);
681 
682  // Now, we are ready to compute the step delta_lambda
683  const Number delta_lambda_comp = delta_lambda_numerator /
684  delta_lambda_denominator;
685  // Lambda is real-valued
686  const Real delta_lambda = libmesh_real(delta_lambda_comp);
687 
688  // Knowing delta_lambda, we are ready to update delta_u
689  // delta_u = z - delta_lambda*y
690  delta_u->zero();
691  delta_u->add(1., *z);
692  delta_u->add(-delta_lambda, *y);
693  delta_u->close();
694 
695  // Update the system solution and the continuation parameter.
696  solution->add(1., *delta_u);
697  solution->close();
698  *continuation_parameter += delta_lambda;
699 
700  // Did the Newton step actually reduce the residual?
701  rhs_mode = Residual;
702  assembly(true, // Residual
703  false); // Jacobian
704  rhs->close();
705  nonlinear_residual_afterstep = rhs->l2_norm();
706 
707 
708  // In a "normal" Newton step, ||du||/||R|| > 1 since the most recent
709  // step is where you "just were" and the current residual is where
710  // you are now. It can occur that ||du||/||R|| < 1, but these are
711  // likely not good cases to attempt backtracking (?).
712  const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep;
713  if (!quiet)
714  libMesh::out << " norm_du_norm_R=" << norm_du_norm_R << std::endl;
715 
716 
717  // Factor to decrease the stepsize by for backtracking
718  Real newton_stepfactor = 1.;
719 
720  const bool attempt_backtracking =
721  (nonlinear_residual_afterstep > solution_tolerance)
722  && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)
723  && (n_backtrack_steps>0)
724  && (norm_du_norm_R > 1.)
725  ;
726 
727  // If residual is not reduced, do Newton back tracking.
728  if (attempt_backtracking)
729  {
730  if (!quiet)
731  libMesh::out << "Newton step did not reduce residual." << std::endl;
732 
733  // back off the previous step.
734  solution->add(-1., *delta_u);
735  solution->close();
736  *continuation_parameter -= delta_lambda;
737 
738  // Backtracking: start cutting the Newton stepsize by halves until
739  // the new residual is actually smaller...
740  for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step)
741  {
742  newton_stepfactor *= 0.5;
743 
744  if (!quiet)
745  libMesh::out << "Shrinking step size by " << newton_stepfactor << std::endl;
746 
747  // Take fractional step
748  solution->add(newton_stepfactor, *delta_u);
749  solution->close();
750  *continuation_parameter += newton_stepfactor*delta_lambda;
751 
752  rhs_mode = Residual;
753  assembly(true, // Residual
754  false); // Jacobian
755  rhs->close();
756  nonlinear_residual_afterstep = rhs->l2_norm();
757 
758  if (!quiet)
759  libMesh::out << "At shrink step "
760  << backtrack_step
761  << ", nonlinear_residual_afterstep="
762  << nonlinear_residual_afterstep
763  << std::endl;
764 
765  if (nonlinear_residual_afterstep < nonlinear_residual_beforestep)
766  {
767  if (!quiet)
768  libMesh::out << "Backtracking succeeded!" << std::endl;
769 
770  break; // out of backtracking loop
771  }
772 
773  else
774  {
775  // Back off that step
776  solution->add(-newton_stepfactor, *delta_u);
777  solution->close();
778  *continuation_parameter -= newton_stepfactor*delta_lambda;
779  }
780 
781  // Save a copy of the solution from before the Newton step.
782  //UniquePtr<NumericVector<Number>> prior_iterate = solution->clone();
783  }
784  } // end if (attempte_backtracking)
785 
786 
787  // If we tried backtracking but the residual is still not reduced, print message.
788  if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep))
789  {
790  //libMesh::err << "Backtracking failed." << std::endl;
791  libMesh::out << "Backtracking failed." << std::endl;
792 
793  // 1.) Quit, exit program.
794  //libmesh_error_msg("Backtracking failed!");
795 
796  // 2.) Continue with last newton_stepfactor
797  if (newton_step<3)
798  {
799  solution->add(newton_stepfactor, *delta_u);
800  solution->close();
801  *continuation_parameter += newton_stepfactor*delta_lambda;
802  if (!quiet)
803  libMesh::out << "Backtracking could not reduce residual ... continuing anyway!" << std::endl;
804  }
805 
806  // 3.) Break out of Newton iteration loop with newton_converged = false,
807  // reduce the arclength stepsize, and try again.
808  else
809  {
810  break; // out of Newton iteration loop, with newton_converged=false
811  }
812  }
813 
814  // Another type of convergence check: suppose the residual has not been reduced
815  // from its initial value after half of the allowed Newton steps have occurred.
816  // In our experience, this typically means that it isn't going to converge and
817  // we could probably save time by dropping out of the Newton iteration loop and
818  // trying a smaller arcstep.
819  if (this->newton_progress_check)
820  {
821  if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) &&
822  (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations)))
823  {
824  libMesh::out << "Progress check failed: the current residual: "
825  << nonlinear_residual_afterstep
826  << ", is\n"
827  << "larger than the initial residual, and half of the allowed\n"
828  << "number of Newton iterations have elapsed.\n"
829  << "Exiting Newton iterations with converged==false." << std::endl;
830 
831  break; // out of Newton iteration loop, newton_converged = false
832  }
833  }
834 
835  // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value
837  {
838  libMesh::out << "Continuation parameter fell below min-allowable value." << std::endl;
839  break; // out of Newton iteration loop, newton_converged = false
840  }
841 
842  // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value
843  if ( (max_continuation_parameter != 0.0) &&
845  {
846  libMesh::out << "Current continuation parameter value: "
848  << " exceeded max-allowable value."
849  << std::endl;
850  break; // out of Newton iteration loop, newton_converged = false
851  }
852 
853 
854  // Check the convergence of the parameter and the solution. If they are small
855  // enough, we can break out of the Newton iteration loop.
856  const Real norm_delta_u = delta_u->l2_norm();
857  const Real norm_u = solution->l2_norm();
858  libMesh::out << " delta_lambda = " << delta_lambda << std::endl;
859  libMesh::out << " newton_stepfactor*delta_lambda = " << newton_stepfactor*delta_lambda << std::endl;
860  libMesh::out << " lambda_current = " << *continuation_parameter << std::endl;
861  libMesh::out << " ||delta_u|| = " << norm_delta_u << std::endl;
862  libMesh::out << " ||delta_u||/||u|| = " << norm_delta_u / norm_u << std::endl;
863 
864 
865  // Evaluate the residual at the current Newton iterate. We don't want to detect
866  // convergence due to a small Newton step when the residual is still not small.
867  rhs_mode = Residual;
868  assembly(true, // Residual
869  false); // Jacobian
870  rhs->close();
871  const Real norm_residual = rhs->l2_norm();
872  libMesh::out << " ||R||_{L2} = " << norm_residual << std::endl;
873  libMesh::out << " ||R||_{L2}/||u|| = " << norm_residual / norm_u << std::endl;
874 
875 
876  // FIXME: The norm_delta_u tolerance (at least) should be relative.
877  // It doesn't make sense to converge a solution whose size is ~ 10^5 to
878  // a tolerance of 1.e-6. Oh, and we should also probably check the
879  // (relative) size of the residual as well, instead of just the step.
880  if ((std::abs(delta_lambda) < continuation_parameter_tolerance) &&
881  //(norm_delta_u < solution_tolerance) && // This is a *very* strict criterion we can probably skip
882  (norm_residual < solution_tolerance))
883  {
884  if (!quiet)
885  libMesh::out << "Newton iterations converged!" << std::endl;
886 
887  newton_converged = true;
888  break; // out of Newton iterations
889  }
890  } // end nonlinear loop
891 
892  if (!newton_converged)
893  {
894  libMesh::out << "Newton iterations of augmented system did not converge!" << std::endl;
895 
896  // Reduce ds_current, recompute the solution and parameter, and continue to next
897  // arcstep, if there is one.
898  ds_current *= 0.5;
899 
900  // Go back to previous solution and parameter value.
901  *solution = *previous_u;
903 
904  // Compute new predictor with smaller ds
905  apply_predictor();
906  }
907  else
908  {
909  // Set step convergence and break out
910  arcstep_converged=true;
911  break; // out of arclength reduction loop
912  }
913 
914  } // end loop over arclength reductions
915 
916  // Check for convergence of the whole arcstep. If not converged at this
917  // point, we have no choice but to quit.
918  if (!arcstep_converged)
919  libmesh_error_msg("Arcstep failed to converge after max number of reductions! Exiting...");
920 
921  // Print converged solution control parameter and max value.
922  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  libMesh::out.precision(old_precision);
927  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 }
933 
934 
935 
937 {
938  // Solve for the updated tangent du1/ds, d(lambda1)/ds
939  solve_tangent();
940 
941  // Advance the solution and the parameter to the next value.
942  update_solution();
943 }
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
957 {
958  // We shouldn't call this unless the current tangent already makes sense.
960 
961  // Set pointer to underlying Newton solver
962  if (!newton_solver)
963  newton_solver =
964  cast_ptr<NewtonSolver *> (this->time_solver->diff_solver().get());
965 
966  // Assemble the system matrix AND rhs, with rhs = G_{\lambda}
967  this->rhs_mode = G_Lambda;
968 
969  // Assemble Residual and Jacobian
970  this->assembly(true, // Residual
971  true); // Jacobian
972 
973  // Not sure if this is really necessary
974  rhs->close();
975 
976  // Solve G_u*y = G_{\lambda}
977  std::pair<unsigned int, Real> rval =
978  linear_solver->solve(*matrix,
979  *y,
980  *rhs,
981  1.e-12, // relative linear tolerance
982  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  if (!quiet)
988  libMesh::out << "G_u*y = G_{lambda} solver converged at step "
989  << rval.first
990  << " linear tolerance = "
991  << rval.second
992  << "."
993  << std::endl;
994 
995  // Save old solution and parameter tangents for possible use in higher-order
996  // predictor schemes.
998  *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  y->scale(-1.);
1025  const Real ynorm = y->l2_norm();
1026  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  *delta_u = *solution;
1032  delta_u->add(-1., *previous_u);
1033  delta_u->close();
1034 
1035  const Real sgn_dlambda_ds =
1038 
1039  if (sgn_dlambda_ds < 0.)
1040  {
1041  if (!quiet)
1042  libMesh::out << "dlambda_ds is negative." << std::endl;
1043 
1044  dlambda_ds *= -1.;
1045  }
1046 
1047  // Finally, set the new tangent vector, du/ds = dlambda/ds * y.
1048  du_ds->zero();
1049  du_ds->add(dlambda_ds, *y);
1050  du_ds->close();
1051 
1052  if (!quiet)
1053  {
1054  libMesh::out << "d(lambda)/ds = " << dlambda_ds << std::endl;
1055  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  y->scale(-1.);
1060  y->close();
1061 }
1062 
1063 
1064 
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  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  if (!quiet)
1130  libMesh::out << "Setting Normalization Parameter Theta=" << Theta << std::endl;
1131 }
1132 
1133 
1134 
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  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  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  if (Theta_LOCA > 1.e8)
1176  {
1177  Theta_LOCA = 1.e8;
1178 
1179  if (!quiet)
1180  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  if (!quiet)
1187  libMesh::out << "Setting Theta_LOCA=" << Theta_LOCA << std::endl;
1188 }
1189 
1190 
1191 
1193 {
1194  // Set some stream formatting flags
1195  std::streamsize old_precision = libMesh::out.precision();
1196  libMesh::out.precision(16);
1197  libMesh::out.setf(std::ios_base::scientific);
1198 
1199  // We must have a tangent that makes sense before we can update the solution.
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  y_old->close();
1209  y->close();
1210  const Real yoldnorm = y_old->l2_norm();
1211  const Real ynorm = y->l2_norm();
1212  const Number yoldy = y_old->dot(*y);
1213  const Real yold_over_y = yoldnorm/ynorm;
1214 
1215  if (!quiet)
1216  {
1217  libMesh::out << "yoldnorm=" << yoldnorm << std::endl;
1218  libMesh::out << "ynorm=" << ynorm << std::endl;
1219  libMesh::out << "yoldy=" << yoldy << std::endl;
1220  libMesh::out << "yoldnorm/ynorm=" << yoldnorm/ynorm << std::endl;
1221  }
1222 
1223  // Save the current value of ds before updating it
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  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  const Real double_threshold = 0.5;
1267  const Real halve_threshold = 0.5;
1268  if (yold_over_y > double_threshold)
1269  ds_current *= 2.;
1270  else if (yold_over_y < halve_threshold)
1271  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.
1277  {
1279  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  const Real newtonstep_growthfactor =
1288  static_cast<Real>(Nmax) / static_cast<Real>(newton_step+1);
1289 
1290  if (!quiet)
1291  libMesh::out << "newtonstep_growthfactor=" << newtonstep_growthfactor << std::endl;
1292 
1293  ds_current *= newtonstep_growthfactor;
1294  }
1295  }
1296 
1297 
1298  // Don't let the stepsize get above the user's maximum-allowed stepsize.
1299  if (ds_current > ds)
1300  ds_current = ds;
1301 
1302  // Check also for a minimum allowed stepsize.
1303  if (ds_current < ds_min)
1304  {
1305  libMesh::out << "Enforcing minimum-allowed arclength stepsize of " << ds_min << std::endl;
1306  ds_current = ds_min;
1307  }
1308 
1309  if (!quiet)
1310  {
1311  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  set_Theta();
1317 
1318  // Also, recompute the LOCA scaling factor, which attempts to
1319  // maintain a reasonable value of dlambda/ds
1320  set_Theta_LOCA();
1321 
1322  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  *delta_u = *solution;
1327  delta_u->add(-1, *previous_u);
1328  delta_u->close();
1329  const Real normdu = delta_u->l2_norm();
1330  const Real C = (std::log (Theta_LOCA*normdu) /
1332  if (!quiet)
1333  libMesh::out << "C=" << C << std::endl;
1334 
1335  // Save the current value of u and lambda before updating.
1337 
1338  if (!quiet)
1339  {
1340  libMesh::out << "Updating the solution with the tangent guess." << std::endl;
1341  libMesh::out << "||u_old||=" << this->solution->l2_norm() << std::endl;
1342  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  apply_predictor();
1349 
1350  if (!quiet)
1351  {
1352  libMesh::out << "||u_new||=" << this->solution->l2_norm() << std::endl;
1353  libMesh::out << "lambda_new=" << *continuation_parameter << std::endl;
1354  }
1355 
1356  // Unset previous stream flags
1357  libMesh::out.precision(old_precision);
1358  libMesh::out.unsetf(std::ios_base::scientific);
1359 }
1360 
1361 
1362 
1363 
1365 {
1366  // Save the old solution vector
1367  *previous_u = *solution;
1368 
1369  // Save the old value of lambda
1371 }
1372 
1373 
1374 
1376 {
1377  if (predictor == Euler)
1378  {
1379  // 1.) Euler Predictor
1380  // Predict next the solution
1381  solution->add(ds_current, *du_ds);
1382  solution->close();
1383 
1384  // Predict next parameter value
1386  }
1387 
1388 
1389  else if (predictor == AB2)
1390  {
1391  // 2.) 2nd-order explicit AB predictor
1392  libmesh_assert_not_equal_to (previous_ds, 0.0);
1393  const Real stepratio = ds_current/previous_ds;
1394 
1395  // Build up next solution value.
1396  solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
1397  solution->add(-0.5*ds_current*stepratio , *previous_du_ds);
1398  solution->close();
1399 
1400  // Next parameter value
1402  0.5*ds_current*((2.+stepratio)*dlambda_ds -
1403  stepratio*previous_dlambda_ds);
1404  }
1405 
1406  else
1407  libmesh_error_msg("Unknown predictor!");
1408 }
1409 
1410 } // namespace libMesh
T libmesh_real(T a)
double abs(double a)
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
This is the EquationSystems class.
Real Theta
Arclength normalization parameter.
Real previous_dlambda_ds
The old parameter tangent value.
Real dlambda_ds
The most recent value of the derivative of the continuation parameter with respect to s...
NumericVector< Number > * previous_u
The solution at the previous value of the continuation variable.
virtual void solve() libmesh_override
Perform a standard "solve" of the system, without doing continuation.
UniquePtr< LinearSolver< Number > > linear_solver
We maintain our own linear solver interface, for solving custom systems of equations and/or things wh...
virtual Real l2_norm() const =0
Real initial_newton_tolerance
How much to try to reduce the residual by at the first (inexact) Newton step.
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Definition: diff_solver.h:156
Real old_continuation_parameter
The system also keeps track of the old value of the continuation parameter.
NumericVector< Number > * previous_du_ds
The value of du_ds from the previous solution.
Real * continuation_parameter
The continuation parameter must be a member variable of the derived class, and the "continuation_para...
NumericVector< Number > * delta_u
Temporary vector "delta u" ...
const class libmesh_nullptr_t libmesh_nullptr
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
NumericVector< Number > * rhs
The system matrix.
Real solution_tolerance
How tightly should the Newton iterations attempt to converge ||delta_u|| Defaults to 1...
The libMesh namespace provides an interface to certain functionality in the library.
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
Definition: diff_system.C:69
This class provides a specific system class.
Definition: fem_system.h:53
const std::string & name() const
Definition: system.h:1998
void continuation_solve()
Perform a continuation solve of the system.
virtual void solve() libmesh_override
Invokes the solver associated with the system.
Definition: diff_system.C:152
virtual void zero()=0
Set all entries to zero.
void set_Theta_LOCA()
A centralized function for setting the other normalization parameter, i.e.
Real Theta_LOCA
Another normalization parameter, which is described in the LOCA manual.
libmesh_assert(j)
Real ds_min
The minimum-allowed steplength, defaults to 1.e-8.
void advance_arcstep()
Call this function after a continuation solve to compute the tangent and get the next initial guess...
unsigned int n_arclength_reductions
Number of arclength reductions to try when Newton fails to reduce the residual.
bool tangent_initialized
False until initialize_tangent() is called.
virtual void scale(const T factor)=0
Scale each element of the vector by the given factor.
This base class can be inherited from to provide interfaces to linear solvers from different packages...
Definition: linear_solver.h:58
NumericVector< Number > * du_ds
Extra work vectors used by the continuation algorithm.
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
Real max_continuation_parameter
The maximum-allowable value of the continuation parameter.
NumericVector< Number > * z
Temporary vector "z" ...
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:148
NumericVector< Number > * y_old
Temporary vector "y_old" ...
NumericVector< Number > * y
Temporary vector "y" ...
Real ds
The initial arclength stepsize, selected by the user.
void apply_predictor()
Applies the predictor to the current solution to get a guess for the next solution.
void unsetf(std::ios_base::fmtflags mask)
Clear the associated flags.
Real newton_stepgrowth_aggressiveness
A measure of how rapidly one should attempt to grow the arclength stepsize based on the number of New...
bool quiet
If quiet==false, the System prints extra information about what it is doing.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
Real previous_ds
The previous arcstep length used.
void set_Theta()
A centralized function for setting the normalization parameter theta.
ContinuationSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
Real min_continuation_parameter
The minimum-allowable value of the continuation parameter.
double pow(double a, int b)
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
Real ds_current
Value of stepsize currently in use.
NewtonSolver * newton_solver
A pointer to the underlying Newton solver used by the DiffSystem.
void initialize_tangent()
Before starting arclength continuation, we need at least 2 prior solutions (both solution and u_previ...
void save_current_solution()
Stores the current solution and continuation parameter (as "previous_u" and "old_continuation_paramet...
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool on_command_line(const std::string &arg)
Definition: libmesh.C:921
std::ios_base::fmtflags setf(std::ios_base::fmtflags fmtfl)
Set the associated flags.
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override
Prepares matrix or rhs for matrix assembly.
Definition: fem_system.C:852
virtual ~ContinuationSystem()
Destructor.
virtual T dot(const NumericVector< T > &v) const =0
bool newton_progress_check
True by default, the Newton progress check allows the Newton loop to exit if half the allowed iterati...
unsigned int n_backtrack_steps
Another scaling parameter suggested by the LOCA people.
Real continuation_parameter_tolerance
How tightly should the Newton iterations attempt to converge delta_lambda.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
void solve_tangent()
Special solve algorithm for solving the tangent system.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
void update_solution()
This function (which assumes the most recent tangent vector has been computed) updates the solution a...
long double min(long double a, double b)
std::streamsize precision() const
Get the associated write precision.
Second-order explicit Adams-Bashforth predictor.
First-order Euler predictor.
unsigned int newton_step
Loop counter for nonlinear (Newton) iteration loop.