CppNoddy  0.90
Public Member Functions | List of all members
CppNoddy::ODE_BVP< _Type, _Xtype > Class Template Reference

A templated object for real/complex vector system of first-order ordinary differential equations. More...

#include <ODE_BVP.h>

Inheritance diagram for CppNoddy::ODE_BVP< _Type, _Xtype >:
CppNoddy::ArcLength_base< _Type > CppNoddy::Uncopyable

Public Member Functions

 ODE_BVP (Equation_1matrix< _Type, _Xtype > *ptr_to_equation, const DenseVector< _Xtype > &nodes, Residual< _Type > *ptr_to_left_residual, Residual< _Type > *ptr_to_right_residual)
 The class is defined by a vector function for the system. More...
 
virtual ~ODE_BVP ()
 Destructor. More...
 
virtual void actions_before_linear_solve (BandedMatrix< _Type > &a, DenseVector< _Type > &b)
 A virtual method that is called prior to the linear solve stage of the solve2() method. More...
 
void solve2 ()
 Formulate and solve the ODE using Newton iteration and a second-order finite difference scheme. More...
 
std::pair< unsigned, unsigned > adapt (const double &adapt_tol)
 Adapt the computational mesh ONCE. More...
 
void adapt_until (const double &adapt_tol)
 Adaptively solve the system until no refinements or unrefinements are applied. More...
 
void set_monitor_det (bool flag)
 Set the flag that determines if the determinant will be monitored The default is to monitor. More...
 
void init_arc (_Type *p, const double &length, const double &max_length)
 Initialise the class ready for arc length continuation. More...
 
double arclength_solve (const double &step)
 Arc-length solve the system. More...
 
OneD_Node_Mesh< _Type, _Xtype > & solution ()
 
double & tolerance ()
 Access method to the tolerance. More...
 
int & max_itns ()
 Access method to the maximum number of iterations. More...
 
template<>
std::pair< unsigned, unsigned > adapt (const double &adapt_tol)
 
template<>
std::pair< unsigned, unsigned > adapt (const double &adapt_tol)
 
template<>
void adapt_until (const double &adapt_tol)
 
template<>
void adapt_until (const double &adapt_tol)
 
- Public Member Functions inherited from CppNoddy::ArcLength_base< _Type >
 ArcLength_base ()
 
virtual ~ArcLength_base ()
 
void init_arc (DenseVector< _Type > x, _Type *p, const double &length, const double &max_length)
 Initialise the class ready for arc-length continuation. More...
 
double & ds ()
 Return a handle to the arclength step. More...
 
double & arcstep_multiplier ()
 Used to set the multiplication constant used when increasing or decreasing the arclength step. More...
 
bool & rescale_theta ()
 Handle to the RESCALE_THETA flag. More...
 
double & theta ()
 Set the arclength theta parameter. More...
 
double & desired_arc_proportion ()
 Handle to the desired proportion of the parameter to be used in the arc length solver. More...
 

Additional Inherited Members

- Protected Member Functions inherited from CppNoddy::ArcLength_base< _Type >
void update (const DenseVector< _Type > &x)
 A method called by arclength_solve and init_arc which stores the current converged state and parameter and hence computes the derivatives w.r.t the arc-length. More...
 
double arclength_residual (const DenseVector< _Type > &x) const
 The extra constraint that is to be used to replace the unknown arc length. More...
 
DenseVector< _Type > Jac_arclength_residual (DenseVector< _Type > &x) const
 The derivative of the arclength_residual function with respect to each of the state variables. More...
 
void update_theta (const DenseVector< _Type > &x)
 Automatically update the Keller THETA such that the proportion of the arclength obtained from the parameter is the desired value. More...
 
- Protected Attributes inherited from CppNoddy::ArcLength_base< _Type >
_Type * p_PARAM
 pointer to the parameter in arclength solves More...
 
DenseVector< _Type > LAST_X
 state variable at the last computed solution More...
 
DenseVector< _Type > X_DERIV_S
 derivative of the state variable w.r.t. arc length More...
 
_Type LAST_PARAM
 parameter value at the last computed solution More...
 
_Type PARAM_DERIV_S
 derivative of the parameter w.r.t arc length More...
 
double DS
 size of the arc length step More...
 
double MAX_DS
 maximum arc length step to be taken More...
 
double ARCSTEP_MULTIPLIER
 step change multiplier More...
 
bool INITIALISED
 for the arc-length solver - to show it has been initialised More...
 
double THETA
 the arclength theta More...
 

Detailed Description

