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

#include <OneD_TVDLF_Mesh.h>

Public Member Functions

 OneD_TVDLF_Mesh (const DenseVector< double > &X, OneD_Hyperbolic_System *ptr, fn_ptr init_ptr)
 Constructor for the Finite Volume Mesh using linear elements.
 ~OneD_TVDLF_Mesh ()
 Empty desctructor.
DenseVector< double > get_mid_node_vector () const
 Get the nodal positions in the middle of each element.
DenseVector< double > get_face_pos_vector () const
 Get the positions of the element faces.
void set_limiter (unsigned id)
 Set the limiter type to be applied in the slope values.
double update (const double &CFL, const double &max_dt=std::numeric_limits< long double >::max())
 Update all the elements in this mesh to a new time step.
void update_to (const double &CFL, const double &t_end)
 Update all the elements in this mesh to a USER SPECIFIED time step.
double update_to_red (const double &CFL, const double &max_dt)
double update_to_black (const double &CFL, const double &max_dt)
const double & get_time () const
 Get a const reference to the time value for the current mesh.
DenseVector< double > integrate () const
 Integrate the concentration values across the entire mesh.
OneD_Node_Mesh< double > get_soln (std::string location="centre", std::string colour="black")
 Get a OneD_Mesh<double> object containing the one dimensional data in the usual format.
OneD_Node_Mesh< double > get_slope ()

Detailed Description

Definition at line 20 of file OneD_TVDLF_Mesh.h.

Constructor & Destructor Documentation

CppNoddy::OneD_TVDLF_Mesh::OneD_TVDLF_Mesh ( const DenseVector< double > &  X,
OneD_Hyperbolic_System ptr,
fn_ptr  init_ptr 
)

Constructor for the Finite Volume Mesh using linear elements.

Parameters
XA vector of nodal locations at which the element FACES will positioned
ptrA pointer to the hyperbolic system applied to this mesh
init_ptrA pointer to a function that defines the initial conditions
XA vector of nodal locations at which the element FACES will positioned

Definition at line 18 of file OneD_TVDLF_Mesh.cpp.

References CppNoddy::Example::s, and CppNoddy::DenseVector< _Type >::size().

