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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo