CppNoddy  0.85
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Groups Pages
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.
virtual ~ODE_BVP ()
 Destructor.
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.
void solve2 ()
 Formulate and solve the ODE using Newton iteration and a second-order finite difference scheme.
std::pair< unsigned, unsigned > adapt (const double &adapt_tol)
 Adapt the computational mesh ONCE.
void adapt_until (const double &adapt_tol)
 Adaptively solve the system until no refinements or unrefinements are applied.
void set_monitor_det (bool flag)
 Set the flag that determines if the determinant will be monitored The default is to monitor.
void init_arc (_Type *p, const double &length, const double &max_length)
 Initialise the class ready for arc length continuation.
double arclength_solve (const double &step)
 Arc-length solve the system.
OneD_Node_Mesh< _Type, _Xtype > & solution ()
double & tolerance ()
 Access method to the tolerance.
int & max_itns ()
 Access method to the maximum number of iterations.
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.
virtual void solve (DenseVector< _Type > &x)=0
 Compute a solution for that state & parameter variables that are an arc-length 'ds' from the current state.
double & ds ()
 Return a handle to the arclength step.
double & arcstep_multiplier ()
 Used to set the multiplication constant used when increasing or decreasing the arclength step.
bool & rescale_theta ()
 Handle to the RESCALE_THETA flag.
double & theta ()
 Set the arclength theta parameter.
double & desired_arc_proportion ()
 Handle to the desired proportion of the parameter to be used in the arc length solver.

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.
double arclength_residual (const DenseVector< _Type > &x) const
 The extra constraint that is to be used to replace the unknown 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.
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.
- Protected Attributes inherited from CppNoddy::ArcLength_base< _Type >
_Type * p_PARAM
 pointer to the parameter in arclength solves
DenseVector< _Type > LAST_X
 state variable at the last computed solution
DenseVector< _Type > X_DERIV_S
 derivative of the state variable w.r.t. arc length
_Type LAST_PARAM
 parameter value at the last computed solution
_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
bool INITIALISED
 for the arc-length solver - to show it has been initialised
double THETA
 the arclength theta

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.

:
ArcLength_base<_Type> (),
MAX_ITERATIONS( 12 ),
TOL( 1.e-8 ),
LAST_DET_SIGN( 0 ),
p_EQUATION( ptr_to_equation ),
p_LEFT_RESIDUAL( ptr_to_left_residual ),
p_RIGHT_RESIDUAL( ptr_to_right_residual ),
MONITOR_DET( true )
{
// set up the solution mesh using these nodes
SOLUTION = OneD_Node_Mesh<_Type, _Xtype>( nodes, p_EQUATION -> get_order() );
if ( ( p_LEFT_RESIDUAL -> get_number_of_vars() != p_EQUATION -> get_order() ) ||
( p_RIGHT_RESIDUAL -> get_number_of_vars() != p_EQUATION -> get_order() ) ||
( p_LEFT_RESIDUAL -> get_order() + p_RIGHT_RESIDUAL -> get_order() != p_EQUATION -> get_order() ) )
{
std::cout << "order " << p_EQUATION -> get_order() << "\n";
std::cout << "left nvars " << p_LEFT_RESIDUAL -> get_number_of_vars() << "\n";
std::cout << "right nvars " << p_RIGHT_RESIDUAL -> get_number_of_vars() << "\n";
std::cout << "left order " << p_LEFT_RESIDUAL -> get_order() << "\n";
std::cout << "right order " << p_RIGHT_RESIDUAL -> get_order() << "\n";
std::string problem;
problem = "It looks like the ODE_BVP equation and boundary conditions are\n";
problem += "not well posed. The number of variables for each boundary condition\n";
problem += "has to be the same as the order of the equation. Also the order of\n";
problem += "both boundary conditions has to sum to the order of the equation.\n";
throw ExceptionRuntime( problem );
}
#ifdef TIME
// timers
T_ASSEMBLE = Timer( "ODE_BVP: Assembling of the matrix (incl. equation updates):" );
T_SOLVE = Timer( "ODE_BVP: Solving of the matrix:" );
T_REFINE = Timer( "ODE_BVP: Refining the mesh:" );
#endif
}
template<typename _Type , typename _Xtype >
CppNoddy::ODE_BVP< _Type, _Xtype >::~ODE_BVP ( )
virtual

