CppNoddy  0.90
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. More...
 
 ~OneD_TVDLF_Mesh ()
 Empty desctructor. More...
 
DenseVector< double > get_mid_node_vector () const
 Get the nodal positions in the middle of each element. More...
 
DenseVector< double > get_face_pos_vector () const
 Get the positions of the element faces. More...
 
void set_limiter (unsigned id)
 Set the limiter type to be applied in the slope values. More...
 
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. More...
 
void update_to (const double &CFL, const double &t_end)
 Update all the elements in this mesh to a USER SPECIFIED time step. More...
 
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. More...
 
DenseVector< double > integrate () const
 Integrate the concentration values across the entire mesh. More...
 
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. More...
 
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().

21  {
22  MESH_TIME = 0.0;
23 #ifdef DEBUG
24  std::cout << "DEBUG: Starting construction of a OneD_TVDLF_Mesh object. \n";
25 #endif
26  unsigned N = X.size();
27  if ( N <= 2 )
28  {
29  std::string problem;
30  problem = " The OneD_TVDLF_Mesh object is trying to construct itself \n";
31  problem += " with just one element! \n";
32  throw ExceptionRuntime( problem );
33  }
34 #ifdef DEBUG
35  std::cout << "\nDEBUG: configuration of the black mesh \n";
36 #endif
37  // set up the fn ptr to the initial conditions fn
38  p_Q_INIT = init_ptr;
39  // store the order of the conservative system here for simplicity
40  ORDER_OF_SYSTEM = ptr -> get_order();
41 
42  // set up the black elements
43  {
44  // first elt
45  BLACK_ELTS.push_back( OneD_TVDLF_Elt( X[0], X[1], ptr, true, -1 ) );
46  for ( std::size_t i = 1; i <= N - 3; ++i )
47  {
48  // interior elts
49  BLACK_ELTS.push_back( OneD_TVDLF_Elt( X[i], X[i+1], ptr ) );
50  }
51  // last elt
52  BLACK_ELTS.push_back( OneD_TVDLF_Elt( X[N-2], X[N-1], ptr, true, 1 ) );
53  }
54 #ifdef DEBUG
55  std::cout << "\nDEBUG: configuration of the red mesh \n";
56 #endif
57  // set up the red elements
58  {
59  // first elt
60  RED_ELTS.push_back( OneD_TVDLF_Elt( X[0], ( X[1] + X[0] ) / 2, ptr, true, -1 ) );
61  // interior elts
62  for ( std::size_t i = 1; i <= N - 2; ++i )
63  {
64  // interior elts
65  RED_ELTS.push_back( OneD_TVDLF_Elt( ( X[i-1] + X[i] ) / 2, ( X[i] + X[i+1] ) / 2, ptr ) );
66  }
67  // last elt
68  RED_ELTS.push_back( OneD_TVDLF_Elt( ( X[N-2] + X[N-1] ) / 2, X[N-1], ptr, true, 1 ) );
69  }
70 #ifdef DEBUG
71  std::cout << "\nDEBUG: linking the two meshes \n";
72 #endif
73  // the only tricky part for this scheme is the mesh interconnections
74  // black to red is easy enough
75  elt_iter er = RED_ELTS.begin();
76  elt_iter eb = BLACK_ELTS.begin();
77  DenseVector<double> s_left( 2, 0. );
78  s_left[ 0 ] = -1.;
79  s_left[ 1 ] = 0.;
80  DenseVector<double> s_right( 2, 0. );
81  s_right[ 0 ] = 0.;
82  s_right[ 1 ] = 1.;
83  DenseVector<double> s_whole( 2, 0. );
84  s_whole[ 0 ] = -1.;
85  s_whole[ 1 ] = 1.;
86  DenseVector<double> s_gen( 2, 0. );
87  // loop over red elts -- and define which black elts contribute
88  // this is straightforward, even for nonuniform meshes.
89  while ( er != RED_ELTS.end() )
90  {
91  if ( er -> get_external_flag() )
92  {
93  // if its an external elt
94  if ( er -> get_external_face_i() < 0 )
95  {
96  // and external face is left
97  er -> add_contribution( &( *eb ), s_left, 1 );
98  ++er;
99  }
100  else
101  {
102  // and external face is right
103  er -> add_contribution( &( *eb ), s_right, -1 );
104  ++er;
105  }
106  }
107  else
108  {
109  //internal elt
110  er -> add_contribution( &( *eb ), s_right, -1 );
111  ++eb;
112  er -> add_contribution( &( *eb ), s_left, 1 );
113  ++er;
114  }
115  }
116  // loop over all black elts and define which red elts contribute
117  // this is more involved if we allow for non-uniform meshes.
118  eb = BLACK_ELTS.begin();
119  er = RED_ELTS.begin();
120  while ( eb != BLACK_ELTS.end() )
121  {
122  if ( eb -> get_external_flag() )
123  {
124  // if its an external elt
125  if ( eb -> get_external_face_i() < 0 )
126  {
127  // and external face is left
128  eb -> add_contribution( &( *er ), s_whole, 0 );
129  ++er;
130  double s = er -> get_s( eb -> get_x( 1.0 ) );
131  s_gen[ 0 ] = -1.0;
132  s_gen[ 1 ] = s;
133  eb -> add_contribution( &( *er ), s_gen, 1 );
134  ++eb;
135  }
136  else
137  {
138  // and external face is right
139  double s = er -> get_s( eb -> get_x( -1.0 ) );
140  s_gen[ 0 ] = s;
141  s_gen[ 1 ] = 1.0;
142  eb -> add_contribution( &( *er ), s_gen, -1 );
143  ++er;
144  eb -> add_contribution( &( *er ), s_whole, 0 );
145  ++eb;
146  }
147  }
148  else
149  {
150  // internal elt
151  double s = er -> get_s( eb -> get_x( -1.0 ) );
152  s_gen[ 0 ] = s;
153  s_gen[ 1 ] = 1.0;
154  eb -> add_contribution( &( *er ), s_gen, -1 );
155  ++er;
156  s = er -> get_s( eb -> get_x( 1.0 ) );
157  s_gen[ 0 ] = -1.0;
158  s_gen[ 1 ] = s;
159  eb -> add_contribution( &( *er ), s_gen, 1 );
160  ++eb;
161  }
162  }
163 
164 #ifdef DEBUG
165  std::cout << "DEBUG: Setting the initial state of the meesh. \n";
166 #endif
167  // set the initial condition for each elt
168  eb = BLACK_ELTS.begin();
169  while ( eb != BLACK_ELTS.end() )
170  {
171  DenseVector<double> Q( ORDER_OF_SYSTEM, 0.0 );
172  p_Q_INIT( eb -> get_x( 0.0 ), Q );
173  eb -> set_Q_mid( Q );
174  ++eb;
175  }
176 
177  //eb = BLACK_ELTS.end();
178  //--eb;
179  //std::cout << "Last elt = " << eb -> get_Q( 0.0 )[ 1 ] << "\n";
180  //std::cout << "size = " << eb -> get_dx() << "\n";
181  //--eb;
182  //std::cout << "Last elt = " << eb -> get_Q( 0.0 )[ 1 ] << "\n";
183  //std::cout << "size = " << eb -> get_dx() << "\n";
184 
185  //DenseVector<double> x( get_mid_node_vector() );
186  //OneD_Mesh<double> soln( x );
187  //soln.set_nvars( ORDER_OF_SYSTEM );
188  //for ( std::size_t i = 0; i < x.size(); ++i )
189  //{
190  //soln.set_nodes_vars( i, BLACK_ELTS[i].get_Q( 0.0 ) );
191  //std::cout << "Get_soln Q = " << soln.get_nodes_vars( i )[ 1 ] << " at x = " << x[i] << "\n";
192  //}
193 
194  // default limiter = 0 (minmod)
195  LIMITER = 0;
196  calc_slopes( &BLACK_ELTS );
197 
198 #ifdef DEBUG
199  std::cout << "DEBUG: Finished construction of a OneD_TVDLF_Mesh object. \n";
200 #endif
201  }
double s
relative rotation rate
std::size_t size() const
A pass-thru definition to get the size of the vector.
Definition: DenseVector.h:365
CppNoddy::OneD_TVDLF_Mesh::~OneD_TVDLF_Mesh ( )