{
MESH_TIME = 0.0;
#ifdef DEBUG
std::cout << "DEBUG: Starting construction of a OneD_TVDLF_Mesh object. \n";
#endif
unsigned N = X.size();
if ( N <= 2 )
{
std::string problem;
problem = " The OneD_TVDLF_Mesh object is trying to construct itself \n";
problem += " with just one element! \n";
throw ExceptionRuntime( problem );
}
#ifdef DEBUG
std::cout << "\nDEBUG: configuration of the black mesh \n";
#endif
// set up the fn ptr to the initial conditions fn
p_Q_INIT = init_ptr;
// store the order of the conservative system here for simplicity
ORDER_OF_SYSTEM = ptr -> get_order();
// set up the black elements
{
// first elt
BLACK_ELTS.push_back( OneD_TVDLF_Elt( X[0], X[1], ptr, true, -1 ) );
for ( std::size_t i = 1; i <= N - 3; ++i )
{
// interior elts
BLACK_ELTS.push_back( OneD_TVDLF_Elt( X[i], X[i+1], ptr ) );
}
// last elt
BLACK_ELTS.push_back( OneD_TVDLF_Elt( X[N-2], X[N-1], ptr, true, 1 ) );
}
#ifdef DEBUG
std::cout << "\nDEBUG: configuration of the red mesh \n";
#endif
// set up the red elements
{
// first elt
RED_ELTS.push_back( OneD_TVDLF_Elt( X[0], ( X[1] + X[0] ) / 2, ptr, true, -1 ) );
// interior elts
for ( std::size_t i = 1; i <= N - 2; ++i )
{
// interior elts
RED_ELTS.push_back( OneD_TVDLF_Elt( ( X[i-1] + X[i] ) / 2, ( X[i] + X[i+1] ) / 2, ptr ) );
}
// last elt
RED_ELTS.push_back( OneD_TVDLF_Elt( ( X[N-2] + X[N-1] ) / 2, X[N-1], ptr, true, 1 ) );
}
#ifdef DEBUG
std::cout << "\nDEBUG: linking the two meshes \n";
#endif
// the only tricky part for this scheme is the mesh interconnections
// black to red is easy enough
elt_iter er = RED_ELTS.begin();
elt_iter eb = BLACK_ELTS.begin();
DenseVector<double> s_left( 2, 0. );
s_left[ 0 ] = -1.;
s_left[ 1 ] = 0.;
DenseVector<double> s_right( 2, 0. );
s_right[ 0 ] = 0.;
s_right[ 1 ] = 1.;
DenseVector<double> s_whole( 2, 0. );
s_whole[ 0 ] = -1.;
s_whole[ 1 ] = 1.;
DenseVector<double> s_gen( 2, 0. );
// loop over red elts -- and define which black elts contribute
// this is straightforward, even for nonuniform meshes.
while ( er != RED_ELTS.end() )
{
if ( er -> get_external_flag() )
{
// if its an external elt
if ( er -> get_external_face_i() < 0 )
{
// and external face is left
er -> add_contribution( &( *eb ), s_left, 1 );
++er;
}
else
{
// and external face is right
er -> add_contribution( &( *eb ), s_right, -1 );
++er;
}
}
else
{
//internal elt
er -> add_contribution( &( *eb ), s_right, -1 );
++eb;
er -> add_contribution( &( *eb ), s_left, 1 );
++er;
}
}
// loop over all black elts and define which red elts contribute
// this is more involved if we allow for non-uniform meshes.
eb = BLACK_ELTS.begin();
er = RED_ELTS.begin();
while ( eb != BLACK_ELTS.end() )
{
if ( eb -> get_external_flag() )
{
// if its an external elt
if ( eb -> get_external_face_i() < 0 )
{
// and external face is left
eb -> add_contribution( &( *er ), s_whole, 0 );
++er;
double s = er -> get_s( eb -> get_x( 1.0 ) );
s_gen[ 0 ] = -1.0;
s_gen[ 1 ] = s;
eb -> add_contribution( &( *er ), s_gen, 1 );
++eb;
}
else
{
// and external face is right
double s = er -> get_s( eb -> get_x( -1.0 ) );
s_gen[ 0 ] = s;
s_gen[ 1 ] = 1.0;
eb -> add_contribution( &( *er ), s_gen, -1 );
++er;
eb -> add_contribution( &( *er ), s_whole, 0 );
++eb;
}
}
else
{
// internal elt
double s = er -> get_s( eb -> get_x( -1.0 ) );
s_gen[ 0 ] = s;
s_gen[ 1 ] = 1.0;
eb -> add_contribution( &( *er ), s_gen, -1 );
++er;
s = er -> get_s( eb -> get_x( 1.0 ) );
s_gen[ 0 ] = -1.0;
s_gen[ 1 ] = s;
eb -> add_contribution( &( *er ), s_gen, 1 );
++eb;
}
}
#ifdef DEBUG
std::cout << "DEBUG: Setting the initial state of the meesh. \n";
#endif
// set the initial condition for each elt
eb = BLACK_ELTS.begin();
while ( eb != BLACK_ELTS.end() )
{
DenseVector<double> Q( ORDER_OF_SYSTEM, 0.0 );
p_Q_INIT( eb -> get_x( 0.0 ), Q );
eb -> set_Q_mid( Q );
++eb;
}
//eb = BLACK_ELTS.end();
//--eb;
//std::cout << "Last elt = " << eb -> get_Q( 0.0 )[ 1 ] << "\n";
//std::cout << "size = " << eb -> get_dx() << "\n";
//--eb;
//std::cout << "Last elt = " << eb -> get_Q( 0.0 )[ 1 ] << "\n";
//std::cout << "size = " << eb -> get_dx() << "\n";
//DenseVector<double> x( get_mid_node_vector() );
//OneD_Mesh<double> soln( x );
//soln.set_nvars( ORDER_OF_SYSTEM );
//for ( std::size_t i = 0; i < x.size(); ++i )
//{
//soln.set_nodes_vars( i, BLACK_ELTS[i].get_Q( 0.0 ) );
//std::cout << "Get_soln Q = " << soln.get_nodes_vars( i )[ 1 ] << " at x = " << x[i] << "\n";
//}
// default limiter = 0 (minmod)
LIMITER = 0;
calc_slopes( &BLACK_ELTS );
#ifdef DEBUG
std::cout << "DEBUG: Finished construction of a OneD_TVDLF_Mesh object. \n";
#endif
}
CppNoddy::OneD_TVDLF_Mesh::~OneD_TVDLF_Mesh ( )

