NeoPZ
TPZDiffusionConsLaw.cpp
Go to the documentation of this file.
1 
5 #include "TPZDiffusionConsLaw.h"
6 #include "TPZCompElDisc.h"
7 #include "pzfmatrix.h"
8 #include "pzvec.h"
9 #include "pzreal.h"
10 using namespace std;
11 
12 REAL TPZDiffusionConsLaw::fGamma = 1.4;
13 REAL TPZDiffusionConsLaw::fDelta = -1.0;
14 REAL TPZDiffusionConsLaw::fCFL = 0.0;
15 std::string TPZDiffusionConsLaw::fArtificialDiffusion = "LS";
16 
17 TPZDiffusionConsLaw::~TPZDiffusionConsLaw(){}
18 
19 TPZDiffusionConsLaw::TPZDiffusionConsLaw(){
20 
21  fDimension = 0;
22  fA.Redim(0,0);
23  fB.Redim(0,0);
24  fC.Redim(0,0);
25 }
26 
27 TPZDiffusionConsLaw::TPZDiffusionConsLaw(TPZVec<REAL> U,REAL gamma,int dim,const std::string &diff) :
28 fA(0,0), fB(0,0), fC(0,0) {
29  fGamma = gamma;
30  fDimension = dim;
31  fArtificialDiffusion = diff;
32  int nstate = 2 + dim;
33  fA.Redim(nstate,nstate);
34  fB.Redim(nstate,nstate);
35  fC.Redim(nstate,nstate);
36  JacobFlux(U,fA,fB,fC);
37 }
38 
39 void TPZDiffusionConsLaw::GradientOfTheFlow(TPZFMatrix<REAL> &DF1,TPZFMatrix<REAL> &DF2,TPZFMatrix<REAL> &DF3){
40 
41  DF1 = fA;
42  DF2 = fB;
43  DF3 = fC;
44 }
45 
46 REAL TPZDiffusionConsLaw::CFL(int degree){
47 
48  if(fCFL != 0.) return fCFL;
49  return (1.0/(2.0*(REAL)degree+1.0));
50 }
51 
52 REAL TPZDiffusionConsLaw::Delta(){
53  return fDelta;
54 }
55 
56 REAL TPZDiffusionConsLaw::DeltaOtimo(){
57 
58  int degree = TPZCompEl::GetgOrder();
59  REAL cfl = CFL(degree);
60  REAL delta = ( (10./3.)*cfl*cfl - (2./3.)*cfl + 1./10. );
61  return delta;
62 }
63 
64 void TPZDiffusionConsLaw::Divergence(TPZVec<REAL> &dphi,TPZFMatrix<REAL> &diverg){
65 
66  fArtificialDiffusion = "LS";
67  std::string type = fArtificialDiffusion;
68  TPZFMatrix<REAL> diffterm(0,0);
69  PointOperator(dphi,diffterm);
70  int i,j,nstate=2+fDimension;
71  diverg.Redim(nstate,1);
72  diverg.Zero();
73  // soma das linhas de diffterm da : div(F) = �1/� [�/�] + �2/� [�/�] + �3/� [�/�]
74  for(i=0;i<nstate;i++) for(j=0;j<nstate;j++) diverg(i,0) += diffterm(i,j);// 5x1 (3D)
75  fArtificialDiffusion = type;
76 }
77 
78 void TPZDiffusionConsLaw::PointOperator(TPZVec<REAL> &dphi,TPZFMatrix<REAL> &diff_term){
79 
80  diff_term.Zero();
81  int size = dphi.NElements();
82  int nstate = 2+fDimension;
83  if(size != 3){
84  //esta forma permite trabalhar com 1, 2 e 3 dimens�s
85  cout << "TPZDiffusionConsLaw::PointOperator error data size";
86  }
87  TPZFMatrix<REAL> Tx(nstate,nstate),Ty(nstate,nstate),Tz(nstate,nstate);
88  TPZFMatrix<REAL> Trx(nstate,nstate),Try(nstate,nstate),Trz(nstate,nstate);
89  diff_term.Redim(nstate,nstate);
90  Tau(Tx,Ty,Tz);
91  Tx.Transpose(&Trx);
92  Ty.Transpose(&Try);
93  Tz.Transpose(&Trz);
94  int i,j;
95  // �/� * (T1)' + �/� * (T2)' + �/� * (T3)'
96  for(i=0;i<nstate;i++){
97  for(j=0;j<nstate;j++){
98  diff_term(i,j) = dphi[0] * Trx(i,j) + dphi[1] * Try(i,j) + dphi[2] * Trz(i,j);
99  }
100  }
101 }
102 
103 void TPZDiffusionConsLaw::Tau(TPZFMatrix<REAL> &Tx,TPZFMatrix<REAL> &Ty,TPZFMatrix<REAL> &Tz){
104 
105  if(!strcmp(fArtificialDiffusion.c_str(),"SUPG")){
106  SUPG(Tx,Ty,Tz);
107  return;
108  }
109  if(!strcmp(fArtificialDiffusion.c_str(),"LS")){
110  LS(Tx,Ty,Tz);
111  return;
112  }
113  if(!strcmp(fArtificialDiffusion.c_str(),"BORNHAUS")){
114  Bornhaus(Tx,Ty,Tz);
115  return;
116  }
117  cout << "\nTPZDiffusionConsLaw::Tau case not implemented\n";
118 }
119 
120 void TPZDiffusionConsLaw::SUPG(TPZFMatrix<REAL> &Tx,TPZFMatrix<REAL> &Ty,TPZFMatrix<REAL> &Tz){
121  cout << "TPZDiffusionConsLaw:: SUPG artificial diffusion SUPG not implemented\n";
122 }
123 
124 void TPZDiffusionConsLaw::LS(TPZFMatrix<REAL> &Tx,TPZFMatrix<REAL> &Ty,TPZFMatrix<REAL> &Tz){
125 
126  fA.Transpose(&Tx);
127  fB.Transpose(&Ty);
128  fC.Transpose(&Tz);
129 }
130 
131 void TPZDiffusionConsLaw::Bornhaus(TPZFMatrix<REAL> &Tx,TPZFMatrix<REAL> &Ty,TPZFMatrix<REAL> &Tz){
132  cout << "TPZDiffusionConsLaw::Bornhaus artificial diffusion Bornhaus not implemented\n";
133 }
134 
135 void TPZDiffusionConsLaw::JacobFlux(TPZVec<REAL> U,TPZFMatrix<REAL> &A,TPZFMatrix<REAL> &B,TPZFMatrix<REAL> &C) {//OK
136 
137  if(fabs(U[0]) < 1.e-6) {
138  cout << "\nTPZDiffusionConsLaw::JacobFlux: Densidade quase nula, o jacobiano falha\n";
139  cout << "Densidade = " << U[0] << endl;
140  U[0] = 1.e-6;
141  cout << "Nova densidade = " << U[0] << endl;
142  //cout << "Programa abortado\n\n";
143  //exit(-1);
144  }
145 
146  int cap = A.Rows();
147  if(cap < 3 || cap > 5){
148  cout << "\n\nTPZDiffusionConsLaw::JacobFlux case not trated\n\n";
149  A.Redim(0,0);
150  B.Redim(0,0);
151  C.Redim(0,0);
152  exit(-1);
153  }
154 
155  REAL u,v,w,e;
156  REAL gamma1 = fGamma-1.;
157  REAL gamma2 = gamma1/2.;
158  REAL gamma3 = fGamma-3;
159 
160  if(cap == 5){
161 
162  A.Redim(5,5);
163  B.Redim(5,5);
164  C.Redim(5,5);
165 
166  u = U[1]/U[0];
167  v = U[2]/U[0];
168  w = U[3]/U[0];
169  e = U[4]/U[0];
170 
171  REAL u2 = u*u;
172  REAL v2 = v*v;
173  REAL w2 = w*w;
174  REAL vel = u2+v2+w2;
175 
176  A(0,0) = 0.;
177  A(0,1) = 1.;
178  A(0,2) = 0.;
179  A(0,3) = 0.;
180  A(0,4) = 0.;
181 
182  A(1,0) = gamma2*vel-u2;
183  A(1,1) = -gamma3*u;
184  A(1,2) = -gamma1*v;
185  A(1,3) = -gamma1*w;
186  A(1,4) = gamma1;
187 
188  A(2,0) = -u*v;
189  A(2,1) = v;
190  A(2,2) = u;
191  A(2,3) = 0.;
192  A(2,4) = 0.;
193 
194  A(3,0) = -u*w;
195  A(3,1) = w;
196  A(3,2) = 0.;
197  A(3,3) = u;
198  A(3,4) = 0.;
199 
200  A(4,0) = -fGamma*e*u + gamma1*u*vel;
201  A(4,1) = fGamma*e - gamma1*u2 - gamma2*vel;
202  A(4,2) = -gamma1*u*v;
203  A(4,3) = -gamma1*u*w;
204  A(4,4) = fGamma*u;
205 
206  B(0,0) = 0.;
207  B(0,1) = 0.;
208  B(0,2) = 1.;
209  B(0,3) = 0.;
210  B(0,4) = 0.;
211 
212  B(1,0) = -u*v;
213  B(1,1) = v;
214  B(1,2) = u;
215  B(1,3) = 0.;
216  B(1,4) = 0.;
217 
218  B(2,0) = gamma2*vel-v2;
219  B(2,1) = -gamma1*u;
220  B(2,2) = -gamma3*v;
221  B(2,3) = -gamma1*w;
222  B(2,4) = gamma1;
223 
224  B(3,0) = -v*w;
225  B(3,1) = 0.;
226  B(3,2) = w;
227  B(3,3) = v;
228  B(3,4) = 0.;
229 
230  B(4,0) = -fGamma*e*v + gamma1*v*vel;
231  B(4,1) = -gamma1*u*v;
232  B(4,2) = fGamma*e - gamma1*v2 - gamma2*vel;
233  B(4,3) = -gamma1*v*w;
234  B(4,4) = fGamma*v;
235 
236  C(0,0) = 0.;
237  C(0,1) = 0.;
238  C(0,2) = 0.;
239  C(0,3) = 1.;
240  C(0,4) = 0.;
241 
242  C(1,0) = -u*w;
243  C(1,1) = w;
244  C(1,2) = 0.;
245  C(1,3) = u;
246  C(1,4) = 0.;
247 
248  C(2,0) = -v*w;
249  C(2,1) = 0.;
250  C(2,2) = w;
251  C(2,3) = v;
252  C(2,4) = 0.;
253 
254  C(3,0) = gamma2*vel-w2;
255  C(3,1) = -gamma1*u;
256  C(3,2) = -gamma1*v;
257  C(3,3) = -gamma3*w;
258  C(3,4) = gamma1;
259 
260  C(4,0) = -fGamma*e*w + gamma1*w*vel;
261  C(4,1) = -gamma1*u*w;
262  C(4,2) = -gamma1*v*w;
263  C(4,3) = fGamma*e - gamma1*w2 - gamma2*vel;
264  C(4,4) = fGamma*w;
265 
266  } else if(cap == 4){
267 
268  A.Redim(4,4);
269  B.Redim(4,4);
270  C.Redim(4,4);
271 
272  u = U[1]/U[0];
273  v = U[2]/U[0];
274  e = U[3]/U[0];
275 
276  REAL u2 = u*u;
277  REAL v2 = v*v;
278  REAL vel = u2+v2;
279 
280  A(0,0) = 0.;
281  A(0,1) = 1.;
282  A(0,2) = 0.;
283  A(0,3) = 0.;
284 
285  A(1,0) = gamma2*vel-u2;
286  A(1,1) = -gamma3*u;
287  A(1,2) = -gamma1*v;
288  A(1,3) = gamma1;
289 
290  A(2,0) = -u*v;
291  A(2,1) = v;
292  A(2,2) = u;
293  A(2,3) = 0.;
294 
295  A(3,0) = -fGamma*e*u + gamma1*u*vel;
296  A(3,1) = fGamma*e - gamma1*u2 - gamma2*vel;
297  A(3,2) = -gamma1*u*v;
298  A(3,3) = fGamma*u;
299 
300  B(0,0) = 0.;
301  B(0,1) = 0.;
302  B(0,2) = 1.;
303  B(0,3) = 0.;
304 
305  B(1,0) = -u*v;
306  B(1,1) = v;
307  B(1,2) = u;
308  B(1,3) = 0.;
309 
310  B(2,0) = gamma2*vel-v2;
311  B(2,1) = -gamma1*u;
312  B(2,2) = -gamma3*v;
313  B(2,3) = gamma1;
314 
315  B(3,0) = -fGamma*e*v + gamma1*v*vel;
316  B(3,1) = -gamma1*u*v;
317  B(3,2) = fGamma*e - gamma1*v2 - gamma2*vel;
318  B(3,3) = fGamma*v;
319 
320  } else if(cap == 3){
321 
322  A.Redim(3,3);
323  B.Redim(3,3);
324  C.Redim(3,3);
325 
326  u = U[1]/U[0];
327  e = U[2]/U[0];
328 
329  REAL u2 = u*u;
330  REAL vel = u2;
331 
332  A(0,0) = 0.;
333  A(0,1) = 1.;
334  A(0,2) = 0.;
335 
336  A(1,0) = gamma2*vel-u2;
337  A(1,1) = -gamma3*u;
338  A(1,2) = gamma1;
339 
340  A(2,0) = -fGamma*e*u + gamma1*u*vel;
341  A(2,1) = fGamma*e - gamma1*u2 - gamma2*vel;
342  A(2,2) = fGamma*u;
343  }
344 }
345 
346 //left = **_f right = **_t
347 void TPZDiffusionConsLaw::Roe_Flux(
348  REAL rho_f, REAL rhou_f, REAL rhov_f, REAL rhow_f,
349  REAL rhoE_f,REAL rho_t, REAL rhou_t, REAL rhov_t, REAL rhow_t,
350  REAL rhoE_t, REAL nx, REAL ny, REAL nz, REAL gam,
351  REAL &flux_rho, REAL &flux_rhou, REAL &flux_rhov,
352  REAL &flux_rhow, REAL &flux_rhoE){
353 
354  REAL alpha1,alpha2,alpha3,alpha4,alpha5,alpha;
355  REAL a1,a2,a3,a4,a5,b1,b2,b3,b4,b5;
356  REAL ep_t, ep_f, p_t, p_f;
357  REAL rhouv_t, rhouv_f, rhouw_t, rhouw_f, rhovw_t, rhovw_f;
358  REAL lambda_f, lambda_t;
359  REAL delta_rho, delta_rhou, delta_rhov, delta_rhow, delta_rhoE;
360  REAL hnx, hny, hnz;
361  REAL tempo11, usc;
362 
363  flux_rho = 0;
364  flux_rhou = 0;
365  flux_rhov = 0;
366  flux_rhow = 0;
367  flux_rhoE = 0;
368 
369  REAL gam1 = gam - 1.0;
370  REAL irho_f = 1.0/rho_f;
371  REAL irho_t = 1.0/rho_t;
372 
373  //
374  //.. Compute the ROE Averages
375  //
376  //.... some useful quantities
377  REAL coef1 = sqrt(rho_f);
378  REAL coef2 = sqrt(rho_t);
379  REAL somme_coef = coef1 + coef2;
380  REAL isomme_coef = 1.0/somme_coef;
381  REAL u_f = rhou_f*irho_f;
382  REAL v_f = rhov_f*irho_f;
383  REAL w_f = rhow_f*irho_f;
384  REAL h_f = (gam * rhoE_f*irho_f) - (.5*gam1) * (u_f * u_f + v_f * v_f + w_f * w_f);
385  REAL u_t = rhou_t*irho_t;
386  REAL v_t = rhov_t*irho_t;
387  REAL w_t = rhow_t*irho_t;
388  REAL h_t = (gam * rhoE_t*irho_t) - (.5*gam1) * (u_t * u_t + v_t * v_t + w_t * w_t);
389 
390  //.... averages
391  //REAL rho_ave = coef1 * coef2;
392  REAL u_ave = (coef1 * u_f + coef2 * u_t) * isomme_coef;
393  REAL v_ave = (coef1 * v_f + coef2 * v_t) * isomme_coef;
394  REAL w_ave = (coef1 * w_f + coef2 * w_t) * isomme_coef;
395  REAL h_ave = (coef1 * h_f + coef2 * h_t) * isomme_coef;
396  //
397  //.. Compute Speed of sound
398  REAL scal = u_ave * nx + v_ave * ny + w_ave * nz;
399  REAL norme = sqrt(nx * nx + ny * ny + nz * nz);
400  REAL inorme = 1.0/norme;
401  REAL u2pv2pw2 = u_ave * u_ave + v_ave * v_ave + w_ave * w_ave;
402  REAL c_speed = gam1 * (h_ave - 0.5 * u2pv2pw2);
403  if(c_speed < 1e-6) c_speed = 1e-6; // avoid division by 0 if critical
404  c_speed = sqrt(c_speed);
405  REAL c_speed2 = c_speed * norme;
406  //
407  //.. Compute the eigenvalues of the Jacobian matrix
408  REAL eig_val1 = scal - c_speed2;
409  REAL eig_val2 = scal;
410  REAL eig_val3 = scal + c_speed2;
411  //
412  //.. Compute the ROE flux
413  //.... In this part many tests upon the eigenvalues
414  //.... are done to simplify calculations
415  //.... Here we use the two formes of the ROE flux :
416  //.... phi(Wl,Wr) = F(Wl) + A-(Wroe)(Wr - Wl)
417  //.... phi(Wl,Wr) = F(Wr) - A+(Wroe)(Wr - Wl)
418  //
419  if(eig_val2 <= 0.0) {
420  p_t = gam1 * (rhoE_t - 0.5 * (rhou_t * rhou_t +
421  rhov_t * rhov_t + rhow_t * rhow_t) * irho_t);
422  ep_t = rhoE_t + p_t;
423  rhouv_t = rhou_t * v_t;
424  rhouw_t = rhou_t * w_t;
425  rhovw_t = rhov_t * w_t;
426  flux_rho = rhou_t * nx + rhov_t * ny + rhow_t * nz;
427  flux_rhou = (rhou_t * u_t + p_t) * nx + rhouv_t * ny + rhouw_t * nz;
428  flux_rhov = rhouv_t * nx + (rhov_t * v_t + p_t) * ny + rhovw_t * nz;
429  flux_rhow = rhouw_t * nx + rhovw_t * ny + (rhow_t * w_t + p_t) * nz;
430  flux_rhoE = ep_t * (u_t * nx + v_t * ny + w_t * nz);
431  //
432  //.... A Entropic modification
433  //
434  p_f = gam1 * (rhoE_f - 0.5 * (rhou_f * rhou_f + rhov_f * rhov_f
435  + rhow_f * rhow_f) * irho_f);
436  lambda_f = u_f * nx + v_f * ny + w_f * nz + norme
437  * sqrt(gam * p_f * irho_f);
438  lambda_t = u_t * nx + v_t * ny + w_t * nz + norme
439  * sqrt(gam * p_t * irho_t);
440  if ((lambda_f < 0.) && (lambda_t > 0.)) {
441  eig_val3 = lambda_t * (eig_val3 - lambda_f) / (lambda_t - lambda_f);
442  }
443  //
444  if (eig_val3 > 0.0) {
445  //.. In this case A+ is obtained by multiplying the last
446  //.. colomne of T-1 with the last row of T with eig_val3 //Cedric
447  delta_rho = rho_t - rho_f; //right - left
448  delta_rhou = rhou_t - rhou_f; //**_t - **_f
449  delta_rhov = rhov_t - rhov_f;
450  delta_rhow = rhow_t - rhow_f;
451  delta_rhoE = rhoE_t - rhoE_f;
452  //
453  scal = scal * inorme;
454  hnx = nx * inorme;
455  hny = ny * inorme;
456  hnz = nz * inorme;
457  usc = 1.0/c_speed;
458  tempo11 = gam1 * usc;
459  //.. Last columne of the matrix T-1
460  a1 = usc;
461  a2 = u_ave * usc + hnx;
462  a3 = v_ave * usc + hny;
463  a4 = w_ave * usc + hnz;
464  a5 = 0.5 * u2pv2pw2 * usc + 2.5 * c_speed + scal;
465  //.. Last row of the matrix T * eig_val3
466  b1 = 0.5 * (0.5 * tempo11 * u2pv2pw2 - scal);
467  b2 = 0.5 * (hnx - tempo11 * u_ave);
468  b3 = 0.5 * (hny - tempo11 * v_ave);
469  b4 = 0.5 * (hnz - tempo11 * w_ave);
470  b5 = 0.5 * tempo11;
471  //
472  alpha1 = b1 * delta_rho;
473  alpha2 = b2 * delta_rhou;
474  alpha3 = b3 * delta_rhov;
475  alpha4 = b4 * delta_rhow;
476  alpha5 = b5 * delta_rhoE;
477  alpha = eig_val3 * (alpha1 + alpha2 + alpha3 + alpha4 + alpha5);
478  //
479  flux_rho -= a1 * alpha;
480  flux_rhou -= a2 * alpha;
481  flux_rhov -= a3 * alpha;
482  flux_rhow -= a4 * alpha;
483  flux_rhoE -= a5 * alpha;
484  }
485  }
486  //
487  if(eig_val2 > 0.0) {
488  p_f = gam1 * (rhoE_f - 0.5 * (rhou_f * rhou_f +
489  rhov_f * rhov_f + rhow_f * rhow_f) * irho_f);
490  ep_f = rhoE_f + p_f;
491  rhouv_f = rhou_f * v_f;
492  rhouw_f = rhou_f * w_f;
493  rhovw_f = rhov_f * w_f;
494  flux_rho = rhou_f * nx + rhov_f * ny + rhow_f * nz;
495  flux_rhou = (rhou_f * u_f + p_f) * nx + rhouv_f * ny + rhouw_f * nz;
496  flux_rhov = rhouv_f * nx + (rhov_f * v_f + p_f) * ny + rhovw_f * nz;
497  flux_rhow = rhouw_f * nx + rhovw_f * ny + (rhow_f * w_f + p_f) * nz;
498  flux_rhoE = ep_f * (u_f * nx + v_f * ny + w_f * nz);
499  //
500  // A Entropic modification
501  //
502  p_t = gam1 * (rhoE_t - 0.5 * (rhou_t * rhou_t +
503  rhov_t * rhov_t + rhow_t * rhow_t) * irho_t);
504  lambda_f = u_f * nx + v_f * ny + w_f * nz - norme
505  * sqrt(gam * p_f * irho_f);
506  lambda_t = u_t * nx + v_t * ny + w_t * nz - norme
507  * sqrt(gam * p_t * irho_t);
508  if ((lambda_f < 0.) && (lambda_t > 0.)) {
509  eig_val1 = lambda_f * (lambda_t - eig_val1) / (lambda_t - lambda_f);
510  }
511  //
512  if (eig_val1 < 0.0) {
513  //.. In this case A+ is obtained by multiplying the first
514  //.. columne of T-1 with the first row of T with eig_val1
515  delta_rho = rho_t - rho_f;
516  delta_rhou = rhou_t - rhou_f;
517  delta_rhov = rhov_t - rhov_f;
518  delta_rhow = rhow_t - rhow_f;
519  delta_rhoE = rhoE_t - rhoE_f;
520  //
521  scal = scal * inorme;
522  hnx = nx * inorme;
523  hny = ny * inorme;
524  hnz = nz * inorme;
525  usc = 1.0/c_speed;
526  tempo11 = gam1 * usc;
527  //.. First colomne of the matrix T-1
528  a1 = usc;
529  a2 = u_ave * usc - hnx;
530  a3 = v_ave * usc - hny;
531  a4 = w_ave * usc - hnz;
532  a5 = 0.5 * u2pv2pw2 * usc + 2.5 * c_speed - scal;
533  //.. First row of the matrix T * eig_val1
534  b1 = 0.5 * (0.5 * tempo11 * u2pv2pw2 + scal);
535  b2 = -0.5 * (hnx + tempo11 * u_ave);
536  b3 = -0.5 * (hny + tempo11 * v_ave);
537  b4 = -0.5 * (hnz + tempo11 * w_ave);
538  b5 = 0.5 * tempo11;
539  //
540  alpha1 = b1 * delta_rho;
541  alpha2 = b2 * delta_rhou;
542  alpha3 = b3 * delta_rhov;
543  alpha4 = b4 * delta_rhow;
544  alpha5 = b5 * delta_rhoE;
545  alpha = eig_val1 * (alpha1 + alpha2 + alpha3 + alpha4 + alpha5);
546  //
547  flux_rho += a1 * alpha;
548  flux_rhou += a2 * alpha;
549  flux_rhov += a3 * alpha;
550  flux_rhow += a4 * alpha;
551  flux_rhoE += a5 * alpha;
552  }
553  }
554 }
555 
556 //left = **_f right = **_t
557 void TPZDiffusionConsLaw::Roe_Flux(REAL rho_f, REAL rhou_f, REAL rhov_f, REAL rhoE_f,
558  REAL rho_t, REAL rhou_t, REAL rhov_t, REAL rhoE_t,
559  REAL nx, REAL ny, REAL gam,
560  REAL &flux_rho, REAL &flux_rhou,REAL &flux_rhov, REAL &flux_rhoE){
561 
562  REAL alpha1,alpha2,alpha3,alpha4,a1,a2,a3,a4,b1,b2,b3,b4,alpha;
563  REAL ep_t, ep_f, p_t, p_f;
564  REAL rhouv_t, rhouv_f;
565  REAL lambda_f, lambda_t;
566  REAL delta_rho, delta_rhou,delta_rhov, delta_rhoE;
567  REAL hnx, hny;
568  REAL tempo11, usc;
569 
570  flux_rho = 0;
571  flux_rhou = 0;
572  flux_rhov = 0;
573  flux_rhoE = 0;
574 
575  REAL gam1 = gam - 1.0;
576  //REAL gam2 = gam * (gam - 1.0);
577  //REAL igam = 1.0 / (gam - 1.0);
578 
579  //
580  //.. Compute the ROE Averages
581  //
582  //.... some useful quantities
583  REAL coef1 = sqrt(rho_f);
584  REAL coef2 = sqrt(rho_t);
585  REAL somme_coef = coef1 + coef2;
586  REAL u_f = rhou_f/rho_f;
587  REAL v_f = rhov_f/rho_f;
588  REAL h_f = (gam * rhoE_f/rho_f) - (gam1 / 2.0) * (u_f * u_f + v_f * v_f);
589  REAL u_t = rhou_t/rho_t;
590  REAL v_t = rhov_t/rho_t;
591  REAL h_t = (gam * rhoE_t/rho_t) - (gam1 / 2.0) * (u_t * u_t + v_t * v_t);
592 
593  //.... averages
594  //REAL rho_ave = coef1 * coef2;
595  REAL u_ave = (coef1 * u_f + coef2 * u_t) / somme_coef;
596  REAL v_ave = (coef1 * v_f + coef2 * v_t) / somme_coef;
597  REAL h_ave = (coef1 * h_f + coef2 * h_t) / somme_coef;
598  //
599  //.. Compute Speed of sound
600  REAL scal = u_ave * nx + v_ave * ny;
601  REAL norme = sqrt(nx * nx + ny * ny);
602  REAL u2pv2 = u_ave * u_ave + v_ave * v_ave;
603  REAL c_speed = gam1 * (h_ave - 0.5 * u2pv2);
604  if(c_speed < 1e-6) c_speed = 1e-6; // avoid division by 0 if critical
605  c_speed = sqrt(c_speed);
606  REAL c_speed2 = c_speed * norme;
607  //
608  //.. Compute the eigenvalues of the Jacobian matrix
609  REAL eig_val1 = scal - c_speed2;
610  REAL eig_val2 = scal;
611  REAL eig_val3 = scal + c_speed2;
612  //
613  //.. Compute the ROE flux
614  //.... In this part many tests upon the eigenvalues
615  //.... are done to simplify calculations
616  //.... Here we use the two formes of the ROE flux :
617  //.... phi(Wl,Wr) = F(Wl) + A-(Wroe)(Wr - Wl)
618  //.... phi(Wl,Wr) = F(Wr) - A+(Wroe)(Wr - Wl)
619  //
620  if(eig_val2 <= 0.0) {
621  p_t = gam1 * (rhoE_t - 0.5 * (rhou_t * rhou_t + rhov_t * rhov_t) / rho_t);
622  ep_t = rhoE_t + p_t;
623  rhouv_t = rhou_t * v_t;
624  flux_rho = rhou_t * nx + rhov_t * ny;
625  flux_rhou = (rhou_t * u_t + p_t) * nx + rhouv_t * ny;
626  flux_rhov = rhouv_t * nx + (rhov_t * v_t + p_t) * ny;
627  flux_rhoE = ep_t * (u_t * nx + v_t * ny);
628  //
629  //.... A Entropic modification
630  //
631  p_f = gam1 * (rhoE_f - 0.5 * (rhou_f * rhou_f + rhov_f * rhov_f) / rho_f);
632  lambda_f = u_f * nx + v_f * ny + norme * sqrt(gam * p_f / rho_f);
633  lambda_t = u_t * nx + v_t * ny + norme
634  * sqrt(gam * p_t / rho_t);
635  if ((lambda_f < 0.) && (lambda_t > 0.)) {
636  eig_val3 = lambda_t * (eig_val3 - lambda_f) / (lambda_t - lambda_f);
637  }
638  //
639  if (eig_val3 > 0.0) {
640  //.. In this case A+ is obtained by multiplying the last
641  //.. colomne of T-1 with the last row of T with eig_val3
642  delta_rho = rho_t - rho_f;
643  delta_rhou = rhou_t - rhou_f;
644  delta_rhov = rhov_t - rhov_f;
645  delta_rhoE = rhoE_t - rhoE_f;
646  //
647  scal = scal / norme;
648  hnx = nx / norme;
649  hny = ny / norme;
650  usc = 1.0/c_speed;
651  tempo11 = gam1 * usc;
652  //.. Last columne of the matrix T-1
653  a1 = usc;
654  a2 = u_ave * usc + hnx;
655  a3 = v_ave * usc + hny;
656  a4 = 0.5 * u2pv2 * usc + 2.5 * c_speed + scal;
657  //.. Last row of the matrix T * eig_val3
658  b1 = 0.5 * eig_val3 * (0.5 * tempo11 * u2pv2 - scal);
659  b2 = 0.5 * eig_val3 * (hnx - tempo11 * u_ave);
660  b3 = 0.5 * eig_val3 * (hny - tempo11 * v_ave);
661  b4 = 0.5 * eig_val3 * tempo11;
662  //
663  alpha1 = a1 * b1 * delta_rho;
664  alpha2 = a1 * b2 * delta_rhou;
665  alpha3 = a1 * b3 * delta_rhov;
666  alpha4 = a1 * b4 * delta_rhoE;
667  alpha = alpha1 + alpha2 + alpha3 + alpha4;
668  //
669  flux_rho -= alpha;
670  flux_rhou -= a2 * b1 * delta_rho + a2 * b2 * delta_rhou +
671  a2 * b3 * delta_rhov + a2 * b4 * delta_rhoE;
672  flux_rhov -= a3 * b1 * delta_rho + a3 * b2 * delta_rhou +
673  a3 * b3 * delta_rhov + a3 * b4 * delta_rhoE;
674  flux_rhoE -= a4 * b1 * delta_rho + a4 * b2 * delta_rhou +
675  a4 * b3 * delta_rhov + a4 * b4 * delta_rhoE;
676  }
677  }
678  //
679  if(eig_val2 > 0.0) {
680  p_f = gam1 * (rhoE_f - REAL(0.5) * (rhou_f * rhou_f +
681  rhov_f * rhov_f) / rho_f);
682  ep_f = rhoE_f + p_f;
683  rhouv_f = rhou_f * v_f;
684  flux_rho = rhou_f * nx + rhov_f * ny;
685  flux_rhou = (rhou_f * u_f + p_f) * nx + rhouv_f * ny;
686  flux_rhov = rhouv_f * nx + (rhov_f * v_f + p_f) * ny;
687  flux_rhoE = ep_f * (u_f * nx + v_f * ny);
688  //
689  // A Entropic modification
690  //
691  p_t = gam1 * (rhoE_t - REAL(0.5) * (rhou_t * rhou_t +
692  rhov_t * rhov_t) / rho_t);
693  lambda_f = u_f * nx + v_f * ny - norme * sqrt(gam * p_f / rho_f);
694  lambda_t = u_t * nx + v_t * ny - norme * sqrt(gam * p_t / rho_t);
695  if ((lambda_f < 0.) && (lambda_t > 0.)) {
696  eig_val1 = lambda_f * (lambda_t - eig_val1) / (lambda_t - lambda_f);
697  }
698  //
699  if (eig_val1 < 0.0) {
700  //.. In this case A+ is obtained by multiplying the first
701  //.. columne of T-1 with the first row of T with eig_val1
702  delta_rho = rho_t - rho_f;
703  delta_rhou = rhou_t - rhou_f;
704  delta_rhov = rhov_t - rhov_f;
705  delta_rhoE = rhoE_t - rhoE_f;
706  //
707  scal = scal / norme;
708  hnx = nx / norme;
709  hny = ny / norme;
710  usc = 1.0/c_speed;
711  tempo11 = gam1 * usc;
712  //.. First colomne of the matrix T-1
713  a1 = usc;
714  a2 = u_ave * usc - hnx;
715  a3 = v_ave * usc - hny;
716  a4 = 0.5 * u2pv2 * usc + 2.5 * c_speed - scal;
717  //.. First row of the matrix T * eig_val1
718  b1 = 0.5 * eig_val1 * (0.5 * tempo11 * u2pv2 + scal);
719  b2 = -0.5 * eig_val1 * (hnx + tempo11 * u_ave);
720  b3 = -0.5 * eig_val1 * (hny + tempo11 * v_ave);
721  b4 = 0.5 * eig_val1 * tempo11;
722  //
723  alpha1 = a1 * b1 * delta_rho;
724  alpha2 = a1 * b2 * delta_rhou;
725  alpha3 = a1 * b3 * delta_rhov;
726  alpha4 = a1 * b4 * delta_rhoE;
727  alpha = alpha1 + alpha2 + alpha3 + alpha4;
728  //
729  flux_rho += alpha;
730  flux_rhou += a2 * b1 * delta_rho + a2 * b2 * delta_rhou +
731  a2 * b3 * delta_rhov + a2 * b4 * delta_rhoE;
732  flux_rhov += a3 * b1 * delta_rho + a3 * b2 * delta_rhou +
733  a3 * b3 * delta_rhov + a3 * b4 * delta_rhoE;
734  flux_rhoE += a4 * b1 * delta_rho + a4 * b2 * delta_rhou +
735  a4 * b3 * delta_rhov + a4 * b4 * delta_rhoE;
736  }
737  }
738 }
TPZDiffusionConsLaw::fArtificialDiffusion
static std::string fArtificialDiffusion
Term that adds stability to the numerical method of approach: SUPG, LS, Bornhaus. ...
Definition: TPZDiffusionConsLaw.h:35
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZDiffusionConsLaw::fA
TPZFMatrix< REAL > fA
Matrix computation derivatives of fluxes: dFx/dx, dFy/dy, dFz/dz.
Definition: TPZDiffusionConsLaw.h:39
TPZDiffusionConsLaw::SUPG
void SUPG(TPZFMatrix< REAL > &Tx, TPZFMatrix< REAL > &Ty, TPZFMatrix< REAL > &Tz)
Definition: TPZDiffusionConsLaw.cpp:120
TPZDiffusionConsLaw::Divergence
void Divergence(TPZVec< REAL > &dphi, TPZFMatrix< REAL > &diverg)
Definition: TPZDiffusionConsLaw.cpp:64
TPZDiffusionConsLaw::PointOperator
void PointOperator(TPZVec< REAL > &dphi, TPZFMatrix< REAL > &diff_term)
Operation product point in the diffusion term.
Definition: TPZDiffusionConsLaw.cpp:78
pzvec.h
Templated vector implementation.
degree
void degree(int root, int adj_num, int adj_row[], int adj[], int mask[], int deg[], int *iccsze, int ls[], int node_num)
Definition: rcm.cpp:875
TPZDiffusionConsLaw::fDelta
static REAL fDelta
Coefficient of size ratio of the element in the diffusion term.
Definition: TPZDiffusionConsLaw.h:29
TPZDiffusionConsLaw::Tau
void Tau(TPZFMatrix< REAL > &Tx, TPZFMatrix< REAL > &Ty, TPZFMatrix< REAL > &Tz)
Definition: TPZDiffusionConsLaw.cpp:103
TPZDiffusionConsLaw::fB
TPZFMatrix< REAL > fB
Definition: TPZDiffusionConsLaw.h:39
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
TPZDiffusionConsLaw::JacobFlux
static void JacobFlux(TPZVec< REAL > U, TPZFMatrix< REAL > &A, TPZFMatrix< REAL > &B, TPZFMatrix< REAL > &C)
Jacobiano of the tensor flux of Euler.
Definition: TPZDiffusionConsLaw.cpp:135
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
TPZDiffusionConsLaw::fC
TPZFMatrix< REAL > fC
Definition: TPZDiffusionConsLaw.h:39
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZDiffusionConsLaw.h
Contains the TPZDiffusionConsLaw class which implements a Euler equation where is introduced a diffus...
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZCompEl::GetgOrder
static int GetgOrder()
Set the default value of the interpolation order.
Definition: pzcompel.h:830
TPZCompElDisc.h
Contains declaration of TPZCompelDisc class which implements a computational element for discontinuou...
TPZDiffusionConsLaw::fDimension
int fDimension
Problem dimension.
Definition: TPZDiffusionConsLaw.h:42
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
gamma
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
Definition: tpzgaussrule.cpp:567
TPZDiffusionConsLaw::Bornhaus
void Bornhaus(TPZFMatrix< REAL > &Tx, TPZFMatrix< REAL > &Ty, TPZFMatrix< REAL > &Tz)
Definition: TPZDiffusionConsLaw.cpp:131
TPZDiffusionConsLaw::~TPZDiffusionConsLaw
~TPZDiffusionConsLaw()
Destructor.
Definition: TPZDiffusionConsLaw.cpp:17
TPZDiffusionConsLaw::LS
void LS(TPZFMatrix< REAL > &Tx, TPZFMatrix< REAL > &Ty, TPZFMatrix< REAL > &Tz)
Definition: TPZDiffusionConsLaw.cpp:124
TPZDiffusionConsLaw::GradientOfTheFlow
void GradientOfTheFlow(TPZFMatrix< REAL > &DF1, TPZFMatrix< REAL > &DF2, TPZFMatrix< REAL > &DF3)
Definition: TPZDiffusionConsLaw.cpp:39
TPZDiffusionConsLaw::fCFL
static REAL fCFL
Parameter that it limits the condition of stability of the numerical approach.
Definition: TPZDiffusionConsLaw.h:32
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
pzreal.h
Contains the declaration of TPZFlopCounter class and TPZCounter struct.
TPZFMatrix::Transpose
void Transpose(TPZMatrix< TVar > *const T) const override
It makes *T the transpose of current matrix.
Definition: pzfmatrix.cpp:1073
TPZDiffusionConsLaw::TPZDiffusionConsLaw
TPZDiffusionConsLaw()
Default constructor.
Definition: TPZDiffusionConsLaw.cpp:19
TPZDiffusionConsLaw::fGamma
static REAL fGamma
Ratio between specific heat is constant and the specific heat the constant volume of a polytropic gas...
Definition: TPZDiffusionConsLaw.h:26
TPZDiffusionConsLaw::Roe_Flux
static void Roe_Flux(REAL rho_f, REAL rhou_f, REAL rhov_f, REAL rhow_f, REAL rhoE_f, REAL rho_t, REAL rhou_t, REAL rhov_t, REAL rhow_t, REAL rhoE_t, REAL nx, REAL ny, REAL nz, REAL gam, REAL &flux_rho, REAL &flux_rhou, REAL &flux_rhov, REAL &flux_rhow, REAL &flux_rhoE)
Flux of Roe (MOUSE program)
Definition: TPZDiffusionConsLaw.cpp:347
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13