template<typename _Type, typename _Xtype = double>
class CppNoddy::ODE_BVP< _Type, _Xtype >

A templated object for real/complex vector system of first-order ordinary differential equations.

The class is double templated, with the first type associated with real/complex data, and second (real/complex) type associated with a problem on the real line or line in the complex plane.

Definition at line 38 of file ODE_BVP.h.

Constructor & Destructor Documentation

template<typename _Type, typename _Xtype>
CppNoddy::ODE_BVP< _Type, _Xtype >::ODE_BVP ( Equation_1matrix< _Type, _Xtype > *  ptr_to_equation,
const DenseVector< _Xtype > &  nodes,
Residual< _Type > *  ptr_to_left_residual,
Residual< _Type > *  ptr_to_right_residual 
)

The class is defined by a vector function for the system.

Parameters
ptr_to_equationA pointer to an Equation_1matrix object.
nodesA vector that defines the nodal positions.
ptr_to_left_residualA pointer to a residual object that defines the LHS boundary conditions.
ptr_to_right_residualA pointer to a residual object that defines the RHS boundary conditions.

Definition at line 32 of file ODE_BVP.cpp.

35  :
36  ArcLength_base<_Type> (),
37  MAX_ITERATIONS( 12 ),
38  TOL( 1.e-8 ),
39  LAST_DET_SIGN( 0 ),
40  p_EQUATION( ptr_to_equation ),
41  p_LEFT_RESIDUAL( ptr_to_left_residual ),
42  p_RIGHT_RESIDUAL( ptr_to_right_residual ),
43  MONITOR_DET( true )
44  {
45  // set up the solution mesh using these nodes
46  SOLUTION = OneD_Node_Mesh<_Type, _Xtype>( nodes, p_EQUATION -> get_order() );
47  if ( ( p_LEFT_RESIDUAL -> get_number_of_vars() != p_EQUATION -> get_order() ) ||
48  ( p_RIGHT_RESIDUAL -> get_number_of_vars() != p_EQUATION -> get_order() ) ||
49  ( p_LEFT_RESIDUAL -> get_order() + p_RIGHT_RESIDUAL -> get_order() != p_EQUATION -> get_order() ) )
50  {
51  std::cout << "order " << p_EQUATION -> get_order() << "\n";
52  std::cout << "left nvars " << p_LEFT_RESIDUAL -> get_number_of_vars() << "\n";
53  std::cout << "right nvars " << p_RIGHT_RESIDUAL -> get_number_of_vars() << "\n";
54  std::cout << "left order " << p_LEFT_RESIDUAL -> get_order() << "\n";
55  std::cout << "right order " << p_RIGHT_RESIDUAL -> get_order() << "\n";
56  std::string problem;
57  problem = "It looks like the ODE_BVP equation and boundary conditions are\n";
58  problem += "not well posed. The number of variables for each boundary condition\n";
59  problem += "has to be the same as the order of the equation. Also the order of\n";
60  problem += "both boundary conditions has to sum to the order of the equation.\n";
61  throw ExceptionRuntime( problem );
62  }
63 #ifdef TIME
64  // timers
65  T_ASSEMBLE = Timer( "ODE_BVP: Assembling of the matrix (incl. equation updates):" );
66  T_SOLVE = Timer( "ODE_BVP: Solving of the matrix:" );
67  T_REFINE = Timer( "ODE_BVP: Refining the mesh:" );
68 #endif
69  }
template<typename _Type , typename _Xtype >
CppNoddy::ODE_BVP< _Type, _Xtype >::~ODE_BVP ( )
virtual

Destructor.

Definition at line 72 of file ODE_BVP.cpp.

73  {
74 #ifdef TIME
75  std::cout << "\n";
76  //T_ASSEMBLE.stop();
77  T_ASSEMBLE.print();
78  //T_SOLVE.stop();
79  T_SOLVE.print();
80  //T_REFINE.stop();
81  T_REFINE.print();
82 #endif
83  }

Member Function Documentation

template<typename _Type, typename _Xtype = double>
virtual void CppNoddy::ODE_BVP< _Type, _Xtype >::actions_before_linear_solve ( BandedMatrix< _Type > &  a,
DenseVector< _Type > &  b 
)
inlinevirtual

A virtual method that is called prior to the linear solve stage of the solve2() method.

This is called prior to any linear solve to allow the user to manually tune the matrix problem directly.

Parameters
aThe Jacbian matrix to be passed to the linear solver
bThe residual vector to be passed to the linear solver

Definition at line 65 of file ODE_BVP.h.

