CppNoddy  0.85
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
CppNoddy::PDE_IBVP< _Type > Class Template Reference

A templated object for real/complex vector system of unsteady equations. More...

#include <PDE_IBVP.h>

Inheritance diagram for CppNoddy::PDE_IBVP< _Type >:
CppNoddy::Uncopyable

Public Member Functions

 PDE_IBVP (Equation_2matrix< _Type > *equation_ptr, const DenseVector< double > &nodes, Residual_with_coords< _Type > *ptr_to_left_residual, Residual_with_coords< _Type > *ptr_to_right_residual)
 The class is defined by a vector function for the system.
 ~PDE_IBVP ()
 Destructor.
void step2 (const double &dt)
 A Crank-Nicolson 'time' stepper.
void assemble_matrix_problem (BandedMatrix< _Type > &a, DenseVector< _Type > &b, const double &dt)
 Assembles the matrix problem for a BVP solve at the current time level.
double & coord ()
 Return a reference to the current value of the 'timelike/parabolic' coordinate.
OneD_Node_Mesh< _Type > & solution ()
double & tolerance ()
 Access method to the tolerance.
int & max_itns ()
 Access method to the maximum number of iterations.

Additional Inherited Members

Detailed Description

template<typename _Type>
class CppNoddy::PDE_IBVP< _Type >

A templated object for real/complex vector system of unsteady equations.

Definition at line 38 of file PDE_IBVP.h.

Constructor & Destructor Documentation

template<typename _Type >
CppNoddy::PDE_IBVP< _Type >::PDE_IBVP ( Equation_2matrix< _Type > *  equation_ptr,
const DenseVector< double > &  nodes,
Residual_with_coords< _Type > *  ptr_to_left_residual,
Residual_with_coords< _Type > *  ptr_to_right_residual 
)

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

Parameters
equation_ptrA pointer to an inherited Equation 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 35 of file PDE_IBVP.cpp.

References CppNoddy::PDE_IBVP< _Type >::coord().

:
TOL( 1.e-8 ),
T( 0.0 ),
MAX_ITERATIONS( 12 ),
p_EQUATION( ptr_to_equation ),
p_LEFT_RESIDUAL( ptr_to_left_residual ),
p_RIGHT_RESIDUAL( ptr_to_right_residual )
{
SOLN = OneD_Node_Mesh<_Type>( nodes, p_EQUATION -> get_order() );
PREV_SOLN = SOLN;
p_EQUATION -> coord( 1 ) = T;
if ( p_EQUATION -> get_order() - p_LEFT_RESIDUAL -> get_order() - p_RIGHT_RESIDUAL -> get_order() != 0 )
{
std::string problem( "\n The PDE_IBVP class has been constructed, but the order of the \n");
problem += "system does not match the number of boundary conditions.\n";
throw ExceptionRuntime( problem );
}
#ifdef TIME
// timers
T_ASSEMBLE = Timer( "Assembling of the matrix (incl. equation updates):" );
T_SOLVE = Timer( "Solving of the matrix:" );
#endif
}
template<typename _Type >
CppNoddy::PDE_IBVP< _Type >::~PDE_IBVP ( )

Destructor.

Definition at line 63 of file PDE_IBVP.cpp.

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

Member Function Documentation

template<typename _Type >
void CppNoddy::PDE_IBVP< _Type >::assemble_matrix_problem ( BandedMatrix< _Type > &  a,
DenseVector< _Type > &  b,
const double &  dt 
)

Assembles the matrix problem for a BVP solve at the current time level.

Parameters
aThe LHS (banded) matrix.
bThe RHS (dense) vector.
dtThe 'time step' to be taken.

Definition at line 152 of file PDE_IBVP.cpp.

References CppNoddy::BandedMatrix< _Type >::assign(), CppNoddy::Utility::dot(), CppNoddy::BandedMatrix< _Type >::get_elt_iter(), CppNoddy::Example::left, CppNoddy::BandedMatrix< _Type >::noffdiag(), and CppNoddy::Example::right.