Empty desctructor.

Definition at line 203 of file OneD_TVDLF_Mesh.cpp.

{}

Member Function Documentation

DenseVector< double > CppNoddy::OneD_TVDLF_Mesh::get_face_pos_vector ( ) const

Get the positions of the element faces.

Returns
An NVector<double> of the spatial points

Definition at line 218 of file OneD_TVDLF_Mesh.cpp.

References CppNoddy::DenseVector< _Type >::push_back().

Referenced by get_soln().

{
DenseVector<double> X;
celt_iter e = BLACK_ELTS.begin();
while ( e != BLACK_ELTS.end() )
{
X.push_back( e -> get_x( -1.0 ) );
++e;
}
--e;
X.push_back( e -> get_x( 1.0 ) );
return X;
}
DenseVector< double > CppNoddy::OneD_TVDLF_Mesh::get_mid_node_vector ( ) const

Get the nodal positions in the middle of each element.

Returns
An NVector<double> of the nodal points

Definition at line 206 of file OneD_TVDLF_Mesh.cpp.

References CppNoddy::DenseVector< _Type >::push_back().

Referenced by get_slope(), and get_soln().

{
DenseVector<double> X;
celt_iter e = BLACK_ELTS.begin();
while ( e != BLACK_ELTS.end() )
{
X.push_back( e -> get_x( 0.0 ) );
++e;
}
return X;
}
OneD_Node_Mesh<double> CppNoddy::OneD_TVDLF_Mesh::get_slope ( )
inline

Definition at line 228 of file OneD_TVDLF_Mesh.h.

References get_mid_node_vector(), CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::set_nodes_vars(), and CppNoddy::DenseVector< _Type >::size().

Referenced by update_to_black(), and update_to_red().

{
vector_of_elts* elts( get_elts_from_colour( "black" ) );
DenseVector<double> X( get_mid_node_vector() );
// the variables are the slope for each variable
OneD_Node_Mesh<double> slope_mesh( X, ORDER_OF_SYSTEM );
for ( std::size_t i = 0; i < X.size(); ++i )
{
slope_mesh.set_nodes_vars( i, ( *elts )[i].get_slope( ) );
}
return slope_mesh;
}
OneD_Node_Mesh< double > CppNoddy::OneD_TVDLF_Mesh::get_soln ( std::string  location = "centre",
std::string  colour = "black" 
)

Get a OneD_Mesh<double> object containing the one dimensional data in the usual format.

Parameters
locationUse "centre" for mid-elt values and "face_average" for the average values at the (discontinuous) element boundaries
colourWhich mesh to output, unless debugging, this should be "black" otherwise time values will be slightly out
Returns
The required mesh object

Definition at line 270 of file OneD_TVDLF_Mesh.cpp.

References get_face_pos_vector(), get_mid_node_vector(), and CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::set_nodes_vars().

Referenced by main().

{
vector_of_elts* elts( get_elts_from_colour( mesh_colour ) );
OneD_Node_Mesh<double> soln;
if ( location == "centre" )
{
DenseVector<double> X( get_mid_node_vector() );
soln = OneD_Node_Mesh<double>( X, ORDER_OF_SYSTEM );
for ( std::size_t i = 0; i < X.size(); ++i )
{
soln.set_nodes_vars( i, ( *elts )[i].get_Q( 0.0 ) );
}
}
else
if ( location == "face_average" )
{
DenseVector<double> X( get_face_pos_vector() );
soln = OneD_Node_Mesh<double>( X, ORDER_OF_SYSTEM );
std::size_t N = X.size() - 1;
soln.set_nodes_vars( 0, ( *elts )[0].get_Q( -1.0 ) );
for ( std::size_t i = 1; i < N; ++i )
{
soln.set_nodes_vars( i, ( ( *elts )[i-1].get_Q( 1.0 ) + ( *elts )[i].get_Q( -1.0 ) ) / 2 );
}
soln.set_nodes_vars( N, ( *elts )[N-1].get_Q( 1.0 ) );
}
else
{
std::string problem;
problem = " In OneD_TVDLF_Mesh::get_soln you have passed an unrecognised ";
problem += " location for data output. Use 'centre' or 'face_average'. \n";
throw ExceptionRuntime( problem );
}
return soln;
}
const double & CppNoddy::OneD_TVDLF_Mesh::get_time ( ) const