Referenced by CppNoddy::ODE_BVP< _Type, _Xtype >::arclength_solve(), and CppNoddy::ODE_BVP< _Type, _Xtype >::solve2().

66  {}
template<typename _Type , typename _Xtype >
std::pair< unsigned, unsigned > CppNoddy::ODE_BVP< _Type, _Xtype >::adapt ( const double &  adapt_tol)

Adapt the computational mesh ONCE.

We step through each interior node and evaluate the residual at that node. If the residual is less than the convergence tolerance, the node is removed. If the residual is greater than the adapt_tol parameter, then the mesh is adapted with an addition node place either side of the evaluation node.

Parameters
adapt_tolThe residual tolerance at a nodal point that will lead to the mesh adaptation.
Returns
A pair of values indicating the number of refinements and unrefinements.

Definition at line 533 of file ODE_BVP.cpp.

References CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::DenseVector< _Type >::push_back(), and CppNoddy::ArcLength_base< _Type >::update().

Referenced by CppNoddy::ODE_BVP< D_complex >::actions_before_linear_solve(), CppNoddy::ODE_BVP< _Type, _Xtype >::adapt_until(), and main().

534  {
535 #ifdef TIME
536  T_REFINE.start();
537 #endif
538  // the order of the problem
539  unsigned order( p_EQUATION -> get_order() );
540  // get the number of nodes in the mesh
541  // -- this may have been refined by the user since the
542  // last call.
543  unsigned N( SOLUTION.get_nnodes() );
544  // row counter
545  //std::size_t row( 0 );
546  // local state variable and functions
547  DenseVector<_Type> F_node( order, 0.0 );
548  DenseVector<_Type> R_node( order, 0.0 );
549  // make sure (un)refine flags vector is sized
550  std::vector< bool > refine( N, false );
551  std::vector< bool > unrefine( N, false );
552  // Residual vector at interior nodes
553  DenseVector<double> Res2( N, 0.0 );
554  // reset row counter
555  //row = 0;
556  // inner nodes of the mesh, node = 1,2,...,N-2
557  for ( std::size_t node = 1; node <= N - 2; node += 2 )
558  {
559  // set the current solution at this node
560  for ( unsigned var = 0; var < order; ++var )
561  {
562  //F_node[ var ] = SOLUTION[ var ][ node ];
563  F_node[ var ] = SOLUTION( node, var );
564  }
565  // set the y-pos in the eqn
566  p_EQUATION -> coord(0) = SOLUTION.coord( node );
567  // Update the equation to the nodal position
568  p_EQUATION -> update( F_node );
569  //// evaluate the RHS at the node
570  //p_EQUATION -> get_residual( R_node );
571  // step size centred at the node
572  _Xtype invh = 1. / ( SOLUTION.coord( node + 1 ) - SOLUTION.coord( node - 1 ) );
573  // loop over all the variables
574  DenseVector<_Type> temp( order, 0.0 );
575  for ( unsigned var = 0; var < order; ++var )
576  {
577  // residual
578  //temp[ var ] = p_EQUATION -> residual()[ var ] - ( SOLUTION[ var ][ node + 1 ] - SOLUTION[ var ][ node - 1 ] ) * invh;
579  temp[ var ] = p_EQUATION -> residual()[ var ] - ( SOLUTION( node + 1, var ) - SOLUTION( node - 1, var ) ) * invh;
580  }
581  Res2[ node ] = temp.inf_norm();
582  if ( Res2[ node ] > adapt_tol )
583  {
584  refine[ node ] = true;
585  }
586  else
587  if ( Res2[ node ] < TOL / 10. )
588  {
589  unrefine[ node ] = true;
590  }
591  }
592 
593  std::size_t no_refined( 0 ), no_unrefined( 0 );
594  for ( std::size_t i = 0; i < refine.size(); ++i )
595  {
596  if ( refine[ i ] == true )
597  {
598  no_refined += 1;
599  }
600  if ( unrefine[ i ] == true )
601  {
602  no_unrefined += 1;
603  }
604  }
605 #ifdef DEBUG
606  std::cout << "[ DEBUG ] Refinements = " << no_refined << "\n";
607  std::cout << "[ DEBUG ] Unrefinements = " << no_unrefined << "\n";
608 #endif
609 
610  // make a new refined/unrefined mesh
611  DenseVector<_Xtype> X( SOLUTION.nodes() );
612  DenseVector<_Xtype> newX;
613  newX.push_back( X[ 0 ] );
614  for ( std::size_t i = 1; i < N - 1; ++i )
615  {
616  if ( refine[ i ] )
617  {
618  if ( !refine[ i - 1 ] )
619  {
620  // have not already refined to the left
621  // so refine left AND right with new nodes
622  _Xtype left( X[ i - 1 ] );
623  _Xtype right( X[ i + 1 ] );
624  _Xtype here( X[ i ] );
625  newX.push_back( ( left + here ) / 2. );
626  newX.push_back( here );
627  newX.push_back( ( right + here ) / 2. );
628  }
629  else
630  {
631  // already left refined, so just refine right
632  _Xtype right( X[ i + 1 ] );
633  _Xtype here( X[ i ] );
634  newX.push_back( here );
635  newX.push_back( ( right + here ) / 2. );
636  }
637  }
638  else
639  if ( !unrefine[ i ] )
640  {
641  newX.push_back( X[ i ] );
642  }
643  }
644  newX.push_back( X[ N - 1 ] );
645 
646  // remesh the current solution to this (un)refined mesh
647  SOLUTION.remesh1( newX );
648 #ifdef TIME
649  T_REFINE.stop();
650 #endif
651  // adding nodes will screw up the sign of the determinant of the Jacobian
652  // so here we'll just make it zero
653  LAST_DET_SIGN = 0;
654 
655  std::pair< unsigned, unsigned > feedback;
656  feedback.first = no_refined;
657  feedback.second = no_unrefined;
658  return feedback;
659  }
void update(const DenseVector< _Type > &x)
A method called by arclength_solve and init_arc which stores the current converged state and paramete...
template<>
std::pair< unsigned, unsigned > CppNoddy::ODE_BVP< std::complex< double >, std::complex< double > >::adapt ( const double &  adapt_tol)