Empty desctructor.

Definition at line 203 of file OneD_TVDLF_Mesh.cpp.

204  {}

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

219  {
220  DenseVector<double> X;
221  celt_iter e = BLACK_ELTS.begin();
222  while ( e != BLACK_ELTS.end() )
223  {
224  X.push_back( e -> get_x( -1.0 ) );
225  ++e;
226  }
227  --e;
228  X.push_back( e -> get_x( 1.0 ) );
229  return X;
230  }
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().

207  {
208  DenseVector<double> X;
209  celt_iter e = BLACK_ELTS.begin();
210  while ( e != BLACK_ELTS.end() )
211  {
212  X.push_back( e -> get_x( 0.0 ) );
213  ++e;
214  }
215  return X;
216  }
OneD_Node_Mesh<double> CppNoddy::OneD_TVDLF_Mesh::get_slope ( )
inline

Definition at line 228 of file OneD_TVDLF_Mesh.h.

References CppNoddy::Example::A(), 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().

229  {
230  vector_of_elts* elts( get_elts_from_colour( "black" ) );
231  DenseVector<double> X( get_mid_node_vector() );
232  // the variables are the slope for each variable
233  OneD_Node_Mesh<double> slope_mesh( X, ORDER_OF_SYSTEM );
234  for ( std::size_t i = 0; i < X.size(); ++i )
235  {
236  slope_mesh.set_nodes_vars( i, ( *elts )[i].get_slope( ) );
237  }
238  return slope_mesh;
239  }
DenseVector< double > get_mid_node_vector() const
Get the nodal positions in the middle of each element.
OneD_Node_Mesh< double > get_slope()
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 CppNoddy::Example::A(), get_face_pos_vector(), get_mid_node_vector(), CppNoddy::DenseVector< _Type >::push_back(), CppNoddy::Example::s, CppNoddy::OneD_Node_Mesh< _Type, _Xtype >::set_nodes_vars(), and CppNoddy::DenseVector< _Type >::size().