Get a const reference to the time value for the current mesh.

Returns
time The time level at which the data in the mesh applies.

Definition at line 253 of file OneD_TVDLF_Mesh.cpp.

{
return MESH_TIME;
}
DenseVector< double > CppNoddy::OneD_TVDLF_Mesh::integrate ( ) const

Integrate the concentration values across the entire mesh.

Returns
The values of the integral of each component.

Definition at line 258 of file OneD_TVDLF_Mesh.cpp.

Referenced by main().

{
celt_iter e = BLACK_ELTS.begin();
DenseVector<double> sum( ORDER_OF_SYSTEM, 0.0 );
while ( e != BLACK_ELTS.end() )
{
sum += e -> get_Q( 0.0 ) * e -> get_dx();
++e;
}
return sum;
}
void CppNoddy::OneD_TVDLF_Mesh::set_limiter ( unsigned  id)

Set the limiter type to be applied in the slope values.

0 is no limiter, 1 is Lax-Wendroff, 2 is Beam-Warming, 3 is MC and 4 is Superbee.

Parameters
idThe identifier of the limiter.

Definition at line 232 of file OneD_TVDLF_Mesh.cpp.

Referenced by main().

{
LIMITER = id;
}
double CppNoddy::OneD_TVDLF_Mesh::update ( const double &  CFL,
const double &  max_dt = std::numeric_limits<long double>::max() 
)

Update all the elements in this mesh to a new time step.

Parameters
CFLThe CFL value to be used to determine the time step
max_dtDo not take a time step larger than this irrespective of the CFL value

Definition at line 237 of file OneD_TVDLF_Mesh.cpp.

References update_to_black(), and update_to_red().

Referenced by main(), and update_to().

{
double first_dt = update_to_red( CFL, max_dt / 2.0 );
double second_dt = update_to_black( CFL, max_dt - first_dt );
return first_dt + second_dt;
}
void CppNoddy::OneD_TVDLF_Mesh::update_to ( const double &  CFL,
const double &  t_end 
)

Update all the elements in this mesh to a USER SPECIFIED time step.

Parameters
CFLThe CFL value to be used to determine the time step
t_endThe time level to compute to

Definition at line 244 of file OneD_TVDLF_Mesh.cpp.

References update().

{
do
{
update( CFL, std::abs( t_end - MESH_TIME ) );
}
while ( MESH_TIME < t_end );
}
double CppNoddy::OneD_TVDLF_Mesh::update_to_black ( const double &  CFL,
const double &  max_dt 
)
inline

Definition at line 140 of file OneD_TVDLF_Mesh.h.

References get_slope(), and CppNoddy::Example::source_fn().

Referenced by update().