Definition at line 477 of file ODE_BVP.cpp.

478  {
479  std::string problem;
480  problem = " The ODE_BVP.adapt method has been called for a \n";
481  problem += " problem in the complex plane. This doesn't make sense \n";
482  problem += " as the path is not geometrically defined. \n";
483  throw ExceptionRuntime( problem );
484  }
template<>
std::pair< unsigned, unsigned > CppNoddy::ODE_BVP< double, std::complex< double > >::adapt ( const double &  adapt_tol)

Definition at line 488 of file ODE_BVP.cpp.

489  {
490  std::string problem;
491  problem = " The ODE_BVP.adapt method has been called for a \n";
492  problem += " problem in the complex plane. This doesn't make sense \n";
493  problem += " as the path is not geometrically defined. \n";
494  throw ExceptionRuntime( problem );
495  }
template<typename _Type , typename _Xtype >
void CppNoddy::ODE_BVP< _Type, _Xtype >::adapt_until ( const double &  adapt_tol)

Adaptively solve the system until no refinements or unrefinements are applied.

Parameters
adapt_tolThe residual tolerance at a nodal point that will lead to the mesh adaptation.

Definition at line 520 of file ODE_BVP.cpp.

References CppNoddy::ODE_BVP< _Type, _Xtype >::adapt(), and CppNoddy::ODE_BVP< _Type, _Xtype >::solve2().

Referenced by CppNoddy::ODE_BVP< D_complex >::actions_before_linear_solve().

521  {
522  std::pair< unsigned, unsigned > changes;
523  do
524  {
525  changes = adapt( adapt_tol );
526  solve2();
527  std::cout << "[INFO] Adapting: " << changes.first << " " << changes.second << "\n";
528  }
529  while ( changes.first + changes.second != 0 );
530  }
void solve2()
Formulate and solve the ODE using Newton iteration and a second-order finite difference scheme...
Definition: ODE_BVP.cpp:87
std::pair< unsigned, unsigned > adapt(const double &adapt_tol)
Adapt the computational mesh ONCE.
Definition: ODE_BVP.cpp:533
template<>
void CppNoddy::ODE_BVP< std::complex< double >, std::complex< double > >::adapt_until ( const double &  adapt_tol)

Definition at line 499 of file ODE_BVP.cpp.

500  {
501  std::string problem;
502  problem = " The ODE_BVP.adapt method has been called for a \n";
503  problem += " problem in the complex plane. This doesn't make sense \n";
504  problem += " as the path is not geometrically defined. \n";
505  throw ExceptionRuntime( problem );
506  }
template<>
void CppNoddy::ODE_BVP< double, std::complex< double > >::adapt_until ( const double &  adapt_tol)

Definition at line 510 of file ODE_BVP.cpp.

