CppNoddy  0.90
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. More...
 
 ~PDE_IBVP ()
 Destructor. More...
 
void step2 (const double &dt)
 A Crank-Nicolson 'time' stepper. More...
 
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. More...
 
double & coord ()
 Return a reference to the current value of the 'timelike/parabolic' coordinate. More...
 
OneD_Node_Mesh< _Type > & solution ()
 
OneD_Node_Mesh< _Type > & previous_solution ()
 
double & tolerance ()
 Access method to the tolerance. More...
 
int & max_itns ()
 Access method to the maximum number of iterations. More...
 

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().

38  :
39  TOL( 1.e-8 ),
40  T( 0.0 ),
41  MAX_ITERATIONS( 12 ),
42  p_EQUATION( ptr_to_equation ),
43  p_LEFT_RESIDUAL( ptr_to_left_residual ),
44  p_RIGHT_RESIDUAL( ptr_to_right_residual )
45  {
46  SOLN = OneD_Node_Mesh<_Type>( nodes, p_EQUATION -> get_order() );
47  PREV_SOLN = SOLN;
48  //
49  p_EQUATION -> coord( 1 ) = T;
50  if ( p_EQUATION -> get_order() - p_LEFT_RESIDUAL -> get_order() - p_RIGHT_RESIDUAL -> get_order() != 0 )
51  {
52  std::string problem( "\n The PDE_IBVP class has been constructed, but the order of the \n");
53  problem += "system does not match the number of boundary conditions.\n";
54  throw ExceptionRuntime( problem );
55  }
56 #ifdef TIME
57  // timers
58  T_ASSEMBLE = Timer( "Assembling of the matrix (incl. equation updates):" );
59  T_SOLVE = Timer( "Solving of the matrix:" );
60 #endif
61  }
double & coord()
Return a reference to the current value of the &#39;timelike/parabolic&#39; coordinate.
Definition: PDE_IBVP.h:67
template<typename _Type >
CppNoddy::PDE_IBVP< _Type >::~PDE_IBVP ( )

Destructor.

Definition at line 64 of file PDE_IBVP.cpp.

65  {
66 #ifdef TIME
67  std::cout << "\n";
68  T_ASSEMBLE.stop();
69  T_ASSEMBLE.print();
70  T_SOLVE.stop();
71  T_SOLVE.print();
72 #endif
73  }

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 153 of file PDE_IBVP.cpp.