{
// integrate the red mesh data back onto the black mesh
{
elt_iter eb = BLACK_ELTS.begin();
while ( eb != BLACK_ELTS.end() )
{
eb -> set_Q_mid( eb -> contributed_Q() / eb -> get_dx() );
++eb;
}
}
// step thru the red elements to find a max time step
double second_dt;
{
elt_iter er = RED_ELTS.begin();
second_dt = er -> get_max_dt();
++er;
while ( er != RED_ELTS.end() )
{
second_dt = std::min( second_dt, er -> get_max_dt() );
++er;
}
second_dt *= CFL;
}
if ( second_dt > max_dt )
{
second_dt = max_dt;
}
calc_slopes( &BLACK_ELTS );
// compute the fluxes through the element boundaries
// in the black mesh, using the nodal values of the red mesh
// then update the concentrations in the black mesh elements
{
elt_iter eb = BLACK_ELTS.begin();
DenseVector<double> flux_in_left( ORDER_OF_SYSTEM, 0.0 );
DenseVector<double> flux_out_right( ORDER_OF_SYSTEM, 0.0 );
// start with the left most elt & work out the flux in from the left
eb -> contributed_flux_in_left( second_dt, flux_in_left, MESH_TIME );
while ( eb != BLACK_ELTS.end() )
{
// work out the flux out to the right
eb -> contributed_flux_out_right( second_dt, flux_out_right, MESH_TIME );
// we now have the flux difference
DenseVector<double> deltaQ = ( flux_in_left - flux_out_right ) * second_dt / eb -> get_dx();
// contribution from the source integral
{
double x_mid( eb -> get_x( 0.0 ) );
DenseVector<double> q_mid( eb -> get_Q( 0.0 ) );
DenseVector<double> slope( eb -> get_slope() );
q_mid += ( eb -> get_source_fn( 0.0 )
- eb -> get_Jac_flux_fn( 0.0 ).multiply( slope ) ) * 0.5 * second_dt;
DenseVector<double> r_mid( ORDER_OF_SYSTEM, 0.0 );
eb -> system_ptr -> source_fn( x_mid, q_mid, slope, r_mid );
deltaQ += r_mid * second_dt;
}
eb -> set_Q_mid( eb -> get_Q( 0.0 ) + deltaQ );
// the flux out right in this elt is the flux in left of the next one
flux_in_left = flux_out_right;
++eb;
}
}
// compute the slopes using the specified limiter
// for the black elements
calc_slopes( &BLACK_ELTS );
MESH_TIME += second_dt;
return second_dt;
}
double CppNoddy::OneD_TVDLF_Mesh::update_to_red ( const double &  CFL,
const double &  max_dt 
)
inline

Definition at line 67 of file OneD_TVDLF_Mesh.h.

References get_slope(), and CppNoddy::Example::source_fn().

Referenced by update().

{
// the black mesh slopes are set in the constructor
// and at the end of every 'update' call
// integrate the black mesh data onto the red mesh
{
elt_iter er = RED_ELTS.begin();
while ( er != RED_ELTS.end() )
{
er -> set_Q_mid( er -> contributed_Q() / er -> get_dx() );
++er;
}
}
// step thru the black elements to find a max time step
double first_dt;
{
elt_iter eb = BLACK_ELTS.begin();
first_dt = eb -> get_max_dt();
++eb;
while ( eb != BLACK_ELTS.end() )
{
first_dt = std::min( first_dt, eb -> get_max_dt() );
++eb;
}
first_dt *= CFL;
}
if ( first_dt > max_dt )
{
first_dt = max_dt;
}
calc_slopes( &RED_ELTS );
// compute the fluxes through the element boundaries
// in the red mesh, using the nodal values of the black mesh
// then update the concentrations in the red mesh elements
{
elt_iter er = RED_ELTS.begin();
DenseVector<double> flux_in_left( ORDER_OF_SYSTEM, 0.0 );
DenseVector<double> flux_out_right( ORDER_OF_SYSTEM, 0.0 );
// start with the left most elt & work out the flux in from the left
er -> contributed_flux_in_left( first_dt, flux_in_left, MESH_TIME );
while ( er != RED_ELTS.end() )
{
// work out the flux out to the right
er -> contributed_flux_out_right( first_dt, flux_out_right, MESH_TIME );
// we now have the flux difference
DenseVector<double> deltaQ = ( flux_in_left - flux_out_right ) * first_dt / er -> get_dx();
// contribution from the source integral
{
double x_mid( er -> get_x( 0.0 ) );
DenseVector<double> slope( er -> get_slope() );
DenseVector<double> q_mid( er -> get_Q( 0.0 ) );
q_mid += ( er -> get_source_fn( 0.0 )
- er -> get_Jac_flux_fn( 0.0 ).multiply( slope ) ) * 0.5 * first_dt;
DenseVector<double> r_mid( ORDER_OF_SYSTEM, 0.0 );
er -> system_ptr -> source_fn( x_mid, q_mid, slope, r_mid );
deltaQ += r_mid * first_dt;
}
er -> set_Q_mid( er -> get_Q( 0.0 ) + deltaQ );
// the flux out right in this elt is the flux in left of the next one
flux_in_left = flux_out_right;
++er;
}
}
// compute the slopes using the specified limiter for the red elements
calc_slopes( &RED_ELTS );
MESH_TIME += first_dt;
return first_dt;
}

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

© 2012

R.E. Hewitt