Destructor.

Definition at line 72 of file ODE_BVP.cpp.

{
#ifdef TIME
std::cout << "\n";
T_ASSEMBLE.stop();
T_ASSEMBLE.print();
T_SOLVE.stop();
T_SOLVE.print();
T_REFINE.stop();
T_REFINE.print();
#endif
}

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.

{}
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 535 of file ODE_BVP.cpp.

References CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::Example::left, CppNoddy::DenseVector< _Type >::push_back(), and CppNoddy::Example::right.

Referenced by main().

{
#ifdef TIME
T_REFINE.start();
#endif
// the order of the problem
unsigned order( p_EQUATION -> get_order() );
// get the number of nodes in the mesh
// -- this may have been refined by the user since the
// last call.
unsigned N( SOLUTION.get_nnodes() );
// row counter
//std::size_t row( 0 );
// local state variable and functions
DenseVector<_Type> F_node( order, 0.0 );
DenseVector<_Type> R_node( order, 0.0 );
// make sure (un)refine flags vector is sized
std::vector< bool > refine( N, false );
std::vector< bool > unrefine( N, false );
// Residual vector at interior nodes
DenseVector<double> Res2( N, 0.0 );
// reset row counter
//row = 0;
// inner nodes of the mesh, node = 1,2,...,N-2
for ( std::size_t node = 1; node <= N - 2; node += 2 )
{
// set the current solution at this node
for ( unsigned var = 0; var < order; ++var )
{
//F_node[ var ] = SOLUTION[ var ][ node ];
F_node[ var ] = SOLUTION( node, var );
}
// set the y-pos in the eqn
p_EQUATION -> coord(0) = SOLUTION.coord( node );
// Update the equation to the nodal position
p_EQUATION -> update( F_node );
//// evaluate the RHS at the node
//p_EQUATION -> get_residual( R_node );
// step size centred at the node
_Xtype invh = 1. / ( SOLUTION.coord( node + 1 ) - SOLUTION.coord( node - 1 ) );
// loop over all the variables
DenseVector<_Type> temp( order, 0.0 );
for ( unsigned var = 0; var < order; ++var )
{
// residual
//temp[ var ] = p_EQUATION -> residual()[ var ] - ( SOLUTION[ var ][ node + 1 ] - SOLUTION[ var ][ node - 1 ] ) * invh;
temp[ var ] = p_EQUATION -> residual()[ var ] - ( SOLUTION( node + 1, var ) - SOLUTION( node - 1, var ) ) * invh;
}
Res2[ node ] = temp.inf_norm();
if ( Res2[ node ] > adapt_tol )
{
refine[ node ] = true;
}
else
if ( Res2[ node ] < TOL / 10. )
{
unrefine[ node ] = true;
}
}
std::size_t no_refined( 0 ), no_unrefined( 0 );
for ( std::size_t i = 0; i < refine.size(); ++i )
{
if ( refine[ i ] == true )
{
no_refined += 1;
}
if ( unrefine[ i ] == true )
{
no_unrefined += 1;
}
}
#ifdef DEBUG
std::cout << "[ DEBUG ] Refinements = " << no_refined << "\n";
std::cout << "[ DEBUG ] Unrefinements = " << no_unrefined << "\n";
#endif
// make a new refined/unrefined mesh
DenseVector<_Xtype> X( SOLUTION.nodes() );
DenseVector<_Xtype> newX;
newX.push_back( X[ 0 ] );
for ( std::size_t i = 1; i < N - 1; ++i )
{
if ( refine[ i ] )
{
if ( !refine[ i - 1 ] )
{
// have not already refined to the left
// so refine left AND right with new nodes
_Xtype left( X[ i - 1 ] );
_Xtype right( X[ i + 1 ] );
_Xtype here( X[ i ] );
newX.push_back( ( left + here ) / 2. );
newX.push_back( here );
newX.push_back( ( right + here ) / 2. );
}
else
{
// already left refined, so just refine right
_Xtype right( X[ i + 1 ] );
_Xtype here( X[ i ] );
newX.push_back( here );
newX.push_back( ( right + here ) / 2. );
}
}
else
if ( !unrefine[ i ] )
{
newX.push_back( X[ i ] );
}
}
newX.push_back( X[ N - 1 ] );
// remesh the current solution to this (un)refined mesh
SOLUTION.remesh1( newX );
#ifdef TIME
T_REFINE.stop();
#endif
// adding nodes will screw up the sign of the determinant of the Jacobian
// so here we'll just make it zero
LAST_DET_SIGN = 0;
std::pair< unsigned, unsigned > feedback;
feedback.first = no_refined;
feedback.second = no_unrefined;
return feedback;
}
template<>
std::pair< unsigned, unsigned > CppNoddy::ODE_BVP< std::complex< double >, std::complex< double > >::adapt ( const double &  adapt_tol)