511  {
512  std::string problem;
513  problem = " The ODE_BVP.adapt method has been called for a \n";
514  problem += " problem in the complex plane. This doesn't make sense \n";
515  problem += " as the path is not geometrically defined. \n";
516  throw ExceptionRuntime( problem );
517  }
template<typename _Type , typename _Xtype >
double CppNoddy::ODE_BVP< _Type, _Xtype >::arclength_solve ( const double &  step)

Arc-length solve the system.

Before this can be called the arc_init method should have been called in order to ensure we know a solution and have derivatives w.r.t. the arc-length parameter.

Todo:
y & z should be solved for simultaneously - but to do this we'd have to extend the LAPACK interface and/or provide the same feature for the native solver.

Definition at line 176 of file ODE_BVP.cpp.

References CppNoddy::ODE_BVP< _Type, _Xtype >::actions_before_linear_solve(), CppNoddy::ArcLength_base< _Type >::arclength_residual(), CppNoddy::ArcLength_base< _Type >::ARCSTEP_MULTIPLIER, CppNoddy::BandedMatrix< _Type >::assign(), CppNoddy::Utility::dot(), CppNoddy::ArcLength_base< _Type >::ds(), CppNoddy::ArcLength_base< _Type >::DS, CppNoddy::LinearSystem_base::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::ArcLength_base< _Type >::Jac_arclength_residual(), CppNoddy::ArcLength_base< _Type >::LAST_PARAM, CppNoddy::ArcLength_base< _Type >::LAST_X, CppNoddy::ArcLength_base< _Type >::MAX_DS, CppNoddy::ArcLength_base< _Type >::p_PARAM, CppNoddy::ArcLength_base< _Type >::PARAM_DERIV_S, CppNoddy::LinearSystem_base::set_monitor_det(), CppNoddy::BandedLinearSystem< _Type >::solve(), CppNoddy::ODE_BVP< _Type, _Xtype >::solve2(), CppNoddy::ArcLength_base< _Type >::THETA, CppNoddy::ArcLength_base< _Type >::update(), CppNoddy::ArcLength_base< _Type >::X_DERIV_S, and CppNoddy::Example::z().

Referenced by CppNoddy::ODE_BVP< D_complex >::actions_before_linear_solve(), and main().

