SPI
SLEPc PETSc Interface is used for rapid development and intuitive matrix operations similar to MatLab or NumPy style of operations
SPIbaseflow.cpp
Go to the documentation of this file.
1 #include "SPIbaseflow.hpp"
2 
3 namespace SPI{
4  // constructors
6  SPIbaseflow::SPIbaseflow(
7  std::string _name
8  ){ this->name = _name; }
10  SPIbaseflow::SPIbaseflow(
11  SPIVec U,
12  SPIVec V,
13  SPIVec Ux,
14  SPIVec Uxx,
15  SPIVec Uy,
16  SPIVec Uxy,
17  SPIVec Vx,
18  SPIVec Vxx,
19  SPIVec Vy,
20  SPIVec W,
21  SPIVec Wx,
22  SPIVec Wy,
23  SPIVec Wxy,
24  SPIVec P,
25  std::string _name
26  ){
27  this->name = _name;
28  this->U = U; this->U.name = "U";
29  this->V = V; this->V.name = "V";
30  this->Ux = Ux; this->Ux.name = "Ux";
31  this->Uxx = Uxx; this->Uxx.name = "Uxx";
32  this->Uy = Uy; this->Uy.name = "Uy";
33  this->Uxy = Uxy; this->Uxy.name = "Uxy";
34  this->Vx = Vx; this->Vx.name = "Vx";
35  this->Vxx = Vxx; this->Vxx.name = "Vxx";
36  this->Vy = Vy; this->Vy.name = "Vy";
37  this->W = W; this->W.name = "W";
38  this->Wx = Wx; this->Wx.name = "Wx";
39  this->Wy = Wy; this->Wy.name = "Wy";
40  this->Wxy = Wxy; this->Wxy.name = "Wxy";
41  this->P = P; this->P.name = "P";
42  this->flag_init=PETSC_TRUE;
43  }
45  PetscInt SPIbaseflow::print(){
46  // print baseflow values if initialized
47  PetscPrintf(PETSC_COMM_WORLD,("---------------- "+name+"---done-------\n\n").c_str());
48  if(this->flag_init){
49  this->U.print();
50  this->V.print();
51  this->Ux.print();
52  this->Uxx.print();
53  this->Uy.print();
54  this->Uxy.print();
55  this->Vx.print();
56  this->Vxx.print();
57  this->Vy.print();
58  this->W.print();
59  this->Wx.print();
60  this->Wy.print();
61  this->Wxy.print();
62  this->P.print();
63  }else{
64  SPI::printf(this->name+std::string(" not initialized"));
65  }
66  PetscPrintf(PETSC_COMM_WORLD,("---------------- "+name+"---done-------\n\n").c_str());
67  return 0;
68  }
69 
71  SPIbaseflow::~SPIbaseflow(){
72  if(this->flag_init){
73  this->U.~SPIVec();
74  this->V.~SPIVec();
75  this->Ux.~SPIVec();
76  this->Uxx.~SPIVec();
77  this->Uy.~SPIVec();
78  this->Uxy.~SPIVec();
79  this->Vx.~SPIVec();
80  this->Vxx.~SPIVec();
81  this->Vy.~SPIVec();
82  this->W.~SPIVec();
83  this->Wx.~SPIVec();
84  this->Wy.~SPIVec();
85  this->Wxy.~SPIVec();
86  this->P.~SPIVec();
87  }
88  }
89  /* \brief set Blasius boundary layer flow \return SPIbaseflow at the current conditions */
90  SPIbaseflow blasius(
91  SPIparams &params,
92  SPIgrid1D &grid
93  ){ // if base flow is Blasius Flat-Plate
94  PetscInt multiply_nypts = 8; //data.multiply_nypts_for_bblf;
95  PetscScalar dy,jfloat;
96  PetscScalar multiply_nypts_float=multiply_nypts;
97  PetscScalar eta[multiply_nypts*grid.ny];
98  PetscScalar deta[multiply_nypts*grid.ny-1];
99  // TODO: fairly inefficient way to calculate on all cores, then save it in parallel
100  for(int i=0; i<grid.ny; ++i){
101  dy = grid.y(i+1,PETSC_TRUE)-grid.y(i,PETSC_TRUE);
102  //std::cout<<"y["<<i<<"] = "<<data.y[i]<<std::endl;
103  for(int j=0; j<multiply_nypts; j++){
104  PetscInt ij = i*multiply_nypts+j;
105  jfloat=(PetscScalar)j;
106  eta[ij] = (grid.y(i,PETSC_TRUE)+dy*jfloat/multiply_nypts_float)*PetscSqrtScalar(1./(2.*params.nu*params.x)); // asume U_inf = 1
107  }
108  }
109  for(int i=0; i<grid.ny*multiply_nypts-1; ++i){ // set spacing of eta
110  deta[i] = eta[i+1] - eta[i];
111  }
112  // initialize variable
113  PetscScalar **fs = new PetscScalar*[grid.ny*multiply_nypts];// ny pts to solve
114  for(int i=0; i<grid.ny*multiply_nypts; ++i) fs[i] = new PetscScalar[3]; // 3 variables
115  PetscScalar k1[3],k2[3],k3[3],k4[3];// RK4 intermediate variables
116  PetscScalar temp[3],temp2[3];// temporary array to hold values
117  // set initial conditions (ICs)
118  //fs[0][0] = 0.469600; // f''[0] = 0.469600
119  //fs[0][0] = 0.33205733621519630*sqrt(2.); // f''[0] = 0.469600
120  fs[0][0] = 0.332057336215195*std::sqrt(2.);
121  fs[0][1] = 0.; // f'[0] = 0
122  fs[0][2] = 0.; // f[0] = 0
123 
124  // march through each eta value using _bblf function evaluation
125  for(int i=0; i<grid.ny*multiply_nypts-1; ++i){ // for each eta value
126  _bblf(fs[i],temp); // calculate f(i)
127  for(int j=0; j<3; ++j){
128  k1[j] = deta[i]*temp[j];
129  temp2[j] = fs[i][j] + k1[j]/2.;
130  }
131  _bblf(temp2,temp); // calculate f(i+k1/2);
132  for(int j=0; j<3; ++j){
133  k2[j] = deta[i]*temp[j];
134  temp2[j] = fs[i][j] + k2[j]/2.;
135  }
136  _bblf(temp2,temp); // calculate f(i+k2/2);
137  for(int j=0; j<3; ++j){
138  k3[j] = deta[i]*temp[j];
139  temp2[j] = fs[i][j] + k3[j];
140  }
141  _bblf(temp2,temp); // calculate f(i+k3);
142  for(int j=0; j<3; ++j){
143  k4[j] = deta[i]*temp[j];
144  }
145  for(int j=0; j<3; ++j){ // f[i+1] = f[i] + (k1+2k2 + 2k3 + k4)/6
146  fs[i+1][j] = fs[i][j] + (k1[j] + (k2[j]*2.) + (k3[j]*2.) + k4[j])/6.;
147  //std::cout<<"fs["<<i+1<<"]["<<j<<"] = "<<fs[i+1][j]<<std::endl;
148  //std::cout<<" k1,k2,k3,k4 = "<<k1[j]<<", "<<k2[j]<<", "<<k3[j]<<", "<<k4[j]<<std::endl;
149  }
150  }
151  // print eta and fp, fpp to compare to textbook values
152  //for(int i=0; i<grid.ny; ++i){
153  //PetscInt j=0;
154  //PetscInt ij = i*multiply_nypts+j;
155  //printf("eta = %g and f = %g and f' = %g and f'' = %g",PetscRealPart(eta[ij]),PetscRealPart(fs[ij][2]),PetscRealPart(fs[ij][1]),PetscRealPart(fs[ij][0]));
156  //}
157 
158  // save values
159  PetscScalar *fpp = new PetscScalar[grid.ny*multiply_nypts];
160  PetscScalar *fp = new PetscScalar[grid.ny*multiply_nypts];
161  PetscScalar *f = new PetscScalar[grid.ny*multiply_nypts];
162  PetscScalar *Utmp = new PetscScalar[grid.ny*multiply_nypts];
163  PetscScalar *Uxtmp = new PetscScalar[grid.ny*multiply_nypts];
164  PetscScalar *Uytmp = new PetscScalar[grid.ny*multiply_nypts];
165  PetscScalar *Vtmp = new PetscScalar[grid.ny*multiply_nypts];
166  PetscScalar *Vxtmp = new PetscScalar[grid.ny*multiply_nypts];
167  PetscScalar *Vytmp = new PetscScalar[grid.ny*multiply_nypts];
168  for (int i=0; i<grid.ny*multiply_nypts; ++i){
169  fpp[i] = fs[i][0];
170  fp [i] = fs[i][1];
171  f [i] = fs[i][2];
172 
173  // save U
174  Utmp[i] = fp[i]; // U =f' assuming Uinf=1
175  Uxtmp[i]= fpp[i]*(-eta[i]/(2.*params.x)); // Ux =f''(-eta/(2x))
176  Uytmp[i]= fpp[i]*PetscSqrtScalar(1./(2.*params.nu*params.x)); // Uy =f''sqrt(1./(2 nu x))
177  // save V
178  Vtmp[i] = PetscSqrtScalar(params.nu/(2.*params.x))*(eta[i]*fp[i] - f[i]); // V = sqrt(nu/2x)(eta f' - f)
179  Vxtmp[i]= PetscSqrtScalar(params.nu/(8.*PetscPowScalar(params.x,3)))*
180  (
181  - eta[i]*fp[i]
182  + f[i]
183  - PetscPowScalar(eta[i],2)*fpp[i]
184  ); // Vx = sqrt(nu/8x^3)(-eta fp + f - eta^2 f'')
185  Vytmp[i]= (1./(2.*params.x))*eta[i]*fpp[i]; // Vy = (1/2x) eta f''
186 
187 
188  }
189 
190  // Checks
191  // Check divergence
192  PetscScalar *divergence = new PetscScalar[grid.ny*multiply_nypts];
193  for(int i=0; i<grid.ny*multiply_nypts; ++i) divergence[i] = Uxtmp[i] + Vytmp[i];
194  PetscScalar sum_divergence = 0.;
195  for(int i=0; i<grid.ny*multiply_nypts; ++i) sum_divergence += PetscAbsScalar(divergence[i]);
196  printfc("the sum of the divergence of Blasius boundary flow is %g+%gi",sum_divergence);
197 
198  // save values for SPE in grid structure
199  SPIVec U(grid.ny), Uy(grid.ny), Ux(grid.ny);
200  SPIVec V(grid.ny), Vy(grid.ny), Vx(grid.ny);
201  SPIVec O(zeros(grid.ny));;
202  for(int i=0; i<grid.ny; i++){
203  U(i,Utmp[i*multiply_nypts]);
204  Uy(i,Uytmp[i*multiply_nypts]);
205  Ux(i,Uxtmp[i*multiply_nypts]);
206  V(i,Vtmp[i*multiply_nypts]);
207  Vy(i,Vytmp[i*multiply_nypts]);
208  Vx(i,Vxtmp[i*multiply_nypts]);
209  // set freestream velocities
210  }
211  U();
212  Ux();
213  Uy();
214  V();
215  Vx();
216  Vy();
217 
218  //SPIbaseflow baseflow(U, V, Ux, Uy, grid.Dy*Ux, Vy, O, O, O, O, O); // doesn't project correctly
219  SPIbaseflow baseflow(U, V, Ux, O,grid.Dy*U, grid.Dy*Ux, Vx, O, grid.Dy*V, O, O, O, O, O); // bad for UltraS grid
220  if(grid.ytype==SPI::UltraS){ // fix for UltraS grid
221  SPIgrid1D gridCheby(grid.y,"gridCheby",SPI::Chebyshev);
222  SPIMat Dyp(grid.T*grid.S0invS1inv*grid.Dy*grid.That);
223  baseflow.Uy = Dyp*U;
224  baseflow.Uxy = Dyp*Ux;
225  baseflow.Vy = Dyp*V;
226  }
227 
228  // clear memory
229  for(int i=0; i<grid.ny*multiply_nypts; ++i) delete[] fs[i];
230  delete[] fs;
231  delete[] fpp;
232  delete[] fp;
233  delete[] f;
234  delete[] Utmp;
235  delete[] Uxtmp;
236  delete[] Uytmp;
237  delete[] Vtmp;
238  delete[] Vxtmp;
239  delete[] Vytmp;
240 
241  return baseflow;
242  }
243  int _bblf(
244  const PetscScalar input[3],
245  PetscScalar output[3]){
246  output[0] = -input[2]*input[0]; // f''= \int -f*f'' deta
247  output[1] = input[0]; // f' = \int f'' deta
248  output[2] = input[1]; // f = \int f' deta
249  return 0;
250  }
251 
252  /* \brief calculate baseflow for Plane Poiseuille flow \return SPIbaseflow of the Plane Poiseuille flow */
253  SPIbaseflow channel(
254  SPIparams &params,
255  SPIgrid1D &grid
256  ){
257  SPI::SPIVec U((1.0-((grid.y)^2)),"U");
258  SPI::SPIVec Uy((-2.*grid.y),"Uy");
259  SPI::SPIVec o(U*0.0,"o"); // zero vector
260  SPI::SPIbaseflow channel_flow(U,o,o,o,Uy,o,o,o,o,o,o,o,o,o);
261  return channel_flow;
262  }
263 
264 }
SPI::SPIbaseflow::Uxy
SPIVec Uxy
Definition: SPIbaseflow.hpp:44
SPI::SPIbaseflow::Uxx
SPIVec Uxx
Definition: SPIbaseflow.hpp:42
SPI::SPIbaseflow::V
SPIVec V
Definition: SPIbaseflow.hpp:40
SPI::SPIbaseflow::Vxx
SPIVec Vxx
Definition: SPIbaseflow.hpp:46
SPI::SPIgrid1D
Class to contain various grid parameters.
Definition: SPIgrid.hpp:35
SPI::SPIgrid1D::T
SPIMat T
Chebyshev operator taking it from Chebyshev coefficients to physical space.
Definition: SPIgrid.hpp:58
SPI::SPIbaseflow::Vx
SPIVec Vx
Definition: SPIbaseflow.hpp:45
SPI::SPIgrid1D::y
SPIVec y
grid
Definition: SPIgrid.hpp:47
SPI::SPIgrid1D::Dy
SPIMat Dy
1st derivative operator with respect to y
Definition: SPIgrid.hpp:51
SPI::Chebyshev
@ Chebyshev
Chebyshev collocated grid.
Definition: SPIgrid.hpp:29
SPI::SPIVec
Definition: SPIVec.hpp:13
SPI::SPIbaseflow::name
std::string name
baseflow name
Definition: SPIbaseflow.hpp:35
SPI::printf
PetscInt printf(std::string msg,...)
Definition: SPIprint.cpp:5
SPI::SPIbaseflow::Vy
SPIVec Vy
Definition: SPIbaseflow.hpp:47
SPI
Definition: SPIbaseflow.hpp:16
SPI::SPIbaseflow::print
PetscInt print()
print baseflow values
Definition: SPIbaseflow.cpp:45
SPI::SPIbaseflow::Wy
SPIVec Wy
Definition: SPIbaseflow.hpp:50
SPI::SPIbaseflow::P
SPIVec P
Definition: SPIbaseflow.hpp:52
SPI::SPIVec::print
PetscInt print()
Definition: SPIVec.cpp:470
SPI::UltraS
@ UltraS
UltraSpherical grid and derivatives.
Definition: SPIgrid.hpp:30
SPI::SPIbaseflow::Uy
SPIVec Uy
Definition: SPIbaseflow.hpp:43
SPI::SPIbaseflow::SPIbaseflow
SPIbaseflow(std::string _name="baseflow")
constructor with no arguments
Definition: SPIbaseflow.cpp:6
SPI::SPIbaseflow::Wx
SPIVec Wx
Definition: SPIbaseflow.hpp:49
SPI::SPIgrid1D::That
SPIMat That
Chebyshev operator taking it from physical space to Chebyshev coefficients.
Definition: SPIgrid.hpp:59
SPI::SPIparams::x
PetscScalar x
current streamwise position
Definition: SPIparams.hpp:18
SPI::SPIVec::name
std::string name
Vec name (important for hdf5 i/o)
Definition: SPIVec.hpp:26
SPI::blasius
SPIbaseflow blasius(SPIparams &params, SPIgrid1D &grid)
Definition: SPIbaseflow.cpp:90
SPI::SPIgrid1D::S0invS1inv
SPIMat S0invS1inv
[in] inverse of S0^-1 * S1^-1
Definition: SPIgrid.hpp:56
SPI::SPIparams::nu
PetscScalar nu
kinematic viscosity (typically 1/Re)
Definition: SPIparams.hpp:20
SPI::SPIbaseflow::Ux
SPIVec Ux
Definition: SPIbaseflow.hpp:41
SPI::printfc
PetscInt printfc(std::string msg, const PetscScalar val)
Definition: SPIprint.cpp:29
SPI::_bblf
int _bblf(const PetscScalar input[3], PetscScalar output[3])
Definition: SPIbaseflow.cpp:243
SPI::SPIbaseflow::Wxy
SPIVec Wxy
Definition: SPIbaseflow.hpp:51
SPI::SPIgrid1D::ny
PetscInt ny
Definition: SPIgrid.hpp:39
SPI::channel
SPIbaseflow channel(SPIparams &params, SPIgrid1D &grid)
Definition: SPIbaseflow.cpp:253
SPI::SPIbaseflow::W
SPIVec W
Definition: SPIbaseflow.hpp:48
SPIbaseflow.hpp
SPI::SPIbaseflow::U
SPIVec U
Definition: SPIbaseflow.hpp:39
SPI::SPIbaseflow::flag_init
PetscBool flag_init
flag if it has been initialized
Definition: SPIbaseflow.hpp:37
SPI::SPIparams
Definition: SPIparams.hpp:8
SPI::SPIgrid1D::ytype
gridtype ytype
type of grid
Definition: SPIgrid.hpp:48
SPI::sqrt
SPIVec sqrt(const SPIVec &A)
take the atanh of each element in a vector
Definition: SPIVec.cpp:741
SPI::SPIVec::~SPIVec
~SPIVec()
Definition: SPIVec.cpp:483
SPI::SPIMat
Definition: SPIMat.hpp:17
SPI::SPIbaseflow::~SPIbaseflow
~SPIbaseflow()
destructor of saved SPIVec of baseflow
Definition: SPIbaseflow.cpp:71
SPI::SPIbaseflow
Definition: SPIbaseflow.hpp:17
SPI::zeros
SPIMat zeros(const PetscInt m, const PetscInt n)
create, form, and return zeros matrix of size mxn
Definition: SPIMat.cpp:708
Generated by   doxygen 1.8.17