Definition at line 479 of file ODE_BVP.cpp.

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

Definition at line 490 of file ODE_BVP.cpp.

{
std::string problem;
problem = " The ODE_BVP.adapt method has been called for a \n";
problem += " problem in the complex plane. This doesn't make sense \n";
problem += " as the path is not geometrically defined. \n";
throw ExceptionRuntime( problem );
}
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 522 of file ODE_BVP.cpp.

{
std::pair< unsigned, unsigned > changes;
do
{
changes = adapt( adapt_tol );
solve2();
std::cout << "[INFO] Adapting: " << changes.first << " " << changes.second << "\n";
}
while ( changes.first + changes.second != 0 );
}
template<>
void CppNoddy::ODE_BVP< std::complex< double >, std::complex< double > >::adapt_until ( const double &  adapt_tol)

Definition at line 501 of file ODE_BVP.cpp.

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

Definition at line 512 of file ODE_BVP.cpp.

{
std::string problem;
problem = " The ODE_BVP.adapt method has been called for a \n";
problem += " problem in the complex plane. This doesn't make sense \n";
problem += " as the path is not geometrically defined. \n";
throw ExceptionRuntime( problem );
}
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::Utility::dot(), CppNoddy::LinearSystem_base::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::LinearSystem_base::set_monitor_det(), CppNoddy::BandedLinearSystem< _Type >::solve(), and CppNoddy::Example::z().

Referenced by main().