Referenced by main(), and update_to_black().

271  {
272  vector_of_elts* elts( get_elts_from_colour( mesh_colour ) );
273  OneD_Node_Mesh<double> soln;
274  if ( location == "centre" )
275  {
276  DenseVector<double> X( get_mid_node_vector() );
277  soln = OneD_Node_Mesh<double>( X, ORDER_OF_SYSTEM );
278  for ( std::size_t i = 0; i < X.size(); ++i )
279  {
280  soln.set_nodes_vars( i, ( *elts )[i].get_Q( 0.0 ) );
281  }
282  }
283  else
284  if ( location == "face_average" )
285  {
286  DenseVector<double> X( get_face_pos_vector() );
287  soln = OneD_Node_Mesh<double>( X, ORDER_OF_SYSTEM );
288  std::size_t N = X.size() - 1;
289  soln.set_nodes_vars( 0, ( *elts )[0].get_Q( -1.0 ) );
290  for ( std::size_t i = 1; i < N; ++i )
291  {
292  soln.set_nodes_vars( i, ( ( *elts )[i-1].get_Q( 1.0 ) + ( *elts )[i].get_Q( -1.0 ) ) / 2 );
293  }
294  soln.set_nodes_vars( N, ( *elts )[N-1].get_Q( 1.0 ) );
295  }
296  else
297  {
298  std::string problem;
299  problem = " In OneD_TVDLF_Mesh::get_soln you have passed an unrecognised ";
300  problem += " location for data output. Use 'centre' or 'face_average'. \n";
301  throw ExceptionRuntime( problem );
302  }
303 
304  return soln;
305  }
DenseVector< double > get_face_pos_vector() const
Get the positions of the element faces.
DenseVector< double > get_mid_node_vector() const
Get the nodal positions in the middle of each element.
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.