{
// clear the Jacobian matrix
a.assign( 0.0 );
// inverse of the time step
const double inv_dt( 1. / dt );
// the order of the problem
const unsigned order( p_EQUATION -> get_order() );
// number of spatial nodes
const unsigned ny( SOLN.get_nnodes() );
// row counter
std::size_t row( 0 );
// a matrix that is used in the Jacobian of the mass matrix terms
DenseMatrix<_Type> h0( order, order, 0.0 );
DenseMatrix<_Type> h1( order, order, 0.0 );
// local state variable and functions
DenseVector<_Type> F_midpt( order, 0.0 );
DenseVector<_Type> O_midpt( order, 0.0 );
DenseVector<_Type> state( order, 0.0 );
DenseVector<_Type> state_dt( order, 0.0 );
DenseVector<_Type> state_dy( order, 0.0 );
// BCn equation is evaluated at the next time step
p_LEFT_RESIDUAL -> coord( 0 ) = T + dt;
// update the BC residuals for the current iteration
p_LEFT_RESIDUAL -> update( SOLN.get_nodes_vars( 0 ) );
// add the (linearised) LHS BCs to the matrix problem
for ( unsigned i = 0; i < p_LEFT_RESIDUAL -> get_order(); ++i )
{
// loop thru variables at LHS of the domain
for ( unsigned var = 0; var < order; ++var )
{
a( row, var ) = p_LEFT_RESIDUAL -> jacobian()( i, var );
}
b[ row ] = - p_LEFT_RESIDUAL -> residual()[ i ];
++row;
}
// inner nodes of the mesh, node = 0,1,2,...,N-2
for ( std::size_t node = 0; node <= ny - 2; ++node )
{
const std::size_t l_node = node;
const std::size_t r_node = node + 1;
// inverse of step size
const double inv_dy = 1. / ( SOLN.coord( r_node ) - SOLN.coord( l_node ) );
// set the current solution at this node by 2nd order evaluation at mid point
for ( unsigned var = 0; var < order; ++var )
{
const _Type F_midpt = ( SOLN( l_node, var ) + SOLN( r_node, var ) ) / 2.;
const _Type O_midpt = ( PREV_SOLN( l_node, var ) + PREV_SOLN( r_node, var ) ) / 2.;
state_dy[ var ] = ( SOLN( r_node, var ) - SOLN( l_node, var )
+ PREV_SOLN( r_node, var ) - PREV_SOLN( l_node, var ) ) * inv_dy / 2.;
state[ var ] = ( F_midpt + O_midpt ) / 2.;
state_dt[ var ] = ( F_midpt - O_midpt ) * inv_dt;
}
// set the equation's y & t values to be mid points
// y
p_EQUATION -> coord(0) = 0.5 * ( SOLN.coord( l_node ) + SOLN.coord( r_node ) );
// t
p_EQUATION -> coord(1) = T + dt / 2;
// Update the equation to the mid point position
p_EQUATION -> update( state );
// evaluate the Jacobian of mass contribution multiplied by state_dy
p_EQUATION -> get_jacobian_of_matrix0_mult_vector( state, state_dy, h0 );
// evaluate the Jacobian of mass contribution multiplied by state_dt
p_EQUATION -> get_jacobian_of_matrix1_mult_vector( state, state_dt, h1 );
// loop over all the variables
//
// to avoid repeated mapping arithmetic within operator() of the
// BandedMatrix class we'll access the matrix with the iterator.
//
// il = position for (0, lnode*order)
typename BandedMatrix<_Type>::elt_iter l_iter( a.get_elt_iter( 0, l_node * order ) );
// ir = position for (0, rnode*order)
typename BandedMatrix<_Type>::elt_iter r_iter( a.get_elt_iter( 0, r_node * order ) );
//
for ( unsigned var = 0; var < order; ++var )
{
// add the matrix mult terms to the linearised problem
for ( unsigned i = 0; i < order; ++i ) // dummy index
{
// add the Jacobian terms
// indirect access: a( row, order * l_node + i ) -= jac_midpt( var, i ) * 0.5;
const std::size_t offset( i * a.noffdiag() * 3 + row );
const typename BandedMatrix<_Type>::elt_iter left( l_iter + offset );
const typename BandedMatrix<_Type>::elt_iter right( r_iter + offset );
//
*left -= p_EQUATION -> jacobian()( var, i ) * 0.5;
// indirect access: a( row, order * r_node + i ) -= jac_midpt( var, i ) * 0.5;
*right -= p_EQUATION -> jacobian()( var, i ) * 0.5;
// add the Jacobian of mass terms for dt terms
// indirect access: a( row, order * l_node + i ) += h1( var, i ) * 0.5;
*left += h0( var, i ) * 0.5;
// indirect access: a( row, order * r_node + i ) += h1( var, i ) * 0.5;
*right += h0( var, i ) * 0.5;
// add the Jacobian of mass terms for dt terms
// indirect access: a( row, order * l_node + i ) += h2( var, i ) * 0.5;
*left += h1( var, i ) * 0.5;
// indirect access: a( row, order * r_node + i ) += h2( var, i ) * 0.5;
*right += h1( var, i ) * 0.5;
// add the mass matrix terms
// indirect access: a( row, order * l_node + i ) -= mass_midpt( var, i ) * inv_dy;
*left -= p_EQUATION -> matrix0()( var, i ) * inv_dy;
// indirect access: a( row, order * r_node + i ) += mass_midpt( var, i ) * inv_dy;
*right += p_EQUATION -> matrix0()( var, i ) * inv_dy;
// indirect access: a( row, order * l_node + i ) += mass_midpt( var, i ) * inv_dt;
*left += p_EQUATION -> matrix1()( var, i ) * inv_dt;
// indirect access: a( row, order * r_node + i ) += mass_midpt( var, i ) * inv_dt;
*right += p_EQUATION -> matrix1()( var, i ) * inv_dt;
}
// RHS
b[ row ] = p_EQUATION -> residual()[ var ];
b[ row ] -= Utility::dot( p_EQUATION -> matrix0()[ var ], state_dy );
b[ row ] -= Utility::dot( p_EQUATION -> matrix1()[ var ], state_dt );
b[ row ] *= 2;
// increment the row
row += 1;
}
}
// BCn equation is evaluated at the next step time point as
// they cannot depend on d/dt terms
p_RIGHT_RESIDUAL -> coord( 0 ) = T + dt;
// update the BC residuals for the current iteration
p_RIGHT_RESIDUAL -> update( SOLN.get_nodes_vars( ny - 1 ) );
// add the (linearised) RHS BCs to the matrix problem
for ( unsigned i = 0; i < p_RIGHT_RESIDUAL -> get_order(); ++i )
{
// loop thru variables at RHS of the domain
for ( unsigned var = 0; var < order; ++var )
{
a( row, order * ( ny - 1 ) + var ) = p_RIGHT_RESIDUAL -> jacobian()( i, var );
}
b[ row ] = - p_RIGHT_RESIDUAL -> residual()[ i ];
++row;
}
#ifdef PARANOID
if ( row != ny * order )
{
std::string problem( "\n The ODE_BVP has an incorrect number of boundary conditions. \n" );
throw ExceptionRuntime( problem );
}
#endif
}
template<typename _Type>
double& CppNoddy::PDE_IBVP< _Type >::coord ( )
inline