{
this -> ds() = step;
// order of the equation
unsigned order( p_EQUATION -> get_order() );
// the number of nodes in the BVP
unsigned n( SOLUTION.get_nnodes() );
// Banded Jacobian
BandedMatrix<_Type> Jac( n * order, 2 * order - 1, 0.0 );
// residuals over all nodes vectors
DenseVector<_Type> Res1( n * order, 0.0 );
DenseVector<_Type> Res2( n * order, 0.0 );
// RHS vectors for the linear solver(s)
DenseVector<_Type> y( n * order, 0.0 );
DenseVector<_Type> z( n * order, 0.0 );
#ifdef LAPACK
BandedLinearSystem<_Type> system1( &Jac, &y, "lapack" );
BandedLinearSystem<_Type> system2( &Jac, &z, "lapack" );
#else
BandedLinearSystem<_Type> system1( &Jac, &y, "native" );
BandedLinearSystem<_Type> system2( &Jac, &z, "native" );
#endif
// we use the member data 'monitor_det' to set the linear system determinant monitor
system1.set_monitor_det( MONITOR_DET );
// make backups in case we can't find a converged solution
DenseVector<_Type> backup_state( SOLUTION.vars_as_vector() );
_Type backup_parameter( *( this -> p_PARAM ) );
// we can generate a 1st-order guess for the next state and parameter
DenseVector<_Type> x( this -> LAST_X + this -> X_DERIV_S * this -> DS );
*( this -> p_PARAM ) = this -> LAST_PARAM + this -> PARAM_DERIV_S * this -> DS;
// the class keeps the solution in a OneD_mesh object, so we update it here
SOLUTION.set_vars_from_vector( x );
// determinant monitor
int det_sign( 0 );
// a flag to indicate if we have made a successful step
bool step_succeeded( false );
// iteration counter
int counter = 0;
// loop until converged or too many iterations
do
{
/// \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.
++counter;
// extra constraint corresponding to the new unknow parameter (arclength)
double E1 = this -> arclength_residual( x );
// construct the Jacobian/residual matrix problem
assemble_matrix_problem( Jac, Res1 );
//BandedMatrix<_Type> Jac_copy( Jac );
y = Res1;
#ifdef DEBUG
//std::cout << " [ DEBUG ] : max_residual = " << Res1.inf_norm() << "\n";
#endif
if ( Res1.inf_norm() < TOL && counter > 1 )
{
step_succeeded = true;
break;
}
try
{
system1.solve();
det_sign = system1.get_det_sign();
}
catch ( ExceptionExternal )
{
step_succeeded = false;
break;
}
// derivs w.r.t parameter
const double delta( 1.e-8 );
*( this -> p_PARAM ) += delta;
// we actually just want the Res2 & e2 vectors, so we still
// use the original Jacobian j ... but it was overwritten by the
// previous solve.
assemble_matrix_problem( Jac, Res2 );
//Jac = Jac_copy;
double E2 = this -> arclength_residual( x );
*( this -> p_PARAM ) -= delta;
DenseVector<_Type> dRes_dp( ( Res2 - Res1 ) / delta );
double dE_dp = ( E2 - E1 ) / delta;
z = dRes_dp;
try
{
system2.solve();
}
catch ( ExceptionExternal )
{
step_succeeded = false;
break;
}
// this is the extra (full) row in the augmented matrix problem
// bottom row formed from dE/dx_j
DenseVector<_Type> Jac_E( this -> Jac_arclength_residual( x ) );
// given the solutions y & z, the bordering algo provides the
// solution to the full sparse system via
_Type delta_p = - ( E1 + Utility::dot( Jac_E, y ) ) /
( dE_dp + Utility::dot( Jac_E, z ) );
DenseVector<_Type> delta_x = y + z * delta_p;
#ifdef DEBUG
std::cout << " [ DEBUG ] : Arclength_solve, dp = " << delta_p
<< " dx.inf_norm() = " << delta_x.inf_norm()
<< " theta = " << this -> THETA << "\n";
#endif
// update the state variables and the parameter with corrections
x += delta_x;
*( this -> p_PARAM ) += delta_p;
// the corrections must be put back into the mesh container
SOLUTION.set_vars_from_vector( x );
// converged?
if ( delta_x.inf_norm() < TOL )
{
step_succeeded = true;
break;
}
// too many iterations?
if ( counter > MAX_ITERATIONS )
{
step_succeeded = false;
break;
}
}
while ( true );
//
if ( !step_succeeded )
{
// if not a successful step then restore things
SOLUTION.set_vars_from_vector( backup_state );
*( this -> p_PARAM ) = backup_parameter;
// reduce our step length
this -> DS /= this -> ARCSTEP_MULTIPLIER;
#ifdef DEBUG
std::cout << "[ DEBUG ] : REJECTING STEP \n";
std::cout << "[ DEBUG ] : I decreased DS to " << this -> DS << "\n";
#endif
}
else
{
// update the variables needed for arc-length continuation
this -> update( SOLUTION.vars_as_vector() );
if ( LAST_DET_SIGN * det_sign < 0 )
{
// a change in the sign of the determinant of the Jacobian
// has been found
LAST_DET_SIGN = det_sign;
std::string problem;
problem = "[ INFO ] : Determinant monitor has changed signs in the Newton class.\n";
problem += "[ INFO ] : Bifurcation detected.\n";
throw ExceptionBifurcation( problem );
}
#ifdef DEBUG
std::cout << " [ DEBUG ] : Number of iterations = " << counter << "\n";
std::cout << " [ DEBUG ] : Parameter p = " << *( this -> p_PARAM )
<< "; arclength DS = " << this -> DS << "\n";
std::cout << " [ DEBUG ] : Arclength THETA = " << this -> THETA << "\n";
#endif
if ( counter > 8 || std::abs( this -> DS ) > this -> MAX_DS )
{
// converging too slowly, so decrease DS
this -> DS /= this -> ARCSTEP_MULTIPLIER;
#ifdef DEBUG
std::cout << " [ DEBUG ] : I decreased DS to " << this -> DS << "\n";
#endif
}
if ( counter < 4 )
{
if ( std::abs( this -> DS * this -> ARCSTEP_MULTIPLIER ) < this -> MAX_DS )
{
// converging too quickly, so increase DS
this -> DS *= this -> ARCSTEP_MULTIPLIER;
#ifdef DEBUG
std::cout << " [ DEBUG ] : I increased DS to " << this -> DS << "\n";
#endif
}
}
}
return this -> DS;
}
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 main().