177  {
178  this -> ds() = step;
179  // order of the equation
180  unsigned order( p_EQUATION -> get_order() );
181  // the number of nodes in the BVP
182  unsigned n( SOLUTION.get_nnodes() );
183  // Banded Jacobian
184  BandedMatrix<_Type> Jac( n * order, 2 * order - 1, 0.0 );
185  // residuals over all nodes vectors
186  DenseVector<_Type> Res1( n * order, 0.0 );
187  DenseVector<_Type> Res2( n * order, 0.0 );
188  // RHS vectors for the linear solver(s)
189  DenseVector<_Type> y( n * order, 0.0 );
190  DenseVector<_Type> z( n * order, 0.0 );
191 #ifdef LAPACK
192  BandedLinearSystem<_Type> system1( &Jac, &y, "lapack" );
193  BandedLinearSystem<_Type> system2( &Jac, &z, "lapack" );
194 #else
195  BandedLinearSystem<_Type> system1( &Jac, &y, "native" );
196  BandedLinearSystem<_Type> system2( &Jac, &z, "native" );
197 #endif
198  // we use the member data 'monitor_det' to set the linear system determinant monitor
199  system1.set_monitor_det( MONITOR_DET );
200  // make backups in case we can't find a converged solution
201  DenseVector<_Type> backup_state( SOLUTION.vars_as_vector() );
202  _Type backup_parameter( *( this -> p_PARAM ) );
203  // we can generate a 1st-order guess for the next state and parameter
204  DenseVector<_Type> x( this -> LAST_X + this -> X_DERIV_S * this -> DS );
205  *( this -> p_PARAM ) = this -> LAST_PARAM + this -> PARAM_DERIV_S * this -> DS;
206  // the class keeps the solution in a OneD_mesh object, so we update it here
207  SOLUTION.set_vars_from_vector( x );
208  // determinant monitor
209  int det_sign( 0 );
210  // a flag to indicate if we have made a successful step
211  bool step_succeeded( false );
212  // iteration counter
213  int counter = 0;
214  // loop until converged or too many iterations
215  do
216  {
217  /// \todo y & z should be solved for simultaneously - but to do this
218  /// we'd have to extend the LAPACK interface and/or provide the
219  /// same feature for the native solver.
220  ++counter;
221  // extra constraint corresponding to the new unknow parameter (arclength)
222  double E1 = this -> arclength_residual( x );
223  // construct the Jacobian/residual matrix problem
224  assemble_matrix_problem( Jac, Res1 );
225  actions_before_linear_solve( Jac, Res1 );
226  //BandedMatrix<_Type> Jac_copy( Jac );
227  y = Res1;
228 #ifdef DEBUG
229  //std::cout << " [ DEBUG ] : max_residual = " << Res1.inf_norm() << "\n";
230 #endif
231  if ( Res1.inf_norm() < TOL && counter > 1 )
232  {
233  step_succeeded = true;
234  break;
235  }
236  try
237  {
238  system1.solve();
239  det_sign = system1.get_det_sign();
240  }
241  catch ( ExceptionExternal )
242  {
243  step_succeeded = false;
244  break;
245  }
246  // derivs w.r.t parameter
247  const double delta( 1.e-8 );
248  *( this -> p_PARAM ) += delta;
249  // we actually just want the Res2 & e2 vectors, so we still
250  // use the original Jacobian j ... but it was overwritten by the
251  // previous solve.
252  assemble_matrix_problem( Jac, Res2 );
253  //Jac = Jac_copy;
254  double E2 = this -> arclength_residual( x );
255  actions_before_linear_solve( Jac, Res2 );
256  *( this -> p_PARAM ) -= delta;
257  DenseVector<_Type> dRes_dp( ( Res2 - Res1 ) / delta );
258  double dE_dp = ( E2 - E1 ) / delta;
259  z = dRes_dp;
260  try
261  {
262  system2.solve();
263  }
264  catch ( ExceptionExternal )
265  {
266  step_succeeded = false;
267  break;
268  }
269  // this is the extra (full) row in the augmented matrix problem
270  // bottom row formed from dE/dx_j
271  DenseVector<_Type> Jac_E( this -> Jac_arclength_residual( x ) );
272  // given the solutions y & z, the bordering algo provides the
273  // solution to the full sparse system via
274  _Type delta_p = - ( E1 + Utility::dot( Jac_E, y ) ) /
275  ( dE_dp + Utility::dot( Jac_E, z ) );
276  DenseVector<_Type> delta_x = y + z * delta_p;
277 #ifdef DEBUG
278  std::cout << " [ DEBUG ] : Arclength_solve, dp = " << delta_p
279  << " dx.inf_norm() = " << delta_x.inf_norm()
280  << " theta = " << this -> THETA << "\n";
281 #endif
282  // update the state variables and the parameter with corrections
283  x += delta_x;
284  *( this -> p_PARAM ) += delta_p;
285  // the corrections must be put back into the mesh container
286  SOLUTION.set_vars_from_vector( x );
287  // converged?
288  if ( delta_x.inf_norm() < TOL )
289  {
290  step_succeeded = true;
291  break;
292  }
293  // too many iterations?
294  if ( counter > MAX_ITERATIONS )
295  {
296  step_succeeded = false;
297  break;
298  }
299  }
300  while ( true );
301  //
302  if ( !step_succeeded )
303  {
304  // if not a successful step then restore things
305  SOLUTION.set_vars_from_vector( backup_state );
306  *( this -> p_PARAM ) = backup_parameter;
307  // reduce our step length
308  this -> DS /= this -> ARCSTEP_MULTIPLIER;
309 #ifdef DEBUG
310  std::cout << "[ DEBUG ] : REJECTING STEP \n";
311  std::cout << "[ DEBUG ] : I decreased DS to " << this -> DS << "\n";
312 #endif
313  }
314  else
315  {
316  // update the variables needed for arc-length continuation
317  this -> update( SOLUTION.vars_as_vector() );
318  if ( LAST_DET_SIGN * det_sign < 0 )
319  {
320  // a change in the sign of the determinant of the Jacobian
321  // has been found
322  LAST_DET_SIGN = det_sign;
323  std::string problem;
324  problem = "[ INFO ] : Determinant monitor has changed signs in the Newton class.\n";
325  problem += "[ INFO ] : Bifurcation detected.\n";
326  throw ExceptionBifurcation( problem );
327  }
328 
329 #ifdef DEBUG
330  std::cout << " [ DEBUG ] : Number of iterations = " << counter << "\n";
331  std::cout << " [ DEBUG ] : Parameter p = " << *( this -> p_PARAM )
332  << "; arclength DS = " << this -> DS << "\n";
333  std::cout << " [ DEBUG ] : Arclength THETA = " << this -> THETA << "\n";
334 #endif
335  if ( counter > 8 || std::abs( this -> DS ) > this -> MAX_DS )
336  {
337  // converging too slowly, so decrease DS
338  this -> DS /= this -> ARCSTEP_MULTIPLIER;
339 #ifdef DEBUG
340  std::cout << " [ DEBUG ] : I decreased DS to " << this -> DS << "\n";
341 #endif
342  }
343  if ( counter < 4 )
344  {
345  if ( std::abs( this -> DS * this -> ARCSTEP_MULTIPLIER ) < this -> MAX_DS )
346  {
347  // converging too quickly, so increase DS
348  this -> DS *= this -> ARCSTEP_MULTIPLIER;
349 #ifdef DEBUG
350  std::cout << " [ DEBUG ] : I increased DS to " << this -> DS << "\n";
351 #endif
352  }
353  }
354  }
355  return this -> DS;
356  }
double z(const double &x)
Topography shape.
_Type dot(const DenseVector< _Type > &X, const DenseVector< _Type > &Y)
Templated dot product.
Definition: Utility.h:381
DenseVector< _Type > LAST_X
state variable at the last computed solution
virtual void actions_before_linear_solve(BandedMatrix< _Type > &a, DenseVector< _Type > &b)
A virtual method that is called prior to the linear solve stage of the solve2() method.
Definition: ODE_BVP.h:65
double arclength_residual(const DenseVector< _Type > &x) const
The extra constraint that is to be used to replace the unknown arc length.
_Type PARAM_DERIV_S
derivative of the parameter w.r.t arc length
double DS
size of the arc length step
double MAX_DS
maximum arc length step to be taken
double ARCSTEP_MULTIPLIER
step change multiplier
_Type LAST_PARAM
parameter value at the last computed solution
double THETA
the arclength theta
_Type * p_PARAM
pointer to the parameter in arclength solves
double & ds()
Return a handle to the arclength step.
void update(const DenseVector< _Type > &x)
A method called by arclength_solve and init_arc which stores the current converged state and paramete...
DenseVector< _Type > X_DERIV_S
derivative of the state variable w.r.t. arc length
DenseVector< _Type > Jac_arclength_residual(DenseVector< _Type > &x) const
The derivative of the arclength_residual function with respect to each of the state variables...
template<typename _Type, typename _Xtype >
void CppNoddy::ODE_BVP< _Type, _Xtype >::init_arc ( _Type *  p,
const double &  length,
const double &  max_length 
)

