assembly_handler.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// OOOMPH-LIB includes
27#include "assembly_handler.h"
28#include "elements.h"
29#include "problem.h"
30#include "mesh.h"
31
32namespace oomph
33{
34 /// ////////////////////////////////////////////////////////////////////
35 // Non-inline functions for the AssemblyHandler class
36 /// ///////////////////////////////////////////////////////////////////
37
38 //===================================================================
39 /// Get the number of elemental degrees of freedom. Direct call
40 /// to the function in the element.
41 //===================================================================
42 unsigned AssemblyHandler::ndof(GeneralisedElement* const& elem_pt)
43 {
44 return elem_pt->ndof();
45 }
46
47 //===================================================================
48 /// Get the vector of dofs in the element elem_pt at time level t
49 /// Direct call to the function in the element.
50 //===================================================================
52 const unsigned& t,
53 Vector<double>& dof)
54 {
55 return elem_pt->dof_vector(t, dof);
56 }
57
58 //===================================================================
59 /// Get the vector of pointers to dofs in the element elem_pt
60 /// Direct call to the function in the element.
61 //===================================================================
63 Vector<double*>& dof_pt)
64 {
65 return elem_pt->dof_pt_vector(dof_pt);
66 }
67
68 /// Return the t-th level of storage associated with the i-th
69 /// (local) dof stored in the problem
70 double& AssemblyHandler::local_problem_dof(Problem* const& problem_pt,
71 const unsigned& t,
72 const unsigned& i)
73 {
74 return *(problem_pt->dof_pt(i) + t);
75 }
76
77
78 //==================================================================
79 /// Get the global equation number of the local unknown. Direct call
80 /// to the function in the element.
81 //==================================================================
82 unsigned long AssemblyHandler::eqn_number(GeneralisedElement* const& elem_pt,
83 const unsigned& ieqn_local)
84 {
85 return elem_pt->eqn_number(ieqn_local);
86 }
87
88 //==================================================================
89 /// Get the residuals by calling the underlying element's residuals
90 /// directly.
91 //=================================================================
93 Vector<double>& residuals)
94 {
95 elem_pt->get_residuals(residuals);
96 }
97
98 //=======================================================================
99 /// Calculate the elemental Jacobian matrix "d equation / d variable" by
100 /// calling the element's get_jacobian function.
101 //======================================================================
103 Vector<double>& residuals,
104 DenseMatrix<double>& jacobian)
105 {
106 elem_pt->get_jacobian(residuals, jacobian);
107 }
108
109 //=======================================================================
110 /// Calculate all desired vectors and matrices that are required by
111 /// the element by calling the get_jacobian function
112 //=======================================================================
114 GeneralisedElement* const& elem_pt,
117 {
118 get_jacobian(elem_pt, vec[0], matrix[0]);
119 }
120
121 //=======================================================================
122 /// Calculate the derivative of the residuals with respect to
123 /// a parameter, by calling the elemental function
124 //======================================================================
126 GeneralisedElement* const& elem_pt,
127 double* const& parameter_pt,
128 Vector<double>& dres_dparam)
129 {
130 elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
131 }
132
133
134 //=====================================================================
135 /// Calculate the derivative of the residuals and jacobian
136 /// with respect to a parameter by calling the elemental function
137 //========================================================================
139 GeneralisedElement* const& elem_pt,
140 double* const& parameter_pt,
141 Vector<double>& dres_dparam,
142 DenseMatrix<double>& djac_dparam)
143 {
144 elem_pt->get_djacobian_dparameter(parameter_pt, dres_dparam, djac_dparam);
145 }
146
147 /// Calculate the product of the Hessian (derivative of Jacobian with
148 /// respect to all variables) an eigenvector, Y, and
149 /// other specified vectors, C
150 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}
151 /// At the moment the dof pointer is passed in to enable
152 /// easy calculation of finite difference default
154 GeneralisedElement* const& elem_pt,
155 Vector<double> const& Y,
156 DenseMatrix<double> const& C,
157 DenseMatrix<double>& product)
158 {
159 elem_pt->get_hessian_vector_products(Y, C, product);
160 }
161
162
163 //=======================================================================
164 /// Return the eigenfunction(s) associated with the bifurcation that
165 /// has been detected in bifurcation tracking problems. Default
166 /// Broken implementation
167 //=======================================================================
169 {
170 std::ostringstream error_stream;
171 error_stream << "There is no bifurcation parameter associated with the "
172 "current assembly handler.\n"
173 << "Eigenfunction are only calculated by the Fold, PitchFork "
174 "and Hopf Handlers"
175 << "\n";
176
177 throw OomphLibError(
178 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
179 return 0;
180 }
181
182
183 //=======================================================================
184 /// Return the eigenfunction(s) associated with the bifurcation that
185 /// has been detected in bifurcation tracking problems. Default
186 /// Broken implementation
187 //=======================================================================
189 {
190 std::ostringstream error_stream;
191 error_stream << "There is no eigenfunction associated with the current "
192 "assembly handler.\n"
193 << "Eigenfunction are only calculated by the Fold, PitchFork "
194 "and Hopf Handlers"
195 << "\n";
196
197 throw OomphLibError(
198 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
199 }
200
201 /// ==========================================================================
202 /// Compute the inner products of the given vector of pairs of
203 /// history values over the element. The values of the index in the pair
204 /// may be the same.
205 //===========================================================================
207 GeneralisedElement* const& elem_pt,
208 Vector<std::pair<unsigned, unsigned>> const& history_index,
209 Vector<double>& inner_product)
210 {
211 elem_pt->get_inner_products(history_index, inner_product);
212 }
213
214 //==========================================================================
215 /// Compute the vectors that when taken as a dot product with
216 /// other history values give the inner product over the element.
217 /// In other words if we call get_inner_product_vectors({0,1},output)
218 /// output[0] will be a vector such that dofs.output[0] gives the inner
219 /// product of current dofs with themselves.
220 //==========================================================================
222 GeneralisedElement* const& elem_pt,
223 Vector<unsigned> const& history_index,
224 Vector<Vector<double>>& inner_product_vector)
225 {
226 elem_pt->get_inner_product_vectors(history_index, inner_product_vector);
227 }
228
229
230 /// ////////////////////////////////////////////////////////////////////
231 // Non-inline functions for the ExplicitTimeStepHandler class
232 /// ///////////////////////////////////////////////////////////////////
233
234
235 //===================================================================
236 /// Get the number of elemental degrees of freedom. Direct call
237 /// to the function in the element.
238 //===================================================================
240 {
241 return elem_pt->ndof();
242 }
243
244 //==================================================================
245 /// Get the global equation number of the local unknown. Direct call
246 /// to the function in the element.
247 //==================================================================
249 GeneralisedElement* const& elem_pt, const unsigned& ieqn_local)
250 {
251 return elem_pt->eqn_number(ieqn_local);
252 }
253
254 //==================================================================
255 /// Call the element's residuals
256 //=================================================================
258 GeneralisedElement* const& elem_pt, Vector<double>& residuals)
259 {
260 elem_pt->get_residuals(residuals);
261 }
262
263 //=======================================================================
264 /// Replace get jacobian with the call to get the mass matrix
265 //======================================================================
267 Vector<double>& residuals,
268 DenseMatrix<double>& jacobian)
269 {
270 elem_pt->get_mass_matrix(residuals, jacobian);
271 }
272
273
274 //=======================================================================
275 /// Calculate all desired vectors and matrices that are required by
276 /// the problem by calling those of the underlying element.
277 //=======================================================================
279 GeneralisedElement* const& elem_pt,
282 {
283#ifdef PARANOID
284 // Check dimension
285 if (matrix.size() != 1)
286 {
287 throw OomphLibError("ExplicitTimeSteps should return one matrix",
288 OOMPH_CURRENT_FUNCTION,
289 OOMPH_EXCEPTION_LOCATION);
290 }
291#endif
292 // Get the residuals and <mass matrix
293 elem_pt->get_mass_matrix(vec[0], matrix[0]);
294 }
295
296
297 /// ////////////////////////////////////////////////////////////////////
298 // Non-inline functions for the EigenProblemHandler class
299 /// ///////////////////////////////////////////////////////////////////
300
301
302 //===================================================================
303 /// Get the number of elemental degrees of freedom. Direct call
304 /// to the function in the element.
305 //===================================================================
307 {
308 return elem_pt->ndof();
309 }
310
311 //==================================================================
312 /// Get the global equation number of the local unknown. Direct call
313 /// to the function in the element.
314 //==================================================================
316 GeneralisedElement* const& elem_pt, const unsigned& ieqn_local)
317 {
318 return elem_pt->eqn_number(ieqn_local);
319 }
320
321 //==================================================================
322 /// Cannot call get_residuals for an eigenproblem, so throw an error
323 //=================================================================
325 Vector<double>& residuals)
326 {
327 throw OomphLibError(
328 "An eigenproblem does not have a get_residuals function",
329 OOMPH_CURRENT_FUNCTION,
330 OOMPH_EXCEPTION_LOCATION);
331 }
332
333 //=======================================================================
334 /// Cannot call get_jacobian for an eigenproblem, so throw an error
335 //======================================================================
337 Vector<double>& residuals,
338 DenseMatrix<double>& jacobian)
339 {
340 throw OomphLibError("An eigenproblem does not have a get_jacobian function",
341 OOMPH_CURRENT_FUNCTION,
342 OOMPH_EXCEPTION_LOCATION);
343 }
344
345
346 //=======================================================================
347 /// Calculate all desired vectors and matrices that are required by
348 /// the problem by calling those of the underlying element.
349 //=======================================================================
351 GeneralisedElement* const& elem_pt,
354 {
355#ifdef PARANOID
356 // Check dimension
357 if (matrix.size() != 2)
358 {
359 throw OomphLibError("EigenProblems should return two matrices",
360 OOMPH_CURRENT_FUNCTION,
361 OOMPH_EXCEPTION_LOCATION);
362 }
363#endif
364 // Find the number of variables
365 unsigned n_var = elem_pt->ndof();
366 // Assign a dummy residuals vector
367 Vector<double> dummy(n_var);
368 // Get the jacobian and mass matrices
369 elem_pt->get_jacobian_and_mass_matrix(dummy, matrix[0], matrix[1]);
370
371 // If we have a non-zero shift, then shift the A matrix
372 if (Sigma_real != 0.0)
373 {
374 // Set the shifted matrix
375 for (unsigned i = 0; i < n_var; i++)
376 {
377 for (unsigned j = 0; j < n_var; j++)
378 {
379 matrix[0](i, j) -= Sigma_real * matrix[1](i, j);
380 }
381 }
382 }
383 }
384
385
386 //======================================================================
387 /// Clean up the memory that may have been allocated by the solver
388 //=====================================================================
390 {
391 if (Alpha_pt != 0)
392 {
393 delete Alpha_pt;
394 }
395 if (E_pt != 0)
396 {
397 delete E_pt;
398 }
399 }
400
401 //===================================================================
402 /// Use a block factorisation to solve the augmented system
403 /// associated with a PitchFork bifurcation.
404 //===================================================================
406 DoubleVector& result)
407 {
408 // if the result is setup then it should not be distributed
409#ifdef PARANOID
410 if (result.built())
411 {
412 if (result.distributed())
413 {
414 throw OomphLibError("The result vector must not be distributed",
415 OOMPH_CURRENT_FUNCTION,
416 OOMPH_EXCEPTION_LOCATION);
417 }
418 }
419#endif
420
421 // Locally cache the pointer to the handler.
422 FoldHandler* handler_pt =
423 static_cast<FoldHandler*>(problem_pt->assembly_handler_pt());
424
425 // Switch the handler to "block solver" mode
426 handler_pt->solve_augmented_block_system();
427
428 // We need to find out the number of dofs in the problem
429 unsigned n_dof = problem_pt->ndof();
430
431 // create the linear algebra distribution for this solver
432 // currently only global (non-distributed) distributions are allowed
433 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
434 this->build_distribution(dist);
435
436 // if the result vector is not setup then rebuild with distribution = global
437 if (!result.built())
438 {
439 result.build(this->distribution_pt(), 0.0);
440 }
441
442 // Setup storage for temporary vectors
443 DoubleVector a(this->distribution_pt(), 0.0),
444 b(this->distribution_pt(), 0.0);
445
446 // Allocate storage for Alpha which can be used in the resolve
447 if (Alpha_pt != 0)
448 {
449 delete Alpha_pt;
450 }
451 Alpha_pt = new DoubleVector(this->distribution_pt(), 0.0);
452
453 // We are going to do resolves using the underlying linear solver
455
456 // Solve the first system Aa = R
457 Linear_solver_pt->solve(problem_pt, a);
458
459 // The vector in the top right-hand block is the jacobian multiplied
460 // by the null vector
461
462 // Get the present null vector from the handler
463 DoubleVector y(this->distribution_pt(), 0.0);
464 for (unsigned n = 0; n < (n_dof - 1); ++n)
465 {
466 y[n] = handler_pt->Y[n];
467 }
468 // For simplicity later, add a zero at the end
469 y[n_dof - 1] = 0.0;
470
471 // Loop over the elements and assemble the product
472 // Local storage for the product terms
473 DoubleVector Jy(this->distribution_pt(), 0.0);
474
475 // Calculate the product of the jacobian matrices, etc
476 unsigned long n_element = problem_pt->mesh_pt()->nelement();
477 for (unsigned long e = 0; e < n_element; e++)
478 {
479 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
480 // Loop over the ndofs in each element
481 unsigned n_var = elem_pt->ndof();
482 // Get the jacobian matrix
483 DenseMatrix<double> jac(n_var);
484 // storage for the residual
485 Vector<double> res(n_var);
486 // Get unperturbed raw jacobian
487 elem_pt->get_jacobian(res, jac);
488
489 // Multiply the dofs
490 for (unsigned n = 0; n < n_var; n++)
491 {
492 unsigned eqn_number = elem_pt->eqn_number(n);
493 for (unsigned m = 0; m < n_var; m++)
494 {
495 unsigned unknown = elem_pt->eqn_number(m);
496 Jy[eqn_number] += jac(n, m) * y[unknown];
497 }
498 }
499 }
500 // The final entry of the vector will be zero
501
502 // Now resolve to find alpha
504
505 // The vector that multiplie the product matrix is actually y - alpha
506 DoubleVector y_minus_alpha(this->distribution_pt(), 0.0);
507 for (unsigned n = 0; n < n_dof; n++)
508 {
509 y_minus_alpha[n] = y[n] - (*Alpha_pt)[n];
510 }
511
512 // We can now construct our multipliers
513 // Prepare to scale
514 double dof_length = 0.0, a_length = 0.0, alpha_length = 0.0;
515 for (unsigned n = 0; n < n_dof; n++)
516 {
517 if (std::fabs(problem_pt->dof(n)) > dof_length)
518 {
519 dof_length = std::fabs(problem_pt->dof(n));
520 }
521 if (std::fabs(a[n]) > a_length)
522 {
523 a_length = std::fabs(a[n]);
524 }
525 if (std::fabs(y_minus_alpha[n]) > alpha_length)
526 {
527 alpha_length = std::fabs(y_minus_alpha[n]);
528 }
529 }
530
531 double a_mult = dof_length / a_length;
532 double alpha_mult = dof_length / alpha_length;
533 const double FD_step = 1.0e-8;
534 a_mult += FD_step;
535 alpha_mult += FD_step;
536 a_mult *= FD_step;
537 alpha_mult *= FD_step;
538
539 // Local storage for the product terms
540 DoubleVector Jprod_a(this->distribution_pt(), 0.0),
541 Jprod_alpha(this->distribution_pt(), 0.0);
542
543 // Calculate the product of the jacobian matrices, etc
544 for (unsigned long e = 0; e < n_element; e++)
545 {
546 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
547 // Loop over the ndofs in each element
548 unsigned n_var = handler_pt->ndof(elem_pt);
549 // Get the jacobian matrices
550 DenseMatrix<double> jac(n_var), jac_a(n_var), jac_alpha(n_var);
551 // elemental residual storage
552 Vector<double> res_elemental(n_var);
553 // Get unperturbed jacobian
554 handler_pt->get_jacobian(elem_pt, res_elemental, jac);
555
556 // Backup the dofs
557 Vector<double> dof_bac(n_var);
558 // Perturb the dofs
559 for (unsigned n = 0; n < n_var; n++)
560 {
561 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
562 dof_bac[n] = problem_pt->dof(eqn_number);
563 // Pertub by vector a
564 problem_pt->dof(eqn_number) += a_mult * a[eqn_number];
565 }
566
568
569 // Now get the new jacobian
570 handler_pt->get_jacobian(elem_pt, res_elemental, jac_a);
571
572 // Perturb the dofs
573 for (unsigned n = 0; n < n_var; n++)
574 {
575 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
576 problem_pt->dof(eqn_number) = dof_bac[n];
577 // Pertub by vector a
578 problem_pt->dof(eqn_number) += alpha_mult * y_minus_alpha[eqn_number];
579 }
580
582
583 // Now get the new jacobian
584 handler_pt->get_jacobian(elem_pt, res_elemental, jac_alpha);
585
586 // Reset the dofs
587 for (unsigned n = 0; n < n_var; n++)
588 {
589 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
590 problem_pt->dof(eqn_number) = dof_bac[n];
591 }
592
594
595 // OK, now work out the products
596 // Note the (n_var-1), we are only interested in the non-augmented
597 // jacobian
598 for (unsigned n = 0; n < (n_var - 1); n++)
599 {
600 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
601 double prod_a = 0.0, prod_alpha = 0.0;
602 for (unsigned m = 0; m < (n_var - 1); m++)
603 {
604 unsigned unknown = handler_pt->eqn_number(elem_pt, m);
605 prod_a += (jac_a(n, m) - jac(n, m)) * y[unknown];
606 prod_alpha += (jac_alpha(n, m) - jac(n, m)) * y[unknown];
607 }
608 Jprod_a[eqn_number] += prod_a / a_mult;
609 Jprod_alpha[eqn_number] += prod_alpha / alpha_mult;
610 }
611 }
612
613 Jprod_alpha[n_dof - 1] = 0.0;
614 Jprod_a[n_dof - 1] = 0.0;
615
616 // OK, now we can formulate the next vectors
617 // The assumption here is that the result has been set to the
618 // residuals.
619 for (unsigned n = 0; n < n_dof - 1; n++)
620 {
621 b[n] = result[n_dof + n] - Jprod_a[n];
622 }
623 // The final residual is the entry corresponding to the original parameter
624 b[n_dof - 1] = result[n_dof - 1];
625
626 // Allocate storage for E which can be used in the resolve
627 if (E_pt != 0)
628 {
629 delete E_pt;
630 }
631 E_pt = new DoubleVector(this->distribution_pt(), 0.0);
632 DoubleVector f(this->distribution_pt(), 0.0);
634 Linear_solver_pt->resolve(Jprod_alpha, *E_pt);
635
636 // Calculate the final entry in the vector e
637 const double e_final = (*E_pt)[n_dof - 1];
638 // Calculate the final entry in the vector d
639 const double d_final = f[n_dof - 1] / e_final;
640 // Assemble the final corrections
641 for (unsigned n = 0; n < n_dof - 1; n++)
642 {
643 result[n] = a[n] - (*Alpha_pt)[n] * d_final + d_final * y[n];
644 result[n_dof + n] = f[n] - (*E_pt)[n] * d_final;
645 }
646 // The result corresponding to the parameter
647 result[n_dof - 1] = a[n_dof - 1] - (*Alpha_pt)[n_dof - 1] * d_final;
648
649 // The sign of the jacobian is the sign of the final entry in e
650 problem_pt->sign_of_jacobian() =
651 static_cast<int>(std::fabs(e_final) / e_final);
652
653 // Switch things to our block solver
654 handler_pt->solve_full_system();
655
656 // If we are not storing things, clear the results
657 if (!Enable_resolve)
658 {
659 // We no longer need to store the matrix
661 delete Alpha_pt;
662 Alpha_pt = 0;
663 delete E_pt;
664 E_pt = 0;
665 }
666 // Otherwise, also store the problem pointer
667 else
668 {
669 Problem_pt = problem_pt;
670 }
671 }
672
673
674 //======================================================================
675 // Hack the re-solve to use the block factorisation
676 //======================================================================
678 DoubleVector& result)
679 {
680 // Check that the factors have been stored
681 if (Alpha_pt == 0)
682 {
683 throw OomphLibError("The required vectors have not been stored",
684 OOMPH_CURRENT_FUNCTION,
685 OOMPH_EXCEPTION_LOCATION);
686 }
687
688 // Get the pointer to the problem
689 Problem* const problem_pt = Problem_pt;
690
691 FoldHandler* handler_pt =
692 static_cast<FoldHandler*>(problem_pt->assembly_handler_pt());
693
694 // Switch things to our block solver
695 handler_pt->solve_augmented_block_system();
696
697 // We need to find out the number of dofs in the problem
698 unsigned n_dof = problem_pt->ndof();
699
700#ifdef PARANOID
701 // if the result is setup then it should not be distributed
702 if (result.built())
703 {
704 if (result.distributed())
705 {
706 throw OomphLibError("The result vector must not be distributed",
707 OOMPH_CURRENT_FUNCTION,
708 OOMPH_EXCEPTION_LOCATION);
709 }
710 }
711 // the rhs must be setup
712 if (!rhs.built())
713 {
714 throw OomphLibError("The rhs vector must be setup",
715 OOMPH_CURRENT_FUNCTION,
716 OOMPH_EXCEPTION_LOCATION);
717 }
718#endif
719
720 // if the result vector is not setup then rebuild with distribution = global
721 if (!result.built())
722 {
723 result.build(rhs.distribution_pt(), 0.0);
724 }
725
726 // Setup storage
727 DoubleVector a(this->distribution_pt(), 0.0),
728 b(this->distribution_pt(), 0.0);
729
730 // Set the values of the a vector
731 for (unsigned n = 0; n < (n_dof - 1); n++)
732 {
733 a[n] = rhs[n];
734 }
735 // The entry associated with the additional parameter is zero
736 a[n_dof - 1] = 0.0;
737
739
740 // Copy rhs vector into local storage so it doesn't get overwritten
741 // if the linear solver decides to initialise the solution vector, say,
742 // which it's quite entitled to do!
743 DoubleVector input_a(a);
744
745 Linear_solver_pt->resolve(input_a, a);
746
747 // We can now construct our multipliers
748 // Prepare to scale
749 double dof_length = 0.0, a_length = 0.0;
750 for (unsigned n = 0; n < n_dof; n++)
751 {
752 if (std::fabs(problem_pt->dof(n)) > dof_length)
753 {
754 dof_length = std::fabs(problem_pt->dof(n));
755 }
756
757 if (std::fabs(a[n]) > a_length)
758 {
759 a_length = std::fabs(a[n]);
760 }
761 }
762 double a_mult = dof_length / a_length;
763 const double FD_step = 1.0e-8;
764 a_mult += FD_step;
765 a_mult *= FD_step;
766
767 DoubleVector Jprod_a(this->distribution_pt(), 0.0);
768
769 unsigned long n_element = problem_pt->mesh_pt()->nelement();
770 for (unsigned long e = 0; e < n_element; e++)
771 {
772 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
773 // Loop over the ndofs in each element
774 unsigned n_var = handler_pt->ndof(elem_pt);
775 // Get some jacobian matrices
776 DenseMatrix<double> jac(n_var), jac_a(n_var);
777 // the elemental residual
778 Vector<double> res_elemental(n_var);
779 // Get unperturbed jacobian
780 handler_pt->get_jacobian(elem_pt, res_elemental, jac);
781
782 // Backup the dofs
783 Vector<double> dof_bac(n_var);
784 // Perturb the dofs
785 for (unsigned n = 0; n < n_var; n++)
786 {
787 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
788 dof_bac[n] = problem_pt->dof(eqn_number);
789 // Pertub by vector a
790 problem_pt->dof(eqn_number) += a_mult * a[eqn_number];
791 }
792
794
795 // Now get the new jacobian
796 handler_pt->get_jacobian(elem_pt, res_elemental, jac_a);
797
798 // Reset the dofs
799 for (unsigned n = 0; n < n_var; n++)
800 {
801 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
802 problem_pt->dof(eqn_number) = dof_bac[n];
803 }
804
806
807 // OK, now work out the products
808 for (unsigned n = 0; n < (n_var - 1); n++)
809 {
810 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
811 double prod_a = 0.0;
812 for (unsigned m = 0; m < (n_var - 1); m++)
813 {
814 unsigned unknown = handler_pt->eqn_number(elem_pt, m);
815 prod_a += (jac_a(n, m) - jac(n, m)) * handler_pt->Y[unknown];
816 }
817 Jprod_a[eqn_number] += prod_a / a_mult;
818 }
819 }
820
821 Jprod_a[n_dof - 1] = 0.0;
822
823 // OK, now we can formulate the next vectors
824 for (unsigned n = 0; n < n_dof - 1; n++)
825 {
826 b[n] = rhs[n_dof + n] - Jprod_a[n];
827 }
828 // The residuals for the final entry should be those associated
829 // with the parameter
830 b[n_dof - 1] = rhs[n_dof - 1];
831
832 DoubleVector f(this->distribution_pt(), 0.0);
833
835
836 // Calculate the final entry in the vector d
837 const double d_final = f[n_dof - 1] / (*E_pt)[n_dof - 1];
838 // Assemble the final corrections
839 for (unsigned n = 0; n < n_dof - 1; n++)
840 {
841 result[n] = a[n] - (*Alpha_pt)[n] * d_final + d_final * handler_pt->Y[n];
842 result[n_dof + n] = f[n] - (*E_pt)[n] * d_final;
843 }
844 // The result corresponding to the paramater
845 result[n_dof - 1] = a[n_dof - 1] - (*Alpha_pt)[n_dof - 1] * d_final;
846
848
849 // Switch things to our block solver
850 handler_pt->solve_full_system();
851 }
852
853
854 /// ////////////////////////////////////////////////////////////////////
855 // Non-inline functions for the FoldHandler class
856 /// ///////////////////////////////////////////////////////////////////
857
858
859 //========================================================================
860 /// Constructor: Initialise the fold handler,
861 /// by setting initial guesses for Y, Phi and calculating Count.
862 /// If the system changes, a new handler must be constructed.
863 //========================================================================
865 double* const& parameter_pt)
866 : Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
867 {
868 // Set the problem pointer
869 Problem_pt = problem_pt;
870 // Set the number of degrees of freedom
871 Ndof = problem_pt->ndof();
872
873 // create the linear algebra distribution for this solver
874 // currently only global (non-distributed) distributions are allowed
876 new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
877
878 // Resize the vectors of additional dofs and constants
879 Phi.resize(Ndof);
880 Y.resize(Ndof);
881 Count.resize(Ndof, 0);
882
883 // Loop over all the elements in the problem
884 unsigned n_element = problem_pt->mesh_pt()->nelement();
885 for (unsigned e = 0; e < n_element; e++)
886 {
887 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
888 // Loop over the local freedoms in the element
889 unsigned n_var = elem_pt->ndof();
890 for (unsigned n = 0; n < n_var; n++)
891 {
892 // Increase the associated global equation number counter
893 ++Count[elem_pt->eqn_number(n)];
894 }
895 }
896
897 // Calculate the value Phi by
898 // solving the system JPhi = dF/dlambda
899
900 // Locally cache the linear solver
901 LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();
902
903 // Save the status before entry to this routine
904 bool enable_resolve = linear_solver_pt->is_resolve_enabled();
905
906 // We need to do a resolve
907 linear_solver_pt->enable_resolve();
908
909 // Storage for the solution
910 DoubleVector x(dist_pt, 0.0);
911
912 // Solve the standard problem, we only want to make sure that
913 // we factorise the matrix, if it has not been factorised. We shall
914 // ignore the return value of x.
915 linear_solver_pt->solve(problem_pt, x);
916
917 // Get the vector dresiduals/dparameter
918 problem_pt->get_derivative_wrt_global_parameter(parameter_pt, x);
919
920 // Copy rhs vector into local storage so it doesn't get overwritten
921 // if the linear solver decides to initialise the solution vector, say,
922 // which it's quite entitled to do!
923 DoubleVector input_x(x);
924
925 // Now resolve the system with the new RHS and overwrite the solution
926 linear_solver_pt->resolve(input_x, x);
927
928 // Restore the storage status of the linear solver
929 if (enable_resolve)
930 {
931 linear_solver_pt->enable_resolve();
932 }
933 else
934 {
935 linear_solver_pt->disable_resolve();
936 }
937
938 // Add the global parameter as an unknown to the problem
939 problem_pt->Dof_pt.push_back(parameter_pt);
940
941
942 // Normalise the initial guesses for phi
943 double length = 0.0;
944 for (unsigned n = 0; n < Ndof; n++)
945 {
946 length += x[n] * x[n];
947 }
948 length = sqrt(length);
949
950 // Now add the null space components to the problem unknowns
951 // and initialise them and Phi to the same normalised values
952 for (unsigned n = 0; n < Ndof; n++)
953 {
954 problem_pt->Dof_pt.push_back(&Y[n]);
955 Y[n] = Phi[n] = -x[n] / length;
956 }
957
958 // delete the dist_pt
959 problem_pt->Dof_distribution_pt->build(
960 problem_pt->communicator_pt(), Ndof * 2 + 1, true);
961 // Remove all previous sparse storage used during Jacobian assembly
963
964 delete dist_pt;
965 }
966
967
968 /// Constructor in which the eigenvector is passed as an initial
969 /// guess
971 double* const& parameter_pt,
972 const DoubleVector& eigenvector)
973 : Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
974 {
975 // Set the problem pointer
976 Problem_pt = problem_pt;
977 // Set the number of degrees of freedom
978 Ndof = problem_pt->ndof();
979
980 // create the linear algebra distribution for this solver
981 // currently only global (non-distributed) distributions are allowed
983 new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
984
985 // Resize the vectors of additional dofs and constants
986 Phi.resize(Ndof);
987 Y.resize(Ndof);
988 Count.resize(Ndof, 0);
989
990 // Loop over all the elements in the problem
991 unsigned n_element = problem_pt->mesh_pt()->nelement();
992 for (unsigned e = 0; e < n_element; e++)
993 {
994 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
995 // Loop over the local freedoms in the element
996 unsigned n_var = elem_pt->ndof();
997 for (unsigned n = 0; n < n_var; n++)
998 {
999 // Increase the associated global equation number counter
1000 ++Count[elem_pt->eqn_number(n)];
1001 }
1002 }
1003
1004 // Add the global parameter as an unknown to the problem
1005 problem_pt->Dof_pt.push_back(parameter_pt);
1006
1007
1008 // Normalise the initial guesses for the eigenvecor
1009 double length = 0.0;
1010 for (unsigned n = 0; n < Ndof; n++)
1011 {
1012 length += eigenvector[n] * eigenvector[n];
1013 }
1014 length = sqrt(length);
1015
1016 // Now add the null space components to the problem unknowns
1017 // and initialise them and Phi to the same normalised values
1018 for (unsigned n = 0; n < Ndof; n++)
1019 {
1020 problem_pt->Dof_pt.push_back(&Y[n]);
1021 Y[n] = Phi[n] = eigenvector[n] / length;
1022 }
1023
1024 // delete the dist_pt
1025 problem_pt->Dof_distribution_pt->build(
1026 problem_pt->communicator_pt(), Ndof * 2 + 1, true);
1027 // Remove all previous sparse storage used during Jacobian assembly
1029
1030 delete dist_pt;
1031 }
1032
1033
1034 /// Constructor in which the eigenvector and normalisation
1035 /// vector are passed as an initial guess
1037 double* const& parameter_pt,
1038 const DoubleVector& eigenvector,
1039 const DoubleVector& normalisation)
1040 : Solve_which_system(Full_augmented), Parameter_pt(parameter_pt)
1041 {
1042 // Set the problem pointer
1043 Problem_pt = problem_pt;
1044 // Set the number of degrees of freedom
1045 Ndof = problem_pt->ndof();
1046
1047 // create the linear algebra distribution for this solver
1048 // currently only global (non-distributed) distributions are allowed
1049 LinearAlgebraDistribution* dist_pt =
1050 new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
1051
1052 // Resize the vectors of additional dofs and constants
1053 Phi.resize(Ndof);
1054 Y.resize(Ndof);
1055 Count.resize(Ndof, 0);
1056
1057 // Loop over all the elements in the problem
1058 unsigned n_element = problem_pt->mesh_pt()->nelement();
1059 for (unsigned e = 0; e < n_element; e++)
1060 {
1061 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
1062 // Loop over the local freedoms in the element
1063 unsigned n_var = elem_pt->ndof();
1064 for (unsigned n = 0; n < n_var; n++)
1065 {
1066 // Increase the associated global equation number counter
1067 ++Count[elem_pt->eqn_number(n)];
1068 }
1069 }
1070
1071 // Add the global parameter as an unknown to the problem
1072 problem_pt->Dof_pt.push_back(parameter_pt);
1073
1074
1075 // Normalise the initial guesses for the eigenvecor
1076 double length = 0.0;
1077 for (unsigned n = 0; n < Ndof; n++)
1078 {
1079 length += eigenvector[n] * normalisation[n];
1080 }
1081 length = sqrt(length);
1082
1083 // Now add the null space components to the problem unknowns
1084 // and initialise them and Phi to the same normalised values
1085 for (unsigned n = 0; n < Ndof; n++)
1086 {
1087 problem_pt->Dof_pt.push_back(&Y[n]);
1088 Y[n] = eigenvector[n] / length;
1089 Phi[n] = normalisation[n];
1090 }
1091
1092 // delete the dist_pt
1093 problem_pt->Dof_distribution_pt->build(
1094 problem_pt->communicator_pt(), Ndof * 2 + 1, true);
1095 // Remove all previous sparse storage used during Jacobian assembly
1097
1098 delete dist_pt;
1099 }
1100
1101
1102 //=================================================================
1103 /// The number of unknowns is 2n+1
1104 //================================================================
1105 unsigned FoldHandler::ndof(GeneralisedElement* const& elem_pt)
1106 {
1107 unsigned raw_ndof = elem_pt->ndof();
1108 // Return different values depending on the type of block decomposition
1109 switch (Solve_which_system)
1110 {
1111 case Full_augmented:
1112 return (2 * raw_ndof + 1);
1113 break;
1114
1115 case Block_augmented_J:
1116 return (raw_ndof + 1);
1117 break;
1118
1119 case Block_J:
1120 return raw_ndof;
1121 break;
1122
1123 default:
1124 std::ostringstream error_stream;
1125 error_stream
1126 << "The Solve_which_system flag can only take values 0, 1, 2"
1127 << " not " << Solve_which_system << "\n";
1128 throw OomphLibError(
1129 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1130 }
1131 }
1132
1133 //=====================================================================
1134 /// Return the global equation number associated with local equation
1135 /// number ieqn_local. We choose to number the unknowns according
1136 /// to the augmented system.
1137 //=======================================================================
1138 unsigned long FoldHandler::eqn_number(GeneralisedElement* const& elem_pt,
1139 const unsigned& ieqn_local)
1140 {
1141 // Find the number of non-augmented dofs in the element
1142 unsigned raw_ndof = elem_pt->ndof();
1143 // Storage for the global eqn number
1144 unsigned long global_eqn = 0;
1145 // If we are a "standard" unknown, just return the standard equation number
1146 if (ieqn_local < raw_ndof)
1147 {
1148 global_eqn = elem_pt->eqn_number(ieqn_local);
1149 }
1150 // Otherwise if we are at an unknown corresponding to the bifurcation
1151 // parameter return the global equation number of the parameter
1152 else if (ieqn_local == raw_ndof)
1153 {
1154 global_eqn = Ndof;
1155 }
1156 // Otherwise we are in the unknown corresponding to a null vector
1157 // return the global unknown Ndof + 1 + global unknown of "standard"
1158 // unknown.
1159 else
1160 {
1161 global_eqn = Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof);
1162 }
1163
1164 // Return the global equation number
1165 return global_eqn;
1166 }
1167
1168 //====================================================================
1169 /// Formulate the augmented system
1170 //====================================================================
1172 Vector<double>& residuals)
1173 {
1174 // Need to get raw residuals and jacobian
1175 unsigned raw_ndof = elem_pt->ndof();
1176
1177 // Find out which system we are solving
1178 switch (Solve_which_system)
1179 {
1180 // If we are solving the standard system
1181 case Block_J:
1182 {
1183 // Get the basic residuals
1184 elem_pt->get_residuals(residuals);
1185 }
1186 break;
1187
1188 // If we are solving the augmented-by-one system
1189 case Block_augmented_J:
1190 {
1191 // Get the basic residuals
1192 elem_pt->get_residuals(residuals);
1193
1194 // Zero the final residual
1195 residuals[raw_ndof] = 0.0;
1196 }
1197 break;
1198
1199 // If we are solving the full augmented system
1200 case Full_augmented:
1201 {
1202 DenseMatrix<double> jacobian(raw_ndof);
1203 // Get the basic residuals and jacobian initially
1204 elem_pt->get_jacobian(residuals, jacobian);
1205
1206 // The normalisation equation must be initialised to -1.0/number of
1207 // elements
1208 residuals[raw_ndof] = -1.0 / Problem_pt->mesh_pt()->nelement();
1209
1210 // Now assemble the equations Jy = 0
1211 for (unsigned i = 0; i < raw_ndof; i++)
1212 {
1213 residuals[raw_ndof + 1 + i] = 0.0;
1214 for (unsigned j = 0; j < raw_ndof; j++)
1215 {
1216 residuals[raw_ndof + 1 + i] +=
1217 jacobian(i, j) * Y[elem_pt->eqn_number(j)];
1218 }
1219 // Add the contribution to phi.y=1
1220 // Need to divide by the number of elements that contribute to this
1221 // unknown so that we do get phi.y exactly.
1222 unsigned global_eqn = elem_pt->eqn_number(i);
1223 residuals[raw_ndof] +=
1224 (Phi[global_eqn] * Y[global_eqn]) / Count[global_eqn];
1225 }
1226 }
1227 break;
1228
1229 default:
1230 std::ostringstream error_stream;
1231 error_stream
1232 << "The Solve_which_system flag can only take values 0, 1, 2"
1233 << " not " << Solve_which_system << "\n";
1234 throw OomphLibError(
1235 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1236 }
1237 }
1238
1239
1240 //=================================================================
1241 /// Calculate the elemental Jacobian matrix "d equation
1242 /// / d variable" for the augmented system
1243 //=================================================================
1245 Vector<double>& residuals,
1246 DenseMatrix<double>& jacobian)
1247 {
1248 // Find the number of augmented dofs
1249 unsigned augmented_ndof = ndof(elem_pt);
1250 // Find the non-augmented dofs
1251 unsigned raw_ndof = elem_pt->ndof();
1252
1253 // Which system are we solving
1254 switch (Solve_which_system)
1255 {
1256 // If we are solving the original system
1257 case Block_J:
1258 {
1259 // Just get the raw jacobian and residuals
1260 elem_pt->get_jacobian(residuals, jacobian);
1261 }
1262 break;
1263
1264 // If we are solving the augmented-by-one system
1265 case Block_augmented_J:
1266 {
1267 // Get the full residuals, we need them
1268 get_residuals(elem_pt, residuals);
1269
1270 // Need to get the raw jacobian (and raw residuals)
1271 Vector<double> newres(augmented_ndof);
1272 elem_pt->get_jacobian(newres, jacobian);
1273
1274 // Now do finite differencing stuff
1275 const double FD_step = 1.0e-8;
1276 // Fill in the first lot of finite differences
1277 {
1278 // increase the global parameter
1279 double* unknown_pt = Problem_pt->Dof_pt[Ndof];
1280 double init = *unknown_pt;
1281 *unknown_pt += FD_step;
1282
1283 // Now do any possible updates
1285
1286 // Get the new (modified) residuals
1287 get_residuals(elem_pt, newres);
1288
1289 // The final column is given by the difference
1290 // between the residuals
1291 for (unsigned n = 0; n < raw_ndof; n++)
1292 {
1293 jacobian(n, augmented_ndof - 1) =
1294 (newres[n] - residuals[n]) / FD_step;
1295 }
1296 // Reset the global parameter
1297 *unknown_pt = init;
1298
1299 // Now do any possible updates
1301 }
1302
1303 // Fill in the bottom row
1304 for (unsigned n = 0; n < raw_ndof; n++)
1305 {
1306 unsigned local_eqn = elem_pt->eqn_number(n);
1307 jacobian(augmented_ndof - 1, n) = Phi[local_eqn] / Count[local_eqn];
1308 }
1309 }
1310 break;
1311
1312 // Otherwise solving the full system
1313 case Full_augmented:
1314 {
1315 // Get the residuals for the augmented system
1316 get_residuals(elem_pt, residuals);
1317
1318 // Need to get the raw residuals and jacobian
1319 Vector<double> newres(raw_ndof);
1320 DenseMatrix<double> newjac(raw_ndof);
1321 elem_pt->get_jacobian(newres, jacobian);
1322
1323 // Fill in the jacobian on the diagonal sub-block of
1324 // the null-space equations
1325 for (unsigned n = 0; n < raw_ndof; n++)
1326 {
1327 for (unsigned m = 0; m < raw_ndof; m++)
1328 {
1329 jacobian(raw_ndof + 1 + n, raw_ndof + 1 + m) = jacobian(n, m);
1330 }
1331 }
1332
1333 // Now finite difference wrt the global unknown
1334 const double FD_step = 1.0e-8;
1335 // Fill in the first lot of finite differences
1336 {
1337 // increase the global parameter
1338 double* unknown_pt = Problem_pt->Dof_pt[Ndof];
1339 double init = *unknown_pt;
1340 *unknown_pt += FD_step;
1341 // Need to update the function
1343
1344 // Get the new raw residuals and jacobian
1345 elem_pt->get_jacobian(newres, newjac);
1346
1347 // The end of the first row is given by the difference
1348 // between the residuals
1349 for (unsigned n = 0; n < raw_ndof; n++)
1350 {
1351 jacobian(n, raw_ndof) = (newres[n] - residuals[n]) / FD_step;
1352 // The end of the second row is given by the difference multiplied
1353 // by the product
1354 for (unsigned l = 0; l < raw_ndof; l++)
1355 {
1356 jacobian(raw_ndof + 1 + n, raw_ndof) +=
1357 (newjac(n, l) - jacobian(n, l)) * Y[elem_pt->eqn_number(l)] /
1358 FD_step;
1359 }
1360 }
1361 // Reset the global parameter
1362 *unknown_pt = init;
1363
1364 // Need to update the function
1366 }
1367
1368 // Now fill in the first column of the second rows
1369 {
1370 for (unsigned n = 0; n < raw_ndof; n++)
1371 {
1372 unsigned long global_eqn = eqn_number(elem_pt, n);
1373 // Increase the first lot
1374 double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
1375 double init = *unknown_pt;
1376 *unknown_pt += FD_step;
1378
1379 // Get the new jacobian
1380 elem_pt->get_jacobian(newres, newjac);
1381
1382 // Work out the differences
1383 for (unsigned k = 0; k < raw_ndof; k++)
1384 {
1385 // jacobian(raw_ndof+k,n) = 0.0;
1386 for (unsigned l = 0; l < raw_ndof; l++)
1387 {
1388 jacobian(raw_ndof + 1 + k, n) +=
1389 (newjac(k, l) - jacobian(k, l)) * Y[elem_pt->eqn_number(l)] /
1390 FD_step;
1391 }
1392 }
1393 *unknown_pt = init;
1395 }
1396 }
1397
1398 // Fill in the row corresponding to the parameter
1399 for (unsigned n = 0; n < raw_ndof; n++)
1400 {
1401 unsigned global_eqn = elem_pt->eqn_number(n);
1402 jacobian(raw_ndof, raw_ndof + 1 + n) =
1403 Phi[global_eqn] / Count[global_eqn];
1404 }
1405
1406 // //Loop over the dofs
1407 // for(unsigned n=0;n<augmented_ndof;n++)
1408 // {
1409 // unsigned long global_eqn = eqn_number(elem_pt,n);
1410 // double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
1411 // double init = *unknown_pt;
1412 // *unknown_pt += FD_step;
1413
1414 // //Get the new residuals
1415 // get_residuals(elem_pt,newres);
1416
1417 // for(unsigned m=0;m<augmented_ndof;m++)
1418 // {
1419 // jacobian(m,n) = (newres[m] - residuals[m])/FD_step;
1420 // }
1421 // //Reset the unknown
1422 // *unknown_pt = init;
1423 // }
1424 }
1425 break;
1426
1427 default:
1428 std::ostringstream error_stream;
1429 error_stream
1430 << "The Solve_which_system flag can only take values 0, 1, 2"
1431 << " not " << Solve_which_system << "\n";
1432 throw OomphLibError(
1433 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1434 }
1435 }
1436
1437
1438 //====================================================================
1439 /// Formulate the derivatives of the augmented system with respect
1440 /// to a parameter
1441 //====================================================================
1443 GeneralisedElement* const& elem_pt,
1444 double* const& parameter_pt,
1445 Vector<double>& dres_dparam)
1446 {
1447 // Need to get raw residuals and jacobian
1448 unsigned raw_ndof = elem_pt->ndof();
1449
1450 // Find out which system we are solving
1451 switch (Solve_which_system)
1452 {
1453 // If we are solving the standard system
1454 case Block_J:
1455 {
1456 // Get the basic residual derivatives
1457 elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
1458 }
1459 break;
1460
1461 // If we are solving the augmented-by-one system
1462 case Block_augmented_J:
1463 {
1464 // Get the basic residual derivatives
1465 elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
1466
1467 // Zero the final derivative
1468 dres_dparam[raw_ndof] = 0.0;
1469 }
1470 break;
1471
1472 // If we are solving the full augmented system
1473 case Full_augmented:
1474 {
1475 DenseMatrix<double> djac_dparam(raw_ndof);
1476 // Get the basic residuals and jacobian derivatives initially
1477 elem_pt->get_djacobian_dparameter(
1478 parameter_pt, dres_dparam, djac_dparam);
1479
1480 // The normalisation equation does not depend on the parameter
1481 dres_dparam[raw_ndof] = 0.0;
1482
1483 // Now assemble the equations dJy/dparameter = 0
1484 for (unsigned i = 0; i < raw_ndof; i++)
1485 {
1486 dres_dparam[raw_ndof + 1 + i] = 0.0;
1487 for (unsigned j = 0; j < raw_ndof; j++)
1488 {
1489 dres_dparam[raw_ndof + 1 + i] +=
1490 djac_dparam(i, j) * Y[elem_pt->eqn_number(j)];
1491 }
1492 }
1493 }
1494 break;
1495
1496 default:
1497 std::ostringstream error_stream;
1498 error_stream
1499 << "The Solve_which_system flag can only take values 0, 1, 2"
1500 << " not " << Solve_which_system << "\n";
1501 throw OomphLibError(
1502 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1503 }
1504 }
1505
1506
1507 //========================================================================
1508 /// Overload the derivative of the residuals and jacobian
1509 /// with respect to a parameter so that it breaks because it should not
1510 /// be required
1511 //========================================================================
1513 double* const& parameter_pt,
1514 Vector<double>& dres_dparam,
1515 DenseMatrix<double>& djac_dparam)
1516 {
1517 std::ostringstream error_stream;
1518 error_stream
1519 << "This function has not been implemented because it is not required\n";
1520 error_stream << "in standard problems.\n";
1521 error_stream
1522 << "If you find that you need it, you will have to implement it!\n\n";
1523
1524 throw OomphLibError(
1525 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1526 }
1527
1528
1529 //=====================================================================
1530 /// Overload the hessian vector product function so that
1531 /// it breaks because it should not be required
1532 //========================================================================
1534 GeneralisedElement* const& elem_pt,
1535 Vector<double> const& Y,
1536 DenseMatrix<double> const& C,
1537 DenseMatrix<double>& product)
1538 {
1539 std::ostringstream error_stream;
1540 error_stream
1541 << "This function has not been implemented because it is not required\n";
1542 error_stream << "in standard problems.\n";
1543 error_stream
1544 << "If you find that you need it, you will have to implement it!\n\n";
1545
1546 throw OomphLibError(
1547 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1548 }
1549
1550
1551 //==========================================================================
1552 /// Return the eigenfunction(s) associated with the bifurcation that
1553 /// has been detected in bifurcation tracking problems
1554 //==========================================================================
1556 {
1557 // There is only one (real) null vector
1558 eigenfunction.resize(1);
1560 // Rebuild the vector
1561 eigenfunction[0].build(&dist, 0.0);
1562 // Set the value to be the null vector
1563 for (unsigned n = 0; n < Ndof; n++)
1564 {
1565 eigenfunction[0][n] = Y[n];
1566 }
1567 }
1568
1569
1570 //=======================================================================
1571 /// Destructor return the problem to its original (non-augmented) state
1572 //=======================================================================
1574 {
1575 // If we are using the block solver reset the problem's linear solver
1576 // to the original one
1577 AugmentedBlockFoldLinearSolver* block_fold_solver_pt =
1578 dynamic_cast<AugmentedBlockFoldLinearSolver*>(
1580
1581 if (block_fold_solver_pt)
1582 {
1583 // Reset the problem's linear solver
1584 Problem_pt->linear_solver_pt() = block_fold_solver_pt->linear_solver_pt();
1585 // Delete the block solver
1586 delete block_fold_solver_pt;
1587 }
1588
1589 // Resize the number of dofs
1590 Problem_pt->Dof_pt.resize(Ndof);
1592 Problem_pt->communicator_pt(), Ndof, false);
1593 // Remove all previous sparse storage used during Jacobian assembly
1595 }
1596
1597 //====================================================================
1598 /// Set to solve the augmented-by-one block system.
1599 //===================================================================
1601 {
1602 // Only bother to do anything if we haven't already set the flag
1604 {
1605 // If we were solving the system with the original jacobian add the
1606 // parameter
1608 {
1609 Problem_pt->Dof_pt.push_back(Parameter_pt);
1610 }
1611
1612 // Restrict the problem to the standard variables and
1613 // the bifurcation parameter only
1614 Problem_pt->Dof_pt.resize(Ndof + 1);
1616 Problem_pt->communicator_pt(), Ndof + 1, false);
1617 // Remove all previous sparse storage used during Jacobian assembly
1619 // Set the solve flag
1621 }
1622 }
1623
1624
1625 //====================================================================
1626 /// Set to solve the block system. The boolean flag specifies
1627 /// whether the block decomposition should use exactly the same jacobian
1628 //===================================================================
1630 {
1631 // Only bother to do anything if we haven't already set the flag
1633 {
1634 // Restrict the problem to the standard variables
1635 Problem_pt->Dof_pt.resize(Ndof);
1637 Problem_pt->communicator_pt(), Ndof, false);
1638 // Remove all previous sparse storage used during Jacobian assembly
1640
1641 // Set the solve flag
1643 }
1644 }
1645
1646 //=================================================================
1647 /// Set to Solve non-block system
1648 //=================================================================
1650 {
1651 // Only do something if we are not solving the full system
1653 {
1654 // If we were solving the problem without any augmentation,
1655 // add the parameter again
1657 {
1658 Problem_pt->Dof_pt.push_back(Parameter_pt);
1659 }
1660
1661 // Always add the additional unknowns back into the problem
1662 for (unsigned n = 0; n < Ndof; n++)
1663 {
1664 Problem_pt->Dof_pt.push_back(&Y[n]);
1665 }
1666
1667 // update the Dof distribution pt
1669 Problem_pt->communicator_pt(), Ndof * 2 + 1, false);
1670 // Remove all previous sparse storage used during Jacobian assembly
1672
1674 }
1675 }
1676
1677
1678 //======================================================================
1679 /// Clean up the memory that may have been allocated by the solver
1680 //=====================================================================
1682 {
1683 if (B_pt != 0)
1684 {
1685 delete B_pt;
1686 }
1687 if (C_pt != 0)
1688 {
1689 delete C_pt;
1690 }
1691 if (D_pt != 0)
1692 {
1693 delete D_pt;
1694 }
1695 if (dJy_dparam_pt != 0)
1696 {
1697 delete dJy_dparam_pt;
1698 }
1699 }
1700
1701 //===================================================================
1702 /// Use a block factorisation to solve the augmented system
1703 /// associated with a PitchFork bifurcation.
1704 //===================================================================
1706 DoubleVector& result)
1707 {
1708 std::cout << "Block pitchfork solve" << std::endl;
1709 // Locally cache the pointer to the handler.
1710 PitchForkHandler* handler_pt =
1711 static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
1712
1713 // Get the augmented distribution from the handler and use it as
1714 // the distribution of the linear solver
1715 LinearAlgebraDistribution aug_dist =
1717 this->build_distribution(aug_dist);
1718
1719 // If the result vector is not setup then die
1720 if (!result.built())
1721 {
1722 throw OomphLibError("Result vector must be built\n",
1723 OOMPH_CURRENT_FUNCTION,
1724 OOMPH_EXCEPTION_LOCATION);
1725 }
1726
1727 // Store the distribution of the result, which is probably not in
1728 // the natural distribution of the augmented solver
1729 LinearAlgebraDistribution result_dist(result.distribution_pt());
1730
1731 // Locally cache a pointer to the parameter
1732 double* const parameter_pt = handler_pt->Parameter_pt;
1733
1734 // Firstly get the derivatives of the full augmented system with
1735 // respect to the parameter
1736
1737 // Allocate storage for the derivatives of the residuals with respect
1738 // to the global parameter using the augmented distribution
1739 DoubleVector dRdparam;
1740 // Then get the appropriate derivatives which will come back in
1741 // some distribution
1742 problem_pt->get_derivative_wrt_global_parameter(parameter_pt, dRdparam);
1743 // Redistribute into the augmented distribution
1744 dRdparam.redistribute(&aug_dist);
1745
1746 // Switch the handler to "block solver" mode (sort out distribution)
1747 handler_pt->solve_block_system();
1748
1749 // Temporary vector for the result (I shouldn't have to set this up)
1750 DoubleVector x1;
1751
1752 // We are going to do resolves using the underlying linear solver
1754 // Solve the first (standard) system Jx1 = R
1755 Linear_solver_pt->solve(problem_pt, x1);
1756
1757 // Allocate storage for B, C and D which can be used in the resolve
1758 // and the derivatives of the jacobian/eigenvector product with
1759 // respect to the parameter
1760 if (B_pt != 0)
1761 {
1762 delete B_pt;
1763 }
1765 if (C_pt != 0)
1766 {
1767 delete C_pt;
1768 }
1770 if (D_pt != 0)
1771 {
1772 delete D_pt;
1773 }
1775 // Need this to be in the distribution of the dofs
1776 if (dJy_dparam_pt != 0)
1777 {
1778 delete dJy_dparam_pt;
1779 }
1780 dJy_dparam_pt = new DoubleVector(handler_pt->Dof_distribution_pt, 0.0);
1781
1782 // Get the symmetry vector from the handler should have the
1783 // Dof distribution
1784 DoubleVector psi = handler_pt->Psi;
1785
1786 // Temporary vector for the rhs that is dR/dparam (augmented distribution)
1787 // DoubleVector F(Linear_solver_pt->distribution_pt(),0.0);
1788 DoubleVector F(handler_pt->Dof_distribution_pt);
1789 const unsigned n_dof_local = F.nrow_local();
1790
1791 // f.nrow local copied from dRdparam
1792 for (unsigned n = 0; n < n_dof_local; n++)
1793 {
1794 F[n] = dRdparam[n];
1795 }
1796 // Fill in the rhs that is dJy/dparam //dJy_dparam nrow local
1797
1798 // In standard cases the offset here will be one
1799 unsigned offset = 1;
1800#ifdef OOMPH_HAS_MPI
1801 // If we are distributed and not on the first processor
1802 // then there is no offset and the eigenfunction is
1803 // directly after the standard dofs
1804 int my_rank = problem_pt->communicator_pt()->my_rank();
1805 if ((problem_pt->distributed()) && (my_rank != 0))
1806 {
1807 offset = 0;
1808 }
1809#endif
1810 for (unsigned n = 0; n < n_dof_local; n++)
1811 {
1812 (*dJy_dparam_pt)[n] = dRdparam[n_dof_local + offset + n];
1813 }
1814
1815 // Now resolve to find c and d
1816 // First, redistribute F and psi into the linear algebra distribution
1817 F.redistribute(Linear_solver_pt->distribution_pt());
1819
1822
1823 // We can now construct various dot products
1824 double psi_d = psi.dot(*D_pt);
1825 double psi_c = psi.dot(*C_pt);
1826 double psi_x1 = psi.dot(x1);
1827
1828 // Calculate another intermediate constant
1829 const double Psi = psi_d / psi_c;
1830
1831 // Construct the second intermediate value,
1832 // assumption is that result has been
1833 // set to the current value of the residuals
1834 // First, redistribute into the Natural distribution of the augmented system
1835 result.redistribute(&aug_dist);
1836 // The required parameter is that corresponding to the dof, which
1837 // is only stored on the root processor
1838 double parameter_residual = result[n_dof_local];
1839#ifdef OOMPH_HAS_MPI
1840 // Broadcast to all others, if we have a distributed problem
1841 if (problem_pt->distributed())
1842 {
1843 MPI_Bcast(&parameter_residual,
1844 1,
1845 MPI_DOUBLE,
1846 0,
1847 problem_pt->communicator_pt()->mpi_comm());
1848 }
1849#endif
1850 // Now construct the value
1851 double x2 = (psi_x1 - parameter_residual) / psi_c;
1852
1853 // Now construct the vectors that multiply the jacobian terms
1855 D_and_X1[0].build(Linear_solver_pt->distribution_pt(), 0.0);
1856 D_and_X1[1].build(Linear_solver_pt->distribution_pt(), 0.0);
1857 // Fill in the appropriate terms
1858 // Get the number of local dofs from the Linear_solver_pt distribution
1859 const unsigned n_dof_local_linear_solver =
1861 for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
1862 {
1863 const double C_ = (*C_pt)[n];
1864 D_and_X1[0][n] = (*D_pt)[n] - Psi * C_;
1865 D_and_X1[1][n] = x1[n] - x2 * C_;
1866 }
1867
1868 // Local storage for the result of the product terms
1869 Vector<DoubleVectorWithHaloEntries> Jprod_D_and_X1(2);
1870
1871 // Redistribute the inputs to have the Dof distribution pt
1872 D_and_X1[0].redistribute(handler_pt->Dof_distribution_pt);
1873 D_and_X1[1].redistribute(handler_pt->Dof_distribution_pt);
1874
1875 // Set up the halo scheme
1876#ifdef OOMPH_HAS_MPI
1877 D_and_X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
1878 D_and_X1[1].build_halo_scheme(problem_pt->Halo_scheme_pt);
1879#endif
1880
1881 // Get the products from the new problem function
1882 problem_pt->get_hessian_vector_products(
1883 handler_pt->Y, D_and_X1, Jprod_D_and_X1);
1884
1885 // OK, now we can formulate the next vectors
1886 //(again assuming result contains residuals)
1887 // Need to redistribute F to the dof distribution
1888 // F.redistribute(handler_pt->Dof_distribution_pt);
1889 DoubleVector G(handler_pt->Dof_distribution_pt);
1890
1891 for (unsigned n = 0; n < n_dof_local; n++)
1892 {
1893 G[n] = result[n_dof_local + offset + n] - Jprod_D_and_X1[1][n] -
1894 x2 * dRdparam[n_dof_local + offset + n];
1895 Jprod_D_and_X1[0][n] *= -1.0;
1896 Jprod_D_and_X1[0][n] -= Psi * dRdparam[n_dof_local + offset + n];
1897 }
1898
1899
1900 // Then redistribute back to the linear solver distribution
1902 Jprod_D_and_X1[0].redistribute(Linear_solver_pt->distribution_pt());
1903 Jprod_D_and_X1[1].redistribute(Linear_solver_pt->distribution_pt());
1904
1905 // Linear solve to get B
1906 Linear_solver_pt->resolve(Jprod_D_and_X1[0], *B_pt);
1907 // Liner solve to get x3
1909 Linear_solver_pt->resolve(G, x3);
1910
1911 // Construst a couple of additional products
1912 double l_x3 = psi.dot(x3);
1913 double l_b = psi.dot(*B_pt);
1914
1915 // get the last intermediate variable
1916 // The required parameter is that corresponding to the dof, which
1917 // is only stored on the root processor
1918 double sigma_residual = result[2 * (n_dof_local + offset) - 1];
1919#ifdef OOMPH_HAS_MPI
1920 // Broadcast to all others, if we have a distributed problem
1921 if (problem_pt->distributed())
1922 {
1923 MPI_Bcast(&sigma_residual,
1924 1,
1925 MPI_DOUBLE,
1926 0,
1927 problem_pt->communicator_pt()->mpi_comm());
1928 }
1929#endif
1930
1931
1932 const double delta_sigma = (l_x3 - sigma_residual) / l_b;
1933 const double delta_lambda = x2 - delta_sigma * Psi;
1934
1935 // Need to do some more rearrangements here because result is global
1936 // but the other vectors are not!
1937
1938 // Create temporary DoubleVectors to hold the results
1941
1942 for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
1943 {
1944 res1[n] = x1[n] - delta_lambda * (*C_pt)[n] - delta_sigma * (*D_pt)[n];
1945 res2[n] = x3[n] - delta_sigma * (*B_pt)[n];
1946 }
1947
1948 // Now redistribute these into the Dof distribution pointer
1949 res1.redistribute(handler_pt->Dof_distribution_pt);
1950 res2.redistribute(handler_pt->Dof_distribution_pt);
1951
1952 // Now can copy over into results into the result vector
1953 for (unsigned n = 0; n < n_dof_local; n++)
1954 {
1955 result[n] = res1[n];
1956 result[n_dof_local + offset + n] = res2[n];
1957 }
1958
1959 // Add the final contributions to the residuals
1960 // only on the root processor if we have a distributed problem
1961#ifdef OOMPH_HAS_MPI
1962 if ((!problem_pt->distributed()) || (my_rank == 0))
1963#endif
1964 {
1965 result[n_dof_local] = delta_lambda;
1966 result[2 * n_dof_local + 1] = delta_sigma;
1967 }
1968
1969
1970 // The sign of the determinant is given by the sign of
1971 // the product psi_c and l_b
1972 // NOT CHECKED YET!
1973 problem_pt->sign_of_jacobian() =
1974 static_cast<int>(std::fabs(psi_c * l_b) / (psi_c * l_b));
1975
1976 // Redistribute the result into its incoming distribution
1977 result.redistribute(&result_dist);
1978
1979 // Switch things to our block solver
1980 handler_pt->solve_full_system();
1981
1982 // If we are not storing things, clear the results
1983 if (!Enable_resolve)
1984 {
1985 // We no longer need to store the matrix
1987 delete B_pt;
1988 B_pt = 0;
1989 delete C_pt;
1990 C_pt = 0;
1991 delete D_pt;
1992 D_pt = 0;
1993 delete dJy_dparam_pt;
1994 dJy_dparam_pt = 0;
1995 }
1996 // Otherwise also store the pointer to the problem
1997 else
1998 {
1999 Problem_pt = problem_pt;
2000 }
2001 }
2002
2003
2004 //==============================================================
2005 // Hack the re-solve to use the block factorisation
2006 //==============================================================
2008 DoubleVector& result)
2009 {
2010 std::cout << "Block pitchfork resolve" << std::endl;
2011 // Check that the factors have been stored
2012 if (B_pt == 0)
2013 {
2014 throw OomphLibError("The required vectors have not been stored",
2015 OOMPH_CURRENT_FUNCTION,
2016 OOMPH_EXCEPTION_LOCATION);
2017 }
2018
2019
2020 // Cache pointer to the problem
2021 Problem* const problem_pt = Problem_pt;
2022
2023 // Locally cache pointer to the handler
2024 PitchForkHandler* handler_pt =
2025 static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
2026
2027 // Get the augmented distribution from the handler
2028 LinearAlgebraDistribution aug_dist =
2030 // this->build_distribution(aug_dist);
2031
2032 // Find the number of dofs of the augmented system
2033 // const unsigned n_aug_dof = problem_pt->ndof();
2034
2035 // Storage for the result distribution
2036 LinearAlgebraDistribution result_dist;
2037
2038 // if the result vector is not setup then rebuild with distribution
2039 // = natural distribution of augmented solver
2040 if (!result.built())
2041 {
2042 result.build(&aug_dist, 0.0);
2043 }
2044 // Otherwise store the incoming distribution and redistribute
2045 else
2046 {
2047 result_dist.build(result.distribution_pt());
2048 result.redistribute(&aug_dist);
2049 }
2050
2051
2052 // Locally cache a pointer to the parameter
2053 // double* const parameter_pt = handler_pt->Parameter_pt;
2054
2055 // Switch things to our block solver
2056 handler_pt->solve_block_system();
2057 // We need to find out the number of dofs
2058 // unsigned n_dof = problem_pt->ndof();
2059
2060 // create the linear algebra distribution for this solver
2061 // currently only global (non-distributed) distributions are allowed
2062 // LinearAlgebraDistribution
2063 // dist(problem_pt->communicator_pt(),n_dof,false);
2064 // this->build_distribution(dist);
2065
2066 // if the result vector is not setup then rebuild with distribution = global
2067 // if (!result.built())
2068 // {
2069 // result.build(this->distribution_pt(),0.0);
2070 // }
2071
2072
2073 // Copy the rhs into local storage
2074 DoubleVector rhs_local = rhs;
2075 // and redistribute into the augmented distribution
2076 rhs_local.redistribute(&aug_dist);
2077
2078 // Setup storage for output
2081
2082 // Create input for RHS with the natural distribution
2083 DoubleVector input_x1(handler_pt->Dof_distribution_pt);
2084 const unsigned n_dof_local = input_x1.nrow_local();
2085
2086 // Set the values of the a vector
2087 for (unsigned n = 0; n < n_dof_local; n++)
2088 {
2089 input_x1[n] = rhs_local[n];
2090 }
2091 // Need to redistribute into the linear algebra distribution
2093
2094 // We want to resolve, of course
2096 // Now solve for the first vector
2097 Linear_solver_pt->resolve(input_x1, x1);
2098
2099 // Get the symmetry vector from the handler
2100 DoubleVector psi = handler_pt->Psi;
2101 // redistribute local copy
2103
2104 // We can now construct various dot products
2105 double psi_d = psi.dot(*D_pt);
2106 double psi_c = psi.dot(*C_pt);
2107 double psi_x1 = psi.dot(x1);
2108
2109 // Calculate another intermediate constant
2110 const double Psi = psi_d / psi_c;
2111
2112 // Construct the second intermediate value,
2113 // assumption is that rhs has been set to the current value of the residuals
2114 double parameter_residual = rhs_local[n_dof_local];
2115#ifdef OOMPH_HAS_MPI
2116 // Broadcast to all others, if we have a distributed problem
2117 if (problem_pt->distributed())
2118 {
2119 MPI_Bcast(&parameter_residual,
2120 1,
2121 MPI_DOUBLE,
2122 0,
2123 problem_pt->communicator_pt()->mpi_comm());
2124 }
2125#endif
2126
2127 double x2 = (psi_x1 - parameter_residual) / psi_c;
2128
2129 // Now construct the vectors that multiply the jacobian terms
2130 // Vector<double> X1(n_dof/*+1*/);
2132 X1[0].build(Linear_solver_pt->distribution_pt(), 0.0);
2133
2134 const unsigned n_dof_local_linear_solver =
2136 for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
2137 {
2138 X1[0][n] = x1[n] - x2 * (*C_pt)[n];
2139 }
2140
2141 // Local storage for the product term
2143
2144 // Redistribute the inputs to have the Dof distribution pt
2145 X1[0].redistribute(handler_pt->Dof_distribution_pt);
2146
2147 // Set up the halo scheme
2148#ifdef OOMPH_HAS_MPI
2149 X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
2150#endif
2151
2152 // Get the product from the problem
2153 problem_pt->get_hessian_vector_products(handler_pt->Y, X1, Jprod_X1);
2154
2155 // In standard cases the offset here will be one
2156 unsigned offset = 1;
2157#ifdef OOMPH_HAS_MPI
2158 // If we are distributed and not on the first processor
2159 // then there is no offset and the eigenfunction is
2160 // directly after the standard dofs
2161 int my_rank = problem_pt->communicator_pt()->my_rank();
2162 if ((problem_pt->distributed()) && (my_rank != 0))
2163 {
2164 offset = 0;
2165 }
2166#endif
2167
2168 // OK, now we can formulate the next vectors
2169 //(again assuming result contains residuals)
2170 // Local storage for the product terms
2171 DoubleVector Mod_Jprod_X1(handler_pt->Dof_distribution_pt, 0.0);
2172
2173 for (unsigned n = 0; n < n_dof_local; n++)
2174 {
2175 Mod_Jprod_X1[n] = rhs_local[n_dof_local + offset + n] - Jprod_X1[0][n] -
2176 x2 * (*dJy_dparam_pt)[n];
2177 }
2178
2179 // Redistribute back to the linear solver distribution
2181
2182 // Liner solve to get x3
2183 Linear_solver_pt->resolve(Mod_Jprod_X1, x3);
2184
2185 // Construst a couple of additional products
2186 double l_x3 = psi.dot(x3);
2187 double l_b = psi.dot(*B_pt);
2188
2189 // get the last intermediate variable
2190 // The required parameter is that corresponding to the dof, which
2191 // is only stored on the root processor
2192 double sigma_residual = rhs_local[2 * (n_dof_local + offset) - 1];
2193#ifdef OOMPH_HAS_MPI
2194 // Broadcast to all others, if we have a distributed problem
2195 if (problem_pt->distributed())
2196 {
2197 MPI_Bcast(&sigma_residual,
2198 1,
2199 MPI_DOUBLE,
2200 0,
2201 problem_pt->communicator_pt()->mpi_comm());
2202 }
2203#endif
2204
2205 // get the last intermediate variable
2206 const double delta_sigma = (l_x3 - sigma_residual) / l_b;
2207 const double delta_lambda = x2 - delta_sigma * Psi;
2208
2209 // Create temporary DoubleVectors to hold the results
2212
2213 for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
2214 {
2215 res1[n] = x1[n] - delta_lambda * (*C_pt)[n] - delta_sigma * (*D_pt)[n];
2216 res2[n] = x3[n] - delta_sigma * (*B_pt)[n];
2217 }
2218
2219 // Now redistribute these into the Dof distribution pointer
2220 res1.redistribute(handler_pt->Dof_distribution_pt);
2221 res2.redistribute(handler_pt->Dof_distribution_pt);
2222
2223 // Now can copy over into results into the result vector
2224 for (unsigned n = 0; n < n_dof_local; n++)
2225 {
2226 result[n] = res1[n];
2227 result[n_dof_local + offset + n] = res2[n];
2228 }
2229
2230 // Add the final contributions to the residuals
2231 // only on the root processor if we have a distributed problem
2232#ifdef OOMPH_HAS_MPI
2233 if ((!problem_pt->distributed()) || (my_rank == 0))
2234#endif
2235 {
2236 result[n_dof_local] = delta_lambda;
2237 result[2 * n_dof_local + 1] = delta_sigma;
2238 }
2239
2240 // Redistribute the result into its incoming distribution (if it had one)
2241 if (result_dist.built())
2242 {
2243 result.redistribute(&result_dist);
2244 }
2245
2247
2248 // Switch things to our block solver
2249 handler_pt->solve_full_system();
2250 }
2251
2252
2253 //--------------------------------------------------------------
2254
2255
2256 //======================================================================
2257 /// Clean up the memory that may have been allocated by the solver
2258 //=====================================================================
2260 {
2261 if (Alpha_pt != 0)
2262 {
2263 delete Alpha_pt;
2264 }
2265 if (E_pt != 0)
2266 {
2267 delete E_pt;
2268 }
2269 }
2270
2271 //===================================================================
2272 /// Use a block factorisation to solve the augmented system
2273 /// associated with a PitchFork bifurcation.
2274 //===================================================================
2276 DoubleVector& result)
2277 {
2278 std::cout << "Augmented pitchfork solve" << std::endl;
2279 // if the result is setup then it should not be distributed
2280#ifdef PARANOID
2281 if (result.built())
2282 {
2283 if (result.distributed())
2284 {
2285 throw OomphLibError("The result vector must not be distributed",
2286 OOMPH_CURRENT_FUNCTION,
2287 OOMPH_EXCEPTION_LOCATION);
2288 }
2289 }
2290#endif
2291
2292
2293 // Locally cache the pointer to the handler.
2294 PitchForkHandler* handler_pt =
2295 static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
2296
2297 // Switch the handler to "block solver" mode
2298 handler_pt->solve_augmented_block_system();
2299
2300 // We need to find out the number of dofs in the problem
2301 unsigned n_dof = problem_pt->ndof();
2302
2303 // create the linear algebra distribution for this solver
2304 // currently only global (non-distributed) distributions are allowed
2305 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
2306 this->build_distribution(dist);
2307
2308 // if the result vector is not setup then rebuild with distribution = global
2309 if (!result.built())
2310 {
2311 result.build(this->distribution_pt(), 0.0);
2312 }
2313
2314 // Setup storage for temporary vectors
2315 DoubleVector a(this->distribution_pt(), 0.0),
2316 b(this->distribution_pt(), 0.0);
2317
2318 // Allocate storage for Alpha which can be used in the resolve
2319 if (Alpha_pt != 0)
2320 {
2321 delete Alpha_pt;
2322 }
2323 Alpha_pt = new DoubleVector(this->distribution_pt(), 0.0);
2324
2325 // We are going to do resolves using the underlying linear solver
2327 // Solve the first system Aa = R
2328 Linear_solver_pt->solve(problem_pt, a);
2329
2330 // Get the symmetry vector from the handler
2331 DoubleVector psi(this->distribution_pt(), 0.0);
2332 for (unsigned n = 0; n < (n_dof - 1); ++n)
2333 {
2334 psi[n] = handler_pt->Psi[n];
2335 }
2336 // Set the final entry to zero
2337 psi[n_dof - 1] = 0.0;
2338
2339 // Now resolve to find alpha
2341
2342 // We can now construct our multipliers
2343 // Prepare to scale
2344 double dof_length = 0.0, a_length = 0.0, alpha_length = 0.0;
2345 for (unsigned n = 0; n < n_dof; n++)
2346 {
2347 if (std::fabs(problem_pt->dof(n)) > dof_length)
2348 {
2349 dof_length = std::fabs(problem_pt->dof(n));
2350 }
2351 if (std::fabs(a[n]) > a_length)
2352 {
2353 a_length = std::fabs(a[n]);
2354 }
2355 if (std::fabs((*Alpha_pt)[n]) > alpha_length)
2356 {
2357 alpha_length = std::fabs((*Alpha_pt)[n]);
2358 }
2359 }
2360
2361 double a_mult = dof_length / a_length;
2362 double alpha_mult = dof_length / alpha_length;
2363 const double FD_step = 1.0e-8;
2364 a_mult += FD_step;
2365 alpha_mult += FD_step;
2366 a_mult *= FD_step;
2367 alpha_mult *= FD_step;
2368
2369 // Local storage for the product terms
2370 DoubleVector Jprod_a(this->distribution_pt(), 0.0),
2371 Jprod_alpha(this->distribution_pt(), 0.0);
2372
2373 // Calculate the product of the jacobian matrices, etc
2374 unsigned long n_element = problem_pt->mesh_pt()->nelement();
2375 for (unsigned long e = 0; e < n_element; e++)
2376 {
2377 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
2378 // Loop over the ndofs in each element
2379 unsigned n_var = handler_pt->ndof(elem_pt);
2380 // Get the jacobian matrices
2381 DenseMatrix<double> jac(n_var), jac_a(n_var), jac_alpha(n_var);
2382 // the elemental residual
2383 Vector<double> res_elemental(n_var);
2384 // Get unperturbed jacobian
2385 handler_pt->get_jacobian(elem_pt, res_elemental, jac);
2386
2387 // Backup the dofs
2388 Vector<double> dof_bac(n_var);
2389 // Perturb the dofs
2390 for (unsigned n = 0; n < n_var; n++)
2391 {
2392 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
2393 dof_bac[n] = problem_pt->dof(eqn_number);
2394 // Pertub by vector a
2395 problem_pt->dof(eqn_number) += a_mult * a[eqn_number];
2396 }
2397
2399
2400 // Now get the new jacobian
2401 handler_pt->get_jacobian(elem_pt, res_elemental, jac_a);
2402
2403 // Perturb the dofs
2404 for (unsigned n = 0; n < n_var; n++)
2405 {
2406 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
2407 problem_pt->dof(eqn_number) = dof_bac[n];
2408 // Pertub by vector a
2409 problem_pt->dof(eqn_number) += alpha_mult * (*Alpha_pt)[eqn_number];
2410 }
2411
2413
2414 // Now get the new jacobian
2415 handler_pt->get_jacobian(elem_pt, res_elemental, jac_alpha);
2416
2417 // Reset the dofs
2418 for (unsigned n = 0; n < n_var; n++)
2419 {
2420 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
2421 problem_pt->dof(eqn_number) = dof_bac[n];
2422 }
2423
2425
2426 // OK, now work out the products
2427 for (unsigned n = 0; n < (n_var - 1); n++)
2428 {
2429 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
2430 double prod_a = 0.0, prod_alpha = 0.0;
2431 for (unsigned m = 0; m < (n_var - 1); m++)
2432 {
2433 unsigned unknown = handler_pt->eqn_number(elem_pt, m);
2434 prod_a += (jac_a(n, m) - jac(n, m)) * handler_pt->Y[unknown];
2435 prod_alpha += (jac_alpha(n, m) - jac(n, m)) * handler_pt->Y[unknown];
2436 }
2437 Jprod_a[eqn_number] += prod_a / a_mult;
2438 Jprod_alpha[eqn_number] += prod_alpha / alpha_mult;
2439 }
2440 }
2441
2442 Jprod_alpha[n_dof - 1] = 0.0;
2443 Jprod_a[n_dof - 1] = 0.0;
2444
2445 // OK, now we can formulate the next vectors
2446 // The assumption here is that the result has been set to the
2447 // residuals.
2448 for (unsigned n = 0; n < n_dof - 1; n++)
2449 {
2450 b[n] = result[n_dof + n] - Jprod_a[n];
2451 }
2452 b[n_dof - 1] = result[2 * n_dof - 1];
2453
2454
2455 // Allocate storage for E which can be used in the resolve
2456 if (E_pt != 0)
2457 {
2458 delete E_pt;
2459 }
2460 E_pt = new DoubleVector(this->distribution_pt(), 0.0);
2461
2462 DoubleVector f(this->distribution_pt(), 0.0);
2463
2465 Linear_solver_pt->resolve(Jprod_alpha, *E_pt);
2466
2467 // Calculate the final entry in the vector e
2468 const double e_final = (*E_pt)[n_dof - 1];
2469 // Calculate the final entry in the vector d
2470 const double d_final = -f[n_dof - 1] / e_final;
2471 // Assemble the final corrections
2472 for (unsigned n = 0; n < n_dof - 1; n++)
2473 {
2474 result[n] = a[n] - (*Alpha_pt)[n] * d_final;
2475 result[n_dof + n] = f[n] + (*E_pt)[n] * d_final;
2476 }
2477
2478 result[n_dof - 1] = a[n_dof - 1] - (*Alpha_pt)[n_dof - 1] * d_final;
2479 result[2 * n_dof - 1] = d_final;
2480
2481 // The sign of the jacobian is the sign of the final entry in e
2482 problem_pt->sign_of_jacobian() =
2483 static_cast<int>(std::fabs(e_final) / e_final);
2484
2485
2486 // Switch things to our block solver
2487 handler_pt->solve_full_system();
2488
2489 // If we are not storing things, clear the results
2490 if (!Enable_resolve)
2491 {
2492 // We no longer need to store the matrix
2494 delete Alpha_pt;
2495 Alpha_pt = 0;
2496 delete E_pt;
2497 E_pt = 0;
2498 }
2499 // Otherwise also store the pointer to the problem
2500 else
2501 {
2502 Problem_pt = problem_pt;
2503 }
2504 }
2505
2506
2507 //==============================================================
2508 // Hack the re-solve to use the block factorisation
2509 //==============================================================
2511 DoubleVector& result)
2512 {
2513 std::cout << "Augmented pitchfork resolve" << std::endl;
2514 // Check that the factors have been stored
2515 if (Alpha_pt == 0)
2516 {
2517 throw OomphLibError("The required vectors have not been stored",
2518 OOMPH_CURRENT_FUNCTION,
2519 OOMPH_EXCEPTION_LOCATION);
2520 }
2521
2522 // Get the pointer to the problem
2523 Problem* const problem_pt = Problem_pt;
2524
2525 PitchForkHandler* handler_pt =
2526 static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
2527
2528 // Switch things to our block solver
2529 handler_pt->solve_augmented_block_system();
2530 // We need to find out the number of dofs
2531 unsigned n_dof = problem_pt->ndof();
2532
2533 // create the linear algebra distribution for this solver
2534 // currently only global (non-distributed) distributions are allowed
2535 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
2536 this->build_distribution(dist);
2537
2538 // if the result vector is not setup then rebuild with distribution = global
2539 if (!result.built())
2540 {
2541 result.build(this->distribution_pt(), 0.0);
2542 }
2543
2544
2545 // Setup storage
2546 DoubleVector a(this->distribution_pt(), 0.0),
2547 b(this->distribution_pt(), 0.0);
2548
2549 // Set the values of the a vector
2550 for (unsigned n = 0; n < n_dof; n++)
2551 {
2552 a[n] = rhs[n];
2553 }
2554
2556
2557 // Copy rhs vector into local storage so it doesn't get overwritten
2558 // if the linear solver decides to initialise the solution vector, say,
2559 // which it's quite entitled to do!
2560 DoubleVector input_a(a);
2561
2562 Linear_solver_pt->resolve(input_a, a);
2563
2564 // We can now construct our multipliers
2565 // Prepare to scale
2566 double dof_length = 0.0, a_length = 0.0;
2567 for (unsigned n = 0; n < n_dof; n++)
2568 {
2569 if (std::fabs(problem_pt->dof(n)) > dof_length)
2570 {
2571 dof_length = std::fabs(problem_pt->dof(n));
2572 }
2573
2574 if (std::fabs(a[n]) > a_length)
2575 {
2576 a_length = std::fabs(a[n]);
2577 }
2578 }
2579 double a_mult = dof_length / a_length;
2580 const double FD_step = 1.0e-8;
2581 a_mult += FD_step;
2582 a_mult *= FD_step;
2583
2584 DoubleVector Jprod_a(this->distribution_pt(), 0.0);
2585
2586 unsigned long n_element = problem_pt->mesh_pt()->nelement();
2587 for (unsigned long e = 0; e < n_element; e++)
2588 {
2589 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
2590 // Loop over the ndofs in each element
2591 unsigned n_var = handler_pt->ndof(elem_pt);
2592 // Get some jacobian matrices
2593 DenseMatrix<double> jac(n_var), jac_a(n_var);
2594 // the elemental residual
2595 Vector<double> res_elemental(n_var);
2596 // Get unperturbed jacobian
2597 handler_pt->get_jacobian(elem_pt, res_elemental, jac);
2598
2599 // Backup the dofs
2600 DoubleVector dof_bac(this->distribution_pt(), 0.0);
2601 // Perturb the dofs
2602 for (unsigned n = 0; n < n_var; n++)
2603 {
2604 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
2605 dof_bac[n] = problem_pt->dof(eqn_number);
2606 // Pertub by vector a
2607 problem_pt->dof(eqn_number) += a_mult * a[eqn_number];
2608 }
2609
2611
2612 // Now get the new jacobian
2613 handler_pt->get_jacobian(elem_pt, res_elemental, jac_a);
2614
2615 // Reset the dofs
2616 for (unsigned n = 0; n < n_var; n++)
2617 {
2618 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
2619 problem_pt->dof(eqn_number) = dof_bac[n];
2620 }
2621
2623
2624 // OK, now work out the products
2625 for (unsigned n = 0; n < (n_var - 1); n++)
2626 {
2627 unsigned eqn_number = handler_pt->eqn_number(elem_pt, n);
2628 double prod_a = 0.0;
2629 for (unsigned m = 0; m < (n_var - 1); m++)
2630 {
2631 unsigned unknown = handler_pt->eqn_number(elem_pt, m);
2632 prod_a += (jac_a(n, m) - jac(n, m)) * handler_pt->Y[unknown];
2633 }
2634 Jprod_a[eqn_number] += prod_a / a_mult;
2635 }
2636 }
2637
2638 Jprod_a[n_dof - 1] = 0.0;
2639
2640 // OK, now we can formulate the next vectors
2641 for (unsigned n = 0; n < n_dof - 1; n++)
2642 {
2643 b[n] = rhs[n_dof + n] - Jprod_a[n];
2644 }
2645 b[n_dof - 1] = rhs[2 * n_dof - 1];
2646
2647 DoubleVector f(this->distribution_pt(), 0.0);
2648
2650
2651 // Calculate the final entry in the vector d
2652 const double d_final = -f[n_dof - 1] / (*E_pt)[n_dof - 1];
2653 // Assemble the final corrections
2654 for (unsigned n = 0; n < n_dof - 1; n++)
2655 {
2656 result[n] = a[n] - (*Alpha_pt)[n] * d_final;
2657 result[n_dof + n] = f[n] + (*E_pt)[n] * d_final;
2658 }
2659
2660 result[n_dof - 1] = a[n_dof - 1] - (*Alpha_pt)[n_dof - 1] * d_final;
2661 result[2 * n_dof - 1] = d_final;
2662
2664
2665 // Switch things to our block solver
2666 handler_pt->solve_full_system();
2667 }
2668
2669 /// ////////////////////////////////////////////////////////////////////
2670 // Non-inline functions for the PitchForkHandler class
2671 /// ///////////////////////////////////////////////////////////////////
2672
2673 //==================================================================
2674 /// Constructor: Initialise the PitchForkHandler by setting intial
2675 /// guesses for Sigma, Y, specifying C and Psi and calculating count.
2676 //==================================================================
2678 Problem* const& problem_pt,
2679 AssemblyHandler* const& assembly_handler_pt,
2680 double* const& parameter_pt,
2681 const DoubleVector& symmetry_vector)
2682 : Solve_which_system(Full_augmented), Sigma(0.0), Parameter_pt(parameter_pt)
2683 {
2684 // Set the problem pointer
2685 Problem_pt = problem_pt;
2686 // Set the assembly handler
2687 Assembly_handler_pt = assembly_handler_pt;
2688 // Set the number of degrees of freedom
2689 Ndof = Problem_pt->ndof();
2690 // Backup the distribution
2692#ifdef OOMPH_HAS_MPI
2693 // Set the distribution flag
2695#endif
2696
2697 // Find the elements in the problem
2698 unsigned n_element = Problem_pt->mesh_pt()->nelement();
2699
2700#ifdef OOMPH_HAS_MPI
2701
2702 // Work out all the global equations to which this processor
2703 // contributes and store the halo data in the problem
2705
2706#endif
2707
2708
2709 // Now use the dof distribution for all double vectors
2714#ifdef OOMPH_HAS_MPI
2716#endif
2717
2718 // Loop over the elements and count the entries
2719 // and number of (non-halo) elements
2720 unsigned n_non_halo_element_local = 0;
2721 for (unsigned e = 0; e < n_element; e++)
2722 {
2724#ifdef OOMPH_HAS_MPI
2725 // Ignore halo elements
2726 if (!elem_pt->is_halo())
2727 {
2728#endif
2729 // Increment the number of non halo elements
2730 ++n_non_halo_element_local;
2731 // Now count the number of times the element contributes to a value
2732 unsigned n_var = assembly_handler_pt->ndof(elem_pt);
2733 for (unsigned n = 0; n < n_var; n++)
2734 {
2735 ++Count.global_value(assembly_handler_pt->eqn_number(elem_pt, n));
2736 }
2737#ifdef OOMPH_HAS_MPI
2738 }
2739#endif
2740 }
2741
2742 // Add together all the counts
2743#ifdef OOMPH_HAS_MPI
2745
2746 // If distributed, find the total number of elements in the problem
2747 if (Distributed)
2748 {
2749 // Need to gather the total number of non halo elements
2750 MPI_Allreduce(&n_non_halo_element_local,
2751 &Nelement,
2752 1,
2753 MPI_UNSIGNED,
2754 MPI_SUM,
2755 Problem_pt->communicator_pt()->mpi_comm());
2756 }
2757 // Otherwise the total number is the same on each processor
2758 else
2759#endif
2760 {
2761 Nelement = n_non_halo_element_local;
2762 }
2763
2764#ifdef OOMPH_HAS_MPI
2765 // Only add the parameter to the first processor, if distributed
2766 int my_rank = Problem_pt->communicator_pt()->my_rank();
2767 if ((!Distributed) || (my_rank == 0))
2768#endif
2769 {
2770 // Add the parameter to the problem
2771 Problem_pt->Dof_pt.push_back(parameter_pt);
2772 }
2773
2774 // Find length of the symmetry vector
2775 double length = symmetry_vector.norm();
2776
2777 // Add the unknowns for the null vector to the problem
2778 // Normalise the symmetry vector and initialise the null vector to the
2779 // symmetry vector and set the constant vector, c, to the
2780 // normalised symmetry vector.
2781
2782 // Need to redistribute the symmetry vector into the (natural)
2783 // distribution of the Dofs
2784#ifdef OOMPH_HAS_MPI
2785 // Check that the symmetry vector has the correct distribution
2786 if (!(*symmetry_vector.distribution_pt() == *Dof_distribution_pt))
2787 {
2788 throw OomphLibError(
2789 "The symmetry vector must have the same distribution as the dofs\n",
2790 OOMPH_CURRENT_FUNCTION,
2791 OOMPH_EXCEPTION_LOCATION);
2792 }
2793#endif
2794
2795 // Only loop over the local unknowns
2796 const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
2797 for (unsigned n = 0; n < n_dof_local; n++)
2798 {
2799 Problem_pt->Dof_pt.push_back(&Y[n]);
2800 // Psi[n] = symmetry_vector[n];
2801 Psi[n] = Y[n] = C[n] = symmetry_vector[n] / length;
2802 }
2803
2804
2805#ifdef OOMPH_HAS_MPI
2806 // Set up the required halo schemes (which also synchronises)
2810
2811
2812 if ((!Distributed) || (my_rank == 0))
2813#endif
2814 // Add the slack parameter to the problem on the first processor
2815 // if distributed
2816 {
2817 Problem_pt->Dof_pt.push_back(&Sigma);
2818 }
2819
2820#ifdef OOMPH_HAS_MPI
2821 if (Distributed)
2822 {
2823 unsigned augmented_first_row = 0;
2824 unsigned augmented_n_row_local = 0;
2825
2826 // Set up the translation scheme on every processor
2827 Global_eqn_number.resize(2 * Ndof + 2);
2828 int n_proc = Problem_pt->communicator_pt()->nproc();
2829 unsigned global_eqn_count = 0;
2830 for (int d = 0; d < n_proc; d++)
2831 {
2832 // Find out the first row of the current processor
2833 if (my_rank == d)
2834 {
2835 augmented_first_row = global_eqn_count;
2836 }
2837
2838 const unsigned n_row_local = Dof_distribution_pt->nrow_local(d);
2839 const unsigned first_row = Dof_distribution_pt->first_row(d);
2840 // Add the basic equations
2841 for (unsigned n = 0; n < n_row_local; n++)
2842 {
2843 Global_eqn_number[first_row + n] = global_eqn_count;
2844 ++global_eqn_count;
2845 }
2846 // If on the first processor add the pointer to the parameter
2847 if (d == 0)
2848 {
2849 Global_eqn_number[Ndof] = global_eqn_count;
2850 ++global_eqn_count;
2851 }
2852 // Add the eigenfunction
2853 for (unsigned n = 0; n < n_row_local; n++)
2854 {
2855 Global_eqn_number[Ndof + 1 + first_row + n] = global_eqn_count;
2856 ++global_eqn_count;
2857 }
2858 // Finally add the slack parameter
2859 if (d == 0)
2860 {
2861 Global_eqn_number[2 * Ndof + 1] = global_eqn_count;
2862 ++global_eqn_count;
2863 }
2864 // Find out the number of rows of the current processor
2865 if (my_rank == d)
2866 {
2867 augmented_n_row_local = global_eqn_count - augmented_first_row;
2868 }
2869 }
2870
2871 // Make a new linear algebra distribution
2874 augmented_first_row,
2875 augmented_n_row_local);
2876 }
2877 else
2878#endif
2879 {
2881 Problem_pt->communicator_pt(), 2 * Ndof + 2, false);
2882 }
2883
2884 // resize
2885 // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
2886 // 2*Ndof+2,false);
2887
2889
2890 // Remove all previous sparse storage used during Jacobian assembly
2892 }
2893
2894 //=====================================================================
2895 /// Destructor return the problem to its original (non-augmented) state
2896 //=====================================================================
2898 {
2899 // If we are using the block solver reset the problem's linear solver
2900 // to the original one
2901 BlockPitchForkLinearSolver* block_pitchfork_solver_pt =
2903
2904 if (block_pitchfork_solver_pt)
2905 {
2906 // Reset the problem's linear solver
2908 block_pitchfork_solver_pt->linear_solver_pt();
2909 // Delete the block solver
2910 delete block_pitchfork_solver_pt;
2911 }
2912
2913 // If we are using the augmented
2914 // block solver reset the problem's linear solver
2915 // to the original one
2916 AugmentedBlockPitchForkLinearSolver* augmented_block_pitchfork_solver_pt =
2919
2920 if (augmented_block_pitchfork_solver_pt)
2921 {
2922 // Reset the problem's linear solver
2924 augmented_block_pitchfork_solver_pt->linear_solver_pt();
2925 // Delete the block solver
2926 delete augmented_block_pitchfork_solver_pt;
2927 }
2928
2929 // Now return the problem to its original size
2931
2932 // Problem_pt->Dof_pt.resize(Ndof);
2933 // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
2934 // Ndof,false);
2935
2937 // Remove all previous sparse storage used during Jacobian assembly
2939 }
2940
2941 //================================================================
2942 /// Get the number of elemental degrees of freedom
2943 //================================================================
2945 {
2946 unsigned raw_ndof = elem_pt->ndof();
2947
2948 // Return different values depending on the type of block decomposition
2949 switch (Solve_which_system)
2950 {
2951 case Full_augmented:
2952 return (2 * raw_ndof + 2);
2953 break;
2954
2955 case Block_augmented_J:
2956 return (raw_ndof + 1);
2957 break;
2958
2959 case Block_J:
2960 return raw_ndof;
2961 break;
2962
2963 default:
2964 std::ostringstream error_stream;
2965 error_stream
2966 << "The Solve_which_system flag can only take values 0, 1, 2"
2967 << " not " << Solve_which_system << "\n";
2968 throw OomphLibError(
2969 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2970 }
2971 }
2972
2973 /// Get the global equation number of the local unknown
2975 const unsigned& ieqn_local)
2976 {
2977 // Get the raw value
2978 unsigned raw_ndof = elem_pt->ndof();
2979 unsigned long global_eqn = 0;
2980
2981 // Which system are we solving
2982 switch (Solve_which_system)
2983 {
2984 // In the block case, it's just the standard numbering
2985 case Block_J:
2986 global_eqn = elem_pt->eqn_number(ieqn_local);
2987 break;
2988
2989 // Block augmented not done properly yet
2990 case Block_augmented_J:
2991 // Not done the distributed case
2992#ifdef OOMPH_HAS_MPI
2993 if (Distributed)
2994 {
2995 throw OomphLibError(
2996 "Block Augmented solver not implemented for distributed case\n",
2997 OOMPH_CURRENT_FUNCTION,
2998 OOMPH_EXCEPTION_LOCATION);
2999 }
3000#endif
3001
3002 // Full case
3003 case Full_augmented:
3004 // The usual equations
3005 if (ieqn_local < raw_ndof)
3006 {
3007 global_eqn = this->global_eqn_number(elem_pt->eqn_number(ieqn_local));
3008 }
3009 // The bifurcation parameter equation
3010 else if (ieqn_local == raw_ndof)
3011 {
3012 global_eqn = this->global_eqn_number(Ndof);
3013 }
3014 // If we are assembling the full system we also have
3015 // The components of the null vector
3016 else if (ieqn_local < (2 * raw_ndof + 1))
3017 {
3018 global_eqn = this->global_eqn_number(
3019 Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof));
3020 }
3021 // The slack parameter
3022 else
3023 {
3024 global_eqn = this->global_eqn_number(2 * Ndof + 1);
3025 }
3026 break;
3027 } // End of switch
3028 return global_eqn;
3029 }
3030
3031 //==============================================================
3032 /// Get the residuals
3033 //==============================================================
3035 Vector<double>& residuals)
3036 {
3037 // Need to get raw residuals and jacobian
3038 unsigned raw_ndof = elem_pt->ndof();
3039
3040 // Find out which system we are solving
3041 switch (Solve_which_system)
3042 {
3043 // If we are solving the original system
3044 case Block_J:
3045 {
3046 // get the basic residuals
3047 elem_pt->get_residuals(residuals);
3048 // Now multiply to fill in the residuals for the final term
3049 for (unsigned i = 0; i < raw_ndof; i++)
3050 {
3051 unsigned local_eqn = elem_pt->eqn_number(i);
3052 // Add the slack parameter to the final residuals
3053 residuals[i] +=
3054 Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
3055 }
3056 }
3057 break;
3058
3059 // If we are solving the augmented-by-one system
3060 case Block_augmented_J:
3061 {
3062 // Get the basic residuals
3063 elem_pt->get_residuals(residuals);
3064
3065 // Zero the final residual
3066 residuals[raw_ndof] = 0.0;
3067 // Now multiply to fill in the residuals for the final term
3068 for (unsigned i = 0; i < raw_ndof; i++)
3069 {
3070 unsigned local_eqn = elem_pt->eqn_number(i);
3071 // Add the slack parameter to the final residuals
3072 residuals[i] +=
3073 Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
3074 // Final term that specifies the symmetry
3075 residuals[raw_ndof] += ((*Problem_pt->global_dof_pt(local_eqn)) *
3076 Psi.global_value(local_eqn)) /
3077 Count.global_value(local_eqn);
3078 }
3079 }
3080 break;
3081
3082 // Otherwise we are solving the fully augemented system
3083 case Full_augmented:
3084 {
3085 DenseMatrix<double> jacobian(raw_ndof);
3086
3087 // Get the basic residuals and jacobian
3088 elem_pt->get_jacobian(residuals, jacobian);
3089
3090 // Initialise the final residuals
3091 residuals[raw_ndof] = 0.0;
3092 residuals[2 * raw_ndof + 1] = -1.0 / Nelement;
3093
3094 // Now multiply to fill in the residuals associated
3095 // with the null vector condition
3096 for (unsigned i = 0; i < raw_ndof; i++)
3097 {
3098 unsigned local_eqn = elem_pt->eqn_number(i);
3099 residuals[raw_ndof + 1 + i] = 0.0;
3100 for (unsigned j = 0; j < raw_ndof; j++)
3101 {
3102 unsigned local_unknown = elem_pt->eqn_number(j);
3103 residuals[raw_ndof + 1 + i] +=
3104 jacobian(i, j) * Y.global_value(local_unknown);
3105 }
3106
3107 // Add the slack parameter to the governing equations
3108 residuals[i] +=
3109 Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
3110
3111 // Specify the symmetry
3112 residuals[raw_ndof] += ((*Problem_pt->global_dof_pt(local_eqn)) *
3113 Psi.global_value(local_eqn)) /
3114 Count.global_value(local_eqn);
3115 // Specify the normalisation
3116 residuals[2 * raw_ndof + 1] +=
3117 (Y.global_value(local_eqn) * C.global_value(local_eqn)) /
3118 Count.global_value(local_eqn);
3119 }
3120 }
3121 break;
3122
3123 default:
3124 std::ostringstream error_stream;
3125 error_stream
3126 << "The Solve_which_system flag can only take values 0, 1, 2"
3127 << " not " << Solve_which_system << "\n";
3128 throw OomphLibError(
3129 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3130 }
3131 }
3132
3133 //======================================================================
3134 /// Calculate the elemental Jacobian matrix "d equation
3135 /// / d variable".
3136 //======================================================================
3138 Vector<double>& residuals,
3139 DenseMatrix<double>& jacobian)
3140 {
3141 unsigned augmented_ndof = ndof(elem_pt);
3142 unsigned raw_ndof = elem_pt->ndof();
3143
3144 // Which system are we solving
3145 switch (Solve_which_system)
3146 {
3147 // If we are solving the original system
3148 case Block_J:
3149 {
3150 // get the raw residuals and jacobian
3151 elem_pt->get_jacobian(residuals, jacobian);
3152
3153 // Now multiply to fill in the residuals for the final term
3154 for (unsigned i = 0; i < raw_ndof; i++)
3155 {
3156 unsigned local_eqn = elem_pt->eqn_number(i);
3157 // Add the slack parameter to the final residuals
3158 residuals[i] +=
3159 Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
3160 }
3161 }
3162 break;
3163
3164 // If we are solving the augmented-by-one system
3165 case Block_augmented_J:
3166 {
3167 // Get the full residuals, we need them
3168 get_residuals(elem_pt, residuals);
3169
3170 // Need to get the raw jacobian (and raw residuals)
3171 Vector<double> newres(augmented_ndof);
3172 elem_pt->get_jacobian(newres, jacobian);
3173
3174 // Now do finite differencing stuff
3175 const double FD_step = 1.0e-8;
3176 // Fill in the first lot of finite differences
3177 {
3178 // increase the global parameter
3179 double* unknown_pt = Parameter_pt;
3180 double init = *unknown_pt;
3181 *unknown_pt += FD_step;
3182
3183 // Not do any possible updates
3185
3186 // Get the new (modified) residuals
3187 get_residuals(elem_pt, newres);
3188
3189 // The final column is given by the difference
3190 // between the residuals
3191 for (unsigned n = 0; n < raw_ndof; n++)
3192 {
3193 jacobian(n, augmented_ndof - 1) =
3194 (newres[n] - residuals[n]) / FD_step;
3195 }
3196 // Reset the global parameter
3197 *unknown_pt = init;
3198
3199 // Now do any possible updates
3201 }
3202
3203 // Fill in the bottom row
3204 for (unsigned n = 0; n < raw_ndof; n++)
3205 {
3206 unsigned local_eqn = elem_pt->eqn_number(n);
3207 jacobian(augmented_ndof - 1, n) =
3208 Psi.global_value(local_eqn) / Count.global_value(local_eqn);
3209 }
3210 }
3211 break;
3212
3213 // Otherwise solving the full system
3214 case Full_augmented:
3215 {
3216 /// ALICE:
3218 // Get the basic jacobian and residuals
3219 elem_pt->get_jacobian(residuals, jacobian);
3220 // get the proper residuals
3221 get_residuals(elem_pt, residuals);
3222
3223 // Now fill in the next diagonal jacobian entry
3224 for (unsigned n = 0; n < raw_ndof; n++)
3225 {
3226 for (unsigned m = 0; m < raw_ndof; m++)
3227 {
3228 jacobian(raw_ndof + 1 + n, raw_ndof + 1 + m) = jacobian(n, m);
3229 }
3230 unsigned local_eqn = elem_pt->eqn_number(n);
3231 // Add in the sigma contribution
3232 jacobian(n, 2 * raw_ndof + 1) =
3233 Psi.global_value(local_eqn) / Count.global_value(local_eqn);
3234 // Symmetry constraint
3235 jacobian(raw_ndof, n) =
3236 Psi.global_value(local_eqn) / Count.global_value(local_eqn);
3237 // Non-zero constraint
3238 jacobian(2 * raw_ndof + 1, raw_ndof + 1 + n) =
3239 C.global_value(local_eqn) / Count.global_value(local_eqn);
3240 }
3241
3242 // Finite difference the remaining blocks
3243 const double FD_step = 1.0e-8;
3244
3245 Vector<double> newres_p(augmented_ndof);
3246
3247 // Loop over the ndofs only
3248 for (unsigned n = 0; n < raw_ndof; ++n)
3249 {
3250 // Get the original (unaugmented) global equation number
3251 unsigned long global_eqn =
3252 Assembly_handler_pt->eqn_number(elem_pt, n);
3253 double* unknown_pt = Problem_pt->global_dof_pt(global_eqn);
3254 double init = *unknown_pt;
3255 *unknown_pt += FD_step;
3257 // Get the new residuals
3258 get_residuals(elem_pt, newres_p);
3259 // Fill in the entries in the block d(Jy)/dx
3260 for (unsigned m = 0; m < raw_ndof; m++)
3261 {
3262 jacobian(raw_ndof + 1 + m, n) =
3263 (newres_p[raw_ndof + 1 + m] - residuals[raw_ndof + 1 + m]) /
3264 (FD_step);
3265 }
3266
3267 // Reset the unknown
3268 *unknown_pt = init;
3270 }
3271
3272 {
3273 // Now increase the global parameter
3274 double* unknown_pt = Parameter_pt;
3275 double init = *unknown_pt;
3276 *unknown_pt += FD_step;
3277
3279 // Get the new residuals
3280 get_residuals(elem_pt, newres_p);
3281
3282 // Add in the first block d/dx
3283 for (unsigned m = 0; m < raw_ndof; m++)
3284 {
3285 jacobian(m, raw_ndof) = (newres_p[m] - residuals[m]) / FD_step;
3286 }
3287 // Add in the second block d/dy
3288 for (unsigned m = raw_ndof + 1; m < augmented_ndof - 1; m++)
3289 {
3290 jacobian(m, raw_ndof) = (newres_p[m] - residuals[m]) / FD_step;
3291 }
3292
3293 // Reset the unknown
3294 *unknown_pt = init;
3296 }
3297 }
3298 break;
3299
3300 default:
3301 std::ostringstream error_stream;
3302 error_stream
3303 << "The Solve_which_system flag can only take values 0, 1, 2"
3304 << " not " << Solve_which_system << "\n";
3305 throw OomphLibError(
3306 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3307 }
3308 /// ALICE
3310 }
3311
3312
3313 //==============================================================
3314 /// Get the derivatives of the residuals with respect to a parameter
3315 //==============================================================
3317 GeneralisedElement* const& elem_pt,
3318 double* const& parameter_pt,
3319 Vector<double>& dres_dparam)
3320 {
3321 // Need to get raw residuals and jacobian
3322 unsigned raw_ndof = elem_pt->ndof();
3324 // Find out which system we are solving
3325 switch (Solve_which_system)
3326 {
3327 // If we are solving the original system
3328 case Block_J:
3329 {
3330 // get the basic residual derivatives
3331 elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
3332 // Slack parameter term does not depened explicitly on the parameter
3333 // so is not added
3334 }
3335 break;
3336
3337 // If we are solving the augmented-by-one system
3338 case Block_augmented_J:
3339 {
3340 // Get the basic residual derivatives
3341 elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
3342
3343 // Zero the final residuals derivative
3344 dres_dparam[raw_ndof] = 0.0;
3345
3346 // Other terms must not depend on the parameter
3347 }
3348 break;
3349
3350 // Otherwise we are solving the fully augemented system
3351 case Full_augmented:
3352 {
3353 DenseMatrix<double> djac_dparam(raw_ndof);
3354
3355 // Get the basic residuals and jacobian derivatives
3356 elem_pt->get_djacobian_dparameter(
3357 parameter_pt, dres_dparam, djac_dparam);
3358
3359 // The "final" residual derivatives do not depend on the parameter
3360 dres_dparam[raw_ndof] = 0.0;
3361 dres_dparam[2 * raw_ndof + 1] = 0.0;
3362
3363 // Now multiply to fill in the residuals associated
3364 // with the null vector condition
3365 for (unsigned i = 0; i < raw_ndof; i++)
3366 {
3367 dres_dparam[raw_ndof + 1 + i] = 0.0;
3368 for (unsigned j = 0; j < raw_ndof; j++)
3369 {
3370 unsigned local_unknown = elem_pt->eqn_number(j);
3371 dres_dparam[raw_ndof + 1 + i] +=
3372 djac_dparam(i, j) * Y.global_value(local_unknown);
3373 }
3374 }
3375 }
3376 break;
3377
3378 default:
3379 std::ostringstream error_stream;
3380 error_stream
3381 << "The Solve_which_system flag can only take values 0, 1, 2"
3382 << " not " << Solve_which_system << "\n";
3383 throw OomphLibError(
3384 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3385 }
3386 }
3387
3388
3389 //========================================================================
3390 /// Overload the derivative of the residuals and jacobian
3391 /// with respect to a parameter so that it breaks because it should not
3392 /// be required
3393 //========================================================================
3395 GeneralisedElement* const& elem_pt,
3396 double* const& parameter_pt,
3397 Vector<double>& dres_dparam,
3398 DenseMatrix<double>& djac_dparam)
3399 {
3400 std::ostringstream error_stream;
3401 error_stream
3402 << "This function has not been implemented because it is not required\n";
3403 error_stream << "in standard problems.\n";
3404 error_stream
3405 << "If you find that you need it, you will have to implement it!\n\n";
3406
3407 throw OomphLibError(
3408 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3409 }
3410
3411
3412 //=====================================================================
3413 /// Overload the hessian vector product function so that
3414 /// it calls the underlying element's hessian function.
3415 //========================================================================
3417 GeneralisedElement* const& elem_pt,
3418 Vector<double> const& Y,
3419 DenseMatrix<double> const& C,
3420 DenseMatrix<double>& product)
3421 {
3422 elem_pt->get_hessian_vector_products(Y, C, product);
3423 }
3424
3425
3426 //==========================================================================
3427 /// Return the eigenfunction(s) associated with the bifurcation that
3428 /// has been detected in bifurcation tracking problems
3429 //==========================================================================
3431 {
3432 // There is only one (real) null vector
3433 eigenfunction.resize(1);
3434 // Rebuild the vector
3435 eigenfunction[0].build(this->Dof_distribution_pt, 0.0);
3436 // Set the value to be the null vector
3437 const unsigned n_row_local = eigenfunction[0].nrow_local();
3438 for (unsigned n = 0; n < n_row_local; n++)
3439 {
3440 eigenfunction[0][n] = Y[n];
3441 }
3442 }
3443
3444#ifdef OOMPH_HAS_MPI
3445 //=====================================================================
3446 // Synchronise the required data
3447 //====================================================================
3449 {
3450 // Only need to bother if the problem is distributed
3451 if (Distributed)
3452 {
3453 // Need only to synchronise the eigenfunction
3454 Y.synchronise();
3455 // Also need to synchronise the parameter and the slack parameter
3456 double broadcast_data[2];
3457 broadcast_data[0] = *Parameter_pt;
3458 broadcast_data[1] = Sigma;
3459 // Broadcast from the root to all processors
3460 MPI_Bcast(broadcast_data,
3461 2,
3462 MPI_DOUBLE,
3463 0,
3464 Problem_pt->communicator_pt()->mpi_comm());
3465
3466 // Now copy received the values back
3467 *Parameter_pt = broadcast_data[0];
3468 Sigma = broadcast_data[1];
3469 }
3470 }
3471#endif
3472
3473 //====================================================================
3474 /// Set to solve the augmented-by-one block system.
3475 //===================================================================
3477 {
3478 // Only bother to do anything if we haven't already set the flag
3480 {
3481 // If we were solving the system with the original jacobian add the
3482 // parameter back
3484 {
3485 Problem_pt->Dof_pt.push_back(Parameter_pt);
3486 }
3487
3488 // Restrict the problem to the standard variables and
3489 // the bifurcation parameter only
3491 Problem_pt->communicator_pt(), Ndof + 1, false);
3492 Problem_pt->Dof_pt.resize(Ndof + 1);
3493
3494 // Problem_pt->Dof_pt.resize(Ndof+1);
3495 // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
3496 // Ndof+1,false);
3497
3498 // Remove all previous sparse storage used during Jacobian assembly
3500
3501 // Set the solve flag
3503 }
3504 }
3505
3506
3507 //==============================================================
3508 /// Set to solve the block system. The boolean flag is used
3509 /// to specify whether the block decomposition should use exactly
3510 /// the same jacobian as the original system.
3511 //==============================================================
3513 {
3514 // Only bother to do anything if we haven't already set the flag
3516 {
3517 // Restrict the problem to the standard variables
3520
3521 // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
3522 // Ndof,false);
3523
3524 // Remove all previous sparse storage used during Jacobian assembly
3526
3527 // Set the solve flag
3529 }
3530 }
3531
3532 //==============================================================
3533 /// Solve non-block system
3534 //==============================================================
3536 {
3537 // Only do something if we are not solving the full system
3539 {
3540#ifdef OOMPH_HAS_MPI
3541 int my_rank = Problem_pt->communicator_pt()->my_rank();
3542#endif
3543 // If we were solving the problem without any augementation
3544 // add the parameter again
3546 {
3547#ifdef OOMPH_HAS_MPI
3548
3549 if ((!Distributed) || (my_rank == 0))
3550#endif
3551 {
3552 Problem_pt->Dof_pt.push_back(Parameter_pt);
3553 }
3554 }
3555
3556 // Always add the additional unknowns back into the problem
3557 const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
3558 for (unsigned n = 0; n < n_dof_local; n++)
3559 {
3560 Problem_pt->Dof_pt.push_back(&Y[n]);
3561 }
3562
3563
3564#ifdef OOMPH_HAS_MPI
3565 if ((!Distributed) || (my_rank == 0))
3566#endif
3567 // Add the slack parameter to the problem on the first processor
3568 // if distributed
3569 {
3570 Problem_pt->Dof_pt.push_back(&Sigma);
3571 }
3572
3573
3574 // Delete the distribtion created in the augmented block solve
3575 if (Problem_pt->Dof_distribution_pt != this->Dof_distribution_pt)
3576 {
3578 }
3579
3580 // update the Dof distribution
3582 // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
3583 // Ndof*2+2,false);
3584 // Remove all previous sparse storage used during Jacobian assembly
3586
3587 // Set the solve flag
3589 }
3590 }
3591
3592
3593 /// ////////////////////////////////////////////////////////////////////
3594 // Non-inline functions for the BlockHopfLinearSolver class
3595 /// ///////////////////////////////////////////////////////////////////
3596
3597
3598 //======================================================================
3599 /// Clean up the memory that may have been allocated by the solver
3600 //=====================================================================
3602 {
3603 if (A_pt != 0)
3604 {
3605 delete A_pt;
3606 }
3607 if (E_pt != 0)
3608 {
3609 delete E_pt;
3610 }
3611 if (G_pt != 0)
3612 {
3613 delete G_pt;
3614 }
3615 }
3616
3617 //===================================================================
3618 /// Use a block factorisation to solve the augmented system
3619 /// associated with a Hopf bifurcation.
3620 //===================================================================
3621 void BlockHopfLinearSolver::solve(Problem* const& problem_pt,
3622 DoubleVector& result)
3623 {
3624 // if the result is setup then it should not be distributed
3625#ifdef PARANOID
3626 if (result.built())
3627 {
3628 if (result.distributed())
3629 {
3630 throw OomphLibError("The result vector must not be distributed",
3631 OOMPH_CURRENT_FUNCTION,
3632 OOMPH_EXCEPTION_LOCATION);
3633 }
3634 }
3635#endif
3636
3637 // Locally cache the pointer to the handler.
3638 HopfHandler* handler_pt =
3639 static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());
3640
3641 // Find the number of dofs in the augmented problem
3642 unsigned n_dof = problem_pt->ndof();
3643
3644 // create the linear algebra distribution for this solver
3645 // currently only global (non-distributed) distributions are allowed
3646 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
3647 this->build_distribution(dist);
3648
3649 // Firstly, let's calculate the derivative of the residuals wrt
3650 // the parameter
3651 DoubleVector dRdparam(this->distribution_pt(), 0.0);
3652
3653 const double FD_step = 1.0e-8;
3654 {
3655 // Cache pointer to the parameter (second last entry in the vector
3656 double* param_pt = &problem_pt->dof(n_dof - 2);
3657 // Backup the parameter
3658 double old_var = *param_pt;
3659 ;
3660 // Increment the parameter
3661 *param_pt += FD_step;
3663
3664 // Now calculate the new residuals
3665 problem_pt->get_residuals(dRdparam);
3666
3667 // Now calculate the difference assume original residuals in resul
3668 for (unsigned n = 0; n < n_dof; n++)
3669 {
3670 dRdparam[n] = (dRdparam[n] - result[n]) / FD_step;
3671 }
3672
3673 // Reset the parameter
3674 *param_pt = old_var;
3676 }
3677
3678 // Switch the handler to "standard" mode
3679 handler_pt->solve_standard_system();
3680
3681 // Find out the number of dofs in the non-standard problem
3682 n_dof = problem_pt->ndof();
3683
3684 // update the distribution pt
3685 dist.build(problem_pt->communicator_pt(), n_dof, false);
3686 this->build_distribution(dist);
3687
3688 // Setup storage for temporary vector
3689 DoubleVector y1(this->distribution_pt(), 0.0),
3690 alpha(this->distribution_pt(), 0.0);
3691
3692 // Allocate storage for A which can be used in the resolve
3693 if (A_pt != 0)
3694 {
3695 delete A_pt;
3696 }
3697 A_pt = new DoubleVector(this->distribution_pt(), 0.0);
3698
3699 // We are going to do resolves using the underlying linear solver
3701
3702 // Solve the first system Jy1 = R
3703 Linear_solver_pt->solve(problem_pt, y1);
3704
3705 // This should have set the sign of the jacobian, store it
3706 int sign_of_jacobian = problem_pt->sign_of_jacobian();
3707
3708 // Now let's get the appropriate bit of alpha
3709 for (unsigned n = 0; n < n_dof; n++)
3710 {
3711 alpha[n] = dRdparam[n];
3712 }
3713
3714 // Resolve to find A
3715 Linear_solver_pt->resolve(alpha, *A_pt);
3716
3717 // Now set to the complex system
3718 handler_pt->solve_complex_system();
3719
3720 // update the distribution
3721 dist.build(problem_pt->communicator_pt(), n_dof * 2, false);
3722 this->build_distribution(dist);
3723
3724 // Resize the stored vector G
3725 if (G_pt != 0)
3726 {
3727 delete G_pt;
3728 }
3729 G_pt = new DoubleVector(this->distribution_pt(), 0.0);
3730
3731 // Solve the first Zg = (Mz -My)
3732 Linear_solver_pt->solve(problem_pt, *G_pt);
3733
3734 // This should have set the sign of the complex matrix's determinant,
3735 // multiply
3736 sign_of_jacobian *= problem_pt->sign_of_jacobian();
3737
3738 // We can now construct our multipliers
3739 // Prepare to scale
3740 double dof_length = 0.0, a_length = 0.0, y1_length = 0.0;
3741 // Loop over the standard number of dofs
3742 for (unsigned n = 0; n < n_dof; n++)
3743 {
3744 if (std::fabs(problem_pt->dof(n)) > dof_length)
3745 {
3746 dof_length = std::fabs(problem_pt->dof(n));
3747 }
3748 if (std::fabs((*A_pt)[n]) > a_length)
3749 {
3750 a_length = std::fabs((*A_pt)[n]);
3751 }
3752 if (std::fabs(y1[n]) > y1_length)
3753 {
3754 y1_length = std::fabs(y1[n]);
3755 }
3756 }
3757
3758 double a_mult = dof_length / a_length;
3759 double y1_mult = dof_length / y1_length;
3760 a_mult += FD_step;
3761 y1_mult += FD_step;
3762 a_mult *= FD_step;
3763 y1_mult *= FD_step;
3764
3765 // Local storage for the product terms
3766 Vector<double> Jprod_a(2 * n_dof, 0.0), Jprod_y1(2 * n_dof, 0.0);
3767
3768 // Temporary storage
3769 Vector<double> rhs(2 * n_dof);
3770
3771 // find the number of elements
3772 unsigned long n_element = problem_pt->mesh_pt()->nelement();
3773
3774 // Calculate the product of the jacobian matrices, etc
3775 for (unsigned long e = 0; e < n_element; e++)
3776 {
3777 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
3778 // Loop over the ndofs in each element
3779 unsigned n_var = elem_pt->ndof();
3780 // Get the jacobian matrices
3781 DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var);
3782 DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var);
3783 // Get unperturbed jacobian
3784 elem_pt->get_jacobian_and_mass_matrix(rhs, jac, M);
3785
3786 // Backup the dofs
3787 Vector<double> dof_bac(n_var);
3788 // Perturb the original dofs
3789 for (unsigned n = 0; n < n_var; n++)
3790 {
3791 unsigned eqn_number = elem_pt->eqn_number(n);
3792 dof_bac[n] = problem_pt->dof(eqn_number);
3793 // Pertub by vector A
3794 problem_pt->dof(eqn_number) += a_mult * (*A_pt)[eqn_number];
3795 }
3796
3797 // Now get the new jacobian and mass matrix
3798 elem_pt->get_jacobian_and_mass_matrix(rhs, jac_a, M_a);
3799
3800 // Perturb the dofs
3801 for (unsigned n = 0; n < n_var; n++)
3802 {
3803 unsigned eqn_number = elem_pt->eqn_number(n);
3804 problem_pt->dof(eqn_number) = dof_bac[n];
3805 // Pertub by vector y1
3806 problem_pt->dof(eqn_number) += y1_mult * y1[eqn_number];
3807 }
3808
3809 // Now get the new jacobian and mass matrix
3810 elem_pt->get_jacobian_and_mass_matrix(rhs, jac_y1, M_y1);
3811
3812 // Reset the dofs
3813 for (unsigned n = 0; n < n_var; n++)
3814 {
3815 unsigned eqn_number = elem_pt->eqn_number(n);
3816 problem_pt->dof(eqn_number) = dof_bac[n];
3817 }
3818
3819 // OK, now work out the products
3820 for (unsigned n = 0; n < n_var; n++)
3821 {
3822 unsigned eqn_number = elem_pt->eqn_number(n);
3823 double prod_a1 = 0.0, prod_y11 = 0.0;
3824 double prod_a2 = 0.0, prod_y12 = 0.0;
3825 for (unsigned m = 0; m < n_var; m++)
3826 {
3827 const unsigned unknown = elem_pt->eqn_number(m);
3828 const double y = handler_pt->Phi[unknown];
3829 const double z = handler_pt->Psi[unknown];
3830 const double omega = handler_pt->Omega;
3831 // Real part (first line)
3832 prod_a1 +=
3833 (jac_a(n, m) - jac(n, m)) * y + omega * (M_a(n, m) - M(n, m)) * z;
3834 prod_y11 +=
3835 (jac_y1(n, m) - jac(n, m)) * y + omega * (M_y1(n, m) - M(n, m)) * z;
3836 // Imag par (second line)
3837 prod_a2 +=
3838 (jac_a(n, m) - jac(n, m)) * z - omega * (M_a(n, m) - M(n, m)) * y;
3839 prod_y12 +=
3840 (jac_y1(n, m) - jac(n, m)) * z - omega * (M_y1(n, m) - M(n, m)) * y;
3841 }
3842 Jprod_a[eqn_number] += prod_a1 / a_mult;
3843 Jprod_y1[eqn_number] += prod_y11 / y1_mult;
3844 Jprod_a[n_dof + eqn_number] += prod_a2 / a_mult;
3845 Jprod_y1[n_dof + eqn_number] += prod_y12 / y1_mult;
3846 }
3847 }
3848
3849 // The assumption here is still that the result has been set to the
3850 // residuals.
3851 for (unsigned n = 0; n < 2 * n_dof; n++)
3852 {
3853 rhs[n] = result[n_dof + n] - Jprod_y1[n];
3854 }
3855
3856 // Temporary storage
3857 DoubleVector y2(this->distribution_pt(), 0.0);
3858
3859 // DoubleVector for rhs
3860 DoubleVector rhs2(this->distribution_pt(), 0.0);
3861 for (unsigned i = 0; i < 2 * n_dof; i++)
3862 {
3863 rhs2[i] = rhs[i];
3864 }
3865
3866 // Solve it
3867 Linear_solver_pt->resolve(rhs2, y2);
3868
3869 // Assemble the next RHS
3870 for (unsigned n = 0; n < 2 * n_dof; n++)
3871 {
3872 rhs[n] = dRdparam[n_dof + n] - Jprod_a[n];
3873 }
3874
3875 // Resive the storage
3876 if (E_pt != 0)
3877 {
3878 delete E_pt;
3879 }
3880 E_pt = new DoubleVector(this->distribution_pt(), 0.0);
3881
3882 // Solve for the next RHS
3883 for (unsigned i = 0; i < 2 * n_dof; i++)
3884 {
3885 rhs2[i] = rhs[i];
3886 }
3887 Linear_solver_pt->resolve(rhs2, *E_pt);
3888
3889 // We can now calculate the final corrections
3890 // We need to work out a large number of dot products
3891 double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g = 0.0;
3892 double dot_h = 0.0;
3893
3894 for (unsigned n = 0; n < n_dof; n++)
3895 {
3896 // Get the appopriate entry
3897 const double Cn = handler_pt->C[n];
3898 dot_c += Cn * y2[n];
3899 dot_d += Cn * y2[n_dof + n];
3900 dot_e += Cn * (*E_pt)[n];
3901 dot_f += Cn * (*E_pt)[n_dof + n];
3902 dot_g += Cn * (*G_pt)[n];
3903 dot_h += Cn * (*G_pt)[n_dof + n];
3904 }
3905
3906 // Now we should be able to work out the corrections
3907 double denom = dot_e * dot_h - dot_g * dot_f;
3908
3909 // Copy the previous residuals
3910 double R31 = result[3 * n_dof], R32 = result[3 * n_dof + 1];
3911 // Delta parameter
3912 const double delta_param =
3913 ((R32 - dot_d) * dot_g - (R31 - dot_c) * dot_h) / denom;
3914 // Delta frequency
3915 const double delta_w = -((R32 - dot_d) + dot_f * delta_param) / (dot_h);
3916
3917 // Load into the result vector
3918 result[3 * n_dof] = delta_param;
3919 result[3 * n_dof + 1] = delta_w;
3920
3921 // The corrections to the null vector
3922 for (unsigned n = 0; n < 2 * n_dof; n++)
3923 {
3924 result[n_dof + n] =
3925 y2[n] - (*E_pt)[n] * delta_param - (*G_pt)[n] * delta_w;
3926 }
3927
3928 // Finally add the corrections to the unknowns
3929 for (unsigned n = 0; n < n_dof; n++)
3930 {
3931 result[n] = y1[n] - (*A_pt)[n] * delta_param;
3932 }
3933
3934 // The sign of the jacobian is the previous signs multiplied by the
3935 // sign of the denominator
3936 problem_pt->sign_of_jacobian() =
3937 sign_of_jacobian * static_cast<int>(std::fabs(denom) / denom);
3938
3939 // Switch things to our full solver
3940 handler_pt->solve_full_system();
3941
3942 // If we are not storing things, clear the results
3943 if (!Enable_resolve)
3944 {
3945 // We no longer need to store the matrix
3947 delete A_pt;
3948 A_pt = 0;
3949 delete E_pt;
3950 E_pt = 0;
3951 delete G_pt;
3952 G_pt = 0;
3953 }
3954 // Otherwise, also store the problem pointer
3955 else
3956 {
3957 Problem_pt = problem_pt;
3958 }
3959 }
3960
3961 //=====================================================================
3962 /// Solve for two right hand sides, required for (efficient) continuation
3963 /// because otherwise we have to store the inverse of the jacobian
3964 /// and the complex equivalent ... nasty.
3965 //=====================================================================
3967 DoubleVector& result,
3968 const DoubleVector& rhs2,
3969 DoubleVector& result2)
3970 {
3971 // if the result is setup then it should not be distributed
3972#ifdef PARANOID
3973 if (result.built())
3974 {
3975 if (result.distributed())
3976 {
3977 throw OomphLibError("The result vector must not be distributed",
3978 OOMPH_CURRENT_FUNCTION,
3979 OOMPH_EXCEPTION_LOCATION);
3980 }
3981 }
3982 if (result2.built())
3983 {
3984 if (result2.distributed())
3985 {
3986 throw OomphLibError("The result2 vector must not be distributed",
3987 OOMPH_CURRENT_FUNCTION,
3988 OOMPH_EXCEPTION_LOCATION);
3989 }
3990 }
3991#endif
3992
3993 // Locally cache the pointer to the handler.
3994 HopfHandler* handler_pt =
3995 static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());
3996
3997 // Find the number of dofs in the augmented problem
3998 unsigned n_dof = problem_pt->ndof();
3999
4000 // create the linear algebra distribution for this solver
4001 // currently only global (non-distributed) distributions are allowed
4002 LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
4003 this->build_distribution(dist);
4004
4005 // if the result vector is not setup then rebuild with distribution = global
4006 if (!result.built())
4007 {
4008 result.build(this->distribution_pt(), 0.0);
4009 }
4010 if (!result2.built())
4011 {
4012 result2.build(this->distribution_pt(), 0.0);
4013 }
4014
4015 // Firstly, let's calculate the derivative of the residuals wrt
4016 // the parameter
4017 DoubleVector dRdparam(this->distribution_pt(), 0.0);
4018
4019 const double FD_step = 1.0e-8;
4020 {
4021 // Cache pointer to the parameter (second last entry in the vector
4022 double* param_pt = &problem_pt->dof(n_dof - 2);
4023 // Backup the parameter
4024 double old_var = *param_pt;
4025 ;
4026 // Increment the parameter
4027 *param_pt += FD_step;
4029
4030 // Now calculate the new residuals
4031 problem_pt->get_residuals(dRdparam);
4032
4033 // Now calculate the difference assume original residuals in resul
4034 for (unsigned n = 0; n < n_dof; n++)
4035 {
4036 dRdparam[n] = (dRdparam[n] - result[n]) / FD_step;
4037 }
4038
4039 // Reset the parameter
4040 *param_pt = old_var;
4042 }
4043
4044 // Switch the handler to "standard" mode
4045 handler_pt->solve_standard_system();
4046
4047 // Find out the number of dofs in the non-standard problem
4048 n_dof = problem_pt->ndof();
4049
4050 //
4051 dist.build(problem_pt->communicator_pt(), n_dof, false);
4052 this->build_distribution(dist);
4053
4054 // Setup storage for temporary vector
4055 DoubleVector y1(this->distribution_pt(), 0.0),
4056 alpha(this->distribution_pt(), 0.0),
4057 y1_resolve(this->distribution_pt(), 0.0);
4058
4059 // Allocate storage for A which can be used in the resolve
4060 if (A_pt != 0)
4061 {
4062 delete A_pt;
4063 }
4064 A_pt = new DoubleVector(this->distribution_pt(), 0.0);
4065
4066 // We are going to do resolves using the underlying linear solver
4068
4069 // Solve the first system Jy1 =
4070 Linear_solver_pt->solve(problem_pt, y1);
4071
4072 // This should have set the sign of the jacobian, store it
4073 int sign_of_jacobian = problem_pt->sign_of_jacobian();
4074
4075 // Now let's get the appropriate bit of alpha
4076 for (unsigned n = 0; n < n_dof; n++)
4077 {
4078 alpha[n] = dRdparam[n];
4079 }
4080
4081 // Resolve to find A
4082 Linear_solver_pt->resolve(alpha, *A_pt);
4083
4084 // Get the solution for the second rhs
4085 for (unsigned n = 0; n < n_dof; n++)
4086 {
4087 alpha[n] = rhs2[n];
4088 }
4089
4090 // Resolve to find y1_resolve
4091 Linear_solver_pt->resolve(alpha, y1_resolve);
4092
4093 // Now set to the complex system
4094 handler_pt->solve_complex_system();
4095
4096 // rebuild the Distribution
4097 dist.build(problem_pt->communicator_pt(), n_dof * 2, false);
4098 this->build_distribution(dist);
4099
4100 // Resize the stored vector G
4101 if (G_pt != 0)
4102 {
4103 delete G_pt;
4104 }
4105 G_pt = new DoubleVector(this->distribution_pt(), 0.0);
4106
4107 // Solve the first Zg = (Mz -My)
4108 Linear_solver_pt->solve(problem_pt, *G_pt);
4109
4110 // This should have set the sign of the complex matrix's determinant,
4111 // multiply
4112 sign_of_jacobian *= problem_pt->sign_of_jacobian();
4113
4114 // We can now construct our multipliers
4115 // Prepare to scale
4116 double dof_length = 0.0, a_length = 0.0, y1_length = 0.0,
4117 y1_resolve_length = 0.0;
4118 // Loop over the standard number of dofs
4119 for (unsigned n = 0; n < n_dof; n++)
4120 {
4121 if (std::fabs(problem_pt->dof(n)) > dof_length)
4122 {
4123 dof_length = std::fabs(problem_pt->dof(n));
4124 }
4125 if (std::fabs((*A_pt)[n]) > a_length)
4126 {
4127 a_length = std::fabs((*A_pt)[n]);
4128 }
4129 if (std::fabs(y1[n]) > y1_length)
4130 {
4131 y1_length = std::fabs(y1[n]);
4132 }
4133 if (std::fabs(y1_resolve[n]) > y1_resolve_length)
4134 {
4135 y1_resolve_length = std::fabs(y1[n]);
4136 }
4137 }
4138
4139
4140 double a_mult = dof_length / a_length;
4141 double y1_mult = dof_length / y1_length;
4142 double y1_resolve_mult = dof_length / y1_resolve_length;
4143 a_mult += FD_step;
4144 y1_mult += FD_step;
4145 y1_resolve_mult += FD_step;
4146 a_mult *= FD_step;
4147 y1_mult *= FD_step;
4148 y1_resolve_mult *= FD_step;
4149
4150 // Local storage for the product terms
4151 Vector<double> Jprod_a(2 * n_dof, 0.0), Jprod_y1(2 * n_dof, 0.0);
4152 Vector<double> Jprod_y1_resolve(2 * n_dof, 0.0);
4153
4154 // Temporary storage
4155 Vector<double> rhs(2 * n_dof);
4156
4157 // find the number of elements
4158 unsigned long n_element = problem_pt->mesh_pt()->nelement();
4159
4160 // Calculate the product of the jacobian matrices, etc
4161 for (unsigned long e = 0; e < n_element; e++)
4162 {
4163 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
4164 // Loop over the ndofs in each element
4165 unsigned n_var = elem_pt->ndof();
4166 // Get the jacobian matrices
4167 DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var),
4168 jac_y1_resolve(n_var);
4169 DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var),
4170 M_y1_resolve(n_var);
4171 // Get unperturbed jacobian
4172 elem_pt->get_jacobian_and_mass_matrix(rhs, jac, M);
4173
4174 // Backup the dofs
4175 Vector<double> dof_bac(n_var);
4176 // Perturb the original dofs
4177 for (unsigned n = 0; n < n_var; n++)
4178 {
4179 unsigned eqn_number = elem_pt->eqn_number(n);
4180 dof_bac[n] = problem_pt->dof(eqn_number);
4181 // Pertub by vector A
4182 problem_pt->dof(eqn_number) += a_mult * (*A_pt)[eqn_number];
4183 }
4184
4185 // Now get the new jacobian and mass matrix
4186 elem_pt->get_jacobian_and_mass_matrix(rhs, jac_a, M_a);
4187
4188 // Perturb the dofs
4189 for (unsigned n = 0; n < n_var; n++)
4190 {
4191 unsigned eqn_number = elem_pt->eqn_number(n);
4192 problem_pt->dof(eqn_number) = dof_bac[n];
4193 // Pertub by vector y1
4194 problem_pt->dof(eqn_number) += y1_mult * y1[eqn_number];
4195 }
4196
4197 // Now get the new jacobian and mass matrix
4198 elem_pt->get_jacobian_and_mass_matrix(rhs, jac_y1, M_y1);
4199
4200 // Perturb the dofs
4201 for (unsigned n = 0; n < n_var; n++)
4202 {
4203 unsigned eqn_number = elem_pt->eqn_number(n);
4204 problem_pt->dof(eqn_number) = dof_bac[n];
4205 // Pertub by vector y1
4206 problem_pt->dof(eqn_number) += y1_resolve_mult * y1_resolve[eqn_number];
4207 }
4208
4209 // Now get the new jacobian and mass matrix
4210 elem_pt->get_jacobian_and_mass_matrix(rhs, jac_y1_resolve, M_y1_resolve);
4211
4212 // Reset the dofs
4213 for (unsigned n = 0; n < n_var; n++)
4214 {
4215 unsigned eqn_number = elem_pt->eqn_number(n);
4216 problem_pt->dof(eqn_number) = dof_bac[n];
4217 }
4218
4219 // OK, now work out the products
4220 for (unsigned n = 0; n < n_var; n++)
4221 {
4222 unsigned eqn_number = elem_pt->eqn_number(n);
4223 double prod_a1 = 0.0, prod_y11 = 0.0, prod_y1_resolve1 = 0.0;
4224 double prod_a2 = 0.0, prod_y12 = 0.0, prod_y1_resolve2 = 0.0;
4225 for (unsigned m = 0; m < n_var; m++)
4226 {
4227 const unsigned unknown = elem_pt->eqn_number(m);
4228 const double y = handler_pt->Phi[unknown];
4229 const double z = handler_pt->Psi[unknown];
4230 const double omega = handler_pt->Omega;
4231 // Real part (first line)
4232 prod_a1 +=
4233 (jac_a(n, m) - jac(n, m)) * y + omega * (M_a(n, m) - M(n, m)) * z;
4234 prod_y11 +=
4235 (jac_y1(n, m) - jac(n, m)) * y + omega * (M_y1(n, m) - M(n, m)) * z;
4236 prod_y1_resolve1 += (jac_y1_resolve(n, m) - jac(n, m)) * y +
4237 omega * (M_y1_resolve(n, m) - M(n, m)) * z;
4238 // Imag par (second line)
4239 prod_a2 +=
4240 (jac_a(n, m) - jac(n, m)) * z - omega * (M_a(n, m) - M(n, m)) * y;
4241 prod_y12 +=
4242 (jac_y1(n, m) - jac(n, m)) * z - omega * (M_y1(n, m) - M(n, m)) * y;
4243 prod_y1_resolve2 += (jac_y1_resolve(n, m) - jac(n, m)) * z -
4244 omega * (M_y1_resolve(n, m) - M(n, m)) * y;
4245 }
4246 Jprod_a[eqn_number] += prod_a1 / a_mult;
4247 Jprod_y1[eqn_number] += prod_y11 / y1_mult;
4248 Jprod_y1_resolve[eqn_number] += prod_y1_resolve1 / y1_resolve_mult;
4249 Jprod_a[n_dof + eqn_number] += prod_a2 / a_mult;
4250 Jprod_y1[n_dof + eqn_number] += prod_y12 / y1_mult;
4251 Jprod_y1_resolve[n_dof + eqn_number] +=
4252 prod_y1_resolve2 / y1_resolve_mult;
4253 }
4254 }
4255
4256 // The assumption here is still that the result has been set to the
4257 // residuals.
4258 for (unsigned n = 0; n < 2 * n_dof; n++)
4259 {
4260 rhs[n] = result[n_dof + n] - Jprod_y1[n];
4261 }
4262
4263 // Temporary storage
4264 DoubleVector y2(this->distribution_pt(), 0.0);
4265
4266 // Solve it
4267 DoubleVector temp_rhs(this->distribution_pt(), 0.0);
4268 for (unsigned i = 0; i < 2 * n_dof; i++)
4269 {
4270 temp_rhs[i] = rhs[i];
4271 }
4272 Linear_solver_pt->resolve(temp_rhs, y2);
4273
4274 // Assemble the next RHS
4275 for (unsigned n = 0; n < 2 * n_dof; n++)
4276 {
4277 rhs[n] = dRdparam[n_dof + n] - Jprod_a[n];
4278 }
4279
4280 // Resive the storage
4281 if (E_pt != 0)
4282 {
4283 delete E_pt;
4284 }
4285 E_pt = new DoubleVector(this->distribution_pt(), 0.0);
4286
4287 // Solve for the next RHS
4288 for (unsigned i = 0; i < 2 * n_dof; i++)
4289 {
4290 temp_rhs[i] = rhs[i];
4291 }
4292 Linear_solver_pt->resolve(temp_rhs, *E_pt);
4293
4294 // Assemble the next RHS
4295 for (unsigned n = 0; n < 2 * n_dof; n++)
4296 {
4297 rhs[n] = rhs2[n_dof + n] - Jprod_y1_resolve[n];
4298 }
4299
4300 DoubleVector y2_resolve(this->distribution_pt(), 0.0);
4301 for (unsigned i = 0; i < 2 * n_dof; i++)
4302 {
4303 temp_rhs[i] = rhs[i];
4304 }
4305 Linear_solver_pt->resolve(temp_rhs, y2_resolve);
4306
4307
4308 // We can now calculate the final corrections
4309 // We need to work out a large number of dot products
4310 double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g = 0.0;
4311 double dot_h = 0.0;
4312
4313 double dot_c_resolve = 0.0, dot_d_resolve = 0.0;
4314
4315 for (unsigned n = 0; n < n_dof; n++)
4316 {
4317 // Get the appopriate entry
4318 const double Cn = handler_pt->C[n];
4319 dot_c += Cn * y2[n];
4320 dot_d += Cn * y2[n_dof + n];
4321 dot_c_resolve += Cn * y2_resolve[n];
4322 dot_d_resolve += Cn * y2_resolve[n_dof + n];
4323 dot_e += Cn * (*E_pt)[n];
4324 dot_f += Cn * (*E_pt)[n_dof + n];
4325 dot_g += Cn * (*G_pt)[n];
4326 dot_h += Cn * (*G_pt)[n_dof + n];
4327 }
4328
4329 // Now we should be able to work out the corrections
4330 double denom = dot_e * dot_h - dot_g * dot_f;
4331
4332 // Copy the previous residuals
4333 double R31 = result[3 * n_dof], R32 = result[3 * n_dof + 1];
4334 // Delta parameter
4335 const double delta_param =
4336 ((R32 - dot_d) * dot_g - (R31 - dot_c) * dot_h) / denom;
4337 // Delta frequency
4338 const double delta_w = -((R32 - dot_d) + dot_f * delta_param) / (dot_h);
4339
4340 // Corrections
4341 double R31_resolve = rhs2[3 * n_dof], R32_resolve = rhs2[3 * n_dof + 1];
4342 // Delta parameter
4343 const double delta_param_resolve = ((R32_resolve - dot_d_resolve) * dot_g -
4344 (R31_resolve - dot_c_resolve) * dot_h) /
4345 denom;
4346 // Delta frequency
4347 const double delta_w_resolve =
4348 -((R32_resolve - dot_d_resolve) + dot_f * delta_param_resolve) / (dot_h);
4349
4350
4351 // Load into the result vector
4352 result[3 * n_dof] = delta_param;
4353 result[3 * n_dof + 1] = delta_w;
4354
4355 // The corrections to the null vector
4356 for (unsigned n = 0; n < 2 * n_dof; n++)
4357 {
4358 result[n_dof + n] =
4359 y2[n] - (*E_pt)[n] * delta_param - (*G_pt)[n] * delta_w;
4360 }
4361
4362 // Finally add the corrections to the unknowns
4363 for (unsigned n = 0; n < n_dof; n++)
4364 {
4365 result[n] = y1[n] - (*A_pt)[n] * delta_param;
4366 }
4367
4368 // Load into the result vector
4369 result2[3 * n_dof] = delta_param_resolve;
4370 result2[3 * n_dof + 1] = delta_w_resolve;
4371
4372 // The corrections to the null vector
4373 for (unsigned n = 0; n < 2 * n_dof; n++)
4374 {
4375 result2[n_dof + n] = y2_resolve[n] - (*E_pt)[n] * delta_param_resolve -
4376 (*G_pt)[n] * delta_w_resolve;
4377 }
4378
4379 // Finally add the corrections to the unknowns
4380 for (unsigned n = 0; n < n_dof; n++)
4381 {
4382 result2[n] = y1_resolve[n] - (*A_pt)[n] * delta_param_resolve;
4383 }
4384
4385
4386 // The sign of the jacobian is the previous signs multiplied by the
4387 // sign of the denominator
4388 problem_pt->sign_of_jacobian() =
4389 sign_of_jacobian * static_cast<int>(std::fabs(denom) / denom);
4390
4391 // Switch things to our full solver
4392 handler_pt->solve_full_system();
4393
4394 // If we are not storing things, clear the results
4395 if (!Enable_resolve)
4396 {
4397 // We no longer need to store the matrix
4399 delete A_pt;
4400 A_pt = 0;
4401 delete E_pt;
4402 E_pt = 0;
4403 delete G_pt;
4404 G_pt = 0;
4405 }
4406 // Otherwise, also store the problem pointer
4407 else
4408 {
4409 Problem_pt = problem_pt;
4410 }
4411 }
4412
4413
4414 //======================================================================
4415 // Hack the re-solve to use the block factorisation
4416 //======================================================================
4418 DoubleVector& result)
4419 {
4420 throw OomphLibError("resolve() is not implemented for this solver",
4421 OOMPH_CURRENT_FUNCTION,
4422 OOMPH_EXCEPTION_LOCATION);
4423 }
4424
4425
4426 /// ////////////////////////////////////////////////////////////////////
4427 // Non-inline functions for the HopfHandler class
4428 /// ///////////////////////////////////////////////////////////////////
4429
4430 //====================================================================
4431 /// Constructor: Initialise the hopf handler,
4432 /// by setting initial guesses for Phi, Psi and calculating Count.
4433 /// If the system changes, a new handler must be constructed.
4434 //===================================================================
4436 double* const& parameter_pt)
4437 : Solve_which_system(0), Parameter_pt(parameter_pt), Omega(0.0)
4438 {
4439 // Set the problem pointer
4440 Problem_pt = problem_pt;
4441 // Set the number of non-augmented degrees of freedom
4442 Ndof = problem_pt->ndof();
4443
4444 // create the linear algebra distribution for this solver
4445 // currently only global (non-distributed) distributions are allowed
4446 LinearAlgebraDistribution* dist_pt =
4447 new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
4448
4449 // Resize the vectors of additional dofs
4450 Phi.resize(Ndof);
4451 Psi.resize(Ndof);
4452 C.resize(Ndof);
4453 Count.resize(Ndof, 0);
4454
4455 // Loop over all the elements in the problem
4456 unsigned n_element = problem_pt->mesh_pt()->nelement();
4457 for (unsigned e = 0; e < n_element; e++)
4458 {
4459 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
4460 // Loop over the local freedoms in an element
4461 unsigned n_var = elem_pt->ndof();
4462 for (unsigned n = 0; n < n_var; n++)
4463 {
4464 // Increase the associated global equation number counter
4465 ++Count[elem_pt->eqn_number(n)];
4466 }
4467 }
4468
4469 // Calculate the value Phi by
4470 // solving the system JPhi = dF/dlambda
4471
4472 // Locally cache the linear solver
4473 LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();
4474
4475 // Save the status before entry to this routine
4476 bool enable_resolve = linear_solver_pt->is_resolve_enabled();
4477
4478 // We need to do a resolve
4479 linear_solver_pt->enable_resolve();
4480
4481 // Storage for the solution
4482 DoubleVector x(dist_pt, 0.0);
4483
4484 // Solve the standard problem, we only want to make sure that
4485 // we factorise the matrix, if it has not been factorised. We shall
4486 // ignore the return value of x.
4487 linear_solver_pt->solve(problem_pt, x);
4488
4489 // Get the vector dresiduals/dparameter
4490 problem_pt->get_derivative_wrt_global_parameter(parameter_pt, x);
4491
4492 // Copy rhs vector into local storage so it doesn't get overwritten
4493 // if the linear solver decides to initialise the solution vector, say,
4494 // which it's quite entitled to do!
4495 DoubleVector input_x(x);
4496
4497 // Now resolve the system with the new RHS and overwrite the solution
4498 linear_solver_pt->resolve(input_x, x);
4499
4500 // Restore the storage status of the linear solver
4501 if (enable_resolve)
4502 {
4503 linear_solver_pt->enable_resolve();
4504 }
4505 else
4506 {
4507 linear_solver_pt->disable_resolve();
4508 }
4509
4510 // Normalise the solution x
4511 double length = 0.0;
4512 for (unsigned n = 0; n < Ndof; n++)
4513 {
4514 length += x[n] * x[n];
4515 }
4516 length = sqrt(length);
4517
4518 // Now add the real part of the null space components to the problem
4519 // unknowns and initialise it all
4520 // This is dumb at the moment ... fix with eigensolver?
4521 for (unsigned n = 0; n < Ndof; n++)
4522 {
4523 problem_pt->Dof_pt.push_back(&Phi[n]);
4524 C[n] = Phi[n] = -x[n] / length;
4525 }
4526
4527 // Set the imaginary part so that the appropriate residual is
4528 // zero initially (eigensolvers)
4529 for (unsigned n = 0; n < Ndof; n += 2)
4530 {
4531 // Make sure that we are not at the end of an array of odd length
4532 if (n != Ndof - 1)
4533 {
4534 Psi[n] = C[n + 1];
4535 Psi[n + 1] = -C[n];
4536 }
4537 // If it's odd set the final entry to zero
4538 else
4539 {
4540 Psi[n] = 0.0;
4541 }
4542 }
4543
4544 // Next add the imaginary parts of the null space components to the problem
4545 for (unsigned n = 0; n < Ndof; n++)
4546 {
4547 problem_pt->Dof_pt.push_back(&Psi[n]);
4548 }
4549 // Now add the parameter
4550 problem_pt->Dof_pt.push_back(parameter_pt);
4551 // Finally add the frequency
4552 problem_pt->Dof_pt.push_back(&Omega);
4553
4554 // rebuild the dof dist
4556 Problem_pt->communicator_pt(), Ndof * 3 + 2, false);
4557 // Remove all previous sparse storage used during Jacobian assembly
4559
4560 // delete the dist_pt
4561 delete dist_pt;
4562 }
4563
4564 //====================================================================
4565 /// Constructor: Initialise the hopf handler,
4566 /// by setting initial guesses for Phi, Psi, Omega and calculating Count.
4567 /// If the system changes, a new handler must be constructed.
4568 //===================================================================
4570 double* const& parameter_pt,
4571 const double& omega,
4572 const DoubleVector& phi,
4573 const DoubleVector& psi)
4574 : Solve_which_system(0), Parameter_pt(parameter_pt), Omega(omega)
4575 {
4576 // Set the problem pointer
4577 Problem_pt = problem_pt;
4578 // Set the number of non-augmented degrees of freedom
4579 Ndof = problem_pt->ndof();
4580
4581 // Resize the vectors of additional dofs
4582 Phi.resize(Ndof);
4583 Psi.resize(Ndof);
4584 C.resize(Ndof);
4585 Count.resize(Ndof, 0);
4586
4587 // Loop over all the elements in the problem
4588 unsigned n_element = problem_pt->mesh_pt()->nelement();
4589 for (unsigned e = 0; e < n_element; e++)
4590 {
4591 GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
4592 // Loop over the local freedoms in an element
4593 unsigned n_var = elem_pt->ndof();
4594 for (unsigned n = 0; n < n_var; n++)
4595 {
4596 // Increase the associated global equation number counter
4597 ++Count[elem_pt->eqn_number(n)];
4598 }
4599 }
4600
4601 // Normalise the guess for phi
4602 double length = 0.0;
4603 for (unsigned n = 0; n < Ndof; n++)
4604 {
4605 length += phi[n] * phi[n];
4606 }
4607 length = sqrt(length);
4608
4609 // Now add the real part of the null space components to the problem
4610 // unknowns and initialise it all
4611 for (unsigned n = 0; n < Ndof; n++)
4612 {
4613 problem_pt->Dof_pt.push_back(&Phi[n]);
4614 C[n] = Phi[n] = phi[n] / length;
4615 Psi[n] = psi[n] / length;
4616 }
4617
4618 // Next add the imaginary parts of the null space components to the problem
4619 for (unsigned n = 0; n < Ndof; n++)
4620 {
4621 problem_pt->Dof_pt.push_back(&Psi[n]);
4622 }
4623
4624 // Now add the parameter
4625 problem_pt->Dof_pt.push_back(parameter_pt);
4626 // Finally add the frequency
4627 problem_pt->Dof_pt.push_back(&Omega);
4628
4629 // rebuild the Dof_distribution_pt
4631 Problem_pt->communicator_pt(), Ndof * 3 + 2, false);
4632 // Remove all previous sparse storage used during Jacobian assembly
4634 }
4635
4636
4637 //=======================================================================
4638 /// Destructor return the problem to its original (non-augmented) state
4639 //=======================================================================
4641 {
4642 // If we are using the block solver reset the problem's linear solver
4643 // to the original one
4644 BlockHopfLinearSolver* block_hopf_solver_pt =
4646 if (block_hopf_solver_pt)
4647 {
4648 // Reset the problem's linear solver
4649 Problem_pt->linear_solver_pt() = block_hopf_solver_pt->linear_solver_pt();
4650 // Delete the block solver
4651 delete block_hopf_solver_pt;
4652 }
4653 // Now return the problem to its original size
4654 Problem_pt->Dof_pt.resize(Ndof);
4656 Problem_pt->communicator_pt(), Ndof, false);
4657 // Remove all previous sparse storage used during Jacobian assembly
4659 }
4660
4661
4662 //=============================================================
4663 /// Get the number of elemental degrees of freedom
4664 //=============================================================
4665 unsigned HopfHandler::ndof(GeneralisedElement* const& elem_pt)
4666 {
4667 unsigned raw_ndof = elem_pt->ndof();
4668 switch (Solve_which_system)
4669 {
4670 // Full augmented system
4671 case 0:
4672 return (3 * raw_ndof + 2);
4673 break;
4674 // Standard non-augmented system
4675 case 1:
4676 return raw_ndof;
4677 break;
4678 // Complex system
4679 case 2:
4680 return (2 * raw_ndof);
4681 break;
4682
4683 default:
4684 throw OomphLibError("Solve_which_system can only be 0,1 or 2",
4685 OOMPH_CURRENT_FUNCTION,
4686 OOMPH_EXCEPTION_LOCATION);
4687 }
4688 }
4689
4690 //=============================================================
4691 /// Get the global equation number of the local unknown
4692 //============================================================
4693 unsigned long HopfHandler::eqn_number(GeneralisedElement* const& elem_pt,
4694 const unsigned& ieqn_local)
4695 {
4696 // Get the raw value
4697 unsigned raw_ndof = elem_pt->ndof();
4698 unsigned long global_eqn;
4699 if (ieqn_local < raw_ndof)
4700 {
4701 global_eqn = elem_pt->eqn_number(ieqn_local);
4702 }
4703 else if (ieqn_local < 2 * raw_ndof)
4704 {
4705 global_eqn = Ndof + elem_pt->eqn_number(ieqn_local - raw_ndof);
4706 }
4707 else if (ieqn_local < 3 * raw_ndof)
4708 {
4709 global_eqn = 2 * Ndof + elem_pt->eqn_number(ieqn_local - 2 * raw_ndof);
4710 }
4711 else if (ieqn_local == 3 * raw_ndof)
4712 {
4713 global_eqn = 3 * Ndof;
4714 }
4715 else
4716 {
4717 global_eqn = 3 * Ndof + 1;
4718 }
4719 return global_eqn;
4720 }
4721
4722 //==================================================================
4723 /// Get the residuals
4724 //=================================================================
4726 Vector<double>& residuals)
4727 {
4728 // Should only call get residuals for the full system
4729 if (Solve_which_system == 0)
4730 {
4731 // Need to get raw residuals and jacobian
4732 unsigned raw_ndof = elem_pt->ndof();
4733
4734 DenseMatrix<double> jacobian(raw_ndof), M(raw_ndof);
4735 // Get the basic residuals, jacobian and mass matrix
4736 elem_pt->get_jacobian_and_mass_matrix(residuals, jacobian, M);
4737
4738 // Initialise the pen-ultimate residual
4739 residuals[3 * raw_ndof] =
4740 -1.0 / (double)(Problem_pt->mesh_pt()->nelement());
4741 residuals[3 * raw_ndof + 1] = 0.0;
4742
4743 // Now multiply to fill in the residuals
4744 for (unsigned i = 0; i < raw_ndof; i++)
4745 {
4746 residuals[raw_ndof + i] = 0.0;
4747 residuals[2 * raw_ndof + i] = 0.0;
4748 for (unsigned j = 0; j < raw_ndof; j++)
4749 {
4750 unsigned global_unknown = elem_pt->eqn_number(j);
4751 // Real part
4752 residuals[raw_ndof + i] += jacobian(i, j) * Phi[global_unknown] +
4753 Omega * M(i, j) * Psi[global_unknown];
4754 // Imaginary part
4755 residuals[2 * raw_ndof + i] += jacobian(i, j) * Psi[global_unknown] -
4756 Omega * M(i, j) * Phi[global_unknown];
4757 }
4758 // Get the global equation number
4759 unsigned global_eqn = elem_pt->eqn_number(i);
4760
4761 // Real part
4762 residuals[3 * raw_ndof] +=
4763 (Phi[global_eqn] * C[global_eqn]) / Count[global_eqn];
4764 // Imaginary part
4765 residuals[3 * raw_ndof + 1] +=
4766 (Psi[global_eqn] * C[global_eqn]) / Count[global_eqn];
4767 }
4768 }
4769 else
4770 {
4771 throw OomphLibError("Solve_which_system can only be 0",
4772 OOMPH_CURRENT_FUNCTION,
4773 OOMPH_EXCEPTION_LOCATION);
4774 }
4775 }
4776
4777
4778 //===============================================================
4779 /// Calculate the elemental Jacobian matrix "d equation
4780 /// / d variable".
4781 //==================================================================
4783 Vector<double>& residuals,
4784 DenseMatrix<double>& jacobian)
4785 {
4786 // The standard case
4787 if (Solve_which_system == 0)
4788 {
4789 unsigned augmented_ndof = ndof(elem_pt);
4790 unsigned raw_ndof = elem_pt->ndof();
4791
4792 // Get the basic residuals and jacobian
4793 DenseMatrix<double> M(raw_ndof);
4794 elem_pt->get_jacobian_and_mass_matrix(residuals, jacobian, M);
4795 // Now fill in the actual residuals
4796 get_residuals(elem_pt, residuals);
4797
4798 // Now the jacobian appears in other entries
4799 for (unsigned n = 0; n < raw_ndof; ++n)
4800 {
4801 for (unsigned m = 0; m < raw_ndof; ++m)
4802 {
4803 jacobian(raw_ndof + n, raw_ndof + m) = jacobian(n, m);
4804 jacobian(raw_ndof + n, 2 * raw_ndof + m) = Omega * M(n, m);
4805 jacobian(2 * raw_ndof + n, 2 * raw_ndof + m) = jacobian(n, m);
4806 jacobian(2 * raw_ndof + n, raw_ndof + m) = -Omega * M(n, m);
4807 unsigned global_eqn = elem_pt->eqn_number(m);
4808 jacobian(raw_ndof + n, 3 * raw_ndof + 1) += M(n, m) * Psi[global_eqn];
4809 jacobian(2 * raw_ndof + n, 3 * raw_ndof + 1) -=
4810 M(n, m) * Phi[global_eqn];
4811 }
4812
4813 unsigned local_eqn = elem_pt->eqn_number(n);
4814 jacobian(3 * raw_ndof, raw_ndof + n) = C[local_eqn] / Count[local_eqn];
4815 jacobian(3 * raw_ndof + 1, 2 * raw_ndof + n) =
4816 C[local_eqn] / Count[local_eqn];
4817 }
4818
4819 const double FD_step = 1.0e-8;
4820
4821 Vector<double> newres_p(augmented_ndof), newres_m(augmented_ndof);
4822
4823 // Loop over the dofs
4824 for (unsigned n = 0; n < raw_ndof; n++)
4825 {
4826 // Just do the x's
4827 unsigned long global_eqn = eqn_number(elem_pt, n);
4828 double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
4829 double init = *unknown_pt;
4830 *unknown_pt += FD_step;
4831
4832 // Get the new residuals
4833 get_residuals(elem_pt, newres_p);
4834
4835 // Reset
4836 *unknown_pt = init;
4837
4838 // Subtract
4839 *unknown_pt -= FD_step;
4840 get_residuals(elem_pt, newres_m);
4841
4842 for (unsigned m = 0; m < raw_ndof; m++)
4843 {
4844 jacobian(raw_ndof + m, n) =
4845 (newres_p[raw_ndof + m] - residuals[raw_ndof + m]) / (FD_step);
4846 jacobian(2 * raw_ndof + m, n) =
4847 (newres_p[2 * raw_ndof + m] - residuals[2 * raw_ndof + m]) /
4848 (FD_step);
4849 }
4850 // Reset the unknown
4851 *unknown_pt = init;
4852 }
4853
4854 {
4855 // Now do the global parameter
4856 double* unknown_pt = Problem_pt->Dof_pt[3 * Ndof];
4857 double init = *unknown_pt;
4858 *unknown_pt += FD_step;
4859
4861 // Get the new residuals
4862 get_residuals(elem_pt, newres_p);
4863
4864 // Reset
4865 *unknown_pt = init;
4866
4867 // Subtract
4868 *unknown_pt -= FD_step;
4869 get_residuals(elem_pt, newres_m);
4870
4871 for (unsigned m = 0; m < augmented_ndof - 2; m++)
4872 {
4873 jacobian(m, 3 * raw_ndof) = (newres_p[m] - residuals[m]) / FD_step;
4874 }
4875 // Reset the unknown
4876 *unknown_pt = init;
4878 }
4879 } // End of standard case
4880 // Normal case
4881 else if (Solve_which_system == 1)
4882 {
4883 // Just get the normal jacobian and residuals
4884 elem_pt->get_jacobian(residuals, jacobian);
4885 }
4886 // Otherwise the augmented complex case
4887 else if (Solve_which_system == 2)
4888 {
4889 unsigned raw_ndof = elem_pt->ndof();
4890
4891 // Get the basic residuals and jacobian
4892 DenseMatrix<double> M(raw_ndof);
4893 elem_pt->get_jacobian_and_mass_matrix(residuals, jacobian, M);
4894
4895 // We now need to fill in the other blocks
4896 for (unsigned n = 0; n < raw_ndof; n++)
4897 {
4898 for (unsigned m = 0; m < raw_ndof; m++)
4899 {
4900 jacobian(n, raw_ndof + m) = Omega * M(n, m);
4901 jacobian(raw_ndof + n, m) = -Omega * M(n, m);
4902 jacobian(raw_ndof + n, raw_ndof + m) = jacobian(n, m);
4903 }
4904 }
4905
4906 // Now overwrite to fill in the residuals
4907 // The decision take is to solve for the mass matrix multiplied
4908 // terms in the residuals because they require no additional
4909 // information to assemble.
4910 for (unsigned n = 0; n < raw_ndof; n++)
4911 {
4912 residuals[n] = 0.0;
4913 residuals[raw_ndof + n] = 0.0;
4914 for (unsigned m = 0; m < raw_ndof; m++)
4915 {
4916 unsigned global_unknown = elem_pt->eqn_number(m);
4917 // Real part
4918 residuals[n] += M(n, m) * Psi[global_unknown];
4919 // Imaginary part
4920 residuals[raw_ndof + n] -= M(n, m) * Phi[global_unknown];
4921 }
4922 }
4923 } // End of complex augmented case
4924 else
4925 {
4926 throw OomphLibError("Solve_which_system can only be 0,1 or 2",
4927 OOMPH_CURRENT_FUNCTION,
4928 OOMPH_EXCEPTION_LOCATION);
4929 }
4930 }
4931
4932
4933 //==================================================================
4934 /// Get the derivatives of the augmented residuals with respect to
4935 /// a parameter
4936 //=================================================================
4938 GeneralisedElement* const& elem_pt,
4939 double* const& parameter_pt,
4940 Vector<double>& dres_dparam)
4941 {
4942 // Should only call get residuals for the full system
4943 if (Solve_which_system == 0)
4944 {
4945 // Need to get raw residuals and jacobian
4946 unsigned raw_ndof = elem_pt->ndof();
4947
4948 DenseMatrix<double> djac_dparam(raw_ndof), dM_dparam(raw_ndof);
4949 // Get the basic residuals, jacobian and mass matrix
4951 parameter_pt, dres_dparam, djac_dparam, dM_dparam);
4952
4953 // Initialise the pen-ultimate residual, which does not
4954 // depend on the parameter
4955 dres_dparam[3 * raw_ndof] = 0.0;
4956 dres_dparam[3 * raw_ndof + 1] = 0.0;
4957
4958 // Now multiply to fill in the residuals
4959 for (unsigned i = 0; i < raw_ndof; i++)
4960 {
4961 dres_dparam[raw_ndof + i] = 0.0;
4962 dres_dparam[2 * raw_ndof + i] = 0.0;
4963 for (unsigned j = 0; j < raw_ndof; j++)
4964 {
4965 unsigned global_unknown = elem_pt->eqn_number(j);
4966 // Real part
4967 dres_dparam[raw_ndof + i] +=
4968 djac_dparam(i, j) * Phi[global_unknown] +
4969 Omega * dM_dparam(i, j) * Psi[global_unknown];
4970 // Imaginary part
4971 dres_dparam[2 * raw_ndof + i] +=
4972 djac_dparam(i, j) * Psi[global_unknown] -
4973 Omega * dM_dparam(i, j) * Phi[global_unknown];
4974 }
4975 }
4976 }
4977 else
4978 {
4979 throw OomphLibError("Solve_which_system can only be 0",
4980 OOMPH_CURRENT_FUNCTION,
4981 OOMPH_EXCEPTION_LOCATION);
4982 }
4983 }
4984
4985
4986 //========================================================================
4987 /// Overload the derivative of the residuals and jacobian
4988 /// with respect to a parameter so that it breaks because it should not
4989 /// be required
4990 //========================================================================
4992 double* const& parameter_pt,
4993 Vector<double>& dres_dparam,
4994 DenseMatrix<double>& djac_dparam)
4995 {
4996 std::ostringstream error_stream;
4997 error_stream
4998 << "This function has not been implemented because it is not required\n";
4999 error_stream << "in standard problems.\n";
5000 error_stream
5001 << "If you find that you need it, you will have to implement it!\n\n";
5002
5003 throw OomphLibError(
5004 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5005 }
5006
5007
5008 //=====================================================================
5009 /// Overload the hessian vector product function so that
5010 /// it breaks because it should not be required
5011 //========================================================================
5013 GeneralisedElement* const& elem_pt,
5014 Vector<double> const& Y,
5015 DenseMatrix<double> const& C,
5016 DenseMatrix<double>& product)
5017 {
5018 std::ostringstream error_stream;
5019 error_stream
5020 << "This function has not been implemented because it is not required\n";
5021 error_stream << "in standard problems.\n";
5022 error_stream
5023 << "If you find that you need it, you will have to implement it!\n\n";
5024
5025 throw OomphLibError(
5026 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5027 }
5028
5029
5030 //==========================================================================
5031 /// Return the eigenfunction(s) associated with the bifurcation that
5032 /// has been detected in bifurcation tracking problems
5033 //==========================================================================
5035 {
5036 // There is a real and imaginary part of the null vector
5037 eigenfunction.resize(2);
5039 // Rebuild the vector
5040 eigenfunction[0].build(&dist, 0.0);
5041 eigenfunction[1].build(&dist, 0.0);
5042 // Set the value to be the null vector
5043 for (unsigned n = 0; n < Ndof; n++)
5044 {
5045 eigenfunction[0][n] = Phi[n];
5046 eigenfunction[1][n] = Psi[n];
5047 }
5048 }
5049
5050
5051 //====================================================================
5052 /// Set to solve the standard (underlying jacobian) system
5053 //===================================================================
5055 {
5056 if (Solve_which_system != 1)
5057 {
5059 // Restrict the problem to the standard variables only
5060 Problem_pt->Dof_pt.resize(Ndof);
5062 Problem_pt->communicator_pt(), Ndof, false);
5063 // Remove all previous sparse storage used during Jacobian assembly
5065 }
5066 }
5067
5068 //====================================================================
5069 /// Set to solve the complex (jacobian and mass matrix) system
5070 //===================================================================
5072 {
5073 // If we were not solving the complex system resize the unknowns
5074 // accordingly
5075 if (Solve_which_system != 2)
5076 {
5078 // Resize to the first Ndofs (will work whichever system we were
5079 // solving before)
5080 Problem_pt->Dof_pt.resize(Ndof);
5081 // Add the first (real) part of the eigenfunction back into the problem
5082 for (unsigned n = 0; n < Ndof; n++)
5083 {
5084 Problem_pt->Dof_pt.push_back(&Phi[n]);
5085 }
5087 Problem_pt->communicator_pt(), Ndof * 2, false);
5088 // Remove all previous sparse storage used during Jacobian assembly
5090 }
5091 }
5092
5093
5094 //=================================================================
5095 /// Set to Solve full system system
5096 //=================================================================
5098 {
5099 // If we are starting from another system
5101 {
5103
5104 // Resize to the first Ndofs (will work whichever system we were
5105 // solving before)
5106 Problem_pt->Dof_pt.resize(Ndof);
5107 // Add the additional unknowns back into the problem
5108 for (unsigned n = 0; n < Ndof; n++)
5109 {
5110 Problem_pt->Dof_pt.push_back(&Phi[n]);
5111 }
5112 for (unsigned n = 0; n < Ndof; n++)
5113 {
5114 Problem_pt->Dof_pt.push_back(&Psi[n]);
5115 }
5116 // Now add the parameter
5117 Problem_pt->Dof_pt.push_back(Parameter_pt);
5118 // Finally add the frequency
5119 Problem_pt->Dof_pt.push_back(&Omega);
5120
5121 //
5123 Problem_pt->communicator_pt(), 3 * Ndof + 2, false);
5124 // Remove all previous sparse storage used during Jacobian assembly
5126 }
5127 }
5128} // namespace oomph
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A class that is used to define the functions used to assemble the elemental contributions to the resi...
virtual unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
virtual void get_inner_product_vectors(GeneralisedElement *const &elem_pt, Vector< unsigned > const &history_index, Vector< Vector< double > > &inner_product_vector)
Compute the vectors that when taken as a dot product with other history values give the inner product...
virtual void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Calculate the product of the Hessian (derivative of Jacobian with respect to all variables) an eigenv...
virtual void get_dresiduals_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
Calculate the derivative of the residuals with respect to a parameter.
virtual double * bifurcation_parameter_pt() const
Return a pointer to the bifurcation parameter in bifurcation tracking problems.
virtual void get_inner_products(GeneralisedElement *const &elem_pt, Vector< std::pair< unsigned, unsigned > > const &history_index, Vector< double > &inner_product)
Compute the inner products of the given vector of pairs of history values over the element.
virtual void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
virtual void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt.
virtual unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Return the global equation number of the local unknown ieqn_local in elem_pt.
virtual void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt.
virtual double & local_problem_dof(Problem *const &problem_pt, const unsigned &t, const unsigned &i)
Return the t-th level of storage associated with the i-th (local) dof stored in the problem.
virtual void dof_vector(GeneralisedElement *const &elem_pt, const unsigned &t, Vector< double > &dof)
Return vector of dofs at time level t in the element elem_pt.
virtual void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable" for elem_pt.
virtual void get_djacobian_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Calculate the derivative of the residuals and jacobian with respect to a parameter.
virtual void dof_pt_vector(GeneralisedElement *const &elem_pt, Vector< double * > &dof_pt)
Return vector of pointers to dofs in the element elem_pt.
A custom linear solver class that is used to solve a block-factorised version of the Fold bifurcation...
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
DoubleVector * E_pt
Pointer to the storage for the vector e.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
~AugmentedBlockFoldLinearSolver()
Destructor: clean up the allocated memory.
DoubleVector * Alpha_pt
Pointer to the storage for the vector alpha.
A custom linear solver class that is used to solve a block-factorised version of the PitchFork bifurc...
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
DoubleVector * Alpha_pt
Pointer to the storage for the vector alpha.
DoubleVector * E_pt
Pointer to the storage for the vector e.
~AugmentedBlockPitchForkLinearSolver()
Destructor: clean up the allocated memory.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
A custom linear solver class that is used to solve a block-factorised version of the Hopf bifurcation...
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
DoubleVector * A_pt
Pointer to the storage for the vector a.
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
~BlockHopfLinearSolver()
Destructor: clean up the allocated memory.
void solve_for_two_rhs(Problem *const &problem_pt, DoubleVector &result, const DoubleVector &rhs2, DoubleVector &result2)
Solve for two right hand sides.
DoubleVector * E_pt
Pointer to the storage for the vector e (0 to n-1)
DoubleVector * G_pt
Pointer to the storage for the vector g (0 to n-1)
A custom linear solver class that is used to solve a block-factorised version of the PitchFork bifurc...
DoubleVector * D_pt
Pointer to the storage for the vector d.
LinearSolver * linear_solver_pt() const
Access function to the original linear solver.
DoubleVector * dJy_dparam_pt
Pointer to the storage for the vector of derivatives with respect to the bifurcation parameter.
~BlockPitchForkLinearSolver()
Destructor: clean up the allocated memory.
DoubleVector * B_pt
Pointer to the storage for the vector b.
Problem * Problem_pt
Pointer to the problem, used in the resolve.
void resolve(const DoubleVector &rhs, DoubleVector &result)
The resolve function also uses the block factorisation.
DoubleVector * C_pt
Pointer to the storage for the vector c.
void solve(Problem *const &problem_pt, DoubleVector &result)
The solve function uses the block factorisation.
LinearSolver * Linear_solver_pt
Pointer to the original linear solver.
bool distributed() const
distribution is serial or distributed
unsigned nrow_local() const
access function for the num of local rows on this processor.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
void build_halo_scheme(DoubleVectorHaloScheme *const &halo_scheme_pt)
Construct the halo scheme and storage for the halo data.
void synchronise()
Synchronise the halo data.
void sum_all_halo_and_haloed_values()
Sum all the data, store in the master (haloed) data and then synchronise.
double & global_value(const unsigned &i)
Direct access to global entry.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
double norm() const
compute the 2 norm of this vector
bool built() const
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt This is deliberately broken in our ei...
double Sigma_real
Storage for the real shift.
void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt.
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable" for elem_pt. Again deliberately bro...
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Return the global equation number of the local unknown ieqn_local in elem_pt.
unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Return the global equation number of the local unknown ieqn_local in elem_pt.
unsigned ndof(GeneralisedElement *const &elem_pt)
Return the number of degrees of freedom in the element elem_pt.
void get_all_vectors_and_matrices(GeneralisedElement *const &elem_pt, Vector< Vector< double > > &vec, Vector< DenseMatrix< double > > &matrix)
Calculate all desired vectors and matrices provided by the element elem_pt.
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable" for elem_pt. Again deliberately bro...
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt This is deliberately broken in our ei...
A class that is used to assemble the augmented system that defines a fold (saddle-node) or limit poin...
Problem * Problem_pt
Pointer to the problem.
Vector< int > Count
A vector that is used to determine how many elements contribute to a particular equation....
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
void solve_full_system()
Solve non-block system.
void get_djacobian_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities....
void solve_augmented_block_system()
Set to solve the augmented block system.
Vector< double > Phi
A constant vector used to ensure that the null vector is not trivial.
unsigned Ndof
Store the number of degrees of freedom in the non-augmented problem.
void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Overload the hessian vector product function so that it breaks.
~FoldHandler()
Destructor, return the problem to its original state before the augmented system was added.
void solve_block_system()
Set to solve the block system.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Get the residuals.
Vector< double > Y
Storage for the null vector.
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
void get_dresiduals_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
Overload the derivatives of the residuals with respect to a parameter to apply to the augmented syste...
FoldHandler(Problem *const &problem_pt, double *const &parameter_pt)
Constructor: initialise the fold handler, by setting initial guesses for Y, Phi and calculating count...
double * Parameter_pt
Storage for the pointer to the parameter.
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
A Generalised Element class.
Definition: elements.h:73
bool is_halo() const
Is this element a halo?
Definition: elements.h:1163
void dof_vector(const unsigned &t, Vector< double > &dof)
Return the vector of dof values at time level t.
Definition: elements.h:841
virtual void get_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Calculate the product of the Hessian (derivative of Jacobian with respect to all variables) an eigenv...
Definition: elements.h:1083
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
virtual void get_inner_product_vectors(Vector< unsigned > const &history_index, Vector< Vector< double > > &inner_product_vector)
Compute the vectors that when taken as a dot product with other history values give the inner product...
Definition: elements.h:1106
virtual void get_inner_products(Vector< std::pair< unsigned, unsigned > > const &history_index, Vector< double > &inner_product)
Return the vector of inner product of the given pairs of history values.
Definition: elements.h:1095
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
Definition: elements.h:980
virtual void get_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &mass_matrix)
Calculate the residuals and the elemental "mass" matrix, the matrix that multiplies the time derivati...
Definition: elements.h:1003
virtual void get_dresiduals_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam)
Calculate the derivatives of the residuals with respect to a parameter.
Definition: elements.h:1034
void dof_pt_vector(Vector< double * > &dof_pt)
Return the vector of pointers to dof values.
Definition: elements.h:866
virtual void get_djacobian_and_dmass_matrix_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
Calculate the derivatives of the elemental Jacobian matrix mass matrix and residuals with respect to ...
Definition: elements.h:1061
virtual void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
Definition: elements.h:990
virtual void get_djacobian_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Calculate the derivatives of the elemental Jacobian matrix and residuals with respect to a parameter.
Definition: elements.h:1046
virtual void get_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Calculate the residuals and jacobian and elemental "mass" matrix, the matrix that multiplies the time...
Definition: elements.h:1016
A class that is used to assemble the augmented system that defines a Hopf bifurcation....
Vector< double > Phi
The real part of the null vector.
Vector< int > Count
A vector that is used to determine how many elements contribute to a particular equation....
double Omega
The critical frequency of the bifurcation.
void get_dresiduals_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
Overload the derivatives of the residuals with respect to a parameter to apply to the augmented syste...
~HopfHandler()
Destructor, return the problem to its original state, before the augmented system was added.
void solve_full_system()
Solve non-block system.
void get_djacobian_dparameter(GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Get the residuals.
void get_hessian_vector_products(GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Overload the hessian vector product function so that it breaks.
unsigned ndof(GeneralisedElement *const &elem_pt)
Get the number of elemental degrees of freedom.
Vector< double > C
A constant vector used to ensure that the null vector is not trivial.
void get_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tra...
Vector< double > Psi
The imaginary part of the null vector.
unsigned long eqn_number(GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
Get the global equation number of the local unknown.
unsigned Solve_which_system
Integer flag to indicate which system should be assembled. There are three possibilities....