Carmen Code
 All Classes Files Functions Variables Macros Pages
Functions
GetBoundaryCells.cpp File Reference

Computes the cells C1, C2, C3, C4 in function of the cells Cell1, Cell2, Cell3, Cell4 taking into account boundary conditions. More...

#include "Carmen.h"
Include dependency graph for GetBoundaryCells.cpp:

Functions

void GetBoundaryCells (Cell &Cell1, Cell &Cell2, Cell &Cell3, Cell &Cell4, Cell &C1, Cell &C2, Cell &C3, Cell &C4, const int AxisNo)
 Transform the 4 cells of the flux Cell1, Cell2, Cell3, Cell4 into C1, C2, C3, C4, to take into account boundary conditions (used in Flux.cpp). More...
 

Detailed Description

Computes the cells C1, C2, C3, C4 in function of the cells Cell1, Cell2, Cell3, Cell4 taking into account boundary conditions.

Function Documentation

void GetBoundaryCells ( Cell &  Cell1,
Cell &  Cell2,
Cell &  Cell3,
Cell &  Cell4,
Cell &  C1,
Cell &  C2,
Cell &  C3,
Cell &  C4,
int  AxisNo 
)

Transform the 4 cells of the flux Cell1, Cell2, Cell3, Cell4 into C1, C2, C3, C4, to take into account boundary conditions (used in Flux.cpp).