Referenced by update_to_black().

254  {
255  return MESH_TIME;
256  }
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(), and update_to_black().

259  {
260  celt_iter e = BLACK_ELTS.begin();
261  DenseVector<double> sum( ORDER_OF_SYSTEM, 0.0 );
262  while ( e != BLACK_ELTS.end() )
263  {
264  sum += e -> get_Q( 0.0 ) * e -> get_dx();
265  ++e;
266  }
267  return sum;
268  }
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().

233  {
234  LIMITER = id;
235  }
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().

238  {
239  double first_dt = update_to_red( CFL, max_dt / 2.0 );
240  double second_dt = update_to_black( CFL, max_dt - first_dt );
241  return first_dt + second_dt;
242  }
double update_to_red(const double &CFL, const double &max_dt)
double update_to_black(const double &CFL, const double &max_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().

245  {
246  do
247  {
248  update( CFL, std::abs( t_end - MESH_TIME ) );
249  }
250  while ( MESH_TIME < t_end );
251  }
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.
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(), get_soln(), get_time(), integrate(), and CppNoddy::Example::source_fn().

Referenced by update().

141  {
142  // integrate the red mesh data back onto the black mesh
143  {
144  elt_iter eb = BLACK_ELTS.begin();
145  while ( eb != BLACK_ELTS.end() )
146  {
147  eb -> set_Q_mid( eb -> contributed_Q() / eb -> get_dx() );
148  ++eb;
149  }
150  }
151  // step thru the red elements to find a max time step
152  double second_dt;
153  {
154  elt_iter er = RED_ELTS.begin();
155  second_dt = er -> get_max_dt();
156  ++er;
157  while ( er != RED_ELTS.end() )
158  {
159  second_dt = std::min( second_dt, er -> get_max_dt() );
160  ++er;
161  }
162  second_dt *= CFL;
163  }
164  if ( second_dt > max_dt )
165  {
166  second_dt = max_dt;
167  }
168 
169  calc_slopes( &BLACK_ELTS );
170  // compute the fluxes through the element boundaries
171  // in the black mesh, using the nodal values of the red mesh
172  // then update the concentrations in the black mesh elements
173  {
174  elt_iter eb = BLACK_ELTS.begin();
175  DenseVector<double> flux_in_left( ORDER_OF_SYSTEM, 0.0 );
176  DenseVector<double> flux_out_right( ORDER_OF_SYSTEM, 0.0 );
177  // start with the left most elt & work out the flux in from the left
178  eb -> contributed_flux_in_left( second_dt, flux_in_left, MESH_TIME );
179  while ( eb != BLACK_ELTS.end() )
180  {
181  // work out the flux out to the right
182  eb -> contributed_flux_out_right( second_dt, flux_out_right, MESH_TIME );
183  // we now have the flux difference
184  DenseVector<double> deltaQ = ( flux_in_left - flux_out_right ) * second_dt / eb -> get_dx();
185  // contribution from the source integral
186  {
187  double x_mid( eb -> get_x( 0.0 ) );
188  DenseVector<double> q_mid( eb -> get_Q( 0.0 ) );
189  DenseVector<double> slope( eb -> get_slope() );
190  q_mid += ( eb -> get_source_fn( 0.0 )
191  - eb -> get_Jac_flux_fn( 0.0 ).multiply( slope ) ) * 0.5 * second_dt;
192  DenseVector<double> r_mid( ORDER_OF_SYSTEM, 0.0 );
193  eb -> system_ptr -> source_fn( x_mid, q_mid, slope, r_mid );
194  deltaQ += r_mid * second_dt;
195  }
196  eb -> set_Q_mid( eb -> get_Q( 0.0 ) + deltaQ );
197  // the flux out right in this elt is the flux in left of the next one
198  flux_in_left = flux_out_right;
199  ++eb;
200  }
201  }
202 
203  // compute the slopes using the specified limiter
204  // for the black elements
205  calc_slopes( &BLACK_ELTS );
206  MESH_TIME += second_dt;
207  return second_dt;
208 
209  }
double source_fn(double &x, double &y)
Define the source function for the Poisson solver.
Definition: Poisson_C.cpp:24
OneD_Node_Mesh< double > get_slope()
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().

68  {
69  // the black mesh slopes are set in the constructor
70  // and at the end of every 'update' call
71 
72  // integrate the black mesh data onto the red mesh
73  {
74  elt_iter er = RED_ELTS.begin();
75  while ( er != RED_ELTS.end() )
76  {
77  er -> set_Q_mid( er -> contributed_Q() / er -> get_dx() );
78  ++er;
79  }
80  }
81 
82  // step thru the black elements to find a max time step
83  double first_dt;
84  {
85  elt_iter eb = BLACK_ELTS.begin();
86  first_dt = eb -> get_max_dt();
87  ++eb;
88  while ( eb != BLACK_ELTS.end() )
89  {
90  first_dt = std::min( first_dt, eb -> get_max_dt() );
91  ++eb;
92  }
93  first_dt *= CFL;
94  }
95  if ( first_dt > max_dt )
96  {
97  first_dt = max_dt;
98  }
99 
100  calc_slopes( &RED_ELTS );
101  // compute the fluxes through the element boundaries
102  // in the red mesh, using the nodal values of the black mesh
103  // then update the concentrations in the red mesh elements
104  {
105  elt_iter er = RED_ELTS.begin();
106  DenseVector<double> flux_in_left( ORDER_OF_SYSTEM, 0.0 );
107  DenseVector<double> flux_out_right( ORDER_OF_SYSTEM, 0.0 );
108  // start with the left most elt & work out the flux in from the left
109  er -> contributed_flux_in_left( first_dt, flux_in_left, MESH_TIME );
110  while ( er != RED_ELTS.end() )
111  {
112  // work out the flux out to the right
113  er -> contributed_flux_out_right( first_dt, flux_out_right, MESH_TIME );
114  // we now have the flux difference
115  DenseVector<double> deltaQ = ( flux_in_left - flux_out_right ) * first_dt / er -> get_dx();
116  // contribution from the source integral
117  {
118  double x_mid( er -> get_x( 0.0 ) );
119  DenseVector<double> slope( er -> get_slope() );
120  DenseVector<double> q_mid( er -> get_Q( 0.0 ) );
121  q_mid += ( er -> get_source_fn( 0.0 )
122  - er -> get_Jac_flux_fn( 0.0 ).multiply( slope ) ) * 0.5 * first_dt;
123  DenseVector<double> r_mid( ORDER_OF_SYSTEM, 0.0 );
124  er -> system_ptr -> source_fn( x_mid, q_mid, slope, r_mid );
125  deltaQ += r_mid * first_dt;
126  }
127  er -> set_Q_mid( er -> get_Q( 0.0 ) + deltaQ );
128  // the flux out right in this elt is the flux in left of the next one
129  flux_in_left = flux_out_right;
130  ++er;
131  }
132  }
133 
134  // compute the slopes using the specified limiter for the red elements
135  calc_slopes( &RED_ELTS );
136  MESH_TIME += first_dt;
137  return first_dt;
138  }
double source_fn(double &x, double &y)
Define the source function for the Poisson solver.
Definition: Poisson_C.cpp:24
OneD_Node_Mesh< double > get_slope()

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

© 2012

R.E. Hewitt