Initialise the class ready for arc length continuation.

The base class requires a vector, so we wrap the base class method here so that the vector can be extracted from the mesh member data.

Parameters
pThe pointer to the parameter
lengthThe initial arc length step to be taken (all in the parameter.
max_lengthThe maximum arc length step to be allowed.

Definition at line 184 of file ODE_BVP.h.

Referenced by CppNoddy::ODE_BVP< D_complex >::actions_before_linear_solve(), and main().

187  {
188  DenseVector<_Type> state( SOLUTION.vars_as_vector() );
189  this -> init_arc( state, p, length, max_length );
190  }
void init_arc(_Type *p, const double &length, const double &max_length)
Initialise the class ready for arc length continuation.
Definition: ODE_BVP.h:184
template<typename _Type, typename _Xtype = double>
int& CppNoddy::ODE_BVP< _Type, _Xtype >::max_itns ( )
inline

Access method to the maximum number of iterations.

Returns
A handle to the private member data MAX_ITERATIONS

Definition at line 123 of file ODE_BVP.h.

124  {
125  return MAX_ITERATIONS;
126  }
template<typename _Type , typename _Xtype >
void CppNoddy::ODE_BVP< _Type, _Xtype >::set_monitor_det ( bool  flag)

Set the flag that determines if the determinant will be monitored The default is to monitor.

Parameters
flagThe boolean value that the flag will be set to

Definition at line 178 of file ODE_BVP.h.

Referenced by CppNoddy::ODE_BVP< D_complex >::actions_before_linear_solve(), and main().

179  {
180  MONITOR_DET = flag;
181  }
template<typename _Type , typename _Xtype >
OneD_Node_Mesh< _Type, _Xtype > & CppNoddy::ODE_BVP< _Type, _Xtype >::solution ( )
inline
Returns
A handle to the solution mesh

Definition at line 193 of file ODE_BVP.h.

Referenced by CppNoddy::ODE_BVP< D_complex >::actions_before_linear_solve(), and main().

194  {
195  return SOLUTION;
196  }
template<typename _Type , typename _Xtype >
void CppNoddy::ODE_BVP< _Type, _Xtype >::solve2 ( )

Formulate and solve the ODE using Newton iteration and a second-order finite difference scheme.

The solution is stored in the publicly accessible 'solution' member data.

Definition at line 87 of file ODE_BVP.cpp.

References CppNoddy::ODE_BVP< _Type, _Xtype >::actions_before_linear_solve(), CppNoddy::LinearSystem_base::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::LinearSystem_base::set_monitor_det(), and CppNoddy::BandedLinearSystem< _Type >::solve().

Referenced by CppNoddy::ODE_BVP< D_complex >::actions_before_linear_solve(), CppNoddy::ODE_BVP< _Type, _Xtype >::adapt_until(), CppNoddy::ODE_BVP< _Type, _Xtype >::arclength_solve(), and main().

88  {
89  // the order of the problem
90  unsigned order( p_EQUATION -> get_order() );
91  // get the number of nodes in the mesh
92  // -- this may have been refined by the user since the
93  // last call.
94  unsigned n( SOLUTION.get_nnodes() );
95  // measure of maximum residual
96  double max_residual( 1.0 );
97  // iteration counter
98  int counter = 0;
99  // determinant monitor
100  int det_sign( LAST_DET_SIGN );
101 
102  // ANY LARGE STORAGE USED IN THE MAIN LOOP IS
103  // DEFINED HERE TO AVOID REPEATED CONSTRUCTION.
104  // Note we blank the A matrix after every iteration.
105  //
106  // Banded LHS matrix - max obove diagonal band width is
107  // from first variable at node i to last variable at node i+1
108  BandedMatrix<_Type> a( n * order, 2 * order - 1, 0.0 );
109  // RHS
110  DenseVector<_Type> b( n * order, 0.0 );
111  // linear solver definition
112 #ifdef LAPACK
113  BandedLinearSystem<_Type> system( &a, &b, "lapack" );
114 #else
115  BandedLinearSystem<_Type> system( &a, &b, "native" );
116 #endif
117  // pass the local monitor_det flag to the linear system
118  system.set_monitor_det( MONITOR_DET );
119 
120  // loop until converged or too many iterations
121  do
122  {
123  // iteration counter
124  ++counter;
125 #ifdef TIME
126  T_ASSEMBLE.start();
127 #endif
128  assemble_matrix_problem( a, b );
130  max_residual = b.inf_norm();
131 #ifdef DEBUG
132  std::cout << " ODE_BVP.solve : Residual_max = " << max_residual << " tol = " << TOL << "\n";
133 #endif
134 #ifdef TIME
135  T_ASSEMBLE.stop();
136  T_SOLVE.start();
137 #endif
138  // linear solver
139  system.solve();
140  // keep the solution in a OneD_Node_Mesh object
141  for ( std::size_t var = 0; var < order; ++var )
142  {
143  for ( std::size_t i = 0; i < n; ++i )
144  {
145  SOLUTION( i, var ) += b[ i * order + var ];
146  }
147  }
148 #ifdef TIME
149  T_SOLVE.stop();
150 #endif
151  }
152  while ( ( max_residual > TOL ) && ( counter < MAX_ITERATIONS ) );
153  if ( counter >= MAX_ITERATIONS )
154  {
155  std::string problem( "\n The ODE_BVP.solve2 method took too many iterations. \n" );
156  throw ExceptionItn( problem, counter, max_residual );
157  }
158 
159  // last_det_sign is set to be 0 in ctor, it must be +/-1 when calculated
160  if ( MONITOR_DET )
161  {
162  det_sign = system.get_det_sign();
163  if ( det_sign * LAST_DET_SIGN < 0 )
164  {
165  LAST_DET_SIGN = det_sign;
166  std::string problem;
167  problem = "[ INFO ] : Determinant monitor has changed signs in ODE_BVP.\n";
168  problem += "[ INFO ] : Bifurcation detected.\n";
169  throw ExceptionBifurcation( problem );
170  }
171  LAST_DET_SIGN = det_sign;
172  }
173  }
virtual void actions_before_linear_solve(BandedMatrix< _Type > &a, DenseVector< _Type > &b)
A virtual method that is called prior to the linear solve stage of the solve2() method.
Definition: ODE_BVP.h:65
template<typename _Type, typename _Xtype = double>
double& CppNoddy::ODE_BVP< _Type, _Xtype >::tolerance ( )
inline

Access method to the tolerance.

Returns
A handle to the private member data TOLERANCE

Definition at line 116 of file ODE_BVP.h.

117  {
118  return TOL;
119  }

The documentation for this class was generated from the following files:

© 2012

R.E. Hewitt