{
DenseVector<_Type> state( SOLUTION.vars_as_vector() );
this -> init_arc( state, p, length, max_length );
}
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.

{
return MAX_ITERATIONS;
}
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 main().

{
MONITOR_DET = flag;
}
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 main().

{
return SOLUTION;
}
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::LinearSystem_base::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::LinearSystem_base::set_monitor_det(), and CppNoddy::BandedLinearSystem< _Type >::solve().

Referenced by main().

{
// the order of the problem
unsigned order( p_EQUATION -> get_order() );
// get the number of nodes in the mesh
// -- this may have been refined by the user since the
// last call.
unsigned n( SOLUTION.get_nnodes() );
// measure of maximum residual
double max_residual( 1.0 );
// iteration counter
int counter = 0;
// determinant monitor
int det_sign( LAST_DET_SIGN );
// ANY LARGE STORAGE USED IN THE MAIN LOOP IS
// DEFINED HERE TO AVOID REPEATED CONSTRUCTION.
// Note we blank the A matrix after every iteration.
//
// Banded LHS matrix - max obove diagonal band width is
// from first variable at node i to last variable at node i+1
BandedMatrix<_Type> a( n * order, 2 * order - 1, 0.0 );
// RHS
DenseVector<_Type> b( n * order, 0.0 );
// linear solver definition
#ifdef LAPACK
BandedLinearSystem<_Type> system( &a, &b, "lapack" );
#else
BandedLinearSystem<_Type> system( &a, &b, "native" );
#endif
// pass the local monitor_det flag to the linear system
system.set_monitor_det( MONITOR_DET );
// loop until converged or too many iterations
do
{
// iteration counter
++counter;
#ifdef TIME
T_ASSEMBLE.stop();
#endif
assemble_matrix_problem( a, b );
max_residual = b.inf_norm();
#ifdef DEBUG
std::cout << " ODE_BVP.solve : Residual_max = " << max_residual << " tol = " << TOL << "\n";
#endif
#ifdef TIME
T_ASSEMBLE.start();
T_SOLVE.start();
#endif
// linear solver
system.solve();
// keep the solution in a OneD_Node_Mesh object
for ( std::size_t var = 0; var < order; ++var )
{
for ( std::size_t i = 0; i < n; ++i )
{
SOLUTION( i, var ) += b[ i * order + var ];
}
}
#ifdef TIME
T_SOLVE.stop();
#endif
}
while ( ( max_residual > TOL ) && ( counter < MAX_ITERATIONS ) );
if ( counter >= MAX_ITERATIONS )
{
std::string problem( "\n The ODE_BVP.solve2 method took too many iterations. \n" );
throw ExceptionItn( problem, counter, max_residual );
}
// last_det_sign is set to be 0 in ctor, it must be +/-1 when calculated
if ( MONITOR_DET )
{
det_sign = system.get_det_sign();
if ( det_sign * LAST_DET_SIGN < 0 )
{
LAST_DET_SIGN = det_sign;
std::string problem;
problem = "[ INFO ] : Determinant monitor has changed signs in ODE_BVP.\n";
problem += "[ INFO ] : Bifurcation detected.\n";
throw ExceptionBifurcation( problem );
}
LAST_DET_SIGN = det_sign;
}
}
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.

{
return TOL;
}

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

© 2012

R.E. Hewitt