problem.h
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// A generic problem class to keep MH happy
27
28// Include guards to prevent multiple inclusion of this header
29#ifndef OOMPH_PROBLEM_CLASS_HEADER
30#define OOMPH_PROBLEM_CLASS_HEADER
31
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38#ifdef OOMPH_HAS_MPI
39#include "mpi.h"
40#endif
41
42// OOMPH-LIB headers
43#include "Vector.h"
44#include "matrices.h"
48#include <complex>
49#include <map>
50
51
52namespace oomph
53{
54 // Forward definition for Data class
55 class Data;
56
57 // Forward definition for Time class
58 class Time;
59
60 // Forward definition for TimeStepper class
61 class TimeStepper;
62
63 // Forward definition for Mesh class
64 class Mesh;
65
66 // Forward definition for RefineableElement class
67 class RefineableElement;
68
69 // Forward definition for PRefineableElement class
70 class PRefineableElement;
71
72 // Forward definition for FiniteElement class
73 class FiniteElement;
74
75 // Forward definition for the GeneralisedElement class
76 class GeneralisedElement;
77
78 // Forward definition for the Linearsolver class
79 class LinearSolver;
80
81 // Forward definition for the Eigensolver class
82 class EigenSolver;
83
84 // Forward definition for the assembly handler
85 class AssemblyHandler;
86
87 // Forward definition for sum of matrices class
88 class SumOfMatrices;
89
90 /// //////////////////////////////////////////////////////////////////
91 /// //////////////////////////////////////////////////////////////////
92 /// //////////////////////////////////////////////////////////////////
93
94
95 //=======================================================================
96 /// The Problem class
97 ///
98 /// The main components of a Problem are:
99 /// - a pointer to the (global) Mesh (which provides ordered
100 /// access to the nodes, elements, etc of all submeshes in the problem --
101 /// if there are any)
102 /// - pointers to the submeshes (if there are any)
103 /// - pointers to any global Data
104 /// - a pointer to the global Time
105 /// - pointers to the Timestepper used in the problem
106 /// - a pointer to the linear solver that will be used in Newton's method
107 /// - pointers to all DOFs in the problem.
108 ///
109 /// Obviously, at this level in the code hierarchy, many things become
110 /// very problem-dependent but when setting up a specific problem
111 /// (as a class that inherits from Problem), the problem constructor
112 /// will/should typically have a structure similar to this:
113 /// \code
114 /// // Set up a timestepper
115 /// TimeStepper* time_stepper_pt=new Steady;
116 ///
117 /// // Add timestepper to problem
118 /// add_time_stepper_pt(time_stepper_pt);
119 ///
120 /// // Build and assign mesh
121 /// Problem::mesh_pt() = new
122 /// SomeMesh<ELEMENT_TYPE>(Nelement,time_stepper_pt);
123 ///
124 /// // Set the boundary conditions for this problem: All nodes are
125 /// // free by default -- just pin the ones that have Dirichlet conditions
126 /// // here (boundary 0 in this example)
127 /// unsigned i_bound = 0
128 /// unsigned num_nod= mesh_pt()->nboundary_node(ibound);
129 /// for (unsigned inod=0;inod<num_nod;inod++)
130 /// {
131 /// mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
132 /// }
133 ///
134 /// // Complete the build of all elements so they are fully functional
135 ///
136 /// // Setup equation numbering scheme
137 /// oomph_info <<"Number of equations: " << assign_eqn_numbers() <<
138 /// std::endl;
139 /// \endcode
140 /// For time-dependent problems, we can then use
141 /// \code assign_initial_values_impulsive (); \endcode
142 /// to generate an initial guess (for an impulsive start) and
143 /// the problem can then be solved, e.g. using the steady
144 /// or unsteady Newton solvers
145 /// \code newton_solve() \endcode
146 /// or
147 /// \code unsteady_newton_solve(...) \endcode
148 ///
149 //=======================================================================
151 {
152 // The handler classes need to be friend
153 friend class FoldHandler;
154 friend class PitchForkHandler;
155 friend class HopfHandler;
156 template<unsigned NNODE_1D>
163
164
165 private:
166 /// The mesh pointer
168
169 /// Vector of pointers to submeshes
171
172 /// Pointer to the linear solver for the problem
174
175 /// Pointer to the linear solver used for explicit time steps (this is
176 /// likely to be different to the linear solver for newton solves because
177 /// explicit time steps only involve inverting a mass matrix. This can be
178 /// done very efficiently by, e.g. CG with a diagonal predconditioner).
180
181 /// Pointer to the eigen solver for the problem
183
184 // Pointer to handler
186
187 /// Pointer to the default linear solver
189
190 /// Pointer to the default eigensolver
192
193 /// Pointer to the default assembly handler
195
196 /// Pointer to global time for the problem
198
199 /// The Vector of time steppers (there could be many
200 /// different ones in multiphysics problems)
202
203 /// Pointer to a single explicit timestepper
205
206 /// Pointer to vector for backup of dofs
208
209 /// Has default set_initial_condition function been called?
210 /// Default: false
212
213 /// Use the globally convergent newton method
215
216 /// Boolean to indicate that empty
217 /// actions_before_read_unstructured_meshes() function has been called.
219
220 /// Boolean to indicate that empty
221 /// actions_after_read_unstructured_meshes() function has been called.
223
224 /// Boolean to indicate whether local dof pointers should be
225 /// stored in the elements
227
228 /// Use values from the time stepper predictor as an initial guess
230
231 protected:
232 /// Vector of pointers to copies of the problem used in adaptive
233 /// bifurcation tracking problems (ALH: TEMPORARY HACK, WILL BE FIXED)
235
236 /// Map used to determine whether the derivatives with respect to
237 /// a parameter should be finite differenced. The default is that
238 /// finite differences should be used
239 std::map<double*, bool> Calculate_dparameter_analytic;
240
241 /// Map used to determine whether the hessian products should be
242 /// computed using finite differences. The default is that finite
243 /// differences will be used
245
246 public:
247 /// Hook for debugging. Can be overloaded in driver code; argument
248 /// allows identification of where we're coming from
249 virtual void debug_hook_fct(const unsigned& i)
250 {
251 oomph_info << "Called empty hook fct with i=" << i << std::endl;
252 }
253
254 /// Function to turn on analytic calculation of the parameter
255 /// derivatives in continuation and bifurcation detection problems
256 inline void set_analytic_dparameter(double* const& parameter_pt)
257 {
258 Calculate_dparameter_analytic[parameter_pt] = true;
259 }
260
261 /// Function to turn off analytic calculation of the parameter
262 /// derivatives in continuation and bifurcation detection problems
263 inline void unset_analytic_dparameter(double* const& parameter_pt)
264 {
265 // Find the iterator to the parameter
266 std::map<double*, bool>::iterator it =
267 Calculate_dparameter_analytic.find(parameter_pt);
268 // If the parameter has been found, erase it
269 if (it != Calculate_dparameter_analytic.end())
270 {
272 }
273 }
274
275 /// Function to determine whether the parameter derivatives
276 /// are calculated analytically
278 double* const& parameter_pt)
279 {
280 // If the entry is found then the iterator returned from the find
281 // command will NOT be equal to the end and the expression will
282 // return true, otherwise it will return false
283 return (Calculate_dparameter_analytic.find(parameter_pt) !=
285 }
286
287 /// Function to turn on analytic calculation of the parameter
288 /// derivatives in continuation and bifurcation detection problems
290 {
292 }
293
294 /// Function to turn off analytic calculation of the parameter
295 /// derivatives in continuation and bifurcation detection problems
297 {
299 }
300
301 /// Function to determine whether the hessian products
302 /// are calculated analytically
304 {
306 }
307
308 /// Set all pinned values to zero.
309 /// Used to set boundary conditions to be homogeneous in the copy
310 /// of the problem used in adaptive bifurcation tracking
311 /// (ALH: TEMPORARY HACK, WILL BE FIXED)
313
314
315 /// Flag to allow suppression of warning messages re reading in
316 /// unstructured meshes during restart.
318
319
320 private:
321 /// Private helper function that actually performs the unsteady
322 /// "doubly" adaptive Newton solve. See actual (non-helper) functions
323 /// for description of parameters.
325 const double& dt,
326 const double& epsilon,
327 const unsigned& max_adapt,
328 const unsigned& suppress_resolve_after_spatial_adapt,
329 const bool& first,
330 const bool& shift = true);
331
332
333 /// Helper function to do compund refinement of (all) refineable
334 /// (sub)mesh(es) uniformly as many times as specified in vector and
335 /// rebuild problem; doc refinement process. Set boolean argument
336 /// to true if you want to prune immediately after refining the meshes
337 /// individually.
338 void refine_uniformly_aux(const Vector<unsigned>& nrefine_for_mesh,
339 DocInfo& doc_info,
340 const bool& prune);
341
342 /// Helper function to do compund p-refinement of (all) p-refineable
343 /// (sub)mesh(es) uniformly as many times as specified in vector and
344 /// rebuild problem; doc refinement process. Set boolean argument
345 /// to true if you want to prune immediately after refining the meshes
346 /// individually.
347 void p_refine_uniformly_aux(const Vector<unsigned>& nrefine_for_mesh,
348 DocInfo& doc_info,
349 const bool& prune);
350
351 /// Helper function to re-setup the Base_mesh enumeration
352 /// (used during load balancing) after pruning
354
355 /// Private helper function that is used to assemble the Jacobian
356 /// matrix in the case when the storage is row or column compressed.
357 /// The boolean Flag indicates
358 /// if we want compressed row format (true) or compressed column.
359 /// This version uses vectors of pairs.
361 Vector<int*>& column_or_row_index,
362 Vector<int*>& row_or_column_start,
363 Vector<double*>& value,
364 Vector<unsigned>& nnz,
365 Vector<double*>& residual,
366 bool compressed_row_flag);
367
368 /// Private helper function that is used to assemble the Jacobian
369 /// matrix in the case when the storage is row or column compressed.
370 /// The boolean Flag indicates
371 /// if we want compressed row format (true) or compressed column.
372 /// This version uses two vectors.
374 Vector<int*>& column_or_row_index,
375 Vector<int*>& row_or_column_start,
376 Vector<double*>& value,
377 Vector<unsigned>& nnz,
378 Vector<double*>& residual,
379 bool compressed_row_flag);
380
381 /// Private helper function that is used to assemble the Jacobian
382 /// matrix in the case when the storage is row or column compressed.
383 /// The boolean Flag indicates
384 /// if we want compressed row format (true) or compressed column.
385 /// This version uses maps
387 Vector<int*>& column_or_row_index,
388 Vector<int*>& row_or_column_start,
389 Vector<double*>& value,
390 Vector<unsigned>& nnz,
391 Vector<double*>& residual,
392 bool compressed_row_flag);
393
394 /// Private helper function that is used to assemble the Jacobian
395 /// matrix in the case when the storage is row or column compressed.
396 /// The boolean Flag indicates
397 /// if we want compressed row format (true) or compressed column.
398 /// This version uses lists
400 Vector<int*>& column_or_row_index,
401 Vector<int*>& row_or_column_start,
402 Vector<double*>& value,
403 Vector<unsigned>& nnz,
404 Vector<double*>& residual,
405 bool compressed_row_flag);
406
407 /// Private helper function that is used to assemble the Jacobian
408 /// matrix in the case when the storage is row or column compressed.
409 /// The boolean Flag indicates
410 /// if we want compressed row format (true) or compressed column.
411 /// This version uses lists
413 Vector<int*>& column_or_row_index,
414 Vector<int*>& row_or_column_start,
415 Vector<double*>& value,
416 Vector<unsigned>& nnz,
417 Vector<double*>& residual,
418 bool compressed_row_flag);
419
420 /// Vector of global data: "Nobody" (i.e. none of the elements etc.)
421 /// is "in charge" of this Data so it would be overlooked when it
422 /// comes to equation-numbering, timestepping etc. Including
423 /// (pointers) to such Data in here, puts the Problem itself
424 /// "in charge" of these tasks.
426
427 /// A function that performs the guts of the continuation
428 /// derivative calculation in arc length continuation problems.
430
431 /// A function that performs the guts of the continuation
432 /// derivative calculation in arc-length continuation problems using
433 /// finite differences
435 double* const& parameter_pt);
436
437 /// A function that is used to adapt a bifurcation-tracking
438 /// problem, which requires separate interpolation of the
439 /// associated eigenfunction. The error measure is chosen to be
440 /// a suitable combination of the errors in the base flow and the
441 /// eigenfunction
442 void bifurcation_adapt_helper(unsigned& n_refined,
443 unsigned& n_unrefined,
444 const unsigned& bifurcation_type,
445 const bool& actually_adapt = true);
446
447 /// A function that is used to document the errors used
448 /// in the adaptive solution of bifurcation problems.
449 void bifurcation_adapt_doc_errors(const unsigned& bifurcation_type);
450
451 // ALH_TEMP_DEVELOPMENT
452 protected:
453 /// The distribution of the DOFs in this problem.
454 /// This object is created in the Problem constructor and setup
455 /// when assign_eqn_numbers(...) is called.
456 /// If this problem is distributed then this distribution will match
457 /// the distribution of the equation numbers.
458 /// If this problem is not distributed then this distribution will
459 /// be uniform over all processors.
461
462 private:
463#ifdef OOMPH_HAS_MPI
464
465 /// Helper method that returns the (unique) global equations to which
466 /// the elements in the range el_lo to el_hi contribute on this
467 /// processor using the given assembly_handler
469 const unsigned& el_lo,
470 const unsigned& el_hi,
471 Vector<unsigned>& my_eqns);
472
473 /// Helper method to assemble CRDoubleMatrices from distributed
474 /// on multiple processors.
476 const LinearAlgebraDistribution* const& dist_pt,
477 Vector<int*>& column_or_row_index,
478 Vector<int*>& row_or_column_start,
479 Vector<double*>& value,
480 Vector<unsigned>& nnz,
481 Vector<double*>& residuals);
482
483 /// A private helper function to
484 /// copy the haloed equation numbers into the halo equation numbers,
485 /// either for the problem's one and only mesh or for all of its
486 /// submeshes. Bools control if we deal with data associated with
487 /// external halo/ed elements/nodes or the "normal" halo/ed ones.
488 void copy_haloed_eqn_numbers_helper(const bool& do_halos,
489 const bool& do_external_halos);
490
491 /// Private helper function to remove repeated data
492 /// in external haloed elements in specified mesh. Bool is true if some data
493 /// was removed -- this usually requires re-running through certain
494 /// parts of the equation numbering procedure.
496 bool& actually_removed_some_data);
497
498
499 /// Consolidate external halo node storage by removing nulled out
500 /// pointers in external halo and haloed schemes for all meshes.
502
503 /// Helper function to re-assign the first and last elements to be
504 /// assembled by each processor during parallel assembly for
505 /// non-distributed problem.
507
508 /// Boolean to switch on assessment of load imbalance in parallel
509 /// assembly of distributed problem
511
512 /// Flag to use "default partition" during load balance.
513 /// Should only be set to true when run in validation mode.
515
516 /// First element to be assembled by given processor for
517 /// non-distributed problem (only kept up to date when default assignment
518 /// is used)
520
521 /// Last element (plus one) to be assembled by given processor
522 /// for non-distributed problem (only kept up to date when default
523 /// assignment is used)
525
526 /// Boolean indicating that the division of elements over processors
527 /// during the assembly process must be re-load-balanced.
528 /// (only used for non-distributed problems)
530
531 /// Map which stores the correspondence between a root element and
532 /// its element number (plus one) within the global mesh at the point
533 /// when it is distributed. NB a root element in this instance is one
534 /// of the elements in the uniformly-refined mesh at the point when
535 /// Problem::distribute() is called, since these elements become roots
536 /// on each of the processors involved in the distribution. Null when
537 /// element doesn't exist following the adjustment of this when pruning.
538 std::map<GeneralisedElement*, unsigned> Base_mesh_element_number_plus_one;
539
540
541 /// Vector to store the correspondence between a root element and
542 /// its element number within the global mesh at the point when it is
543 /// distributed. NB a root element in this instance is one of the elements
544 /// in the uniformly-refined mesh at the point when Problem::distribute()
545 /// is called, since these elements become roots on each of the processors
546 /// involved in the distribution. Null when element doesn't exist
547 /// following the adjustment of this when pruning.
549
550#endif
551
552 protected:
553 /// Vector of pointers to dofs
555
556 /// Counter that records how many elements contribute to each dof.
557 /// Used to determine the (discrete) arc-length automatically.
558 /// It really should be an integer, but is a double so that the
559 /// distribution information can be used for distributed problems
561
562 /// Function that populates the Element_counter_per_dof vector
563 /// with the number of elements that contribute to each dof. For example,
564 /// with linear elements in 1D each dof contains contributions from two
565 /// elements apart from those on the boundary. Returns the total number
566 /// of elements in the problem
568
569
570#ifdef OOMPH_HAS_MPI
571
572 /// Storage for assembly times (used for load balancing)
574
575 /// Pointer to the halo scheme for any global vectors
576 /// that have the Dof_distribution
578
579 /// Storage for the halo degrees of freedom (only required)
580 /// when accessing via the global equation number in distributed problems
582
583 /// Function that is used to setup the halo scheme
585
586#endif
587
588 //--------------------- Newton solver parameters
589
590 /// Relaxation fator for Newton method (only a fractional Newton
591 /// correction is applied if this is less than 1; default 1).
593
594 /// The Tolerance below which the Newton Method is deemed to have
595 /// converged
597
598 /// Maximum number of Newton iterations
600
601 /// Actual number of Newton iterations taken during the most recent
602 /// iteration
604
605 /// Maximum residuals at start and after each newton iteration
607
608 /// Maximum desired residual:
609 /// if the maximum residual exceeds this value, the program will exit
611
612 /// Bool to specify what to do if a Newton solve fails within a
613 /// time adaptive solve. Default (false) is to half the step and try
614 /// again. If true then crash instead.
616
617 /// Is re-use of Jacobian in Newton iteration enabled? Default: false
619
620 /// Has a Jacobian been computed (and can therefore be re-used
621 /// if required)? Default: false
623
624 /// Boolean flag indicating if we're dealing with a linear or
625 /// nonlinear Problem -- if set to false the Newton solver will not check
626 /// the residual before or after the linear solve. Set to true by default;
627 /// can be over-written in specific Problem class.
629
630 /// Protected boolean flag to halt program execution
631 /// during sparse assemble process to assess peak memory
632 /// usage. Initialised to false (obviously!)
634
635 /// Protected boolean flag to provide comprehensive timimings
636 /// during problem distribution. Initialised to false.
638
639 /// Flag to determine which sparse assembly method to use
640 /// By default we use assembly by vectors of pairs.
642
643 /// Enumerated flags to determine which sparse assembly method is used
645 {
651 };
652
653 /// the number of elements to initially allocate for a matrix row
654 /// within the sparse_assembly_with_two_arrays(...) method (if a matrix
655 /// of this size has not been assembled already). If a matrix of this size
656 /// has been assembled then the number of elements in each row in that
657 /// matrix is used as the initial allocation
659
660 /// the number of elements to add to a matrix row when the initial
661 /// allocation is exceeded within the sparse_assemble_with_two_arrays(...)
662 /// method.
664
665 /// the number of elements in each row of a compressed matrix
666 /// in the previous matrix assembly.
668
669 /// A tolerance used to determine whether the entry in a sparse
670 /// matrix is zero. If it is then storage need not be allocated.
672
673 /// Protected helper function that is used to assemble the Jacobian
674 /// matrix in the case when the storage is row or column compressed.
675 /// The boolean Flag indicates
676 /// if we want compressed row format (true) or compressed column.
678 Vector<int*>& column_or_row_index,
679 Vector<int*>& row_or_column_start,
680 Vector<double*>& value,
681 Vector<unsigned>& nnz,
682 Vector<double*>& residual,
683 bool compressed_row_flag);
684
686
687 //---------------------Explicit time-stepping parameters
688
689 /// Is re-use of the mass matrix in explicit timestepping enabled
690 /// Default:false
692
693 /// Has the mass matrix been computed (and can therefore be reused)
694 /// Default: false
696
697 //---------------------Discontinuous control flags
698 /// Is the problem a discontinuous one, i.e. can the
699 /// elemental contributions be treated independently. Default: false
701
702
703 //--------------------- Adaptive time-stepping parameters
704
705 /// Minimum desired dt: if dt falls below this value, exit
707
708 /// Maximum desired dt
710
711 /// Maximum possible increase of dt between time-steps in adaptive
712 /// schemes
714
715 /// Minimum allowed decrease of dt between time-steps in adaptive
716 /// schemes. Lower scaling values will reject the time-step (and retry
717 /// with a smaller dt).
719
720 /// If Minimum_dt_but_still_proceed positive, then dt will not be
721 /// reduced below this value during adaptive timestepping and the
722 /// computation will continue with this value, accepting the larger
723 /// errors that this will incur). Note: This option is disabled by default
724 /// as this value is initialised to -1.0.
726
727
728 //--------------------- Arc-length continuation paramaters
729
730 /// Boolean to control whether arc-length should be scaled
732
733 /// Proportion of the arc-length to taken by the parameter
735
736 /// Value of the scaling parameter required so that the parameter
737 /// occupies the desired proportion of the arc length. NOTE: If you wish
738 /// to change this then make sure to set the value of Scale_arc_length
739 /// to false to ensure the value of this isn't overwritten during the
740 /// arc-length process. Instead of changing this variable, it's better
741 /// to actually update the Desired_proportion_of_arc_length value.
743
744 /// Storage for the sign of the global Jacobian
746
747 /// The direction of the change in parameter that will ensure that
748 /// a branch is followed in one direction only
750
751 /// Storage for the derivative of the global parameter wrt arc-length
753
754 /// Storage for the present value of the global parameter
756
757 /// Boolean to control original or new storage of dof stuff
759
760 /// Storage for the single static continuation timestorage object
762
763 /// Storage for the offset for the continuation derivatives from the
764 /// original dofs (default is 1, but this will be differnet when continuing
765 /// bifurcations and periodic orbits)
767
768 /// Storage for the offset for the current values of dofs from the
769 /// original dofs (default is 2, but this will be differnet when continuing
770 /// bifurcations and periodic orbits)
772
773 /// Storage for the derivative of the problem variables wrt arc-length
775
776 /// Storage for the present values of the variables
778
779 /// Storage for the current step value
781
782 /// The desired number of Newton Steps to reach convergence at each
783 /// step along the arc
785
786 /// Minimum desired value of arc-length
788
789 /// Boolean to control bifurcation detection via determinant of Jacobian
791
792 /// Boolean to control wheter bisection is used to located bifurcation
794
795 /// Boolean to indicate whether a sign change has occured in the Jacobian
797
798 /// Boolean to indicate whether an arc-length step has been taken
800
801 /// Boolean to specify which scheme to use to calculate the continuation
802 /// derivatievs
804
805 public:
806 /// If we have MPI return the "problem has been distributed" flag,
807 /// otherwise it can't be distributed so return false.
808 bool distributed() const
809 {
810#ifdef OOMPH_HAS_MPI
812#else
813 return false;
814#endif
815 }
816
817#ifdef OOMPH_HAS_MPI
818
819
820 /// enum for distribution of distributed jacobians.
821 /// 1 - Automatic - the Problem distribution is employed, unless any
822 /// processor has number of rows equal to 110% of N/P, in which case
823 /// uniform distribution is employed.
824 /// 2 - Problem - the jacobian on processor p only contains rows that
825 /// correspond to equations that are on this processor. (minimises
826 /// communication)
827 /// 3 - Uniform - each processor holds as close to N/P matrix rows as
828 /// possible. (very well load balanced)
832
833 /// accesss function to the distributed matrix distribution method
834 /// 1 - Automatic - the Problem distribution is employed, unless any
835 /// processor has number of rows equal to 110% of N/P, in which case
836 /// uniform distribution is employed.
837 /// 2 - Problem - the jacobian on processor p only contains rows that
838 /// correspond to equations that are on this processor. (minimises
839 /// communication)
840 /// 3 - Uniform - each processor holds as close to N/P matrix rows as
841 /// possible. (very well load balanced)
843 {
845 }
846
847 /// Enable doc of load imbalance in parallel
848 /// assembly of distributed problem
850 {
852 }
853
854 /// Disable doc of load imbalance in parallel
855 /// assembly of distributed problem
857 {
859 }
860
861 /// Return vector of most-recent elemental assembly times
862 /// (used for load balancing). Zero sized if no Jacobian has been
863 /// computed since last re-assignment of equation numbers
865 {
867 }
868
869 /// Clear storage of most-recent elemental assembly times
870 /// (used for load balancing). Next load balancing operation
871 /// will be based on the number of elements associated with a tree root.
873 {
876 }
877
878 private:
879 /// Load balance helper routine: Get data to be sent to other
880 /// processors during load balancing and other information about
881 /// re-distribution.
882 /// - target_domain_for_local_non_halo_element: Input, generated by METIS.
883 /// target_domain_for_local_non_halo_element[e] contains the number
884 /// of the domain [0,1,...,nproc-1] to which non-halo element e on THE
885 /// CURRENT PROCESSOR ONLY has been assigned. The order of the non-halo
886 /// elements is the same as in the Problem's mesh, with the halo
887 /// elements being skipped.
888 /// - send_n: Output, number of data to be sent to each processor
889 /// - send_data: Output, storage for all values to be sent to all processors
890 /// - send_displacement: Output, start location within send_data for data to
891 /// be sent to each processor
892 /// - max_refinement_level_overall: Output, max. refinement level of any
893 /// element
895 const Vector<unsigned>& element_domain_on_this_proc,
896 Vector<int>& send_n,
897 Vector<double>& send_data,
898 Vector<int>& send_displacement,
899 Vector<unsigned>& old_domain_for_base_element,
900 Vector<unsigned>& new_domain_for_base_element,
901 unsigned& max_refinement_level_overall);
902
903 /// Get flat-packed refinement pattern for each root element in
904 /// current mesh (labeled by unique number of root element in unrefined base
905 /// mesh). The vector stored for each root element contains the following
906 /// information:
907 /// - First entry: Number of tree nodes (not just leaves!) in refinement
908 /// tree emanating from this root [Zero if root element is not refineable]
909 /// - Loop over all refinement levels
910 /// - Loop over all tree nodes (not just leaves!)
911 /// - If associated element exists when the mesh has been refined to
912 /// this level (either because it has been refined to this level or
913 /// because it's less refined): 1
914 /// - If the element is to be refined: 1; else: 0
915 /// - else (element doesn't exist when mesh is refined to this level
916 /// (because it's more refined): 0
917 /// .
918 /// .
919 /// .
921 const Vector<unsigned>& old_domain_for_base_element,
922 const Vector<unsigned>& new_domain_for_base_element,
923 const unsigned& max_refinement_level_overall,
924 std::map<unsigned, Vector<unsigned>>&
925 flat_packed_refinement_info_for_root);
926
927
928 /// Load balance helper routine: Send data to other
929 /// processors during load balancing.
930 /// - send_n: Input, number of data to be sent to each processor
931 /// - send_data: Input, storage for all values to be sent to all processors
932 /// - send_displacement: Input, start location within send_data for data to
933 /// be sent to each processor
935 Vector<int>& send_n,
936 Vector<double>& send_data,
937 Vector<int>& send_displacement);
938
939 /// Send refinement information between processors
941 Vector<unsigned>& old_domain_for_base_element,
942 Vector<unsigned>& new_domain_for_base_element,
943 const unsigned& max_refinement_level_overall,
944 std::map<unsigned, Vector<unsigned>>& refinement_info_for_root_local,
945 Vector<Vector<Vector<unsigned>>>& refinement_info_for_root_elements);
946
947 /// The distributed matrix distribution method
948 /// 1 - Automatic - the Problem distribution is employed, unless any
949 /// processor has number of rows equal to 110% of N/P, in which case
950 /// uniform distribution is employed.
951 /// 2 - Problem - the jacobian on processor p only contains rows that
952 /// correspond to equations that are on this processor. (minimises
953 /// communication)
954 /// 3 - Uniform - each processor holds as close to N/P matrix rows as
955 /// possible. (very well load balanced)
957
958 /// The amount of data allocated during the previous parallel sparse
959 /// assemble. Used to optimise the next call to parallel_sparse_assemble()
961
962 protected:
963 /// Has the problem been distributed amongst multiple processors?
965
966 /// Boolean to bypass check of increase in dofs during pruning
968
969 public:
970 /// Enable problem distributed
972 {
974 }
975
976 /// Disable problem distributed
978 {
980 }
981
982 /// Set default first and last elements for parallel assembly
983 /// of non-distributed problem.
985
986 /// Manually set first and last elements for parallel assembly
987 /// of non-distributed problem.
989 Vector<unsigned>& first_el_for_assembly,
990 Vector<unsigned>& last_el_for_assembly)
991 {
992 First_el_for_assembly = first_el_for_assembly;
993 Last_el_plus_one_for_assembly = last_el_for_assembly;
994 unsigned n = last_el_for_assembly.size();
995 for (unsigned i = 0; i < n; i++)
996 {
998 }
999 }
1000
1001
1002 /// Get first and last elements for parallel assembly
1003 /// of non-distributed problem.
1005 Vector<unsigned>& first_el_for_assembly,
1006 Vector<unsigned>& last_el_for_assembly) const
1007 {
1008 first_el_for_assembly = First_el_for_assembly;
1009 last_el_for_assembly = Last_el_plus_one_for_assembly;
1010 unsigned n = last_el_for_assembly.size();
1011 for (unsigned i = 0; i < n; i++)
1012 {
1013 last_el_for_assembly[i] -= 1;
1014 }
1015 }
1016
1017#endif
1018
1019 /// Actions that are to be performed before a mesh adaptation.
1020 /// These might include removing any additional elements, such as traction
1021 /// boundary elements before the adaptation.
1022 virtual void actions_before_adapt() {}
1023
1024 /// Actions that are to be performed after a mesh adaptation.
1025 virtual void actions_after_adapt() {}
1026
1027 protected:
1028 /// Any actions that are to be performed before a complete
1029 /// Newton solve (e.g. adjust boundary conditions). CAREFUL: This
1030 /// step should (and if the FD-based LinearSolver FD_LU is used, must)
1031 /// only update values that are pinned!
1033
1034 /// Any actions that are to be performed after a complete Newton
1035 /// solve, e.g. post processing. CAREFUL: This step should (and if the
1036 /// FD-based LinearSolver FD_LU is used, must) only update values that are
1037 /// pinned!
1039
1040 /// Any actions that are to be performed before
1041 /// the residual is checked in the Newton method, e.g. update
1042 /// any boundary conditions that depend upon
1043 /// variables of the problem; update any `dependent' variables; or
1044 /// perform an update of the nodal positions in SpineMeshes etc.
1045 /// CAREFUL: This
1046 /// step should (and if the FD-based LinearSolver FD_LU is used, must)
1047 /// only update values that are pinned!
1049
1050 /// Any actions that are to be performed before each individual
1051 /// Newton step. Most likely to be used for diagnostic purposes
1052 /// to doc the solution during a non-converging iteration, say.
1054
1055 /// Any actions that are to be performed after each individual
1056 /// Newton step. Most likely to be used for diagnostic purposes
1057 /// to doc the solution during a non-converging iteration, say.
1059
1060 /// Actions that should be performed before each implicit time step.
1061 /// This is needed when one wants
1062 /// to solve a steady problem before timestepping and needs to distinguish
1063 /// between the two cases
1065
1066 /// Actions that should be performed after each implicit time step.
1067 /// This is needed when one wants
1068 /// to solve a steady problem before timestepping and needs to distinguish
1069 /// between the two cases
1071
1072 /// Actions that should be performed after each implicit time step.
1073 /// This is needed if your actions_after_implicit_timestep() modify the
1074 /// solution in a way that affects the error estimate.
1076
1077 /// Actions that should be performed before each explicit time step.
1079
1080 /// Actions that should be performed after each explicit time step.
1082
1083 /// Actions that are to be performed before reading in
1084 /// restart data for problems involving unstructured bulk meshes.
1085 /// If the problem contains such meshes we need
1086 /// to strip out any face elements that are attached to them
1087 /// because restart of unstructured meshes re-creates their elements
1088 /// and nodes from scratch, leading to dangling pointers from the
1089 /// face elements to the old elements and nodes. This function is
1090 /// virtual and (practically) empty
1091 /// but toggles a flag to indicate that it has been called. This is used to
1092 /// issue a warning, prompting the user to consider overloading it
1093 /// if the problem is found to contain unstructured bulk meshes during
1094 /// restarts.
1096 {
1098 }
1099
1100 /// Actions that are to be performed before reading in
1101 /// restart data for problems involving unstructured bulk meshes.
1102 /// Typically used to re-attach FaceElements, say, that were stripped
1103 /// out in actions_before_read_unstructured_meshes(). This function is
1104 /// virtual and (practically) empty
1105 /// but toggles a flag to indicate that it has been called. This is used to
1106 /// issue a warning, prompting the user to consider overloading it
1107 /// if the problem is found to contain unstructured bulk meshes during
1108 /// restarts.
1110 {
1112 }
1113
1114#ifdef OOMPH_HAS_MPI
1115 /// Actions to be performed before a (mesh) distribution
1117
1118 /// Actions to be performed after a (mesh) distribution
1120
1121#endif
1122
1123 /// Actions that are to be performed when the global parameter
1124 /// addressed by parameter_pt
1125 /// has been changed in the function get_derivative_wrt_global_parameter()
1126 /// The default is to call actions_before_newton_solve(),
1127 /// actions_before_newton_convergence_check() and
1128 /// actions_after_newton_solve().
1129 /// This could be amazingly inefficient in certain problems and should
1130 /// be overloaded in such cases. An example would be when a remesh is
1131 /// required in general, but the global parameter does not affect the mesh
1132 /// directly.
1134 double* const& parameter_pt)
1135 {
1136 // Default action should cover all possibilities
1140 }
1141
1142 /// Actions that are to be performed after a change in the parameter
1143 /// that is being varied as part of the solution of a bifurcation detection
1144 /// problem. The default is to call actions_before_newton_solve(),
1145 /// actions_before_newton_convergence_check() and
1146 /// actions_after_newton_solve().
1147 /// This could be amazingly inefficient in certain problems and should
1148 /// be overloaded in such cases. An example would be when a remesh is
1149 /// required in general, but the global parameter does not affect the mesh
1150 /// directly.
1152 {
1153 // Default action should cover all possibilities
1157 }
1158
1159 /// Empty virtual function; provides hook to perform actions after
1160 /// the increase in the arclength parameter (during continuation)
1161 virtual void actions_after_parameter_increase(double* const& parameter_pt)
1162 {
1163 }
1164
1165 /// Access function to the derivative of the i-th (local)
1166 /// dof with respect to
1167 /// the arc length, used in continuation
1168 inline double& dof_derivative(const unsigned& i)
1169 {
1171 {
1172 return Dof_derivative[i];
1173 }
1174 else
1175 {
1176 return (*(Dof_pt[i] + Dof_derivative_offset));
1177 }
1178 }
1179
1180 /// Access function to the current value of the i-th (local)
1181 /// dof at the start of a continuation step
1182 inline double& dof_current(const unsigned& i)
1183 {
1185 {
1186 return Dof_current[i];
1187 }
1188 else
1189 {
1190 return (*(Dof_pt[i] + Dof_current_offset));
1191 }
1192 }
1193
1194 /// Set initial condition (incl previous timesteps).
1195 /// We need to establish this interface
1196 /// because I.C. needs to be reset when problem is adapted during
1197 /// the first timestep.
1199 {
1200 std::ostringstream warn_message;
1201 warn_message
1202 << "Warning: We're using the default (empty) set_initial_condition().\n"
1203 << "If the initial conditions isn't re-assigned after a mesh adaption "
1204 "\n"
1205 << "the initial conditions will represent the interpolation of the \n"
1206 << "initial conditions that were assigned on the original coarse "
1207 "mesh.\n";
1208 OomphLibWarning(warn_message.str(),
1209 "Problem::set_initial_condition()",
1210 OOMPH_EXCEPTION_LOCATION);
1211
1212 // Indicate that this function has been called. This flag is set so
1213 // that the unsteady_newton_solve routine can be instructed whether
1214 // or not to shift the history values. If set_initial_condition() has
1215 // been overloaded than this (default) function won't be called, and
1216 // so this flag will remain false (its default value). If
1217 // set_initial_condition() has not been overloaded then this function
1218 // will be called and the flag set to true.
1220 }
1221
1222 /// Function to calculate a global error norm, used in adaptive
1223 /// timestepping to control the change in timestep. Individual errors for
1224 /// each data object can be obtained via the data timestepper's
1225 /// temporal_error_in_value or temporal_error_in_position functions and
1226 /// should be combined to construct a global norm. For example, in fluids
1227 /// problems a suitable norm is usually the weighted sum of the errors in
1228 /// the velocities; for moving mesh problems is it usually better to use the
1229 /// weighted sum of the errors in position.
1231 {
1232 std::string error_message =
1233 "The global_temporal_error_norm function will be problem-specific:\n";
1234 error_message += "Please write your own in your Problem class";
1235
1236 throw OomphLibError(
1237 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1238 return 0.0;
1239 }
1240
1241 /// The communicator for this problem
1243
1244 public:
1245 /// access function to the oomph-lib communicator
1247 {
1248 return Communicator_pt;
1249 }
1250
1251 /// access function to the oomph-lib communicator, const version
1253 {
1254 return Communicator_pt;
1255 }
1256
1257 /// Function pointer for spatial error estimator
1259 Vector<double>& elemental_error);
1260
1261 /// Function pointer for spatial error estimator with doc
1263 Mesh*& mesh_pt, Vector<double>& elemental_error, DocInfo& doc_info);
1264
1265 /// Constructor: Allocate space for one time stepper
1266 /// and set all pointers to NULL and set defaults for all
1267 /// parameters.
1268 Problem();
1269
1270 /// Broken copy constructor
1271 Problem(const Problem& dummy) = delete;
1272
1273 /// Broken assignment operator
1274 void operator=(const Problem&) = delete;
1275
1276 /// Virtual destructor to clean up memory
1277 virtual ~Problem();
1278
1279 /// Return a pointer to the global mesh
1281 {
1282 return Mesh_pt;
1283 }
1284
1285 /// Return a pointer to the global mesh (const version)
1286 Mesh* const& mesh_pt() const
1287 {
1288 return Mesh_pt;
1289 }
1290
1291 /// Return a pointer to the i-th submesh. If there are
1292 /// no submeshes, the 0-th submesh is the global mesh itself.
1293 Mesh*& mesh_pt(const unsigned& imesh)
1294 {
1295 // If there are no submeshes then the zero-th submesh is the
1296 // mesh itself
1297 if ((Sub_mesh_pt.size() == 0) && (imesh == 0))
1298 {
1299 return Mesh_pt;
1300 }
1301 else
1302 {
1303 return Sub_mesh_pt[imesh];
1304 }
1305 }
1306
1307 /// Return a pointer to the i-th submesh (const version)
1308 Mesh* const& mesh_pt(const unsigned& imesh) const
1309 {
1310 // If there are no submeshes then the zero-th submesh is the
1311 // mesh itself
1312 if ((Sub_mesh_pt.size() == 0) && (imesh == 0))
1313 {
1314 return Mesh_pt;
1315 }
1316 else
1317 {
1318 return Sub_mesh_pt[imesh];
1319 }
1320 }
1321
1322 /// Return number of submeshes
1323 unsigned nsub_mesh() const
1324 {
1325 return Sub_mesh_pt.size();
1326 }
1327
1328 /// Add a submesh to the problem and return its number, i, by which
1329 /// it can be accessed via mesh_pt(i).
1330 unsigned add_sub_mesh(Mesh* const& mesh_pt)
1331 {
1332 Sub_mesh_pt.push_back(mesh_pt);
1333 return Sub_mesh_pt.size() - 1;
1334 }
1335
1336
1337 /// Flush the problem's collection of sub-meshes. Must be followed
1338 /// by call to rebuild_global_mesh().
1340 {
1341 Sub_mesh_pt.clear();
1342 }
1343
1344 /// Build the global mesh by combining the all the submeshes.
1345 /// \b Note: The nodes boundary information refers to the boundary numbers
1346 /// within the submesh!
1347 void build_global_mesh();
1348
1349 /// If one of the submeshes has changed (e.g. by
1350 /// mesh adaptation) we need to update the global mesh.
1351 /// \b Note: The nodes boundary information refers to the
1352 /// boundary numbers within the submesh!
1353 void rebuild_global_mesh();
1354
1355#ifdef OOMPH_HAS_MPI
1356
1357 /// Function to build the Problem's base mesh; this must be supplied
1358 /// by the user if they wish to use the load_balance() functionality,
1359 /// which is only available to problems that have already been distributed.
1360 /// If the problem has multiple meshes, each mesh must be built, added as
1361 /// as a submesh, and a call to build_global_mesh() must be made
1362 /// in this function. On return from this function all meshes must
1363 /// have been refined to the same level that they were in the when
1364 /// Problem::distribute() was first called.
1365 virtual void build_mesh()
1366 {
1367 std::string error_message = "You are likely to have reached this error "
1368 "from Problem::load_balance()\n";
1369 error_message +=
1370 "The build_mesh() function is problem-specific and must be supplied:\n";
1371 error_message +=
1372 "Please write your own in your Problem class, ensuring that\n";
1373 error_message += "any (sub)meshes are built before calling "
1374 "Problem::build_global_mesh().\n";
1375 throw OomphLibError(
1376 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1377 }
1378
1379 /// Balance the load of a (possibly non-uniformly refined) problem
1380 /// that has already been distributed, by re-distributing elements over
1381 /// processors.
1383 {
1384 // Dummy DocInfo
1385 DocInfo doc_info;
1386 doc_info.disable_doc();
1387
1388 // Don't report stats
1389 bool report_stats = false;
1390
1391 // Dummy imposed partioning vector
1392 Vector<unsigned> input_target_domain_for_local_non_halo_element;
1393
1394 // Do it!
1396 doc_info, report_stats, input_target_domain_for_local_non_halo_element);
1397 }
1398
1399 /// Balance the load of a (possibly non-uniformly refined) problem
1400 /// that has already been distributed, by re-distributing elements over
1401 /// processors. Produce explicit stats of load balancing process if
1402 /// boolean, report_stats, is set to true.
1403 void load_balance(const bool& report_stats)
1404 {
1405 // Dummy DocInfo
1406 DocInfo doc_info;
1407 doc_info.disable_doc();
1408
1409 // Dummy imposed partioning vector
1410 Vector<unsigned> input_target_domain_for_local_non_halo_element;
1411
1412 // Do it
1414 doc_info, report_stats, input_target_domain_for_local_non_halo_element);
1415 }
1416
1417
1418 /// Balance the load of a (possibly non-uniformly refined) problem
1419 /// that has already been distributed, by re-distributing elements over
1420 /// processors. Produce explicit stats of load balancing process if
1421 /// boolean, report_stats, is set to true.
1422 void load_balance(DocInfo& doc_info, const bool& report_stats)
1423 {
1424 // Dummy imposed partioning vector
1425 Vector<unsigned> input_target_domain_for_local_non_halo_element;
1426
1427 // Do it
1429 doc_info, report_stats, input_target_domain_for_local_non_halo_element);
1430 }
1431
1432 /// Balance the load of a (possibly non-uniformly refined) problem
1433 /// that has already been distributed, by re-distributing elements over
1434 /// processors. Produce explicit stats of load balancing process if
1435 /// boolean, report_stats, is set to true and doc various bits of
1436 /// data (mainly for debugging) in directory specified by DocInfo object.
1437 /// If final input vector is non-zero-sized it provides an imposed
1438 /// partitioning.
1439 void load_balance(
1440 DocInfo& doc_info,
1441 const bool& report_stats,
1442 const Vector<unsigned>& input_target_domain_for_local_non_halo_element);
1443
1444 /// Set the use of the default partition in the load balance
1446 {
1448 }
1449
1450 /// Do not use the default partition in the load balance
1452 {
1454 }
1455
1456 /// Load balance helper routine: refine each new base (sub)mesh
1457 /// based upon the elements to be refined within each tree at each root
1458 /// on the current processor
1460 Vector<Vector<Vector<unsigned>>>& to_be_refined_on_each_root,
1461 const unsigned& max_level_overall);
1462
1463#endif
1464
1465 /// Return a pointer to the linear solver object
1467 {
1468 return Linear_solver_pt;
1469 }
1470
1471 /// Return a pointer to the linear solver object (const version)
1473 {
1474 return Linear_solver_pt;
1475 }
1476
1477 /// Return a pointer to the linear solver object used for explicit time
1478 /// stepping.
1480 {
1482 }
1483
1484 /// Return a pointer to the linear solver object used for explicit time
1485 /// stepping (const version).
1487 {
1489 }
1490
1491 /// Return a pointer to the eigen solver object
1493 {
1494 return Eigen_solver_pt;
1495 }
1496
1497 /// Return a pointer to the eigen solver object (const version)
1499 {
1500 return Eigen_solver_pt;
1501 }
1502
1503 /// Return a pointer to the global time object
1505 {
1506 return Time_pt;
1507 }
1508
1509 /// Return a pointer to the global time object (const version).
1510 Time* time_pt() const
1511 {
1512 return Time_pt;
1513 }
1514
1515 /// Return the current value of continuous time
1516 double& time();
1517
1518 /// Return the current value of continuous time (const version)
1519 double time() const;
1520
1521
1522 /// Access function for the pointer to the first (presumably only)
1523 /// timestepper
1525 {
1526 if (Time_stepper_pt.size() == 0)
1527 {
1528 throw OomphLibError("No timestepper allocated yet\n",
1529 OOMPH_CURRENT_FUNCTION,
1530 OOMPH_EXCEPTION_LOCATION);
1531 }
1532 return Time_stepper_pt[0];
1533 }
1534
1535 /// Access function for the pointer to the first (presumably only)
1536 /// timestepper
1538 {
1539 if (Time_stepper_pt.size() == 0)
1540 {
1541 throw OomphLibError("No timestepper allocated yet\n",
1542 OOMPH_CURRENT_FUNCTION,
1543 OOMPH_EXCEPTION_LOCATION);
1544 }
1545 return Time_stepper_pt[0];
1546 }
1547
1548 /// Return a pointer to the i-th timestepper
1549 TimeStepper*& time_stepper_pt(const unsigned& i)
1550 {
1551 return Time_stepper_pt[i];
1552 }
1553
1554 /// Return a pointer to the explicit timestepper
1556 {
1558 }
1559
1560 /// Set all problem data to have the same timestepper
1561 /// (timestepper_pt) Return the new number of dofs in the problem
1562 unsigned long set_timestepper_for_all_data(
1564 const bool& preserve_existing_data = false);
1565
1566 /// Shift all values along to prepare for next timestep
1567 virtual void shift_time_values();
1568
1569 /// Return a pointer to the assembly handler object
1571 {
1572 return Assembly_handler_pt;
1573 }
1574
1575 /// Return a pointer to the assembly handler object (const version)
1577 {
1578 return Assembly_handler_pt;
1579 }
1580
1581 /// Access function to min timestep in adaptive timestepping
1582 double& minimum_dt()
1583 {
1584 return Minimum_dt;
1585 }
1586
1587 /// Access function to max timestep in adaptive timestepping
1588 double& maximum_dt()
1589 {
1590 return Maximum_dt;
1591 }
1592
1593 /// Access function to max Newton iterations before giving up.
1595 {
1596 return Max_newton_iterations;
1597 }
1598
1599 /// Access function to Problem_is_nonlinear.
1600 void problem_is_nonlinear(const bool& prob_lin)
1601 {
1602 Problem_is_nonlinear = prob_lin;
1603 }
1604
1605
1606 /// Access function to max residuals in Newton iterations before
1607 /// giving up.
1609 {
1610 return Max_residuals;
1611 }
1612
1613 /// Access function for Time_adaptive_newton_crash_on_solve_fail.
1615 {
1617 }
1618
1619 /// Access function to tolererance of the Newton solver, i.e. the
1620 /// maximum value of the residuals that will be accepted.
1622 {
1624 }
1625
1626 /// Add a timestepper to the problem. The function will automatically
1627 /// create or resize the Time object so that it contains the appropriate
1628 /// number of levels of storage.
1630
1631 /// Set the explicit timestepper for the problem. The function
1632 /// will automatically create or resize the Time object so that it contains
1633 /// the appropriate number of levels of storage
1636
1637
1638 /// Set all timesteps to the same value, dt, and assign
1639 /// weights for all timesteppers in the problem
1640 void initialise_dt(const double& dt);
1641
1642 /// Set the value of the timesteps to be equal to the values passed
1643 /// in a vector and assign weights for all timesteppers in the problem
1644 void initialise_dt(const Vector<double>& dt);
1645
1646 /// Return a pointer to the the i-th global data object
1647 Data*& global_data_pt(const unsigned& i)
1648 {
1649 return Global_data_pt[i];
1650 }
1651
1652 /// Add Data to the Problem's global data -- the Problem will
1653 /// perform equation numbering etc. for such Data
1655 {
1656 Global_data_pt.push_back(global_data_pt);
1657 }
1658
1659
1660 /// Flush the Problem's global data -- resizes container to zero.
1661 /// Data objects are not deleted!
1663 {
1664 Global_data_pt.resize(0);
1665 }
1666
1667 /// Return the pointer to the dof distribution (read-only)
1669 {
1670 return Dof_distribution_pt;
1671 }
1672
1673 /// Return the number of dofs
1674 unsigned long ndof() const
1675 {
1676 return Dof_distribution_pt->nrow();
1677 }
1678
1679 /// Return the number of time steppers
1680 unsigned ntime_stepper() const
1681 {
1682 return Time_stepper_pt.size();
1683 }
1684
1685 /// Return the number of global data values
1686 unsigned nglobal_data() const
1687 {
1688 return Global_data_pt.size();
1689 }
1690
1691 /// Self-test: Check meshes and global data. Return 0 for OK
1692 unsigned self_test();
1693
1694 /// Insist that local dof pointers are set up in each element
1695 /// when equation numbering takes place
1697 {
1699 }
1700
1701 /// Insist that local dof pointers are NOT set up in each element
1702 /// when equation numbering takes place (the default)
1704 {
1706 }
1707
1708 /// Assign all equation numbers for problem: Deals with global
1709 /// data (= data that isn't attached to any elements) and then
1710 /// does the equation numbering for the elements. Virtual so it
1711 /// can be overloaded in MPI problems. Bool argument can be set to false
1712 /// to ignore assigning local equation numbers (found to be necessary in
1713 /// the parallel implementation of locate_zeta between multiple meshes).
1714 unsigned long assign_eqn_numbers(
1715 const bool& assign_local_eqn_numbers = true);
1716
1717 /// Function to describe the dofs in terms of the global
1718 /// equation number, i.e. what type of value (nodal value of
1719 /// a Node; value in a Data object; value of internal Data in an
1720 /// element; etc) is the unknown with a certain global equation number.
1721 /// Output stream defaults to oomph_info.
1722 void describe_dofs(std::ostream& out = *(oomph_info.stream_pt())) const;
1723
1724 /// Indicate that the problem involves discontinuous elements
1725 /// This allows for a more efficiently assembly and inversion of the
1726 /// mass matrix
1728 {
1730 }
1731
1732 /// Disable the use of a discontinuous-element formulation.
1733 /// Note that the methods will all still work even if the elements are
1734 /// discontinuous, we will just be assembling a larger system matrix than
1735 /// necessary.
1737 {
1739 }
1740
1741 /// Return the vector of dofs, i.e. a vector containing the current
1742 /// values of all unknowns.
1743 void get_dofs(DoubleVector& dofs) const;
1744
1745 /// Return vector of the t'th history value of all dofs.
1746 void get_dofs(const unsigned& t, DoubleVector& dofs) const;
1747
1748 /// Set the values of the dofs
1749 void set_dofs(const DoubleVector& dofs);
1750
1751 /// Set the history values of the dofs
1752 void set_dofs(const unsigned& t, DoubleVector& dofs);
1753
1754 /// Set history values of dofs from the type of vector stored in
1755 /// problem::Dof_pt.
1756 void set_dofs(const unsigned& t, Vector<double*>& dof_pt);
1757
1758 /// Add lambda x incremenet_dofs[l] to the l-th dof
1759 void add_to_dofs(const double& lambda, const DoubleVector& increment_dofs);
1760
1761 /// Return a pointer to the dof, indexed by global equation number
1762 /// which may be haloed or stored locally. If it is haloed, a local copy
1763 /// must have been requested on this processor in the Halo_scheme_pt.
1764 inline double* global_dof_pt(const unsigned& i)
1765 {
1766#ifdef OOMPH_HAS_MPI
1767 // If the problem is distributed I have to do something different
1769 {
1770#ifdef PARANOID
1771 if (Halo_scheme_pt == 0)
1772 {
1773 std::ostringstream error_stream;
1774 error_stream
1775 << "Direct access to global dofs in a distributed problem is only\n"
1776 << "possible if the function setup_dof_halo_scheme() has been "
1777 "called.\n"
1778 << "You can access local dofs via Problem::dof(i).\n";
1779 throw OomphLibError(error_stream.str(),
1780 OOMPH_CURRENT_FUNCTION,
1781 OOMPH_EXCEPTION_LOCATION);
1782 }
1783#endif
1784
1785 // Work out the local coordinate of the dof
1786 // based on the original distribution stored in the Halo_scheme
1787 const unsigned first_row_local =
1789 const unsigned n_row_local =
1791
1792 // If we are in range then just call the local value
1793 if ((i >= first_row_local) && (i < first_row_local + n_row_local))
1794 {
1795 return this->Dof_pt[i - first_row_local];
1796 }
1797 // Otherwise the entry is not stored in the local processor
1798 // and we must have haloed it
1799 else
1800 {
1802 }
1803 }
1804 // Otherwise just return the dof
1805 else
1806#endif
1807 {
1808 return this->Dof_pt[i];
1809 }
1810 }
1811
1812 /// i-th dof in the problem
1813 double& dof(const unsigned& i)
1814 {
1815 return *(Dof_pt[i]);
1816 }
1817
1818 /// i-th dof in the problem (const version)
1819 double dof(const unsigned& i) const
1820 {
1821 return *(Dof_pt[i]);
1822 }
1823
1824 /// Pointer to i-th dof in the problem
1825 double*& dof_pt(const unsigned& i)
1826 {
1827 return Dof_pt[i];
1828 }
1829
1830 /// Pointer to i-th dof in the problem (const version)
1831 double* dof_pt(const unsigned& i) const
1832 {
1833 return Dof_pt[i];
1834 }
1835
1836 /// Return the residual vector multiplied by the inverse mass matrix
1837 /// Virtual so that it can be overloaded for mpi problems
1839
1840 /// Get the time derivative of all values (using
1841 /// get_inverse_mass_matrix_times_residuals(..) with all time steppers set
1842 /// to steady) e.g. for use in explicit time steps. The approach used is
1843 /// slighty hacky, beware if you have a residual which is non-linear or
1844 /// implicit in the derivative or if you have overloaded
1845 /// get_jacobian(...).
1846 virtual void get_dvaluesdt(DoubleVector& f);
1847
1848 /// Return the fully-assembled residuals Vector for the problem:
1849 /// Virtual so it can be overloaded in for mpi problems
1850 virtual void get_residuals(DoubleVector& residuals);
1851
1852 /// Return the fully-assembled Jacobian and residuals for the problem
1853 /// Interface for the case when the Jacobian matrix is dense.
1854 /// This is virtual so, if we feel like it (e.g. for testing iterative
1855 /// solvers with specific test matrices, we can bypass the default
1856 /// assembly procedure for the Jacobian and the residual vector.
1857 virtual void get_jacobian(DoubleVector& residuals,
1858 DenseDoubleMatrix& jacobian);
1859
1860 /// Return the fully-assembled Jacobian and residuals for the
1861 /// problem. Interface for the case when the Jacobian is in row-compressed
1862 /// storage format. This is virtual so, if we feel like it (e.g. for testing
1863 /// iterative solvers with specific test matrices), we can bypass the
1864 /// default assembly procedure for the Jacobian and the residual vector.
1865 virtual void get_jacobian(DoubleVector& residuals,
1866 CRDoubleMatrix& jacobian);
1867
1868 /// Return the fully-assembled Jacobian and residuals for the
1869 /// problem. Interface for the case when the Jacobian is in
1870 /// column-compressed storage format. This is virtual so, if we feel like it
1871 /// (e.g. for testing iterative solvers with specific test matrices), we can
1872 /// bypass the default assembly procedure for the Jacobian and the residual
1873 /// vector.
1874 virtual void get_jacobian(DoubleVector& residuals,
1875 CCDoubleMatrix& jacobian);
1876
1877 /// Dummy virtual function that must be overloaded by the problem to
1878 /// specify which matrices should be summed to give the final Jacobian.
1879 virtual void get_jacobian(DoubleVector& residuals, SumOfMatrices& jacobian)
1880 {
1881 std::ostringstream error_stream;
1882 error_stream
1883 << "You must overload this function in your problem class to specify\n"
1884 << "which matrices should be summed to give the final Jacobian.";
1885 throw OomphLibError(
1886 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1887 }
1888
1889 /// Return the fully-assembled Jacobian and residuals, generated by
1890 /// finite differences
1891 void get_fd_jacobian(DoubleVector& residuals,
1892 DenseMatrix<double>& jacobian);
1893
1894 /// Get the derivative of the entire residuals vector wrt a
1895 /// global parameter, used in continuation problems
1896 void get_derivative_wrt_global_parameter(double* const& parameter_pt,
1897 DoubleVector& result);
1898
1899 /// Return the product of the global hessian (derivative of Jacobian
1900 /// matrix with respect to all variables) with
1901 /// an eigenvector, Y, and any number of other specified vectors C
1902 /// (d(J_{ij})/d u_{k}) Y_{j} C_{k}.
1903 /// This function is used in assembling and solving the augmented systems
1904 /// associated with bifurcation tracking.
1905 /// The default implementation is to use finite differences at the global
1906 /// level.
1911
1912
1913 /// Get derivative of an element in the problem wrt a global
1914 /// parameter, used in continuation problems
1915 // void get_derivative_wrt_global_parameter(double* const &parameter_pt,
1916 // GeneralisedElement* const
1917 // &elem_pt, Vector<double>
1918 // &result);
1919
1920 /// Solve an eigenproblem as assembled by EigenElements
1921 /// calculate n_eval eigenvalues and return the corresponding
1922 /// eigenvectors. The boolean flag (default true) specifies whether
1923 /// the steady jacobian should be assembled. If the flag is false
1924 /// then the weighted mass-matrix terms from the timestepper will
1925 /// be included in the jacobian --- this is almost certainly never
1926 /// wanted.
1927 void solve_eigenproblem(const unsigned& n_eval,
1928 Vector<std::complex<double>>& eigenvalue,
1929 Vector<DoubleVector>& eigenvector,
1930 const bool& steady = true);
1931
1932 /// Solve an eigenproblem as assembled by EigenElements,
1933 /// but only return the eigenvalues, not the eigenvectors.
1934 /// The boolean flag (default true) is used to specify whether the
1935 /// weighted mass-matrix terms from the timestepping scheme should
1936 /// be included in the jacobian.
1937 void solve_eigenproblem(const unsigned& n_eval,
1938 Vector<std::complex<double>>& eigenvalue,
1939 const bool& steady = true)
1940 {
1941 // Create temporary storage for the eigenvectors (potentially wasteful)
1942 Vector<DoubleVector> eigenvector;
1943 solve_eigenproblem(n_eval, eigenvalue, eigenvector, steady);
1944 }
1945
1946 /// Get the matrices required by a eigensolver. If the
1947 /// shift parameter is non-zero the second matrix will be shifted
1948 virtual void get_eigenproblem_matrices(CRDoubleMatrix& mass_matrix,
1949 CRDoubleMatrix& main_matrix,
1950 const double& shift = 0.0);
1951
1952 /// Assign the eigenvector passed to the function
1953 /// to the dofs in the problem so that it can be output by
1954 /// the usual routines.
1955 void assign_eigenvector_to_dofs(DoubleVector& eigenvector);
1956
1957 /// Add the eigenvector passed to the function scaled by
1958 /// the constat epsilon to the dofs in the problem so that
1959 /// it can be output by the usual routines
1960 void add_eigenvector_to_dofs(const double& epsilon,
1961 const DoubleVector& eigenvector);
1962
1963
1964 /// Store the current values of the degrees of freedom
1966
1967 /// Restore the stored values of the degrees of freedom
1968 void restore_dof_values();
1969
1970 /// Enable recycling of Jacobian in Newton iteration
1971 /// (if the linear solver allows it).
1972 /// Useful for linear problems with constant Jacobians or nonlinear
1973 /// problems where you're willing to risk the trade-off between
1974 /// faster solve times and degraded Newton convergence rate
1976 {
1979 }
1980
1981 /// Disable recycling of Jacobian in Newton iteration
1983 {
1986 }
1987
1988 /// Is recycling of Jacobian in Newton iteration enabled?
1990 {
1992 }
1993
1995 {
1997 }
1998
1999 /// Use Newton method to solve the problem
2000 void newton_solve();
2001
2002 /// enable globally convergent Newton method
2004 {
2006 }
2007
2008 /// disable globally convergent Newton method
2010 {
2012 }
2013
2014 private:
2015 /// Line search helper for globally convergent Newton method
2017 const Vector<double>& x_old,
2018 const double& half_residual_squared_old,
2019 DoubleVector& gradient,
2020 DoubleVector& newton_dir,
2021 double& half_residual_squared,
2022 const double& stpmax);
2023
2024 public:
2025 /// Adaptive Newton solve: up to max_adapt adaptations of all
2026 /// refineable submeshes are performed to achieve the
2027 /// the error targets specified in the refineable submeshes.
2028 void newton_solve(unsigned const& max_adapt);
2029
2030 /// Solve a steady problem using adaptive Newton's method,
2031 /// but in the context of an overall unstady problem, perhaps to
2032 /// determine an initial condition.
2033 /// This is achieved by setting the weights in the timesteppers to be zero
2034 /// which has the effect of rendering them steady timesteppers.
2035 /// The optional argument max_adapt specifies the max. number of
2036 /// adaptations of all refineable submeshes are performed to
2037 /// achieve the the error targets specified in the refineable submeshes.
2038 void steady_newton_solve(unsigned const& max_adapt = 0);
2039
2040 /// Copy Data values, nodal positions etc from specified problem.
2041 /// Note: This is not a copy constructor. We assume that the current
2042 /// and the "original" problem have both been created by calling
2043 /// the same problem constructor so that all Data objects,
2044 /// time steppers etc. in the two problems are completely independent.
2045 /// This function copies the nodal, internal and global values,
2046 /// and the time parameters from the original problem into
2047 /// "this" one. This functionality is required, e.g. for
2048 /// multigrid computations.
2049 void copy(Problem* orig_problem_pt);
2050
2051 /// Make and return a pointer to the copy of the problem. A virtual
2052 /// function that must be filled in by the user is they wish to perform
2053 /// adaptive refinement in bifurcation tracking or in multigrid problems.
2054 /// ALH: WILL NOT BE NECESSARY IN BIFURCATION TRACKING IN LONG RUN...
2055 virtual Problem* make_copy();
2056
2057 /// Read refinement pattern of all refineable meshes and refine them
2058 /// accordingly, then read all Data and nodal position info from
2059 /// file for restart. Return flag to indicate if the restart was from
2060 /// steady or unsteady solution
2061 virtual void read(std::ifstream& restart_file, bool& unsteady_restart);
2062
2063 /// Read refinement pattern of all refineable meshes and refine them
2064 /// accordingly, then read all Data and nodal position info from
2065 /// file for restart.
2066 virtual void read(std::ifstream& restart_file)
2067 {
2068 bool unsteady_restart;
2069 read(restart_file, unsteady_restart);
2070 }
2071
2072 /// Dump refinement pattern of all refineable meshes and all generic
2073 /// Problem data to file for restart.
2074 virtual void dump(std::ofstream& dump_file) const;
2075
2076 /// Dump refinement pattern of all refineable meshes and all generic
2077 /// Problem data to file for restart.
2078 void dump(const std::string& dump_file_name) const
2079 {
2080 std::ofstream dump_stream(dump_file_name.c_str());
2081#ifdef PARANOID
2082 if (!dump_stream.is_open())
2083 {
2084 std::string err = "Couldn't open file " + dump_file_name;
2085 throw OomphLibError(
2086 err, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
2087 }
2088#endif
2089 dump(dump_stream);
2090 }
2091
2092#ifdef OOMPH_HAS_MPI
2093
2094 /// Get pointers to all possible halo data indexed by global
2095 /// equation number in a map.
2096 void get_all_halo_data(std::map<unsigned, double*>& map_of_halo_data);
2097
2098 /// Classify any non-classified nodes into halo/haloed and
2099 /// synchronise equation numbers. Return the total
2100 /// number of degrees of freedom in the overall problem
2101 long synchronise_eqn_numbers(const bool& assign_local_eqn_numbers = true);
2102
2103 /// Synchronise the degrees of freedom by overwriting
2104 /// the haloed values with their non-halo counterparts held
2105 /// on other processors. Bools control if we deal with data associated with
2106 /// external halo/ed elements/nodes or the "normal" halo/ed ones.
2107 void synchronise_dofs(const bool& do_halos, const bool& do_external_halos);
2108
2109 /// Perform all required synchronisation in solvers
2110 void synchronise_all_dofs();
2111
2112 /// Check the halo/haloed node/element schemes
2113 void check_halo_schemes(DocInfo& doc_info);
2114
2115 /// Check the halo/haloed node/element schemes
2117 {
2118 DocInfo tmp_doc_info;
2119 tmp_doc_info.disable_doc();
2120 check_halo_schemes(tmp_doc_info);
2121 }
2122
2123 /// Distribute the problem and doc, using the specified partition;
2124 /// returns a vector which details the partitioning
2125 Vector<unsigned> distribute(const Vector<unsigned>& element_partition,
2126 DocInfo& doc_info,
2127 const bool& report_stats = false);
2128
2129 /// Distribute the problem; returns a vector which
2130 /// details the partitioning
2132 const bool& report_stats = false);
2133
2134 /// Distribute the problem using the specified partition;
2135 /// returns a vector which details the partitioning
2136 Vector<unsigned> distribute(const Vector<unsigned>& element_partition,
2137 const bool& report_stats = false);
2138
2139 /// Distribute the problem; returns a vector which
2140 /// details the partitioning
2141 Vector<unsigned> distribute(const bool& report_stats = false);
2142
2143 /// Partition the global mesh, return vector specifying the processor
2144 /// number for each element. Virtual so that it can be overloaded by
2145 /// any user; the default is to use METIS to perform the partitioning
2146 /// (with a bit of cleaning up afterwards to sort out "special cases").
2147 virtual void partition_global_mesh(Mesh*& global_mesh_pt,
2148 DocInfo& doc_info,
2149 Vector<unsigned>& element_domain,
2150 const bool& report_stats = false);
2151
2152 /// (Irreversibly) prune halo(ed) elements and nodes, usually
2153 /// after another round of refinement, to get rid of
2154 /// excessively wide halo layers. Note that the current
2155 /// mesh will be now regarded as the base mesh and no unrefinement
2156 /// relative to it will be possible once this function
2157 /// has been called.
2159 const bool& report_stats);
2160
2161 /// (Irreversibly) prune halo(ed) elements and nodes, usually
2162 /// after another round of refinement, to get rid of
2163 /// excessively wide halo layers. Note that the current
2164 /// mesh will be now regarded as the base mesh and no unrefinement
2165 /// relative to it will be possible once this function
2166 /// has been called.
2167 void prune_halo_elements_and_nodes(const bool& report_stats = false)
2168 {
2169 DocInfo doc_info;
2170 doc_info.disable_doc();
2171 prune_halo_elements_and_nodes(doc_info, report_stats);
2172 }
2173
2174 /// Threshold for error throwing in Problem::check_halo_schemes()
2176
2177 /// Access to Problem_has_been_distributed flag
2179 {
2181 }
2182
2183#endif
2184
2185 /// Wrapper function to delete external storage for
2186 /// each submesh of the problem
2188
2189 /// Boolean to indicate if all output is suppressed in
2190 /// Problem::newton_solve(). Defaults to false.
2192
2193
2194 protected:
2195 /// Boolean to indicate whether a Newton step should be taken
2196 /// even if the initial residuals are below the required tolerance
2198
2199 /// What it says: If temporally adaptive Newton solver fails to
2200 /// to converge, reduce timestep by this factor and try again; defaults
2201 /// to 1/2; can be over-written by user in derived problem.
2203
2204
2205 /// Boolean to decide if a timestep is to be rejected if the
2206 /// error estimate post-solve (computed by global_temporal_error_norm())
2207 /// exceeds the tolerance required in the call to
2208 /// adaptive_unsteady_newton_solve(...). Defaults to true.
2210
2211 /// Perform a basic arc-length continuation step using Newton's
2212 /// method. Returns number of Newton steps taken.
2213 unsigned newton_solve_continuation(double* const& parameter_pt);
2214
2215 /// This function performs a basic continuation step using the Newton
2216 /// method. The number of Newton steps taken is returned, to be used in any
2217 /// external step-size control routines.
2218 /// The governing parameter of the problem is passed as a pointer to the
2219 /// routine, as is a vector in which
2220 /// to store the derivatives wrt the parameter, if required.
2221 unsigned newton_solve_continuation(double* const& parameter_pt,
2222 DoubleVector& z);
2223
2224 /// A function to calculate the derivatives wrt the arc-length. This
2225 /// version of the function actually does a linear solve so that the
2226 /// derivatives
2227 /// are calculated "exactly" rather than using the values at the Newton
2228 /// step just before convergence. This is necessary in spatially adaptive
2229 /// problems, in which the number of degrees of freedom changes and so
2230 /// the appropriate derivatives must be calculated for the new variables.
2231 /// This function is called if no Newton steps were taken in the
2232 /// continuation routine ... i.e. the initial residuals were sufficiently
2233 /// small.
2234 void calculate_continuation_derivatives(double* const& parameter_pt);
2235
2236 /// A function to calculate the derivatives with respect to the
2237 /// arc-length required for continuation. The arguments is the solution
2238 /// of the linear system,
2239 /// Jz = dR/dparameter, that gives du/dparameter and the direction
2240 /// output from the newton_solve_continuation function. The derivatives
2241 /// are stored in the ContinuationParameters namespace.
2243
2244 /// A function to calculate the derivatives with respect to the
2245 /// arc-length required for continuation by finite differences, using
2246 /// the previous values of the solution. The derivatives are stored in
2247 /// the ContinuationParameters namespace.
2248 void calculate_continuation_derivatives_fd(double* const& parameter_pt);
2249
2250 /// Return a boolean flag to indicate whether the pointer
2251 /// parameter_pt refers to values stored in a Data object that
2252 /// is contained within the problem
2253 bool does_pointer_correspond_to_problem_data(double* const& parameter_pt);
2254
2255 public:
2256 /// Virtual function that is used to symmetrise the problem so that
2257 /// the current solution exactly satisfies any symmetries within the system.
2258 /// Used when adpativly solving pitchfork detection problems when small
2259 /// asymmetries in the coarse solution can be magnified
2260 /// leading to very inaccurate answers on the fine mesh.
2261 /// This is always problem-specific and must be filled in by the user
2262 /// The default issues a warning
2264
2265 /// Return pointer to the parameter that is used in the
2266 /// bifurcation detection. If we are not tracking a bifurcation then
2267 /// an error will be thrown by the AssemblyHandler
2268 double* bifurcation_parameter_pt() const;
2269
2270
2271 /// Return the eigenfunction calculated as part of a
2272 /// bifurcation tracking process. If we are not tracking a bifurcation
2273 /// then an error will be thrown by the AssemblyHandler
2275
2276
2277 /// Turn on fold tracking using the augmented system specified
2278 /// in the FoldHandler class. After a call to this function subsequent calls
2279 /// of the standard solution methods will converge to a fold (limit) point
2280 /// at a particular value of the variable addressed by parameter_pt.
2281 /// The system may not converge if the initial guess is sufficiently poor
2282 /// or, alternatively, if finite differencing is used to calculate the
2283 /// jacobian matrix in the elements. If the boolean flag block_solver is
2284 /// true (the default) then a block factorisation is used to solve the
2285 /// augmented system which is both faster and uses less memory.
2286 void activate_fold_tracking(double* const& parameter_pt,
2287 const bool& block_solve = true);
2288
2289 /// Activate generic bifurcation tracking for a single (real)
2290 /// eigenvalue where the initial guess for the eigenvector can be specified.
2291 void activate_bifurcation_tracking(double* const& parameter_pt,
2292 const DoubleVector& eigenvector,
2293 const bool& block_solve = true);
2294
2295 /// Activate generic bifurcation tracking for a single (real)
2296 /// eigenvalue where the initial guess for the eigenvector can be specified
2297 /// and the normalisation condition can also be specified.
2298 void activate_bifurcation_tracking(double* const& parameter_pt,
2299 const DoubleVector& eigenvector,
2300 const DoubleVector& normalisation,
2301 const bool& block_solve = true);
2302
2303
2304 /// Turn on pitchfork tracking using the augmented system specified
2305 /// in the PitchForkHandler class.
2306 /// After a call to this function subsequent calls
2307 /// of the standard solution methods will converge to a pitchfork
2308 /// bifurcation
2309 /// at a particular value of the variable addressed by parameter_pt.
2310 /// The symmetry that is to be broken must be specified by supplying
2311 /// a symmetry_vector(ndof). The easiest way to determine such a vector
2312 /// is to solve the associated eigenproblem \f$ Jx = \lambda M x\f$
2313 /// and pass in the eigenvector. This is not always necessary however,
2314 /// if the symmetry is easy to construct.
2315 /// The system may not converge if the initial guess is sufficiently poor
2316 /// or, alternatively, if finite differencing is used to calculate the
2317 /// jacobian matrix in the elements. If the boolean flag block_solver is
2318 /// true (the default) then a block factorisation is used to solve the
2319 /// augmented system which is both faster and requires less memory.
2320 void activate_pitchfork_tracking(double* const& parameter_pt,
2321 const DoubleVector& symmetry_vector,
2322 const bool& block_solve = true);
2323
2324 /// Turn on Hopf bifurcation
2325 /// tracking using the augmented system specified
2326 /// in the HopfHandler class. After a call to this function subsequent calls
2327 /// of the standard solution methods will converge to a Hopf bifuraction
2328 /// at a particular value of the variable addressed by parameter_pt.
2329 /// The system may not converge if the initial guess is sufficiently poor
2330 /// or, alternatively, if finite differencing is used to calculate the
2331 /// jacobian matrix in the elements.
2332 void activate_hopf_tracking(double* const& parameter_pt,
2333 const bool& block_solve = true);
2334
2335 /// Turn on Hopf bifurcation
2336 /// tracking using the augmented system specified
2337 /// in the HopfHandler class. After a call to this function subsequent calls
2338 /// of the standard solution methods will converge to a Hopf bifuraction
2339 /// at a particular value of the variable addressed by parameter_pt.
2340 /// The system may not converge if the initial guess is sufficiently poor
2341 /// or, alternatively, if finite differencing is used to calculate the
2342 /// jacobian matrix in the elements. This interface allows specification
2343 /// of an inital guess for the frequency and real and imaginary parts
2344 /// of the null vector, such as might be obtained from an eigensolve
2345 void activate_hopf_tracking(double* const& parameter_pt,
2346 const double& omega,
2347 const DoubleVector& null_real,
2348 const DoubleVector& null_imag,
2349 const bool& block_solve = true);
2350
2351 /// Deactivate all bifuraction tracking, by reseting
2352 /// the assembly handler to the default
2354 {
2356 }
2357
2358 /// Reset the system to the standard non-augemented state
2360
2361 private:
2362 /// Private helper function that actually contains the guts
2363 /// of the arc-length stepping, parameter_pt is a pointer to the
2364 /// parameter that is traded for the arc-length constraint, ds is
2365 /// the desired arc length and max_adapt is the maximum number of
2366 /// spatial adaptations. The pointer to the parameter may be changed
2367 /// if this is called from the Data-based interface
2368 double arc_length_step_solve_helper(double* const& parameter_pt,
2369 const double& ds,
2370 const unsigned& max_adapt);
2371
2372
2373 // ALH_DEVELOP
2374 protected:
2375 /// Private helper function that is used to set the appropriate
2376 /// pinned values for continuation.
2378
2379 public:
2380 /// Solve a steady problem using arc-length continuation, when the
2381 /// parameter that becomes a variable corresponding to the arc-length
2382 /// constraint equation is an external double:
2383 /// parameter_pt is a pointer to that double,
2384 /// ds is the desired arc_length and max_adapt is the maximum
2385 /// number of spatial adaptations (default zero, no adaptation).
2386 double arc_length_step_solve(double* const& parameter_pt,
2387 const double& ds,
2388 const unsigned& max_adapt = 0);
2389
2390 /// Solve a steady problem using arc-length continuation, when the
2391 /// variable corresponding to the arc-length
2392 /// constraint equation is already stored in data used in the problem:
2393 /// data_pt is a pointer to the appropriate data object,
2394 /// data_index is the index of the value that will be traded for the
2395 /// constriant,
2396 /// ds is the desired arc_length and max_adapt is the maximum
2397 /// number of spatial adaptations (default zero, no adaptation).
2398 /// Note that the value must be pinned in order for this formulation to work
2399 double arc_length_step_solve(Data* const& data_pt,
2400 const unsigned& data_index,
2401 const double& ds,
2402 const unsigned& max_adapt = 0);
2403
2404
2405 /// Reset the "internal" arc-length continuation parameters, so as
2406 /// to allow continuation in another parameter. N.B. The parameters that
2407 /// are reset are the "minimum" that are required, others should perhaps
2408 /// be reset, depending upon the application.
2410 {
2411 Theta_squared = 1.0;
2412 Sign_of_jacobian = 0;
2416 Arc_length_step_taken = false;
2417 Dof_derivative.resize(0);
2418 }
2419
2420 /// Access function for the sign of the global jacobian matrix.
2421 /// This will be set by the linear solver, if possible (direct solver).
2422 /// If not alternative methods must be used to detect bifurcations
2423 /// (solving the associated eigenproblem).
2425 {
2426 return Sign_of_jacobian;
2427 }
2428
2429 /// Take an explicit timestep of size dt and optionally shift
2430 /// any stored values of the time history
2431 void explicit_timestep(const double& dt, const bool& shift_values = true);
2432
2433 /// Advance time by dt and solve by Newton's method.
2434 /// This version always shifts time values
2435 void unsteady_newton_solve(const double& dt);
2436
2437 /// Advance time by dt and solve the system,
2438 /// using Newton's method. The boolean flag is used to control
2439 /// whether the time
2440 /// values should be shifted. If it is true the current data values will
2441 /// be shifted (copied to the locations where there are stored as
2442 /// previous timesteps) before solution.
2443 void unsteady_newton_solve(const double& dt, const bool& shift_values);
2444
2445 /// Unsteady adaptive Newton solve: up to max_adapt adaptations of
2446 /// all refineable submeshes are performed to achieve the the error targets
2447 /// specified in the refineable submeshes. If first==true, the initial
2448 /// conditions are re-assigned after the mesh adaptations. Shifting of time
2449 /// can be suppressed by overwriting the default value of shift (true).
2450 /// [Shifting must be done if first_timestep==true because we're constantly
2451 /// re-assigning the initial conditions; if first_timestep==true and
2452 /// shift==false shifting is performed anyway and a warning is issued.
2453 void unsteady_newton_solve(const double& dt,
2454 const unsigned& max_adapt,
2455 const bool& first,
2456 const bool& shift = true);
2457
2458
2459 /// Unsteady "doubly" adaptive Newton solve: Does temporal
2460 /// adaptation first, i.e. we try to do a timestep with an increment
2461 /// of dt, and adjusting dt until the solution on the given mesh satisfies
2462 /// the temporal error measure with tolerance epsilon. Following
2463 /// this, we do up to max_adapt spatial adaptions (without
2464 /// re-examining the temporal error). If first==true, the initial conditions
2465 /// are re-assigned after the mesh adaptations.
2466 /// Shifting of time can be suppressed by overwriting the
2467 /// default value of shift (true). [Shifting must be done
2468 /// if first_timestep==true because we're constantly re-assigning
2469 /// the initial conditions; if first_timestep==true and shift==false
2470 /// shifting is performed anyway and a warning is issued.
2472 const double& epsilon,
2473 const unsigned& max_adapt,
2474 const bool& first,
2475 const bool& shift = true)
2476 {
2477 // Call helper function with default setting (do re-solve after
2478 // spatial adaptation)
2479 unsigned suppress_resolve_after_spatial_adapt_flag = 0;
2481 dt,
2482 epsilon,
2483 max_adapt,
2484 suppress_resolve_after_spatial_adapt_flag,
2485 first,
2486 shift);
2487 }
2488
2489
2490 /// Unsteady "doubly" adaptive Newton solve: Does temporal
2491 /// adaptation first, i.e. we try to do a timestep with an increment
2492 /// of dt, and adjusting dt until the solution on the given mesh satisfies
2493 /// the temporal error measure with tolerance epsilon. Following
2494 /// this, we do up to max_adapt spatial adaptions (without
2495 /// re-examining the temporal error). If first==true, the initial conditions
2496 /// are re-assigned after the mesh adaptations.
2497 /// Shifting of time can be suppressed by overwriting the
2498 /// default value of shift (true). [Shifting must be done
2499 /// if first_timestep==true because we're constantly re-assigning
2500 /// the initial conditions; if first_timestep==true and shift==false
2501 /// shifting is performed anyway and a warning is issued.
2502 /// Pseudo-Boolean flag suppress_resolve_after_spatial_adapt [0: false;
2503 /// 1: true] does what it says.]
2505 const double& dt,
2506 const double& epsilon,
2507 const unsigned& max_adapt,
2508 const unsigned& suppress_resolve_after_spatial_adapt_flag,
2509 const bool& first,
2510 const bool& shift = true)
2511 {
2512 // Call helper function
2514 dt,
2515 epsilon,
2516 max_adapt,
2517 suppress_resolve_after_spatial_adapt_flag,
2518 first,
2519 shift);
2520 }
2521
2522
2523 /// Attempt to advance timestep by dt_desired. If the solution fails
2524 /// the timestep will be halved until convergence is achieved, or the
2525 /// timestep falls below NewtonSolverParameters::Minimum_time_step. The
2526 /// error control parameter epsilon represents the (approximate) desired
2527 /// magnitude of the global error at each timestep. The routine returns a
2528 /// double that is the suggested next timestep and should be passed as
2529 /// dt_desired the next time the routine is called. This version always
2530 /// shifts the time values.
2531 double adaptive_unsteady_newton_solve(const double& dt_desired,
2532 const double& epsilon);
2533
2534 /// Attempt to advance timestep by dt_desired. If the solution fails
2535 /// the timestep will be halved until convergence is achieved, or the
2536 /// timestep falls below NewtonSolverParameters::Minimum_time_step. The
2537 /// error control parameter epsilon represents the (approximate) desired
2538 /// magnitude of the global error at each timestep. The routine returns a
2539 /// double that is the suggested next timestep and should be passed as
2540 /// dt_desired the next time the routine is called.
2541 /// Once again the boolean flag, shift_values, is used to control whether
2542 /// the time values are shifted.
2543 double adaptive_unsteady_newton_solve(const double& dt_desired,
2544 const double& epsilon,
2545 const bool& shift_values);
2546
2547 /// Initialise data and nodal positions to simulate impulsive
2548 /// start from initial configuration/solution
2550
2551 /// Initialise data and nodal positions to simulate an impulsive
2552 /// start and also set the initial and previous values of dt
2553 void assign_initial_values_impulsive(const double& dt);
2554
2555 /// Calculate predictions
2556 void calculate_predictions();
2557
2558 /// Enable recycling of the mass matrix in explicit timestepping
2559 /// schemes. Useful for timestepping on fixed meshes when you want
2560 /// to avoid the linear solve phase.
2562
2563 /// Turn off recyling of the mass matrix in explicit timestepping
2564 /// schemes
2566
2567 /// Return whether the mass matrix is being reused
2569 {
2571 }
2572
2573 /// Refine refineable sub-meshes, each as many times as
2574 /// specified in the vector and rebuild problem
2575 void refine_uniformly(const Vector<unsigned>& nrefine_for_mesh)
2576 {
2577 DocInfo doc_info;
2578 doc_info.disable_doc();
2579 bool prune = false;
2580 refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2581 }
2582
2583 /// Refine refineable sub-meshes, each as many times as
2584 /// specified in the vector and rebuild problem; doc refinement process
2585 void refine_uniformly(const Vector<unsigned>& nrefine_for_mesh,
2586 DocInfo& doc_info)
2587 {
2588 bool prune = false;
2589 refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2590 }
2591
2592 /// Refine refineable sub-meshes, each as many times as
2593 /// specified in the vector and rebuild problem. Prune after
2594 /// refinements
2595 void refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh)
2596 {
2597 DocInfo doc_info;
2598 doc_info.disable_doc();
2599 bool prune = true;
2600 refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2601 }
2602
2603 /// Refine refineable sub-meshes, each as many times as
2604 /// specified in the vector and rebuild problem; doc refinement process
2605 void refine_uniformly_and_prune(const Vector<unsigned>& nrefine_for_mesh,
2606 DocInfo& doc_info)
2607 {
2608 bool prune = true;
2609 refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2610 }
2611
2612 /// Refine (all) refineable (sub)mesh(es) uniformly and
2613 /// rebuild problem; doc refinement process.
2615 {
2616 // Number of (sub)meshes
2617 unsigned nmesh = std::max(unsigned(1), nsub_mesh());
2618
2619 // Refine each mesh once
2620 Vector<unsigned> nrefine_for_mesh(nmesh, 1);
2621 refine_uniformly(nrefine_for_mesh);
2622 }
2623
2624
2625 /// Refine (all) refineable (sub)mesh(es) uniformly and
2626 /// rebuild problem; doc refinement process.
2628 {
2629 // Number of (sub)meshes
2630 unsigned nmesh = std::max(unsigned(1), nsub_mesh());
2631
2632 // Refine each mesh once
2633 Vector<unsigned> nrefine_for_mesh(nmesh, 1);
2634 refine_uniformly_and_prune(nrefine_for_mesh);
2635 }
2636
2637
2638 /// Refine (all) refineable (sub)mesh(es) uniformly and
2639 /// rebuild problem
2641 {
2642 DocInfo doc_info;
2643 doc_info.disable_doc();
2644 refine_uniformly(doc_info);
2645 }
2646
2647 /// Do uniform refinement for submesh i_mesh with documentation
2648 void refine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
2649
2650 /// Do uniform refinement for submesh i_mesh without documentation
2651 void refine_uniformly(const unsigned& i_mesh)
2652 {
2653 DocInfo doc_info;
2654 doc_info.disable_doc();
2655 refine_uniformly(i_mesh, doc_info);
2656 }
2657
2658 /// p-refine p-refineable sub-meshes, each as many times as
2659 /// specified in the vector and rebuild problem
2660 void p_refine_uniformly(const Vector<unsigned>& nrefine_for_mesh)
2661 {
2662 DocInfo doc_info;
2663 doc_info.disable_doc();
2664 bool prune = false;
2665 p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2666 }
2667
2668 /// p-refine p-refineable sub-meshes, each as many times as
2669 /// specified in the vector and rebuild problem; doc refinement process
2670 void p_refine_uniformly(const Vector<unsigned>& nrefine_for_mesh,
2671 DocInfo& doc_info)
2672 {
2673 bool prune = false;
2674 p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2675 }
2676
2677 /// p-refine p-refineable sub-meshes, each as many times as
2678 /// specified in the vector and rebuild problem. Prune after
2679 /// refinements
2681 {
2682 // Not tested:
2683 throw OomphLibWarning("This functionality has not yet been tested.",
2684 "Problem::p_refine_uniformly_and_prune()",
2685 OOMPH_EXCEPTION_LOCATION);
2686 DocInfo doc_info;
2687 doc_info.disable_doc();
2688 bool prune = true;
2689 p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2690 }
2691
2692 /// p-refine p-refineable sub-meshes, each as many times as
2693 /// specified in the vector and rebuild problem; doc refinement process
2695 DocInfo& doc_info)
2696 {
2697 // Not tested:
2698 throw OomphLibWarning("This functionality has not yet been tested.",
2699 "Problem::p_refine_uniformly_and_prune()",
2700 OOMPH_EXCEPTION_LOCATION);
2701 bool prune = true;
2702 p_refine_uniformly_aux(nrefine_for_mesh, doc_info, prune);
2703 }
2704
2705 /// p-refine (all) p-refineable (sub)mesh(es) uniformly and
2706 /// rebuild problem; doc refinement process.
2708 {
2709 // Number of (sub)meshes
2710 unsigned nmesh = std::max(unsigned(1), nsub_mesh());
2711
2712 // Refine each mesh once
2713 Vector<unsigned> nrefine_for_mesh(nmesh, 1);
2714 p_refine_uniformly(nrefine_for_mesh);
2715 }
2716
2717
2718 /// p-refine (all) p-refineable (sub)mesh(es) uniformly and
2719 /// rebuild problem; doc refinement process.
2721 {
2722 // Not tested:
2723 throw OomphLibWarning("This functionality has not yet been tested.",
2724 "Problem::p_refine_uniformly_and_prune()",
2725 OOMPH_EXCEPTION_LOCATION);
2726 // Number of (sub)meshes
2727 unsigned nmesh = std::max(unsigned(1), nsub_mesh());
2728
2729 // Refine each mesh once
2730 Vector<unsigned> nrefine_for_mesh(nmesh, 1);
2731 p_refine_uniformly_and_prune(nrefine_for_mesh);
2732 }
2733
2734
2735 /// p-refine (all) p-refineable (sub)mesh(es) uniformly and
2736 /// rebuild problem
2738 {
2739 DocInfo doc_info;
2740 doc_info.disable_doc();
2741 p_refine_uniformly(doc_info);
2742 }
2743
2744 /// Do uniform p-refinement for submesh i_mesh with documentation
2745 void p_refine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
2746
2747 /// Do uniform p-refinement for submesh i_mesh without documentation
2748 void p_refine_uniformly(const unsigned& i_mesh)
2749 {
2750 DocInfo doc_info;
2751 doc_info.disable_doc();
2752 p_refine_uniformly(i_mesh, doc_info);
2753 }
2754
2755 /// Refine (one and only!) mesh by splitting the elements identified
2756 /// by their numbers relative to the problems' only mesh, then rebuild
2757 /// the problem.
2759 const Vector<unsigned>& elements_to_be_refined);
2760
2761
2762 /// Refine (one and only!) mesh by splitting the elements identified
2763 /// by their pointers, then rebuild the problem.
2765 const Vector<RefineableElement*>& elements_to_be_refined_pt);
2766
2767 /// Refine specified submesh by splitting the elements identified
2768 /// by their numbers relative to the submesh, then rebuild the problem.
2770 const unsigned& i_mesh, const Vector<unsigned>& elements_to_be_refined);
2771
2772 /// Refine specified submesh by splitting the elements identified
2773 /// by their pointers, then rebuild the problem.
2775 const unsigned& i_mesh,
2776 const Vector<RefineableElement*>& elements_to_be_refined_pt);
2777
2778 /// Refine all submeshes by splitting the elements identified
2779 /// by their numbers relative to each submesh in a Vector of Vectors,
2780 /// then rebuild the problem.
2782 const Vector<Vector<unsigned>>& elements_to_be_refined);
2783
2784 /// Refine all submeshes by splitting the elements identified
2785 /// by their pointers within each submesh in a Vector of Vectors,
2786 /// then rebuild the problem.
2788 const Vector<Vector<RefineableElement*>>& elements_to_be_refined_pt);
2789
2790 /// p-refine (one and only!) mesh by refining the elements identified
2791 /// by their numbers relative to the problems' only mesh, then rebuild
2792 /// the problem.
2794 const Vector<unsigned>& elements_to_be_refined);
2795
2796
2797 /// p-refine (one and only!) mesh by refining the elements identified
2798 /// by their pointers, then rebuild the problem.
2800 const Vector<PRefineableElement*>& elements_to_be_refined_pt);
2801
2802 /// p-refine specified submesh by refining the elements identified
2803 /// by their numbers relative to the submesh, then rebuild the problem.
2805 const unsigned& i_mesh, const Vector<unsigned>& elements_to_be_refined);
2806
2807 /// p-refine specified submesh by refining the elements identified
2808 /// by their pointers, then rebuild the problem.
2810 const unsigned& i_mesh,
2811 const Vector<PRefineableElement*>& elements_to_be_refined_pt);
2812
2813 /// p-refine all submeshes by refining the elements identified
2814 /// by their numbers relative to each submesh in a Vector of Vectors,
2815 /// then rebuild the problem.
2817 const Vector<Vector<unsigned>>& elements_to_be_refined);
2818
2819 /// p-refine all submeshes by refining the elements identified
2820 /// by their pointers within each submesh in a Vector of Vectors,
2821 /// then rebuild the problem.
2823 const Vector<Vector<PRefineableElement*>>& elements_to_be_refined_pt);
2824
2825 /// Refine (all) refineable (sub)mesh(es) uniformly and
2826 /// rebuild problem. Return 0 for success,
2827 /// 1 for failure (if unrefinement has reached the coarsest permitted
2828 /// level)
2829 unsigned unrefine_uniformly();
2830
2831 /// Do uniform refinement for submesh i_mesh without documentation.
2832 /// Return 0 for success, 1 for failure (if unrefinement has reached
2833 /// the coarsest permitted level)
2834 unsigned unrefine_uniformly(const unsigned& i_mesh);
2835
2836 /// p-unrefine (all) p-refineable (sub)mesh(es) uniformly and
2837 /// rebuild problem.
2838 void p_unrefine_uniformly(DocInfo& doc_info);
2839
2840 /// Do uniform p-unrefinement for submesh i_mesh without documentation.
2841 void p_unrefine_uniformly(const unsigned& i_mesh, DocInfo& doc_info);
2842
2843 /// Adapt problem:
2844 /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2845 /// based on their own error estimates and the target errors specified
2846 /// in the mesh(es). Following mesh adaptation,
2847 /// update global mesh, and re-assign equation numbers.
2848 /// Return # of refined/unrefined elements. On return from this
2849 /// function, Problem can immediately be solved again.
2850 void adapt(unsigned& n_refined, unsigned& n_unrefined);
2851
2852 /// Adapt problem:
2853 /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2854 /// based on their own error estimates and the target errors specified
2855 /// in the mesh(es). Following mesh adaptation,
2856 /// update global mesh, and re-assign equation numbers.
2857 /// On return from this function, Problem can immediately be solved again.
2858 /// [Argument-free wrapper]
2859 void adapt()
2860 {
2861 unsigned n_refined, n_unrefined;
2862 adapt(n_refined, n_unrefined);
2863 }
2864
2865 /// p-adapt problem:
2866 /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2867 /// based on their own error estimates and the target errors specified
2868 /// in the mesh(es). Following mesh adaptation,
2869 /// update global mesh, and re-assign equation numbers.
2870 /// Return # of refined/unrefined elements. On return from this
2871 /// function, Problem can immediately be solved again.
2872 void p_adapt(unsigned& n_refined, unsigned& n_unrefined);
2873
2874 /// p-adapt problem:
2875 /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2876 /// based on their own error estimates and the target errors specified
2877 /// in the mesh(es). Following mesh adaptation,
2878 /// update global mesh, and re-assign equation numbers.
2879 /// On return from this function, Problem can immediately be solved again.
2880 /// [Argument-free wrapper]
2881 void p_adapt()
2882 {
2883 unsigned n_refined, n_unrefined;
2884 p_adapt(n_refined, n_unrefined);
2885 }
2886
2887
2888 /// Adapt problem:
2889 /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2890 /// based on the error estimates in elemental_error
2891 /// and the target errors specified
2892 /// in the mesh(es). Following mesh adaptation,
2893 /// update global mesh, and re-assign equation numbers.
2894 /// Return # of refined/unrefined elements. On return from this
2895 /// function, Problem can immediately be solved again.
2897 unsigned& n_refined,
2898 unsigned& n_unrefined,
2899 Vector<Vector<double>>& elemental_error);
2900
2901 /// Adapt problem:
2902 /// Perform mesh adaptation for (all) refineable (sub)mesh(es),
2903 /// based on the error estimates in elemental_error
2904 /// and the target errors specified
2905 /// in the mesh(es). Following mesh adaptation,
2906 /// update global mesh, and re-assign equation numbers.
2907 /// Return # of refined/unrefined elements. On return from this
2908 /// function, Problem can immediately be solved again.
2909 /// [Wrapper without n_refined and n_unrefined arguments]
2911 {
2912 unsigned n_refined, n_unrefined;
2913 adapt_based_on_error_estimates(n_refined, n_unrefined, elemental_error);
2914 }
2915
2916
2917 /// Return the error estimates computed by (all) refineable
2918 /// (sub)mesh(es) in the elemental_error structure, which consists of
2919 /// a vector of vectors of elemental errors, one vector for each (sub)mesh.
2920 void get_all_error_estimates(Vector<Vector<double>>& elemental_error);
2921
2922 /// Get max and min error for all elements in submeshes
2923 void doc_errors(DocInfo& doc_info);
2924
2925 /// Get max and min error for all elements in submeshes
2927 {
2928 DocInfo tmp_doc_info;
2929 tmp_doc_info.disable_doc();
2930 doc_errors(tmp_doc_info);
2931 }
2932
2933 /// Enable the output of information when in the newton solver
2934 /// (Default)
2936 {
2938 }
2939
2940 /// Disable the output of information when in the newton solver
2942 {
2944 }
2945 };
2946
2947
2948 //=======================================================================
2949 /// A class to handle errors in the Newton solver
2950 //=======================================================================
2952 {
2953 public:
2954 /// Error in the linear solver
2956
2957 /// Max. # of iterations performed when the Newton solver died
2958 unsigned iterations;
2959
2960 /// Max. residual when Newton solver died
2961 double maxres;
2962
2963 /// Default constructor, does nothing
2965 {
2966 }
2967
2968 /// Constructor that passes a failure of the linear solver
2969 NewtonSolverError(const bool& Passed_linear_failure)
2970 : linear_solver_error(Passed_linear_failure), iterations(0), maxres(0.0)
2971 {
2972 }
2973
2974 /// Constructor that passes number of iterations and residuals
2975 NewtonSolverError(unsigned Passed_iterations, double Passed_maxres)
2976 : linear_solver_error(false),
2977 iterations(Passed_iterations),
2978 maxres(Passed_maxres)
2979 {
2980 }
2981
2982 /* /// Broken copy constructor */
2983 /* NewtonSolverError(const NewtonSolverError& dummy) = delete; */
2984
2985 /* /// Broken assignment operator */
2986 /* void operator=(const NewtonSolverError&) = delete; */
2987 };
2988
2989
2990} // namespace oomph
2991
2992
2993#endif
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...
A custom linear solver class that is used to solve a block-factorised version of the Fold bifurcation...
A custom linear solver class that is used to solve a block-factorised version of the PitchFork bifurc...
A custom linear solver class that is used to solve a block-factorised version of the Hopf bifurcation...
A custom linear solver class that is used to solve a block-factorised version of the PitchFork bifurc...
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
Definition: matrices.h:2791
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
GeneralisedTimestepper used to store the arclength derivatives and pervious solutions required in con...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
Information for documentation of results: Directory and file number to enable output in the form RESL...
void disable_doc()
Disable documentation.
A class that stores the halo/haloed entries required when using a DoubleVectorWithHaloEntries....
LinearAlgebraDistribution *& distribution_pt()
Return the pointer to the distirbution used to setup the halo information.
unsigned local_index(const unsigned &global_eqn)
Return the local index associated with the global equation.
===================================================================== An extension of DoubleVector th...
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
Base class for all EigenProblem solves. This simply defines standard interfaces so that different sol...
Definition: eigen_solver.h:61
Class for objects than can be advanced in time by an Explicit Timestepper. WARNING: For explicit time...
A Base class for explicit timesteppers.
A class that is used to assemble the augmented system that defines a fold (saddle-node) or limit poin...
A class that is used to assemble the augmented system that defines a Hopf bifurcation....
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
Base class for all linear solvers. This merely defines standard interfaces for linear solvers,...
Definition: linear_solver.h:68
A general mesh class.
Definition: mesh.h:67
A class to handle errors in the Newton solver.
Definition: problem.h:2952
NewtonSolverError(unsigned Passed_iterations, double Passed_maxres)
Constructor that passes number of iterations and residuals.
Definition: problem.h:2975
bool linear_solver_error
Error in the linear solver.
Definition: problem.h:2955
NewtonSolverError()
Default constructor, does nothing.
Definition: problem.h:2964
unsigned iterations
Max. # of iterations performed when the Newton solver died.
Definition: problem.h:2958
NewtonSolverError(const bool &Passed_linear_failure)
Constructor that passes a failure of the linear solver.
Definition: problem.h:2969
double maxres
Max. residual when Newton solver died.
Definition: problem.h:2961
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
std::ostream *& stream_pt()
Access function for the stream pointer.
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
A class that is used to assemble and solve the augmented system of equations associated with calculat...
A class that is used to assemble the augmented system that defines a pitchfork (symmetry-breaking) bi...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
virtual void actions_after_implicit_timestep()
Actions that should be performed after each implicit time step. This is needed when one wants to solv...
Definition: problem.h:1070
bool Always_take_one_newton_step
Boolean to indicate whether a Newton step should be taken even if the initial residuals are below the...
Definition: problem.h:2197
bool Jacobian_reuse_is_enabled
Is re-use of Jacobian in Newton iteration enabled? Default: false.
Definition: problem.h:618
virtual void actions_after_newton_solve()
Any actions that are to be performed after a complete Newton solve, e.g. post processing....
Definition: problem.h:1038
void parallel_sparse_assemble(const LinearAlgebraDistribution *const &dist_pt, Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residuals)
Helper method to assemble CRDoubleMatrices from distributed on multiple processors.
Definition: problem.cc:6556
void describe_dofs(std::ostream &out= *(oomph_info.stream_pt())) const
Function to describe the dofs in terms of the global equation number, i.e. what type of value (nodal ...
Definition: problem.cc:2359
virtual void sparse_assemble_row_or_column_compressed(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Protected helper function that is used to assemble the Jacobian matrix in the case when the storage i...
Definition: problem.cc:4484
double * global_dof_pt(const unsigned &i)
Return a pointer to the dof, indexed by global equation number which may be haloed or stored locally....
Definition: problem.h:1764
bool Store_local_dof_pt_in_elements
Boolean to indicate whether local dof pointers should be stored in the elements.
Definition: problem.h:226
void clear_elemental_assembly_time()
Clear storage of most-recent elemental assembly times (used for load balancing). Next load balancing ...
Definition: problem.h:872
virtual void actions_before_newton_step()
Any actions that are to be performed before each individual Newton step. Most likely to be used for d...
Definition: problem.h:1053
void remove_duplicate_data(Mesh *const &mesh_pt, bool &actually_removed_some_data)
Private helper function to remove repeated data in external haloed elements in specified mesh....
Definition: problem.cc:2568
bool Bifurcation_detection
Boolean to control bifurcation detection via determinant of Jacobian.
Definition: problem.h:790
void get_flat_packed_refinement_pattern_for_load_balancing(const Vector< unsigned > &old_domain_for_base_element, const Vector< unsigned > &new_domain_for_base_element, const unsigned &max_refinement_level_overall, std::map< unsigned, Vector< unsigned > > &flat_packed_refinement_info_for_root)
Get flat-packed refinement pattern for each root element in current mesh (labeled by unique number of...
Definition: problem.cc:19835
void p_refine_uniformly(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem; ...
Definition: problem.h:2670
void adapt()
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on their own error e...
Definition: problem.h:2859
virtual void actions_before_newton_solve()
Any actions that are to be performed before a complete Newton solve (e.g. adjust boundary conditions)...
Definition: problem.h:1032
Vector< unsigned > First_el_for_assembly
First element to be assembled by given processor for non-distributed problem (only kept up to date wh...
Definition: problem.h:519
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Add a submesh to the problem and return its number, i, by which it can be accessed via mesh_pt(i).
Definition: problem.h:1330
virtual void get_jacobian(DoubleVector &residuals, SumOfMatrices &jacobian)
Dummy virtual function that must be overloaded by the problem to specify which matrices should be sum...
Definition: problem.h:1879
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1989
void refine_uniformly_aux(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info, const bool &prune)
Helper function to do compund refinement of (all) refineable (sub)mesh(es) uniformly as many times as...
Definition: problem.cc:15445
void build_global_mesh()
Build the global mesh by combining the all the submeshes. Note: The nodes boundary information refers...
Definition: problem.cc:1493
unsigned unrefine_uniformly()
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem. Return 0 for success,...
Definition: problem.cc:15831
void assign_initial_values_impulsive()
Initialise data and nodal positions to simulate impulsive start from initial configuration/solution.
Definition: problem.cc:11499
void unset_default_partition_in_load_balance()
Do not use the default partition in the load balance.
Definition: problem.h:1451
double Theta_squared
Value of the scaling parameter required so that the parameter occupies the desired proportion of the ...
Definition: problem.h:742
void refine_uniformly(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem; doc ...
Definition: problem.h:2585
virtual void get_eigenproblem_matrices(CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
Get the matrices required by a eigensolver. If the shift parameter is non-zero the second matrix will...
Definition: problem.cc:8368
void refine_uniformly(const Vector< unsigned > &nrefine_for_mesh)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem.
Definition: problem.h:2575
Vector< double > Dof_derivative
Storage for the derivative of the problem variables wrt arc-length.
Definition: problem.h:774
double & dof_current(const unsigned &i)
Access function to the current value of the i-th (local) dof at the start of a continuation step.
Definition: problem.h:1182
virtual void sparse_assemble_row_or_column_compressed_with_lists(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:4930
virtual void sparse_assemble_row_or_column_compressed_with_two_arrays(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:6061
void synchronise_dofs(const bool &do_halos, const bool &do_external_halos)
Synchronise the degrees of freedom by overwriting the haloed values with their non-halo counterparts ...
Definition: problem.cc:16462
virtual void actions_before_distribute()
Actions to be performed before a (mesh) distribution.
Definition: problem.h:1116
LinearSolver *const & linear_solver_pt() const
Return a pointer to the linear solver object (const version)
Definition: problem.h:1472
void load_balance(const bool &report_stats)
Balance the load of a (possibly non-uniformly refined) problem that has already been distributed,...
Definition: problem.h:1403
Distributed_problem_matrix_distribution Dist_problem_matrix_distribution
The distributed matrix distribution method 1 - Automatic - the Problem distribution is employed,...
Definition: problem.h:956
DoubleVectorHaloScheme * Halo_scheme_pt
Pointer to the halo scheme for any global vectors that have the Dof_distribution.
Definition: problem.h:577
virtual void actions_after_change_in_global_parameter(double *const &parameter_pt)
Actions that are to be performed when the global parameter addressed by parameter_pt has been changed...
Definition: problem.h:1133
double dof(const unsigned &i) const
i-th dof in the problem (const version)
Definition: problem.h:1819
void p_refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
p-refine (one and only!) mesh by refining the elements identified by their numbers relative to the pr...
Definition: problem.cc:15159
void setup_dof_halo_scheme()
Function that is used to setup the halo scheme.
Definition: problem.cc:300
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1674
void p_unrefine_uniformly(DocInfo &doc_info)
p-unrefine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem.
Definition: problem.cc:15960
virtual void build_mesh()
Function to build the Problem's base mesh; this must be supplied by the user if they wish to use the ...
Definition: problem.h:1365
void(* SpatialErrorEstimatorWithDocFctPt)(Mesh *&mesh_pt, Vector< double > &elemental_error, DocInfo &doc_info)
Function pointer for spatial error estimator with doc.
Definition: problem.h:1262
static bool Suppress_warning_about_actions_before_read_unstructured_meshes
Flag to allow suppression of warning messages re reading in unstructured meshes during restart.
Definition: problem.h:317
void get_first_and_last_element_for_assembly(Vector< unsigned > &first_el_for_assembly, Vector< unsigned > &last_el_for_assembly) const
Get first and last elements for parallel assembly of non-distributed problem.
Definition: problem.h:1004
OomphCommunicator * communicator_pt()
access function to the oomph-lib communicator
Definition: problem.h:1246
void set_default_partition_in_load_balance()
Set the use of the default partition in the load balance.
Definition: problem.h:1445
bool Use_default_partition_in_load_balance
Flag to use "default partition" during load balance. Should only be set to true when run in validatio...
Definition: problem.h:514
void bifurcation_adapt_helper(unsigned &n_refined, unsigned &n_unrefined, const unsigned &bifurcation_type, const bool &actually_adapt=true)
A function that is used to adapt a bifurcation-tracking problem, which requires separate interpolatio...
Definition: problem.cc:13354
void adapt_based_on_error_estimates(unsigned &n_refined, unsigned &n_unrefined, Vector< Vector< double > > &elemental_error)
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on the error estimat...
Definition: problem.cc:14528
Vector< double > Elemental_assembly_time
Storage for assembly times (used for load balancing)
Definition: problem.h:573
double & dof(const unsigned &i)
i-th dof in the problem
Definition: problem.h:1813
bool Jacobian_has_been_computed
Has a Jacobian been computed (and can therefore be re-used if required)? Default: false.
Definition: problem.h:622
bool Bypass_increase_in_dof_check_during_pruning
Boolean to bypass check of increase in dofs during pruning.
Definition: problem.h:967
virtual void get_dvaluesdt(DoubleVector &f)
Get the time derivative of all values (using get_inverse_mass_matrix_times_residuals(....
Definition: problem.cc:3684
bool & time_adaptive_newton_crash_on_solve_fail()
Access function for Time_adaptive_newton_crash_on_solve_fail.
Definition: problem.h:1614
void enable_store_local_dof_pt_in_elements()
Insist that local dof pointers are set up in each element when equation numbering takes place.
Definition: problem.h:1696
void enable_info_in_newton_solve()
Enable the output of information when in the newton solver (Default)
Definition: problem.h:2935
friend class BlockFoldLinearSolver
Definition: problem.h:158
Problem(const Problem &dummy)=delete
Broken copy constructor.
void p_refine_uniformly(const unsigned &i_mesh)
Do uniform p-refinement for submesh i_mesh without documentation.
Definition: problem.h:2748
AssemblyHandler *const & assembly_handler_pt() const
Return a pointer to the assembly handler object (const version)
Definition: problem.h:1576
void initialise_dt(const double &dt)
Set all timesteps to the same value, dt, and assign weights for all timesteppers in the problem.
Definition: problem.cc:13231
void add_to_dofs(const double &lambda, const DoubleVector &increment_dofs)
Add lambda x incremenet_dofs[l] to the l-th dof.
Definition: problem.cc:3564
void p_refine_uniformly(DocInfo &doc_info)
p-refine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process.
Definition: problem.h:2707
double Timestep_reduction_factor_after_nonconvergence
What it says: If temporally adaptive Newton solver fails to to converge, reduce timestep by this fact...
Definition: problem.h:2202
Vector< double > Dof_current
Storage for the present values of the variables.
Definition: problem.h:777
virtual void debug_hook_fct(const unsigned &i)
Hook for debugging. Can be overloaded in driver code; argument allows identification of where we're c...
Definition: problem.h:249
double Desired_proportion_of_arc_length
Proportion of the arc-length to taken by the parameter.
Definition: problem.h:734
Mesh *& mesh_pt(const unsigned &imesh)
Return a pointer to the i-th submesh. If there are no submeshes, the 0-th submesh is the global mesh ...
Definition: problem.h:1293
virtual void actions_after_implicit_timestep_and_error_estimation()
Actions that should be performed after each implicit time step. This is needed if your actions_after_...
Definition: problem.h:1075
void refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem....
Definition: problem.h:2595
void disable_mass_matrix_reuse()
Turn off recyling of the mass matrix in explicit timestepping schemes.
Definition: problem.cc:11832
EigenSolver * Default_eigen_solver_pt
Pointer to the default eigensolver.
Definition: problem.h:191
void p_refine_uniformly(const Vector< unsigned > &nrefine_for_mesh)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem
Definition: problem.h:2660
unsigned Nnewton_iter_taken
Actual number of Newton iterations taken during the most recent iteration.
Definition: problem.h:603
double * bifurcation_parameter_pt() const
Return pointer to the parameter that is used in the bifurcation detection. If we are not tracking a b...
Definition: problem.cc:10081
void copy_haloed_eqn_numbers_helper(const bool &do_halos, const bool &do_external_halos)
A private helper function to copy the haloed equation numbers into the halo equation numbers,...
Definition: problem.cc:16951
virtual void sparse_assemble_row_or_column_compressed_with_vectors_of_pairs(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:5345
void send_data_to_be_sent_during_load_balancing(Vector< int > &send_n, Vector< double > &send_data, Vector< int > &send_displacement)
Load balance helper routine: Send data to other processors during load balancing.
Definition: problem.cc:18949
bool Time_adaptive_newton_crash_on_solve_fail
Bool to specify what to do if a Newton solve fails within a time adaptive solve. Default (false) is t...
Definition: problem.h:615
Problem()
Constructor: Allocate space for one time stepper and set all pointers to NULL and set defaults for al...
Definition: problem.cc:69
void enable_problem_distributed()
Enable problem distributed.
Definition: problem.h:971
bool is_dparameter_calculated_analytically(double *const &parameter_pt)
Function to determine whether the parameter derivatives are calculated analytically.
Definition: problem.h:277
void flush_sub_meshes()
Flush the problem's collection of sub-meshes. Must be followed by call to rebuild_global_mesh().
Definition: problem.h:1339
EigenSolver *& eigen_solver_pt()
Return a pointer to the eigen solver object.
Definition: problem.h:1492
unsigned Sparse_assembly_method
Flag to determine which sparse assembly method to use By default we use assembly by vectors of pairs.
Definition: problem.h:641
bool & use_predictor_values_as_initial_guess()
Definition: problem.h:1994
void set_analytic_hessian_products()
Function to turn on analytic calculation of the parameter derivatives in continuation and bifurcation...
Definition: problem.h:289
void calculate_predictions()
Calculate predictions.
Definition: problem.cc:11656
const TimeStepper * time_stepper_pt() const
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1537
Vector< unsigned > Last_el_plus_one_for_assembly
Last element (plus one) to be assembled by given processor for non-distributed problem (only kept up ...
Definition: problem.h:524
const OomphCommunicator * communicator_pt() const
access function to the oomph-lib communicator, const version
Definition: problem.h:1252
virtual void actions_after_read_unstructured_meshes()
Actions that are to be performed before reading in restart data for problems involving unstructured b...
Definition: problem.h:1109
void copy(Problem *orig_problem_pt)
Copy Data values, nodal positions etc from specified problem. Note: This is not a copy constructor....
Definition: problem.cc:11865
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3890
void set_explicit_time_stepper_pt(ExplicitTimeStepper *const &explicit_time_stepper_pt)
Set the explicit timestepper for the problem. The function will automatically create or resize the Ti...
Definition: problem.cc:1586
virtual void actions_after_distribute()
Actions to be performed after a (mesh) distribution.
Definition: problem.h:1119
double & minimum_dt()
Access function to min timestep in adaptive timestepping.
Definition: problem.h:1582
virtual void actions_before_implicit_timestep()
Actions that should be performed before each implicit time step. This is needed when one wants to sol...
Definition: problem.h:1064
void p_refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem....
Definition: problem.h:2680
LinearSolver * Linear_solver_pt
Pointer to the linear solver for the problem.
Definition: problem.h:173
virtual void actions_after_explicit_timestep()
Actions that should be performed after each explicit time step.
Definition: problem.h:1081
void disable_jacobian_reuse()
Disable recycling of Jacobian in Newton iteration.
Definition: problem.h:1982
bool Doc_time_in_distribute
Protected boolean flag to provide comprehensive timimings during problem distribution....
Definition: problem.h:637
unsigned Max_newton_iterations
Maximum number of Newton iterations.
Definition: problem.h:599
EigenSolver *const & eigen_solver_pt() const
Return a pointer to the eigen solver object (const version)
Definition: problem.h:1498
Vector< Vector< unsigned > > Sparse_assemble_with_arrays_previous_allocation
the number of elements in each row of a compressed matrix in the previous matrix assembly.
Definition: problem.h:667
virtual void actions_after_parameter_increase(double *const &parameter_pt)
Empty virtual function; provides hook to perform actions after the increase in the arclength paramete...
Definition: problem.h:1161
virtual void read(std::ifstream &restart_file)
Read refinement pattern of all refineable meshes and refine them accordingly, then read all Data and ...
Definition: problem.h:2066
void enable_doc_imbalance_in_parallel_assembly()
Enable doc of load imbalance in parallel assembly of distributed problem.
Definition: problem.h:849
unsigned Dof_current_offset
Storage for the offset for the current values of dofs from the original dofs (default is 2,...
Definition: problem.h:771
void calculate_continuation_derivatives(double *const &parameter_pt)
A function to calculate the derivatives wrt the arc-length. This version of the function actually doe...
Definition: problem.cc:9715
unsigned & max_newton_iterations()
Access function to max Newton iterations before giving up.
Definition: problem.h:1594
void store_current_dof_values()
Store the current values of the degrees of freedom.
Definition: problem.cc:8607
double & dof_derivative(const unsigned &i)
Access function to the derivative of the i-th (local) dof with respect to the arc length,...
Definition: problem.h:1168
bool Problem_has_been_distributed
Has the problem been distributed amongst multiple processors?
Definition: problem.h:964
bool jacobian_reuse_is_enabled()
Is recycling of Jacobian in Newton iteration enabled?
Definition: problem.h:1989
void synchronise_all_dofs()
Perform all required synchronisation in solvers.
Definition: problem.cc:16439
void get_all_error_estimates(Vector< Vector< double > > &elemental_error)
Return the error estimates computed by (all) refineable (sub)mesh(es) in the elemental_error structur...
Definition: problem.cc:14627
virtual void actions_after_change_in_bifurcation_parameter()
Actions that are to be performed after a change in the parameter that is being varied as part of the ...
Definition: problem.h:1151
bool Empty_actions_before_read_unstructured_meshes_has_been_called
Boolean to indicate that empty actions_before_read_unstructured_meshes() function has been called.
Definition: problem.h:218
bool Mass_matrix_has_been_computed
Has the mass matrix been computed (and can therefore be reused) Default: false.
Definition: problem.h:695
unsigned nglobal_data() const
Return the number of global data values.
Definition: problem.h:1686
virtual void actions_before_adapt()
Actions that are to be performed before a mesh adaptation. These might include removing any additiona...
Definition: problem.h:1022
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8783
void(* SpatialErrorEstimatorFctPt)(Mesh *&mesh_pt, Vector< double > &elemental_error)
Function pointer for spatial error estimator.
Definition: problem.h:1258
bool First_jacobian_sign_change
Boolean to indicate whether a sign change has occured in the Jacobian.
Definition: problem.h:796
void reset_arc_length_parameters()
Reset the "internal" arc-length continuation parameters, so as to allow continuation in another param...
Definition: problem.h:2409
void calculate_continuation_derivatives_fd_helper(double *const &parameter_pt)
A function that performs the guts of the continuation derivative calculation in arc-length continuati...
Definition: problem.cc:10004
double Continuation_direction
The direction of the change in parameter that will ensure that a branch is followed in one direction ...
Definition: problem.h:749
unsigned Parallel_sparse_assemble_previous_allocation
The amount of data allocated during the previous parallel sparse assemble. Used to optimise the next ...
Definition: problem.h:960
void enable_mass_matrix_reuse()
Enable recycling of the mass matrix in explicit timestepping schemes. Useful for timestepping on fixe...
Definition: problem.cc:11807
virtual void actions_before_explicit_timestep()
Actions that should be performed before each explicit time step.
Definition: problem.h:1078
void steady_newton_solve(unsigned const &max_adapt=0)
Solve a steady problem using adaptive Newton's method, but in the context of an overall unstady probl...
Definition: problem.cc:9292
void add_global_data(Data *const &global_data_pt)
Add Data to the Problem's global data – the Problem will perform equation numbering etc....
Definition: problem.h:1654
bool Scale_arc_length
Boolean to control whether arc-length should be scaled.
Definition: problem.h:731
double doubly_adaptive_unsteady_newton_solve_helper(const double &dt, const double &epsilon, const unsigned &max_adapt, const unsigned &suppress_resolve_after_spatial_adapt, const bool &first, const bool &shift=true)
Private helper function that actually performs the unsteady "doubly" adaptive Newton solve....
Definition: problem.cc:11379
bool Empty_actions_after_read_unstructured_meshes_has_been_called
Boolean to indicate that empty actions_after_read_unstructured_meshes() function has been called.
Definition: problem.h:222
virtual void get_residuals(DoubleVector &residuals)
Return the fully-assembled residuals Vector for the problem: Virtual so it can be overloaded in for m...
Definition: problem.cc:3714
Vector< GeneralisedElement * > Base_mesh_element_pt
Vector to store the correspondence between a root element and its element number within the global me...
Definition: problem.h:548
void get_fd_jacobian(DoubleVector &residuals, DenseMatrix< double > &jacobian)
Return the fully-assembled Jacobian and residuals, generated by finite differences.
Definition: problem.cc:7823
double doubly_adaptive_unsteady_newton_solve(const double &dt, const double &epsilon, const unsigned &max_adapt, const bool &first, const bool &shift=true)
Unsteady "doubly" adaptive Newton solve: Does temporal adaptation first, i.e. we try to do a timestep...
Definition: problem.h:2471
virtual void sparse_assemble_row_or_column_compressed_with_two_vectors(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:5698
virtual Problem * make_copy()
Make and return a pointer to the copy of the problem. A virtual function that must be filled in by th...
Definition: problem.cc:12011
void disable_problem_distributed()
Disable problem distributed.
Definition: problem.h:977
double Maximum_dt
Maximum desired dt.
Definition: problem.h:709
unsigned long set_timestepper_for_all_data(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data=false)
Set all problem data to have the same timestepper (timestepper_pt) Return the new number of dofs in t...
Definition: problem.cc:11572
void adapt_based_on_error_estimates(Vector< Vector< double > > &elemental_error)
Adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on the error estimat...
Definition: problem.h:2910
double & maximum_dt()
Access function to max timestep in adaptive timestepping.
Definition: problem.h:1588
void activate_pitchfork_tracking(double *const &parameter_pt, const DoubleVector &symmetry_vector, const bool &block_solve=true)
Turn on pitchfork tracking using the augmented system specified in the PitchForkHandler class....
Definition: problem.cc:10188
virtual void dump(std::ofstream &dump_file) const
Dump refinement pattern of all refineable meshes and all generic Problem data to file for restart.
Definition: problem.cc:12029
bool Use_finite_differences_for_continuation_derivatives
Boolean to specify which scheme to use to calculate the continuation derivatievs.
Definition: problem.h:803
void p_refine_uniformly_and_prune(DocInfo &doc_info)
p-refine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process.
Definition: problem.h:2720
bool Arc_length_step_taken
Boolean to indicate whether an arc-length step has been taken.
Definition: problem.h:799
void set_pinned_values_to_zero()
Set all pinned values to zero. Used to set boundary conditions to be homogeneous in the copy of the p...
Definition: problem.cc:4304
bool Default_set_initial_condition_called
Has default set_initial_condition function been called? Default: false.
Definition: problem.h:211
void get_bifurcation_eigenfunction(Vector< DoubleVector > &eigenfunction)
Return the eigenfunction calculated as part of a bifurcation tracking process. If we are not tracking...
Definition: problem.cc:10091
double Relaxation_factor
Relaxation fator for Newton method (only a fractional Newton correction is applied if this is less th...
Definition: problem.h:592
void add_time_stepper_pt(TimeStepper *const &time_stepper_pt)
Add a timestepper to the problem. The function will automatically create or resize the Time object so...
Definition: problem.cc:1545
void refine_selected_elements(const Vector< unsigned > &elements_to_be_refined)
Refine (one and only!) mesh by splitting the elements identified by their numbers relative to the pro...
Definition: problem.cc:14898
LinearAlgebraDistribution *const & dof_distribution_pt() const
Return the pointer to the dof distribution (read-only)
Definition: problem.h:1668
void prune_halo_elements_and_nodes(DocInfo &doc_info, const bool &report_stats)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement,...
Definition: problem.cc:1053
virtual void get_inverse_mass_matrix_times_residuals(DoubleVector &Mres)
Return the residual vector multiplied by the inverse mass matrix Virtual so that it can be overloaded...
Definition: problem.cc:3579
void check_halo_schemes()
Check the halo/haloed node/element schemes.
Definition: problem.h:2116
double Parameter_current
Storage for the present value of the global parameter.
Definition: problem.h:755
void enable_jacobian_reuse()
Enable recycling of Jacobian in Newton iteration (if the linear solver allows it)....
Definition: problem.h:1975
double *& dof_pt(const unsigned &i)
Pointer to i-th dof in the problem.
Definition: problem.h:1825
void assign_eigenvector_to_dofs(DoubleVector &eigenvector)
Assign the eigenvector passed to the function to the dofs in the problem so that it can be output by ...
Definition: problem.cc:8716
Distributed_problem_matrix_distribution
enum for distribution of distributed jacobians. 1 - Automatic - the Problem distribution is employed,...
Definition: problem.h:829
@ Uniform_matrix_distribution
Definition: problem.h:831
@ Default_matrix_distribution
Definition: problem.h:829
@ Problem_matrix_distribution
Definition: problem.h:830
void send_refinement_info_helper(Vector< unsigned > &old_domain_for_base_element, Vector< unsigned > &new_domain_for_base_element, const unsigned &max_refinement_level_overall, std::map< unsigned, Vector< unsigned > > &refinement_info_for_root_local, Vector< Vector< Vector< unsigned > > > &refinement_info_for_root_elements)
Send refinement information between processors.
Definition: problem.cc:18328
unsigned setup_element_count_per_dof()
Function that populates the Element_counter_per_dof vector with the number of elements that contribut...
Definition: problem.cc:230
Data *& global_data_pt(const unsigned &i)
Return a pointer to the the i-th global data object.
Definition: problem.h:1647
void disable_globally_convergent_newton_method()
disable globally convergent Newton method
Definition: problem.h:2009
void problem_is_nonlinear(const bool &prob_lin)
Access function to Problem_is_nonlinear.
Definition: problem.h:1600
void disable_doc_imbalance_in_parallel_assembly()
Disable doc of load imbalance in parallel assembly of distributed problem.
Definition: problem.h:856
OomphCommunicator * Communicator_pt
The communicator for this problem.
Definition: problem.h:1242
double Newton_solver_tolerance
The Tolerance below which the Newton Method is deemed to have converged.
Definition: problem.h:596
void set_consistent_pinned_values_for_continuation()
Private helper function that is used to set the appropriate pinned values for continuation.
Definition: problem.cc:10466
void activate_bifurcation_tracking(double *const &parameter_pt, const DoubleVector &eigenvector, const bool &block_solve=true)
Activate generic bifurcation tracking for a single (real) eigenvalue where the initial guess for the ...
Definition: problem.cc:10129
LinearSolver *& mass_matrix_solver_for_explicit_timestepper_pt()
Return a pointer to the linear solver object used for explicit time stepping.
Definition: problem.h:1479
void prune_halo_elements_and_nodes(const bool &report_stats=false)
(Irreversibly) prune halo(ed) elements and nodes, usually after another round of refinement,...
Definition: problem.h:2167
bool Discontinuous_element_formulation
Is the problem a discontinuous one, i.e. can the elemental contributions be treated independently....
Definition: problem.h:700
bool Doc_imbalance_in_parallel_assembly
Boolean to switch on assessment of load imbalance in parallel assembly of distributed problem.
Definition: problem.h:510
long synchronise_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Classify any non-classified nodes into halo/haloed and synchronise equation numbers....
Definition: problem.cc:16705
void p_refine_uniformly_aux(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info, const bool &prune)
Helper function to do compund p-refinement of (all) p-refineable (sub)mesh(es) uniformly as many time...
Definition: problem.cc:15592
void calculate_continuation_derivatives_helper(const DoubleVector &z)
A function that performs the guts of the continuation derivative calculation in arc length continuati...
Definition: problem.cc:9912
Vector< Problem * > Copy_of_problem_pt
Vector of pointers to copies of the problem used in adaptive bifurcation tracking problems (ALH: TEMP...
Definition: problem.h:234
Vector< unsigned > distribute(const Vector< unsigned > &element_partition, DocInfo &doc_info, const bool &report_stats=false)
Distribute the problem and doc, using the specified partition; returns a vector which details the par...
Definition: problem.cc:404
bool problem_has_been_distributed()
Access to Problem_has_been_distributed flag.
Definition: problem.h:2178
Mesh * Mesh_pt
The mesh pointer.
Definition: problem.h:167
double Parameter_derivative
Storage for the derivative of the global parameter wrt arc-length.
Definition: problem.h:752
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1524
void add_eigenvector_to_dofs(const double &epsilon, const DoubleVector &eigenvector)
Add the eigenvector passed to the function scaled by the constat epsilon to the dofs in the problem s...
Definition: problem.cc:8747
void flush_global_data()
Flush the Problem's global data – resizes container to zero. Data objects are not deleted!
Definition: problem.h:1662
Vector< Data * > Global_data_pt
Vector of global data: "Nobody" (i.e. none of the elements etc.) is "in charge" of this Data so it wo...
Definition: problem.h:425
void recompute_load_balanced_assembly()
Helper function to re-assign the first and last elements to be assembled by each processor during par...
Definition: problem.cc:1689
Vector< double * > Dof_pt
Vector of pointers to dofs.
Definition: problem.h:554
void p_adapt()
p-adapt problem: Perform mesh adaptation for (all) refineable (sub)mesh(es), based on their own error...
Definition: problem.h:2881
double DTSF_max_increase
Maximum possible increase of dt between time-steps in adaptive schemes.
Definition: problem.h:713
bool Must_recompute_load_balance_for_assembly
Boolean indicating that the division of elements over processors during the assembly process must be ...
Definition: problem.h:529
virtual void shift_time_values()
Shift all values along to prepare for next timestep.
Definition: problem.cc:11634
void p_refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
p-refine p-refineable sub-meshes, each as many times as specified in the vector and rebuild problem; ...
Definition: problem.h:2694
double * dof_pt(const unsigned &i) const
Pointer to i-th dof in the problem (const version)
Definition: problem.h:1831
void set_default_first_and_last_element_for_assembly()
Set default first and last elements for parallel assembly of non-distributed problem.
Definition: problem.cc:1613
bool Keep_temporal_error_below_tolerance
Boolean to decide if a timestep is to be rejected if the error estimate post-solve (computed by globa...
Definition: problem.h:2209
bool Shut_up_in_newton_solve
Boolean to indicate if all output is suppressed in Problem::newton_solve(). Defaults to false.
Definition: problem.h:2191
void set_dofs(const DoubleVector &dofs)
Set the values of the dofs.
Definition: problem.cc:3411
void setup_base_mesh_info_after_pruning()
Helper function to re-setup the Base_mesh enumeration (used during load balancing) after pruning.
Definition: problem.cc:20102
Vector< Mesh * > Sub_mesh_pt
Vector of pointers to submeshes.
Definition: problem.h:170
ExplicitTimeStepper *& explicit_time_stepper_pt()
Return a pointer to the explicit timestepper.
Definition: problem.h:1555
Vector< double > Max_res
Maximum residuals at start and after each newton iteration.
Definition: problem.h:606
EigenSolver * Eigen_solver_pt
Pointer to the eigen solver for the problem.
Definition: problem.h:182
bool Bisect_to_find_bifurcation
Boolean to control wheter bisection is used to located bifurcation.
Definition: problem.h:793
void explicit_timestep(const double &dt, const bool &shift_values=true)
Take an explicit timestep of size dt and optionally shift any stored values of the time history.
Definition: problem.cc:10918
void delete_all_external_storage()
Wrapper function to delete external storage for each submesh of the problem.
Definition: problem.cc:16355
double DTSF_min_decrease
Minimum allowed decrease of dt between time-steps in adaptive schemes. Lower scaling values will reje...
Definition: problem.h:718
virtual void symmetrise_eigenfunction_for_adaptive_pitchfork_tracking()
Virtual function that is used to symmetrise the problem so that the current solution exactly satisfie...
Definition: problem.cc:10059
void load_balance(DocInfo &doc_info, const bool &report_stats)
Balance the load of a (possibly non-uniformly refined) problem that has already been distributed,...
Definition: problem.h:1422
double Minimum_dt_but_still_proceed
If Minimum_dt_but_still_proceed positive, then dt will not be reduced below this value during adaptiv...
Definition: problem.h:725
double adaptive_unsteady_newton_solve(const double &dt_desired, const double &epsilon)
Attempt to advance timestep by dt_desired. If the solution fails the timestep will be halved until co...
Definition: problem.cc:11056
unsigned Desired_newton_iterations_ds
The desired number of Newton Steps to reach convergence at each step along the arc.
Definition: problem.h:784
void refine_uniformly_and_prune(DocInfo &doc_info)
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process.
Definition: problem.h:2627
void rebuild_global_mesh()
If one of the submeshes has changed (e.g. by mesh adaptation) we need to update the global mesh....
Definition: problem.cc:1533
void disable_store_local_dof_pt_in_elements()
Insist that local dof pointers are NOT set up in each element when equation numbering takes place (th...
Definition: problem.h:1703
bool Calculate_hessian_products_analytic
Map used to determine whether the hessian products should be computed using finite differences....
Definition: problem.h:244
void set_first_and_last_element_for_assembly(Vector< unsigned > &first_el_for_assembly, Vector< unsigned > &last_el_for_assembly)
Manually set first and last elements for parallel assembly of non-distributed problem.
Definition: problem.h:988
void activate_hopf_tracking(double *const &parameter_pt, const bool &block_solve=true)
Turn on Hopf bifurcation tracking using the augmented system specified in the HopfHandler class....
Definition: problem.cc:10218
Vector< double * > Halo_dof_pt
Storage for the halo degrees of freedom (only required) when accessing via the global equation number...
Definition: problem.h:581
void deactivate_bifurcation_tracking()
Deactivate all bifuraction tracking, by reseting the assembly handler to the default.
Definition: problem.h:2353
Vector< double > elemental_assembly_time()
Return vector of most-recent elemental assembly times (used for load balancing). Zero sized if no Jac...
Definition: problem.h:864
void globally_convergent_line_search(const Vector< double > &x_old, const double &half_residual_squared_old, DoubleVector &gradient, DoubleVector &newton_dir, double &half_residual_squared, const double &stpmax)
Line search helper for globally convergent Newton method.
Definition: problem.cc:9146
void get_hessian_vector_products(DoubleVectorWithHaloEntries const &Y, Vector< DoubleVectorWithHaloEntries > const &C, Vector< DoubleVectorWithHaloEntries > &product)
Return the product of the global hessian (derivative of Jacobian matrix with respect to all variables...
Definition: problem.cc:7985
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1466
virtual double global_temporal_error_norm()
Function to calculate a global error norm, used in adaptive timestepping to control the change in tim...
Definition: problem.h:1230
Assembly_method
Enumerated flags to determine which sparse assembly method is used.
Definition: problem.h:645
@ Perform_assembly_using_two_arrays
Definition: problem.h:650
@ Perform_assembly_using_maps
Definition: problem.h:648
@ Perform_assembly_using_two_vectors
Definition: problem.h:647
@ Perform_assembly_using_vectors_of_pairs
Definition: problem.h:646
@ Perform_assembly_using_lists
Definition: problem.h:649
void enable_globally_convergent_newton_method()
enable globally convergent Newton method
Definition: problem.h:2003
void refine_distributed_base_mesh(Vector< Vector< Vector< unsigned > > > &to_be_refined_on_each_root, const unsigned &max_level_overall)
Load balance helper routine: refine each new base (sub)mesh based upon the elements to be refined wit...
Definition: problem.cc:19990
int Sign_of_jacobian
Storage for the sign of the global Jacobian.
Definition: problem.h:745
std::map< GeneralisedElement *, unsigned > Base_mesh_element_number_plus_one
Map which stores the correspondence between a root element and its element number (plus one) within t...
Definition: problem.h:538
double Max_residuals
Maximum desired residual: if the maximum residual exceeds this value, the program will exit.
Definition: problem.h:610
unsigned nsub_mesh() const
Return number of submeshes.
Definition: problem.h:1323
void unset_analytic_hessian_products()
Function to turn off analytic calculation of the parameter derivatives in continuation and bifurcatio...
Definition: problem.h:296
double & time()
Return the current value of continuous time.
Definition: problem.cc:11531
unsigned self_test()
Self-test: Check meshes and global data. Return 0 for OK.
Definition: problem.cc:13276
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1680
void activate_fold_tracking(double *const &parameter_pt, const bool &block_solve=true)
Turn on fold tracking using the augmented system specified in the FoldHandler class....
Definition: problem.cc:10103
bool Use_globally_convergent_newton_method
Use the globally convergent newton method.
Definition: problem.h:214
LinearSolver * Default_linear_solver_pt
Pointer to the default linear solver.
Definition: problem.h:188
double Minimum_ds
Minimum desired value of arc-length.
Definition: problem.h:787
void get_data_to_be_sent_during_load_balancing(const Vector< unsigned > &element_domain_on_this_proc, Vector< int > &send_n, Vector< double > &send_data, Vector< int > &send_displacement, Vector< unsigned > &old_domain_for_base_element, Vector< unsigned > &new_domain_for_base_element, unsigned &max_refinement_level_overall)
Load balance helper routine: Get data to be sent to other processors during load balancing and other ...
Definition: problem.cc:19279
void restore_dof_values()
Restore the stored values of the degrees of freedom.
Definition: problem.cc:8653
double Numerical_zero_for_sparse_assembly
A tolerance used to determine whether the entry in a sparse matrix is zero. If it is then storage nee...
Definition: problem.h:671
double FD_step_used_in_get_hessian_vector_products
Definition: problem.h:685
virtual void partition_global_mesh(Mesh *&global_mesh_pt, DocInfo &doc_info, Vector< unsigned > &element_domain, const bool &report_stats=false)
Partition the global mesh, return vector specifying the processor number for each element....
Definition: problem.cc:878
unsigned Sparse_assemble_with_arrays_initial_allocation
the number of elements to initially allocate for a matrix row within the sparse_assembly_with_two_arr...
Definition: problem.h:658
void bifurcation_adapt_doc_errors(const unsigned &bifurcation_type)
A function that is used to document the errors used in the adaptive solution of bifurcation problems.
Definition: problem.cc:13648
double Ds_current
Storage for the current step value.
Definition: problem.h:780
virtual void set_initial_condition()
Set initial condition (incl previous timesteps). We need to establish this interface because I....
Definition: problem.h:1198
LinearSolver * Mass_matrix_solver_for_explicit_timestepper_pt
Pointer to the linear solver used for explicit time steps (this is likely to be different to the line...
Definition: problem.h:179
void get_all_halo_data(std::map< unsigned, double * > &map_of_halo_data)
Get pointers to all possible halo data indexed by global equation number in a map.
Definition: problem.cc:16377
void solve_eigenproblem(const unsigned &n_eval, Vector< std::complex< double > > &eigenvalue, const bool &steady=true)
Solve an eigenproblem as assembled by EigenElements, but only return the eigenvalues,...
Definition: problem.h:1937
void load_balance()
Balance the load of a (possibly non-uniformly refined) problem that has already been distributed,...
Definition: problem.h:1382
double Minimum_dt
Minimum desired dt: if dt falls below this value, exit.
Definition: problem.h:706
double arc_length_step_solve(double *const &parameter_pt, const double &ds, const unsigned &max_adapt=0)
Solve a steady problem using arc-length continuation, when the parameter that becomes a variable corr...
Definition: problem.cc:10294
virtual void actions_after_adapt()
Actions that are to be performed after a mesh adaptation.
Definition: problem.h:1025
void p_refine_uniformly()
p-refine (all) p-refineable (sub)mesh(es) uniformly and rebuild problem
Definition: problem.h:2737
bool Use_continuation_timestepper
Boolean to control original or new storage of dof stuff.
Definition: problem.h:758
void calculate_continuation_derivatives_fd(double *const &parameter_pt)
A function to calculate the derivatives with respect to the arc-length required for continuation by f...
Definition: problem.cc:9818
LinearAlgebraDistribution * Dof_distribution_pt
The distribution of the DOFs in this problem. This object is created in the Problem constructor and s...
Definition: problem.h:460
void refine_uniformly()
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem.
Definition: problem.h:2640
static ContinuationStorageScheme Continuation_time_stepper
Storage for the single static continuation timestorage object.
Definition: problem.h:761
bool Problem_is_nonlinear
Boolean flag indicating if we're dealing with a linear or nonlinear Problem – if set to false the New...
Definition: problem.h:628
Distributed_problem_matrix_distribution & distributed_problem_matrix_distribution()
accesss function to the distributed matrix distribution method 1 - Automatic - the Problem distributi...
Definition: problem.h:842
virtual void sparse_assemble_row_or_column_compressed_with_maps(Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
Private helper function that is used to assemble the Jacobian matrix in the case when the storage is ...
Definition: problem.cc:4584
void refine_uniformly_and_prune(const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
Refine refineable sub-meshes, each as many times as specified in the vector and rebuild problem; doc ...
Definition: problem.h:2605
Time * time_pt() const
Return a pointer to the global time object (const version).
Definition: problem.h:1510
void operator=(const Problem &)=delete
Broken assignment operator.
AssemblyHandler * Default_assembly_handler_pt
Pointer to the default assembly handler.
Definition: problem.h:194
void get_my_eqns(AssemblyHandler *const &assembly_handler_pt, const unsigned &el_lo, const unsigned &el_hi, Vector< unsigned > &my_eqns)
Helper method that returns the (unique) global equations to which the elements in the range el_lo to ...
Definition: problem.cc:6509
std::map< double *, bool > Calculate_dparameter_analytic
Map used to determine whether the derivatives with respect to a parameter should be finite difference...
Definition: problem.h:239
virtual void read(std::ifstream &restart_file, bool &unsteady_restart)
Read refinement pattern of all refineable meshes and refine them accordingly, then read all Data and ...
Definition: problem.cc:12251
virtual ~Problem()
Virtual destructor to clean up memory.
Definition: problem.cc:181
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
void get_dofs(DoubleVector &dofs) const
Return the vector of dofs, i.e. a vector containing the current values of all unknowns.
Definition: problem.cc:2479
Time * Time_pt
Pointer to global time for the problem.
Definition: problem.h:197
DoubleVectorWithHaloEntries Element_count_per_dof
Counter that records how many elements contribute to each dof. Used to determine the (discrete) arc-l...
Definition: problem.h:560
void refine_uniformly(DocInfo &doc_info)
Refine (all) refineable (sub)mesh(es) uniformly and rebuild problem; doc refinement process.
Definition: problem.h:2614
double & newton_solver_tolerance()
Access function to tolererance of the Newton solver, i.e. the maximum value of the residuals that wil...
Definition: problem.h:1621
virtual void actions_before_newton_convergence_check()
Any actions that are to be performed before the residual is checked in the Newton method,...
Definition: problem.h:1048
double doubly_adaptive_unsteady_newton_solve(const double &dt, const double &epsilon, const unsigned &max_adapt, const unsigned &suppress_resolve_after_spatial_adapt_flag, const bool &first, const bool &shift=true)
Unsteady "doubly" adaptive Newton solve: Does temporal adaptation first, i.e. we try to do a timestep...
Definition: problem.h:2504
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
TimeStepper *& time_stepper_pt(const unsigned &i)
Return a pointer to the i-th timestepper.
Definition: problem.h:1549
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
Definition: problem.h:1570
Mesh *const & mesh_pt() const
Return a pointer to the global mesh (const version)
Definition: problem.h:1286
unsigned newton_solve_continuation(double *const &parameter_pt)
Perform a basic arc-length continuation step using Newton's method. Returns number of Newton steps ta...
Definition: problem.cc:9376
double Max_permitted_error_for_halo_check
Threshold for error throwing in Problem::check_halo_schemes()
Definition: problem.h:2175
unsigned Sparse_assemble_with_arrays_allocation_increment
the number of elements to add to a matrix row when the initial allocation is exceeded within the spar...
Definition: problem.h:663
void reset_assembly_handler_to_default()
Reset the system to the standard non-augemented state.
Definition: problem.cc:10275
void set_analytic_dparameter(double *const &parameter_pt)
Function to turn on analytic calculation of the parameter derivatives in continuation and bifurcation...
Definition: problem.h:256
virtual void actions_after_newton_step()
Any actions that are to be performed after each individual Newton step. Most likely to be used for di...
Definition: problem.h:1058
bool Pause_at_end_of_sparse_assembly
Protected boolean flag to halt program execution during sparse assemble process to assess peak memory...
Definition: problem.h:633
void remove_null_pointers_from_external_halo_node_storage()
Consolidate external halo node storage by removing nulled out pointers in external halo and haloed sc...
Definition: problem.cc:3178
Vector< double > * Saved_dof_pt
Pointer to vector for backup of dofs.
Definition: problem.h:207
void unsteady_newton_solve(const double &dt)
Advance time by dt and solve by Newton's method. This version always shifts time values.
Definition: problem.cc:10953
AssemblyHandler * Assembly_handler_pt
Definition: problem.h:185
void disable_info_in_newton_solve()
Disable the output of information when in the newton solver.
Definition: problem.h:2941
unsigned Dof_derivative_offset
Storage for the offset for the continuation derivatives from the original dofs (default is 1,...
Definition: problem.h:766
virtual void actions_before_read_unstructured_meshes()
Actions that are to be performed before reading in restart data for problems involving unstructured b...
Definition: problem.h:1095
Mesh *const & mesh_pt(const unsigned &imesh) const
Return a pointer to the i-th submesh (const version)
Definition: problem.h:1308
Vector< TimeStepper * > Time_stepper_pt
The Vector of time steppers (there could be many different ones in multiphysics problems)
Definition: problem.h:201
void doc_errors()
Get max and min error for all elements in submeshes.
Definition: problem.h:2926
double & max_residuals()
Access function to max residuals in Newton iterations before giving up.
Definition: problem.h:1608
int & sign_of_jacobian()
Access function for the sign of the global jacobian matrix. This will be set by the linear solver,...
Definition: problem.h:2424
bool Mass_matrix_reuse_is_enabled
Is re-use of the mass matrix in explicit timestepping enabled Default:false.
Definition: problem.h:691
void get_derivative_wrt_global_parameter(double *const &parameter_pt, DoubleVector &result)
Get the derivative of the entire residuals vector wrt a global parameter, used in continuation proble...
Definition: problem.cc:7887
void enable_discontinuous_formulation()
Indicate that the problem involves discontinuous elements This allows for a more efficiently assembly...
Definition: problem.h:1727
bool distributed() const
If we have MPI return the "problem has been distributed" flag, otherwise it can't be distributed so r...
Definition: problem.h:808
bool are_hessian_products_calculated_analytically()
Function to determine whether the hessian products are calculated analytically.
Definition: problem.h:303
bool does_pointer_correspond_to_problem_data(double *const &parameter_pt)
Return a boolean flag to indicate whether the pointer parameter_pt refers to values stored in a Data ...
Definition: problem.cc:9846
ExplicitTimeStepper * Explicit_time_stepper_pt
Pointer to a single explicit timestepper.
Definition: problem.h:204
void dump(const std::string &dump_file_name) const
Dump refinement pattern of all refineable meshes and all generic Problem data to file for restart.
Definition: problem.h:2078
bool mass_matrix_reuse_is_enabled()
Return whether the mass matrix is being reused.
Definition: problem.h:2568
void disable_discontinuous_formulation()
Disable the use of a discontinuous-element formulation. Note that the methods will all still work eve...
Definition: problem.h:1736
double arc_length_step_solve_helper(double *const &parameter_pt, const double &ds, const unsigned &max_adapt)
Private helper function that actually contains the guts of the arc-length stepping,...
Definition: problem.cc:10515
void refine_uniformly(const unsigned &i_mesh)
Do uniform refinement for submesh i_mesh without documentation.
Definition: problem.h:2651
void solve_eigenproblem(const unsigned &n_eval, Vector< std::complex< double > > &eigenvalue, Vector< DoubleVector > &eigenvector, const bool &steady=true)
Get derivative of an element in the problem wrt a global parameter, used in continuation problems.
Definition: problem.cc:8313
bool Use_predictor_values_as_initial_guess
Use values from the time stepper predictor as an initial guess.
Definition: problem.h:229
void unset_analytic_dparameter(double *const &parameter_pt)
Function to turn off analytic calculation of the parameter derivatives in continuation and bifurcatio...
Definition: problem.h:263
LinearSolver * mass_matrix_solver_for_explicit_timestepper_pt() const
Return a pointer to the linear solver object used for explicit time stepping (const version).
Definition: problem.h:1486
Class for a matrix of the form M = S + G + H + ... where S is the main matrix and G,...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
Class to keep track of discrete/continous time. It is essential to have a single Time object when usi...
Definition: timesteppers.h:63
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...