CppNoddy  0.92
Loading...
Searching...
No Matches
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 19 of file OneD_TVDLF_Mesh.h.

Constructor & Destructor Documentation

◆ OneD_TVDLF_Mesh()

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 17 of file OneD_TVDLF_Mesh.cpp.

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

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

◆ ~OneD_TVDLF_Mesh()

CppNoddy::OneD_TVDLF_Mesh::~OneD_TVDLF_Mesh ( )

Empty desctructor.

Definition at line 183 of file OneD_TVDLF_Mesh.cpp.

184 {}

Member Function Documentation

◆ get_face_pos_vector()

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 196 of file OneD_TVDLF_Mesh.cpp.

196 {
197 DenseVector<double> X;
198 celt_iter e = BLACK_ELTS.begin();
199 while(e != BLACK_ELTS.end()) {
200 X.push_back(e -> get_x(-1.0));
201 ++e;
202 }
203 --e;
204 X.push_back(e -> get_x(1.0));
205 return X;
206 }

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

Referenced by get_soln().

◆ get_mid_node_vector()

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 186 of file OneD_TVDLF_Mesh.cpp.

186 {
187 DenseVector<double> X;
188 celt_iter e = BLACK_ELTS.begin();
189 while(e != BLACK_ELTS.end()) {
190 X.push_back(e -> get_x(0.0));
191 ++e;
192 }
193 return X;
194 }

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

Referenced by get_slope(), and get_soln().

◆ get_slope()

OneD_Node_Mesh< double > CppNoddy::OneD_TVDLF_Mesh::get_slope ( )
inline

Definition at line 216 of file OneD_TVDLF_Mesh.h.

216 {
217 vector_of_elts* elts(get_elts_from_colour("black"));
218 DenseVector<double> X(get_mid_node_vector());
219 // the variables are the slope for each variable
220 OneD_Node_Mesh<double> slope_mesh(X, ORDER_OF_SYSTEM);
221 for(std::size_t i = 0; i < X.size(); ++i) {
222 slope_mesh.set_nodes_vars(i, (*elts)[i].get_slope());
223 }
224 return slope_mesh;
225 }
CppNoddy::OneD_TVDLF_Mesh::get_mid_node_vector
DenseVector< double > get_mid_node_vector() const
Get the nodal positions in the middle of each element.
Definition: OneD_TVDLF_Mesh.cpp:186
CppNoddy::OneD_TVDLF_Mesh::get_slope
OneD_Node_Mesh< double > get_slope()
Definition: OneD_TVDLF_Mesh.h:216

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

Referenced by get_slope(), update_to_black(), and update_to_red().

◆ get_soln()

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 238 of file OneD_TVDLF_Mesh.cpp.

238 {
239 vector_of_elts* elts(get_elts_from_colour(mesh_colour));
240 OneD_Node_Mesh<double> soln;
241 if(location == "centre") {
242 DenseVector<double> X(get_mid_node_vector());
243 soln = OneD_Node_Mesh<double>(X, ORDER_OF_SYSTEM);
244 for(std::size_t i = 0; i < X.size(); ++i) {
245 soln.set_nodes_vars(i, (*elts)[i].get_Q(0.0));
246 }
247 } else if(location == "face_average") {
248 DenseVector<double> X(get_face_pos_vector());
249 soln = OneD_Node_Mesh<double>(X, ORDER_OF_SYSTEM);
250 std::size_t N = X.size() - 1;
251 soln.set_nodes_vars(0, (*elts)[0].get_Q(-1.0));
252 for(std::size_t i = 1; i < N; ++i) {
253 soln.set_nodes_vars(i, ((*elts)[i-1].get_Q(1.0) + (*elts)[i].get_Q(-1.0)) / 2);
254 }
255 soln.set_nodes_vars(N, (*elts)[N-1].get_Q(1.0));
256 } else {
257 std::string problem;
258 problem = " In OneD_TVDLF_Mesh::get_soln you have passed an unrecognised ";
259 problem += " location for data output. Use 'centre' or 'face_average'. \n";
260 throw ExceptionRuntime(problem);
261 }
262
263 return soln;
264 }
CppNoddy::OneD_TVDLF_Mesh::get_face_pos_vector
DenseVector< double > get_face_pos_vector() const
Get the positions of the element faces.
Definition: OneD_TVDLF_Mesh.cpp:196

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

Referenced by main().

◆ get_time()

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 224 of file OneD_TVDLF_Mesh.cpp.

224 {
225 return MESH_TIME;
226 }

◆ integrate()

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 228 of file OneD_TVDLF_Mesh.cpp.

228 {
229 celt_iter e = BLACK_ELTS.begin();
230 DenseVector<double> sum(ORDER_OF_SYSTEM, 0.0);
231 while(e != BLACK_ELTS.end()) {
232 sum += e -> get_Q(0.0) * e -> get_dx();
233 ++e;
234 }
235 return sum;
236 }

Referenced by main().

◆ set_limiter()

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 208 of file OneD_TVDLF_Mesh.cpp.

208 {
209 LIMITER = id;
210 }

Referenced by main().

◆ update()

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 212 of file OneD_TVDLF_Mesh.cpp.

212 {
213 double first_dt = update_to_red(CFL, max_dt / 2.0);
214 double second_dt = update_to_black(CFL, max_dt - first_dt);
215 return first_dt + second_dt;
216 }
CppNoddy::OneD_TVDLF_Mesh::update_to_black
double update_to_black(const double &CFL, const double &max_dt)
Definition: OneD_TVDLF_Mesh.h:133
CppNoddy::OneD_TVDLF_Mesh::update_to_red
double update_to_red(const double &CFL, const double &max_dt)
Definition: OneD_TVDLF_Mesh.h:65

References update_to_black(), and update_to_red().

Referenced by main(), and update_to().

◆ update_to()

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 218 of file OneD_TVDLF_Mesh.cpp.

218 {
219 do {
220 update(CFL, std::abs(t_end - MESH_TIME));
221 } while(MESH_TIME < t_end);
222 }
CppNoddy::OneD_TVDLF_Mesh::update
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.
Definition: OneD_TVDLF_Mesh.cpp:212

References update().

◆ update_to_black()

double CppNoddy::OneD_TVDLF_Mesh::update_to_black ( const double &  CFL,
const double &  max_dt 
)
inline

Definition at line 133 of file OneD_TVDLF_Mesh.h.

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

References get_slope().

Referenced by update().

◆ update_to_red()

double CppNoddy::OneD_TVDLF_Mesh::update_to_red ( const double &  CFL,
const double &  max_dt 
)
inline

Definition at line 65 of file OneD_TVDLF_Mesh.h.

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

References get_slope().

Referenced by update().

The documentation for this class was generated from the following files:
Generated on Tue Jan 1 1980 00:00:00 for CppNoddy by   doxygen 1.9.5

© 2012

R.E. Hewitt