Return a reference to the current value of the 'timelike/parabolic' coordinate.

Returns
A handle to the current time stored in the object

Definition at line 67 of file PDE_IBVP.h.

Referenced by main(), and CppNoddy::PDE_IBVP< _Type >::PDE_IBVP().

{
return T;
}
template<typename _Type>
int& CppNoddy::PDE_IBVP< _Type >::max_itns ( )
inline

Access method to the maximum number of iterations.

Returns
A handle to the private member data MAX_ITERATIONS

Definition at line 84 of file PDE_IBVP.h.

{
return MAX_ITERATIONS;
}
template<typename _Type >
OneD_Node_Mesh< _Type > & CppNoddy::PDE_IBVP< _Type >::solution ( )
inline
Returns
A handle to the solution mesh

Definition at line 117 of file PDE_IBVP.h.

Referenced by main().

{
return SOLN;
}
template<typename _Type >
void CppNoddy::PDE_IBVP< _Type >::step2 ( const double &  dt)

A Crank-Nicolson 'time' stepper.

Definition at line 75 of file PDE_IBVP.cpp.

References CppNoddy::DenseVector< _Type >::inf_norm(), 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 ny( SOLN.get_nnodes() );
// measure of maximum residual
double max_residual( 1.0 );
// iteration counter
int counter = 0;
// store this soln as the 'previous SOLUTION'
PREV_SOLN = SOLN;
// 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( ny * order, 2 * order - 1, 0.0 );
// RHS
DenseVector<_Type> b( ny * order, 0.0 );
// linear solver definition
#ifdef LAPACK
BandedLinearSystem<_Type> system( &a, &b, "lapack" );
#else
BandedLinearSystem<_Type> system( &a, &b, "native" );
#endif
// loop until converged or too many iterations
do
{
// iteration counter
++counter;
#ifdef TIME
T_ASSEMBLE.start();
#endif
max_residual = b.inf_norm();
#ifdef DEBUG
std::cout << " PDE_IBVP.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_GenMesh object
for ( std::size_t var = 0; var < order; ++var )
{
for ( std::size_t i = 0; i < ny; ++i )
{
SOLN( i, var ) += b[ i * order + var ];
}
}
#ifdef TIME
T_SOLVE.stop();
#endif
}
while ( ( max_residual > TOL ) && ( counter < MAX_ITERATIONS ) );
if ( max_residual > TOL )
{
// restore to previous state because this step failed
SOLN = PREV_SOLN;
std::string problem( "\n The PDE_IBVP.step2 method took too many iterations.\n");
problem += "Solution has been restored to the previous accurate state.\n";
throw ExceptionItn( problem, counter, max_residual );
}
#ifdef DEBUG
std::cout << "[DEBUG] time-like variable = " << T << "\n";
#endif
// set the time to the updated level
T += dt;
}
template<typename _Type>
double& CppNoddy::PDE_IBVP< _Type >::tolerance ( )
inline

Access method to the tolerance.

Returns
A handle to the private member data TOL

Definition at line 77 of file PDE_IBVP.h.

{
return TOL;
}

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

© 2012

R.E. Hewitt