References CppNoddy::BandedMatrix< _Type >::assign(), CppNoddy::PDE_IBVP< _Type >::coord(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::coord(), CppNoddy::Utility::dot(), CppNoddy::BandedMatrix< _Type >::get_elt_iter(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::get_nnodes(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::get_nodes_vars(), and CppNoddy::BandedMatrix< _Type >::noffdiag().

Referenced by CppNoddy::PDE_IBVP< _Type >::step2().

154  {
155  // clear the Jacobian matrix
156  a.assign( 0.0 );
157  // inverse of the time step
158  const double inv_dt( 1. / dt );
159  // the order of the problem
160  const unsigned order( p_EQUATION -> get_order() );
161  // number of spatial nodes
162  const unsigned ny( SOLN.get_nnodes() );
163  // row counter
164  std::size_t row( 0 );
165  // a matrix that is used in the Jacobian of the mass matrix terms
166  DenseMatrix<_Type> h0( order, order, 0.0 );
167  DenseMatrix<_Type> h1( order, order, 0.0 );
168  // local state variable and functions
169  DenseVector<_Type> F_midpt( order, 0.0 );
170  DenseVector<_Type> O_midpt( order, 0.0 );
171  DenseVector<_Type> state( order, 0.0 );
172  DenseVector<_Type> state_dt( order, 0.0 );
173  DenseVector<_Type> state_dy( order, 0.0 );
174  // BCn equation is evaluated at the next time step
175  p_LEFT_RESIDUAL -> coord( 0 ) = T + dt;
176  // update the BC residuals for the current iteration
177  p_LEFT_RESIDUAL -> update( SOLN.get_nodes_vars( 0 ) );
178  // add the (linearised) LHS BCs to the matrix problem
179  for ( unsigned i = 0; i < p_LEFT_RESIDUAL -> get_order(); ++i )
180  {
181  // loop thru variables at LHS of the domain
182  for ( unsigned var = 0; var < order; ++var )
183  {
184  a( row, var ) = p_LEFT_RESIDUAL -> jacobian()( i, var );
185  }
186  b[ row ] = - p_LEFT_RESIDUAL -> residual()[ i ];
187  ++row;
188  }
189  // inner nodes of the mesh, node = 0,1,2,...,N-2
190  for ( std::size_t node = 0; node <= ny - 2; ++node )
191  {
192  const std::size_t l_node = node;
193  const std::size_t r_node = node + 1;
194  // inverse of step size
195  const double inv_dy = 1. / ( SOLN.coord( r_node ) - SOLN.coord( l_node ) );
196  // set the current solution at this node by 2nd order evaluation at mid point
197  for ( unsigned var = 0; var < order; ++var )
198  {
199  const _Type F_midpt = ( SOLN( l_node, var ) + SOLN( r_node, var ) ) / 2.;
200  const _Type O_midpt = ( PREV_SOLN( l_node, var ) + PREV_SOLN( r_node, var ) ) / 2.;
201  state_dy[ var ] = ( SOLN( r_node, var ) - SOLN( l_node, var )
202  + PREV_SOLN( r_node, var ) - PREV_SOLN( l_node, var ) ) * inv_dy / 2.;
203  state[ var ] = ( F_midpt + O_midpt ) / 2.;
204  state_dt[ var ] = ( F_midpt - O_midpt ) * inv_dt;
205  }
206  // set the equation's y & t values to be mid points
207  // y
208  p_EQUATION -> coord(0) = 0.5 * ( SOLN.coord( l_node ) + SOLN.coord( r_node ) );
209  // t
210  p_EQUATION -> coord(1) = T + dt / 2;
211  // Update the equation to the mid point position
212  p_EQUATION -> update( state );
213  // evaluate the Jacobian of mass contribution multiplied by state_dy
214  p_EQUATION -> get_jacobian_of_matrix0_mult_vector( state, state_dy, h0 );
215  // evaluate the Jacobian of mass contribution multiplied by state_dt
216  p_EQUATION -> get_jacobian_of_matrix1_mult_vector( state, state_dt, h1 );
217  // loop over all the variables
218  //
219  // to avoid repeated mapping arithmetic within operator() of the
220  // BandedMatrix class we'll access the matrix with the iterator.
221  //
222  // il = position for (0, lnode*order)
223  typename BandedMatrix<_Type>::elt_iter l_iter( a.get_elt_iter( 0, l_node * order ) );
224  // ir = position for (0, rnode*order)
225  typename BandedMatrix<_Type>::elt_iter r_iter( a.get_elt_iter( 0, r_node * order ) );
226  //
227  for ( unsigned var = 0; var < order; ++var )
228  {
229  // add the matrix mult terms to the linearised problem
230  for ( unsigned i = 0; i < order; ++i ) // dummy index
231  {
232  // add the Jacobian terms
233  // indirect access: a( row, order * l_node + i ) -= jac_midpt( var, i ) * 0.5;
234  const std::size_t offset( i * a.noffdiag() * 3 + row );
235  const typename BandedMatrix<_Type>::elt_iter left( l_iter + offset );
236  const typename BandedMatrix<_Type>::elt_iter right( r_iter + offset );
237  //
238  *left -= p_EQUATION -> jacobian()( var, i ) * 0.5;
239  // indirect access: a( row, order * r_node + i ) -= jac_midpt( var, i ) * 0.5;
240  *right -= p_EQUATION -> jacobian()( var, i ) * 0.5;
241  // add the Jacobian of mass terms for dt terms
242  // indirect access: a( row, order * l_node + i ) += h1( var, i ) * 0.5;
243  *left += h0( var, i ) * 0.5;
244  // indirect access: a( row, order * r_node + i ) += h1( var, i ) * 0.5;
245  *right += h0( var, i ) * 0.5;
246  // add the Jacobian of mass terms for dt terms
247  // indirect access: a( row, order * l_node + i ) += h2( var, i ) * 0.5;
248  *left += h1( var, i ) * 0.5;
249  // indirect access: a( row, order * r_node + i ) += h2( var, i ) * 0.5;
250  *right += h1( var, i ) * 0.5;
251  // add the mass matrix terms
252  // indirect access: a( row, order * l_node + i ) -= mass_midpt( var, i ) * inv_dy;
253  *left -= p_EQUATION -> matrix0()( var, i ) * inv_dy;
254  // indirect access: a( row, order * r_node + i ) += mass_midpt( var, i ) * inv_dy;
255  *right += p_EQUATION -> matrix0()( var, i ) * inv_dy;
256  // indirect access: a( row, order * l_node + i ) += mass_midpt( var, i ) * inv_dt;
257  *left += p_EQUATION -> matrix1()( var, i ) * inv_dt;
258  // indirect access: a( row, order * r_node + i ) += mass_midpt( var, i ) * inv_dt;
259  *right += p_EQUATION -> matrix1()( var, i ) * inv_dt;
260  }
261  // RHS
262  b[ row ] = p_EQUATION -> residual()[ var ];
263  b[ row ] -= Utility::dot( p_EQUATION -> matrix0()[ var ], state_dy );
264  b[ row ] -= Utility::dot( p_EQUATION -> matrix1()[ var ], state_dt );
265  b[ row ] *= 2;
266  // increment the row
267  row += 1;
268  }
269  }
270  // BCn equation is evaluated at the next step time point as
271  // they cannot depend on d/dt terms
272  p_RIGHT_RESIDUAL -> coord( 0 ) = T + dt;
273  // update the BC residuals for the current iteration
274  p_RIGHT_RESIDUAL -> update( SOLN.get_nodes_vars( ny - 1 ) );
275  // add the (linearised) RHS BCs to the matrix problem
276  for ( unsigned i = 0; i < p_RIGHT_RESIDUAL -> get_order(); ++i )
277  {
278  // loop thru variables at RHS of the domain
279  for ( unsigned var = 0; var < order; ++var )
280  {
281  a( row, order * ( ny - 1 ) + var ) = p_RIGHT_RESIDUAL -> jacobian()( i, var );
282  }
283  b[ row ] = - p_RIGHT_RESIDUAL -> residual()[ i ];
284  ++row;
285  }
286 #ifdef PARANOID
287  if ( row != ny * order )
288  {
289  std::string problem( "\n The ODE_BVP has an incorrect number of boundary conditions. \n" );
290  throw ExceptionRuntime( problem );
291  }
292 #endif
293  }
_Type dot(const DenseVector< _Type > &X, const DenseVector< _Type > &Y)
Templated dot product.
Definition: Utility.h:381
DenseVector< _Type > get_nodes_vars(const std::size_t &node) const
Get the variables stored at A SPECIFIED node.
std::size_t get_nnodes() const
DenseVector< _Type >::elt_iter elt_iter
Definition: BandedMatrix.h:22
double & coord()
Return a reference to the current value of the &#39;timelike/parabolic&#39; coordinate.
Definition: PDE_IBVP.h:67
const _Xtype & coord(const std::size_t &node) const
Access a nodal position.
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.

References CppNoddy::PDE_IBVP< _Type >::previous_solution(), and CppNoddy::PDE_IBVP< _Type >::solution().

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

68  {
69  return T;
70  }
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 87 of file PDE_IBVP.h.

88  {
89  return MAX_ITERATIONS;
90  }
template<typename _Type >
OneD_Node_Mesh< _Type > & CppNoddy::PDE_IBVP< _Type >::previous_solution ( )
inline
Returns
A handle to the previous step's solution mesh

Definition at line 126 of file PDE_IBVP.h.

Referenced by CppNoddy::PDE_IBVP< _Type >::coord().

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

Definition at line 120 of file PDE_IBVP.h.

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

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

A Crank-Nicolson 'time' stepper.

Definition at line 76 of file PDE_IBVP.cpp.

References CppNoddy::PDE_IBVP< _Type >::assemble_matrix_problem(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::get_nnodes(), CppNoddy::DenseVector< _Type >::inf_norm(), and CppNoddy::BandedLinearSystem< _Type >::solve().

Referenced by main().

77  {
78  // the order of the problem
79  unsigned order( p_EQUATION -> get_order() );
80  // get the number of nodes in the mesh
81  // -- this may have been refined by the user since the last call.
82  unsigned ny( SOLN.get_nnodes() );
83  // measure of maximum residual
84  double max_residual( 1.0 );
85  // iteration counter
86  int counter = 0;
87  // store this soln as the 'previous SOLUTION'
88  PREV_SOLN = SOLN;
89  // ANY LARGE STORAGE USED IN THE MAIN LOOP IS
90  // DEFINED HERE TO AVOID REPEATED CONSTRUCTION.
91  // Note we blank the A matrix after every iteration.
92  //
93  // Banded LHS matrix - max obove diagonal band width is
94  // from first variable at node i to last variable at node i+1
95  BandedMatrix<_Type> a( ny * order, 2 * order - 1, 0.0 );
96  // RHS
97  DenseVector<_Type> b( ny * order, 0.0 );
98  // linear solver definition
99 #ifdef LAPACK
100  BandedLinearSystem<_Type> system( &a, &b, "lapack" );
101 #else
102  BandedLinearSystem<_Type> system( &a, &b, "native" );
103 #endif
104  // loop until converged or too many iterations
105  do
106  {
107  // iteration counter
108  ++counter;
109 #ifdef TIME
110  T_ASSEMBLE.start();
111 #endif
112  assemble_matrix_problem( a, b, dt );
113  max_residual = b.inf_norm();
114 #ifdef DEBUG
115  std::cout << " PDE_IBVP.solve : Residual_max = " << max_residual << " tol = " << TOL << "\n";
116 #endif
117 #ifdef TIME
118  T_ASSEMBLE.stop();
119  T_SOLVE.start();
120 #endif
121  // linear solver
122  system.solve();
123  // keep the solution in a OneD_GenMesh object
124  for ( std::size_t var = 0; var < order; ++var )
125  {
126  for ( std::size_t i = 0; i < ny; ++i )
127  {
128  SOLN( i, var ) += b[ i * order + var ];
129  }
130  }
131 #ifdef TIME
132  T_SOLVE.stop();
133 #endif
134  }
135  while ( ( max_residual > TOL ) && ( counter < MAX_ITERATIONS ) );
136  if ( max_residual > TOL )
137  {
138  // restore to previous state because this step failed
139  SOLN = PREV_SOLN;
140  std::string problem( "\n The PDE_IBVP.step2 method took too many iterations.\n");
141  problem += "Solution has been restored to the previous accurate state.\n";
142  throw ExceptionItn( problem, counter, max_residual );
143  }
144 #ifdef DEBUG
145  std::cout << "[DEBUG] time-like variable = " << T << "\n";
146 #endif
147  // set the time to the updated level
148  T += dt;
149  }
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.
Definition: PDE_IBVP.cpp:153
std::size_t get_nnodes() const
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 80 of file PDE_IBVP.h.

81  {
82  return TOL;
83  }

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

© 2012

R.E. Hewitt