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 >:

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...
std::pair< unsigned, unsigned >	adapt (const double &adapt_tol)
std::pair< unsigned, unsigned >	adapt (const double &adapt_tol)
void	adapt_until (const double &adapt_tol)
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...
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. 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 37 of file ODE_BVP.h.

Constructor & Destructor Documentation

ODE_BVP()

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_equation	A pointer to an Equation_1matrix object.
nodes	A vector that defines the nodal positions.
ptr_to_left_residual	A pointer to a residual object that defines the LHS boundary conditions.
ptr_to_right_residual	A pointer to a residual object that defines the RHS boundary conditions.

Definition at line 31 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
  }

~ODE_BVP()

template<typename _Type , typename _Xtype >

CppNoddy::ODE_BVP< _Type, _Xtype >::~ODE_BVP

virtual

Destructor.

Definition at line 69 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

actions_before_linear_solve()

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

a	The Jacbian matrix to be passed to the linear solver
b	The residual vector to be passed to the linear solver

Definition at line 63 of file ODE_BVP.h.

    {}

adapt() [1/3]

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_tol

The 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 481 of file ODE_BVP.cpp.

                                                                                     {
#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;
  }

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

Referenced by main().

adapt() [2/3]

std::pair< unsigned, unsigned > CppNoddy::ODE_BVP< std::complex< double >, std::complex< double > >::adapt

(

const double &

adapt_tol

)

Definition at line 432 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);
  }

adapt() [3/3]

std::pair< unsigned, unsigned > CppNoddy::ODE_BVP< double, std::complex< double > >::adapt

(

const double &

adapt_tol

)

Definition at line 442 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);
  }

adapt_until() [1/3]

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_tol

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

Definition at line 471 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);
  }

adapt_until() [2/3]

void CppNoddy::ODE_BVP< std::complex< double >, std::complex< double > >::adapt_until

(

const double &

adapt_tol

)

Definition at line 452 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);
  }

adapt_until() [3/3]

void CppNoddy::ODE_BVP< double, std::complex< double > >::adapt_until

(

const double &

adapt_tol

)

Definition at line 462 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);
  }

arclength_solve()

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 164 of file ODE_BVP.cpp.

                                                                   {
    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);
      actions_before_linear_solve(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(const ExceptionExternal &error) {
        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);
      actions_before_linear_solve(Jac, Res2);
      *(this -> p_PARAM)  -= delta;
      DenseVector<_Type> dRes_dp((Res2 - Res1) / delta);
      double dE_dp = (E2 - E1) / delta;
      z = dRes_dp;
      try {
        system2.solve();
      } catch(const ExceptionExternal& error) {
        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) || (delta_x.inf_norm() > 100 )) {
        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;
  }

References CppNoddy::Utility::dot(), CppNoddy::BandedLinearSystem< _Type >::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::BandedLinearSystem< _Type >::set_monitor_det(), and CppNoddy::BandedLinearSystem< _Type >::solve().

Referenced by main().

init_arc()

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

p	The pointer to the parameter
length	The initial arc length step to be taken (all in the parameter.
max_length	The maximum arc length step to be allowed.

Definition at line 179 of file ODE_BVP.h.

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

References p.

Referenced by main().

max_itns()

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 120 of file ODE_BVP.h.

                    {
      return MAX_ITERATIONS;
    }

set_monitor_det()

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

flag	The boolean value that the flag will be set to

Definition at line 174 of file ODE_BVP.h.

                                                        {
    MONITOR_DET = flag;
  }

Referenced by main().

solution()

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 187 of file ODE_BVP.h.

                                                                         {
    return SOLUTION;
  }

Referenced by main().

solve2()

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 83 of file ODE_BVP.cpp.

                                      {
    // 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.start();
#endif
      assemble_matrix_problem(a, b);
      actions_before_linear_solve(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.stop();
      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;
    }
  }

Referenced by main().

tolerance()

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 114 of file ODE_BVP.h.

                        {
      return TOL;
    }

---

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

• include/ODE_BVP.h
• src/ODE_BVP.cpp