Parameters
Cell1second neighbour on the left side
Cell2first neighbour on the left side
Cell3first neighbour on the right side
Cell4second neighbour on the right side
C1Auxiliar cell1
C2Auxiliar cell2
C3Auxiliar cell3
C4Auxiliar cell4
AxisNo...
Returns
void
26 {
27  // --- Local variables ---
28 
29  int InCell1, InCell2, InCell3, InCell4; // Boundary conditions in cells 1, 2, 3, 4
30  real P1, P2, P3, P4; // Pressures in cells 1, 2, 3, 4
31  real T1, T2, T3, T4; // Temperatures in cells 1, 2, 3, 4
32  real rho1, rho2, rho3, rho4; // Densities in cells 1, 2, 3, 4
33  Vector V1(Dimension), V2(Dimension), V3(Dimension), V4(Dimension); // Velocities in cells 1, 2, 3, 4
34  real e1, e2, e3, e4; // Energies in cell 1, 2, 3, 4
35  real Y1=0., Y2=0., Y3=0., Y4=0.; // Partial mass in cell 1, 2, 3, 4
36 
37  int i; // Counter
38 
39  // --- Init C1, C2, C3, C4 ---
40 
41  C1 = Cell1;
42  C2 = Cell2;
43  C3 = Cell3;
44  C4 = Cell4;
45 
46  // --- Depending on the boundary region type, transform C1, C2, C3, C4 ---
47 
48  InCell1 = BoundaryRegion(Cell1.center());
49  InCell2 = BoundaryRegion(Cell2.center());
50  InCell3 = BoundaryRegion(Cell3.center());
51  InCell4 = BoundaryRegion(Cell4.center());
52 
53  // --- Cell2 IN THE BOUNDARY, Cell3 IN THE FLUID -----------------------------------------
54 
55  if (InCell2 !=0 && InCell3 == 0)
56  {
57  switch(InCell2)
58  {
59  // INFLOW
60  case 1:
61  // Dirichlet on temperature
62  T2 = Cell2.temperature();
63  T1 = Cell1.temperature();
64 
65  // Extrapolate pressure
66  P2 = Cell3.oldPressure();
67  P1 = P2;
68  // P1 = 2*P2 - Cell3.pressure();
69 
70  // Compute density
71  rho2 = Gamma*Ma*Ma*P2/T2;
72  rho1 = Gamma*Ma*Ma*P1/T1;
73 
74  // Dirichlet on momentum
75  V2 = (Cell2.density()/rho2)*Cell2.velocity();
76  V1 = (Cell1.density()/rho1)*Cell1.velocity();
77 
78  // Dirichlet on partial mass
79  if (ScalarEqNb == 1)
80  {
81  Y2 = Cell2.average(Dimension+3)/Cell2.density();
82  Y1 = Cell1.average(Dimension+3)/Cell1.density();
83  }
84 
85  // Compute energies
86  e2 = P2/((Gamma-1.)*rho2) + 0.5*N2(V2);
87  e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);
88 
89  // Correct C1 and C2
90  C2.setAverage(1, rho2);
91  C1.setAverage(1, rho1);
92  for (i=1; i<=Dimension; i++)
93  {
94  C2.setAverage(i+1, rho2*V2.value(i));
95  C1.setAverage(i+1, rho1*V1.value(i));
96  }
97  C2.setAverage(Dimension+2,rho2*e2);
98  C1.setAverage(Dimension+2,rho1*e1);
99 
100  if (ScalarEqNb == 1)
101  {
102  C2.setAverage(Dimension+3, rho2*Y2);
103  C1.setAverage(Dimension+3, rho1*Y1);
104  }
105  break;
106 
107  // OUTFLOW : use the old value of the neighbour
108  case 2:
109  C2.setAverage(Cell3.oldAverage());
110  C1.setAverage(Cell3.oldAverage());
111 
112  // Also change the values in the boundary
113  Cell2.setAverage(C2.average());
114  Cell1.setAverage(C1.average());
115  break;
116 
117  // FREE-SLIP WALL : Neuman on all quantities
118  case 3:
119 
120  C2 = Cell3;
121  C1 = Cell4;
122  break;
123 
124  // ADIABATIC WALL
125  case 4:
126 
127  // Dirichlet on velocity
128  V2 = Cell2.velocity();
129  V1 = Cell1.velocity();
130 
131  // Neuman on temperature
132  T2 = Cell3.temperature();
133  T1 = Cell4.temperature();
134 
135  // Neuman on pressure
136  P2 = Cell3.pressure();
137  P1 = Cell4.pressure();
138 
139  // Extrapolate pressure
140  //P2 = 2*Cell3.pressure()-Cell4.pressure();
141  //P1 = P2;
142 
143  // Compute densities
144  rho2 = Gamma*Ma*Ma*P2/T2;
145  rho1 = Gamma*Ma*Ma*P1/T1;
146 
147  // Compute energies
148  e2 = P2/((Gamma-1.)*rho2) + 0.5*N2(V2);
149  e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);
150 
151  // Correct C1 and C2
152  C2.setAverage(1, rho2);
153  C1.setAverage(1, rho1);
154  for (i=1; i<=Dimension; i++)
155  {
156  C2.setAverage(i+1, rho2*V2.value(i));
157  C1.setAverage(i+1, rho1*V1.value(i));
158  }
159  C2.setAverage(Dimension+2,rho2*e2);
160  C1.setAverage(Dimension+2,rho1*e1);
161 
162  // Neuman on partial mass
163  if (ScalarEqNb == 1)
164  {
165  C2.setAverage(Dimension+3, Cell3.average(Dimension+3));
166  C1.setAverage(Dimension+3, Cell4.average(Dimension+3));
167  }
168 
169  break;
170 
171  // ISOTHERMAL WALL
172  case 5:
173 
174  // Dirichlet on velocity
175  V2 = Cell2.velocity();
176  V1 = Cell1.velocity();
177 
178  // Dirichlet on temperature
179  T2 = Cell2.temperature();
180  T1 = Cell1.temperature();
181 
182  // Neuman on pressure
183  P2 = Cell3.pressure();
184  P1 = Cell4.pressure();
185 
186  // Extrapolate pressure
187  //P2 = 2*Cell3.pressure()-Cell4.pressure();
188  //P1 = P2;
189 
190  // Compute densities
191  rho2 = Gamma*Ma*Ma*P2/T2;
192  rho1 = Gamma*Ma*Ma*P1/T1;
193 
194  // Compute energies
195  e2 = P2/((Gamma-1.)*rho2) + 0.5*N2(V2);
196  e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);
197 
198  // Correct C1 and C2
199  C2.setAverage(1, rho2);
200  C1.setAverage(1, rho1);
201  for (i=1; i<=Dimension; i++)
202  {
203  C2.setAverage(i+1, rho2*V2.value(i));
204  C1.setAverage(i+1, rho1*V1.value(i));
205  }
206  C2.setAverage(Dimension+2,rho2*e2);
207  C1.setAverage(Dimension+2,rho1*e1);
208 
209  // Neuman on partial mass
210  if (ScalarEqNb == 1)
211  {
212  C2.setAverage(Dimension+3, Cell3.average(Dimension+3));
213  C1.setAverage(Dimension+3, Cell4.average(Dimension+3));
214  }
215 
216  break;
217  };
218  return;
219  }
220 
221  // --- Cell1 IN THE BOUNDARY, Cell2 IN THE FLUID ------------------------------------------------
222 
223  if (InCell1 != 0 && InCell2 == 0)
224  {
225  switch(InCell1)
226  {
227  // INFLOW
228  case 1:
229  // Dirichlet on temperature
230  T1 = Cell1.temperature();
231 
232  // Extrapolate pressure from old value
233  P1 = Cell2.oldPressure();
234 
235  // Compute density
236  rho1 = Gamma*Ma*Ma*P1/T1;
237 
238  // Dirichlet on momentum
239  V1 = (Cell1.density()/rho1)*Cell1.velocity();
240 
241  // Dirichlet on partial mass
242  if (ScalarEqNb == 1)
243  Y1 = Cell1.average(Dimension+3)/Cell1.density();
244 
245  // Compute energies
246  e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);
247 
248  // Correct C1
249  C1.setAverage(1, rho1);
250  for (i=1; i<=Dimension; i++)
251  C1.setAverage(i+1, rho1*V1.value(i));
252  C1.setAverage(Dimension+2,rho1*e1);
253 
254  if (ScalarEqNb == 1)
255  C1.setAverage(Dimension+3, rho1*Y1);
256  break;
257 
258  // OUTFLOW : Get old value of the neighbour
259  case 2:
260 
261  C1.setAverage(Cell2.oldAverage());
262  break;
263 
264  // FREE-SLIP WALL : Neuman on all quantities
265  case 3:
266 
267  C1 = Cell2;
268  break;
269 
270  // ADIABATIC WALL
271  case 4:
272 
273  // Dirichlet on velocity
274  V1 = Cell1.velocity();
275 
276  // Neuman on temperature
277  T1 = Cell2.temperature();
278 
279  // Neuman on pressure
280  P1 = Cell2.pressure();
281 
282  // Extrapolate pressure
283  //P1 = 2*Cell2.pressure()-Cell3.pressure();
284 
285  // Compute density
286  rho1 = Gamma*Ma*Ma*P1/T1;
287 
288  // Compute energy
289  e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);
290 
291  // Correct C1
292  C1.setAverage(1, rho1);
293  for (i=1; i<=Dimension; i++)
294  C1.setAverage(i+1, rho1*V1.value(i));
295  C1.setAverage(Dimension+2,rho1*e1);
296 
297  // Neuman on partial mass
298  if (ScalarEqNb == 1)
299  C1.setAverage(Dimension+3, Cell2.average(Dimension+3));
300 
301  break;
302 
303  // ISOTHERMAL WALL
304  case 5:
305 
306  // Dirichlet on velocity
307  V1 = Cell1.velocity();
308 
309  // Dirichlet on temperature
310  T1 = Cell1.temperature();
311 
312  // Neuman on pressure
313  P1 = Cell2.pressure();
314 
315  // Extrapolate pressure
316  //P1 = 2*Cell2.pressure()-Cell3.pressure();
317 
318  // Compute density
319  rho1 = Gamma*Ma*Ma*P1/T1;
320 
321  // Compute energies
322  e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);
323 
324  // Correct C1
325  C1.setAverage(1, rho1);
326  for (i=1; i<=Dimension; i++)
327  C1.setAverage(i+1, rho1*V1.value(i));
328  C1.setAverage(Dimension+2,rho1*e1);
329 
330  // Neuman on partial mass
331  if (ScalarEqNb == 1)
332  C1.setAverage(Dimension+3, Cell2.average(Dimension+3));
333 
334  break;
335  };
336  return;
337  }
338  // --- Cell3 IN THE BOUNDARY, Cell2 IN THE FLUID -----------------------------------------
339 
340  if (InCell3 !=0 && InCell2 == 0)
341  {
342  switch(InCell3)
343  {
344  // INFLOW
345  case 1:
346  // Dirichlet on temperature
347  T3 = Cell3.temperature();
348  T4 = Cell4.temperature();
349 
350  // Extrapolate pressure from old value
351  P3 = Cell2.oldPressure();
352  P4 = P3;
353  //P4 = 2*P3 - Cell2.pressure();
354 
355  // Compute densities
356  rho3 = Gamma*Ma*Ma*P3/T3;
357  rho4 = Gamma*Ma*Ma*P4/T4;
358 
359  // Dirichlet on momentum
360  V3 = (Cell3.density()/rho3)*Cell3.velocity();
361  V4 = (Cell4.density()/rho4)*Cell4.velocity();
362 
363  // Dirichlet on partial mass
364  if (ScalarEqNb == 1)
365  {
366  Y3 = Cell3.average(Dimension+3)/Cell3.density();
367  Y4 = Cell4.average(Dimension+3)/Cell4.density();
368  }
369 
370  // Compute energies
371  e3 = P3/((Gamma-1.)*rho3) + 0.5*N2(V3);
372  e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);
373 
374  // Correct C1 and C2
375  C3.setAverage(1, rho3);
376  C4.setAverage(1, rho4);
377  for (i=1; i<=Dimension; i++)
378  {
379  C3.setAverage(i+1, rho3*V3.value(i));
380  C4.setAverage(i+1, rho4*V4.value(i));
381  }
382  C3.setAverage(Dimension+2,rho3*e3);
383  C4.setAverage(Dimension+2,rho4*e4);
384 
385  if (ScalarEqNb == 1)
386  {
387  C3.setAverage(Dimension+3,rho3*Y3);
388  C4.setAverage(Dimension+3,rho4*Y4);
389  }
390  break;
391 
392  // OUTFLOW
393  case 2:
394 
395  C3.setAverage(Cell2.oldAverage());
396  C4.setAverage(Cell2.oldAverage());
397  //C4.setAverage(2*C3.average()-Cell2.average());
398 
399  // Also change the values in the boundary
400  Cell3.setAverage(C3.average());
401  Cell4.setAverage(C4.average());
402  break;
403 
404  // FREE-SLIP WALL : Neuman on all quantities
405  case 3:
406 
407  C3 = Cell2;
408  C4 = Cell1;
409  break;
410 
411  // ADIABATIC WALL
412  case 4:
413 
414  // Dirichlet on velocity
415  V3 = Cell3.velocity();
416  V4 = Cell4.velocity();
417 
418  // Neuman on temperature
419  T3 = Cell2.temperature();
420  T4 = Cell1.temperature();
421 
422  // Neuman on pressure
423  P3 = Cell2.pressure();
424  P4 = Cell1.pressure();
425 
426  // Extrapolate pressure
427  //P3 = 2*Cell2.pressure()-Cell1.pressure();
428  //P4 = P3;
429 
430  // Compute densities
431  rho3 = Gamma*Ma*Ma*P3/T3;
432  rho4 = Gamma*Ma*Ma*P4/T4;
433 
434  // Compute energies
435  e3 = P3/((Gamma-1.)*rho3) + 0.5*N2(V3);
436  e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);
437 
438  // Correct C3 and C4
439  C3.setAverage(1, rho3);
440  C4.setAverage(1, rho4);
441  for (i=1; i<=Dimension; i++)
442  {
443  C3.setAverage(i+1, rho3*V3.value(i));
444  C4.setAverage(i+1, rho4*V4.value(i));
445  }
446  C3.setAverage(Dimension+2,rho3*e3);
447  C4.setAverage(Dimension+2,rho4*e4);
448 
449  // Neuman on partial mass
450  if (ScalarEqNb == 1)
451  {
452  C3.setAverage(Dimension+3, Cell2.average(Dimension+3));
453  C4.setAverage(Dimension+3, Cell1.average(Dimension+3));
454  }
455 
456  break;
457 
458  // ISOTHERMAL WALL
459  case 5:
460 
461  // Dirichlet on velocity
462  V3 = Cell3.velocity();
463  V4 = Cell4.velocity();
464 
465  // Dirichlet on temperature
466  T3 = Cell3.temperature();
467  T4 = Cell4.temperature();
468 
469  // Neuman on pressure
470  P3 = Cell2.pressure();
471  P4 = Cell1.pressure();
472 
473  // Extrapolate pressure
474  //P3 = 2*Cell2.pressure()-Cell1.pressure();
475  //P4 = P3;
476 
477  // Compute densities
478  rho3 = Gamma*Ma*Ma*P3/T3;
479  rho4 = Gamma*Ma*Ma*P4/T4;
480 
481  // Compute energies
482  e3 = P3/((Gamma-1.)*rho3) + 0.5*N2(V3);
483  e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);
484 
485  // Correct C3 and C4
486  C3.setAverage(1, rho3);
487  C4.setAverage(1, rho4);
488  for (i=1; i<=Dimension; i++)
489  {
490  C3.setAverage(i+1, rho3*V3.value(i));
491  C4.setAverage(i+1, rho4*V4.value(i));
492  }
493  C3.setAverage(Dimension+2,rho3*e3);
494  C4.setAverage(Dimension+2,rho4*e4);
495 
496  // Neuman on partial mass
497  if (ScalarEqNb == 1)
498  {
499  C3.setAverage(Dimension+3, Cell2.average(Dimension+3));
500  C4.setAverage(Dimension+3, Cell1.average(Dimension+3));
501  }
502 
503  break;
504  };
505  return;
506  }
507 
508  // --- Cell4 IN THE BOUNDARY, Cell3 IN THE FLUID ------------------------------------------------
509 
510  if (InCell4 != 0 && InCell3 == 0)
511  {
512  switch(InCell4)
513  {
514  // INFLOW
515  case 1:
516  // Dirichlet on temperature
517  T4 = Cell4.temperature();
518 
519  // Extrapolate pressure from old value
520  P4 = Cell3.oldPressure();
521 
522  // Compute density
523  rho4 = Gamma*Ma*Ma*P4/T4;
524 
525  // Dirichlet on momentum
526  V4 = (Cell4.density()/rho4)*Cell4.velocity();
527 
528  // Dirichlet on partial mass
529  if (ScalarEqNb == 1)
530  Y4 = Cell4.average(Dimension+3)/Cell4.density();
531 
532  // Compute energies
533  e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);
534 
535  // Correct C4
536  C4.setAverage(1, rho4);
537  for (i=1; i<=Dimension; i++)
538  C4.setAverage(i+1, rho4*V4.value(i));
539  C4.setAverage(Dimension+2,rho4*e4);
540 
541  if (ScalarEqNb == 1)
542  C4.setAverage(Dimension+3,rho4*Y4);
543  break;
544 
545  // OUTFLOW : Use old cell-average values of the neighbour
546  case 2:
547 
548  C4.setAverage(Cell3.oldAverage());
549  break;
550 
551  // FREE-SLIP WALL : Neuman on all quantities
552  case 3:
553 
554  C4 = Cell3;
555  break;
556 
557  // ADIABATIC WALL
558  case 4:
559 
560  // Dirichlet on velocity
561  V4 = Cell4.velocity();
562 
563  // Neuman on temperature
564  T4 = Cell3.temperature();
565 
566  // Neuman on pressure
567  P4 = Cell3.pressure();
568 
569  // Extrapolate pressure
570  //P4 = 2*Cell3.pressure()-Cell2.pressure();
571 
572  // Compute density
573  rho4 = Gamma*Ma*Ma*P4/T4;
574 
575  // Compute energy
576  e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);
577 
578  // Correct C4
579  C4.setAverage(1, rho4);
580  for (i=1; i<=Dimension; i++)
581  C4.setAverage(i+1, rho4*V4.value(i));
582  C4.setAverage(Dimension+2,rho4*e4);
583 
584  // Neuman on partial mass
585  if (ScalarEqNb == 1)
586  C4.setAverage(Dimension+3, Cell3.average(Dimension+3));
587 
588  break;
589 
590  // ISOTHERMAL WALL
591  case 5:
592 
593  // Dirichlet on velocity
594  V4 = Cell4.velocity();
595 
596  // Dirichlet on temperature
597  T4 = Cell4.temperature();
598 
599  // Neuman on pressure
600  P4 = Cell3.pressure();
601 
602  // Extrapolate pressure
603  //P4 = 2*Cell3.pressure()-Cell2.pressure();
604 
605  // Compute density
606  rho4 = Gamma*Ma*Ma*P4/T4;
607 
608  // Compute energies
609  e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);
610 
611  // Correct C4
612  C4.setAverage(1, rho4);
613  for (i=1; i<=Dimension; i++)
614  C4.setAverage(i+1, rho4*V4.value(i));
615  C4.setAverage(Dimension+2,rho4*e4);
616 
617  // Neuman on partial mass
618  if (ScalarEqNb == 1)
619  C4.setAverage(Dimension+3, Cell3.average(Dimension+3));
620 
621  break;
622  };
623  return;
624  }
625 
626 }
Cell::center
real center(const int AxisNo) const
Returns the component no. AxisNo of the cell-center position.
Definition: Cell.h:1112
Vector
Standard class for every vector in Carmen.
Definition: Vector.h:29
Cell::oldPressure
real oldPressure() const
Identical to the previous one for the values at the instant n-1.
Definition: Cell.cpp:151
Cell::temperature
real temperature() const
Returns the cell-average temperature. Does not work for MHD!
Definition: Cell.cpp:184
N2
real N2(const Vector &V)
Returns the L2-norm of the vector.
Definition: Vector.cpp:1525
Cell::density
real density() const
Returns the cell-average density.
Definition: Cell.h:1262
ScalarEqNb
int ScalarEqNb
Definition: Parameters.cpp:69
Dimension
int Dimension
Definition: Parameters.cpp:74
Cell::velocity
real velocity(const int AxisNo) const
Returns the component no. AxisNo of the cell-average velocity.
Definition: Cell.h:1354
Gamma
real Gamma
Definition: Parameters.cpp:109
Cell::oldAverage
real oldAverage(const int QuantityNo) const
Returns the component no. QuantityNo of the old cell-average values.
Definition: Cell.h:1176
BoundaryRegion
int BoundaryRegion(const Vector &X)
Returns the boundary region type at the position X=(x,y,z). The returned value correspond to: 0 = Fl...
Definition: BoundaryRegion.cpp:63
Cell::average
real average(const int QuantityNo) const
Returns the component no. QuantityNo of the cell-average value.
Definition: Cell.h:1128
Cell::pressure
real pressure() const
Returns the cell-average pressure.
Definition: Cell.cpp:84
Cell::setAverage
void setAverage(const int QuantityNo, const real UserAverage)
Sets the cell-average of the quantity QuantityNo to UserAverage. Example:
Definition: Cell.h:921
Ma
real Ma
Definition: Parameters.cpp:110
real
#define real
Definition: PreProcessor.h:31

Here is the caller graph for this function:

Generated on Wed Sep 6 2017 14:34:41 for Carmen Code by   doxygen 1.8.6