NeoPZ
pzeulerconslaw.cpp
Go to the documentation of this file.
1 
6 #include "pzeulerconslaw.h"
7 #include "pzartdiff.h"
8 #include "pzelmat.h"
9 #include "pzbndcond.h"
10 #include "pzmatrix.h"
11 #include "pzfmatrix.h"
12 #include "pzerror.h"
13 #include "pzreal.h"
14 #include <math.h>
15 #include "pzstring.h"
16 #include "TPZSavable.h"
17 #include "pzerror.h"
18 
19 #include <cmath>
20 using namespace std;
21 
22 #include "pzlog.h"
23 
24 #ifdef LOG4CXX
25 
26 LoggerPtr fluxroe(Logger::getLogger("pz.fluxroe"));
27 LoggerPtr fluxappr(Logger::getLogger("pz.fluxappr"));
28 
29 #endif
30 
31 
32 
33 TPZEulerConsLaw::~TPZEulerConsLaw(){
34 
35 }
36 
37 TPZEulerConsLaw::TPZEulerConsLaw(int nummat,REAL timeStep,
38  REAL gamma,int dim,
39  TPZArtDiffType artdiff) :
40 TPZRegisterClassId(&TPZEulerConsLaw::ClassId),
41 TPZConservationLaw(nummat,timeStep,dim),
42 fArtDiff(artdiff, gamma),
43 fDiff(Explicit_TD),
44 fConvVol(Explicit_TD),
45 fConvFace(Explicit_TD)
46 {
47  fGamma = gamma;
48 }
49 
50 TPZEulerConsLaw::TPZEulerConsLaw() :
51 TPZRegisterClassId(&TPZEulerConsLaw::ClassId),
52 TPZConservationLaw(-1, 0, 3),
53 fArtDiff(LeastSquares_AD, 1.4),
54 fDiff(Explicit_TD),
55 fConvVol(Explicit_TD),
56 fConvFace(Explicit_TD)
57 {
58  fGamma = 1.4;
59 }
60 
61 void TPZEulerConsLaw::SetTimeDiscr(TPZTimeDiscr Diff, TPZTimeDiscr ConvVol, TPZTimeDiscr ConvFace)
62 {
63  fDiff = Diff;
64  fConvVol = ConvVol;
65  fConvFace = ConvFace;
66 }
67 
68 REAL TPZEulerConsLaw::OptimalCFL(int degree)
69 {
70  return 1./((2.0*(REAL)degree) + 1.0);
71 }
72 
73 REAL TPZEulerConsLaw::SetTimeStep(REAL maxveloc,REAL deltax,int degree)
74 {
75  REAL CFL = fCFL;
76  // Notice that fCFL remains 0, so that optimal CFL will
77  // be computed unless CFL is redefined.
78  if(CFL < 0.0) CFL = OptimalCFL(degree);
79 
80  REAL deltaT = CFL*deltax/maxveloc;
81  //cout << "TPZCompMesh::Delta Time : " << deltaT << endl;
82  TPZConservationLaw::SetTimeStep(deltaT);
83 
84  return deltaT;
85 }
86 
87 int TPZEulerConsLaw::NStateVariables(int dim) {
88  return (2 + dim);//U = (rho, rhou, rhov, rhow, rhoe)
89 }
90 
91 int TPZEulerConsLaw::NStateVariables() const {
92  return NStateVariables(Dimension());//U = (rho, rhou, rhov, rhow, rhoe)
93 }
94 
95 STATE TPZEulerConsLaw::Pressure(TPZVec<STATE> &U)
96 {
97  STATE press;
98  TPZEulerConsLaw::Pressure(fGamma, fDim, press, U);
99  return press;
100 }
101 
102 void TPZEulerConsLaw::Print(std::ostream &out) {
103 
104  TPZMaterial::Print(out);
105 
106  TPZConservationLaw::Print(out);
107  out << "Artificial Diffusion: " <<
108  fArtDiff.DiffusionName().Str() << std::endl;
109  out << "Number of State Variables: " << NStateVariables() << std::endl;
110  out << "Number of Fluxes: " << NFluxes() << std::endl;
111 
112  switch(fDiff)
113  {
114  case(Explicit_TD):
115  out << "Explicit Diffusive term\n";
116  break;
117  case(ApproxImplicit_TD):
118  out << "ApproxImplicit Diffusive term\n";
119  break;
120  case(Implicit_TD):
121  out << "Implicit Diffusive term\n";
122  break;
123  default:
124  out << "No Diffusive term\n";
125  }
126 
127  switch(fConvVol)
128  {
129  case(Explicit_TD):
130  out << "Explicit Volume Convective term\n";
131  break;
132  case(Implicit_TD):
133  out << "Implicit Volume Convective term\n";
134  break;
135  default:
136  out << "No Volume Convective term\n";
137  }
138 
139  switch(fConvFace)
140  {
141  case(Explicit_TD):
142  out << "Explicit Face Convective term\n";
143  break;
144  case(Implicit_TD):
145  out << "Implicit Face Convective term\n";
146  break;
147  case(ApproxImplicit_TD):
148  out << "Approximate Implicit Face Convective term\n";
149  break;
150  default:
151  std::cout << "No Face Convective term\n";
152  }
153 
154 
155 }
156 
157 int TPZEulerConsLaw::VariableIndex(const std::string &name) {
158  if( !strcmp(name.c_str(),"density") ) return 1;//rho
159  if( !strcmp(name.c_str(),"velocity") ) return 2;//(u,v,w)
160  if( !strcmp(name.c_str(),"energy") ) return 3;//E
161  if( !strcmp(name.c_str(),"pressure") ) return 4;//p
162  if( !strcmp(name.c_str(),"solution") ) return 5;//(ro,u,v,w,E)
163  if( !strcmp(name.c_str(),"normvelocity") ) return 6;//sqrt(u+v+w)
164  if( !strcmp(name.c_str(),"Mach") ) return 7;//sqrt(u+v+w)/c
165  std::cout << "TPZEulerConsLaw::VariableIndex not defined\n";
166  return TPZMaterial::VariableIndex(name);
167 }
168 
169 int TPZEulerConsLaw::NSolutionVariables(int var){
170 
171  if(var == 1 || var == 3 || var == 4 || var == 6 || var == 7) return 1;
172  if(var == 2) return Dimension();
173  if(var == 5) return NStateVariables();
174 
175  std::cout << "TPZEulerConsLaw::NSolutionVariables not defined\n";
176  return 0;
177 }
178 
179 REAL TPZEulerConsLaw::DeltaX(REAL detJac)
180 {
181  return REAL(2.) * pow(fabs(detJac),(REAL)(1./fDim));
182 }
183 
184 REAL TPZEulerConsLaw::Det(TPZFMatrix<REAL> & Mat)
185 {
186  switch(Mat.Rows())
187  {
188  case 1:
189  return Mat(0,0);
190  // break;
191  case 2:
192  return Mat(0,0) * Mat(1,1) -
193  Mat(1,0) * Mat(0,1);
194  // break;
195  case 3:
196  return Mat(0,0) * Mat(1,1) * Mat(2,2) +
197  Mat(1,0) * Mat(2,1) * Mat(0,2) +
198  Mat(2,0) * Mat(0,1) * Mat(1,2) -
199  Mat(0,2) * Mat(1,1) * Mat(2,0) -
200  Mat(1,2) * Mat(2,1) * Mat(0,0) -
201  Mat(2,2) * Mat(0,1) * Mat(1,0);
202  // break;
203  default:
204  PZError << "TPZEulerConsLaw::Det error: unhandled matrix size: " <<
205  Mat.Rows() << std::endl;
206  }
207 
208  return 0.;
209 }
210 
211 int TPZEulerConsLaw::NFluxes()
212 {
213  return Dimension();
214 }
215 
216 //-----------------Solutions
217 
218 void TPZEulerConsLaw::Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,TPZVec<STATE> &Solout){
219 
220  if(fabs(Sol[0]) < 1.e-10) {
221  PZError << "\nTPZEulerConsLaw::Solution: Density almost null\n"
222  << "Density = " << Sol[0] << std::endl;
223  }
224 
225  if(var == 1) {
226  Solout.Resize(1);
227  Solout[0] = Sol[0];//density
228  return;
229  } else if(var == 2) {
230  int dim = Dimension();
231  Solout.Resize(dim);
232  for(int i=0;i<dim;i++) Solout[i] = Sol[i+1]/Sol[0];//velocity vector
233  return;
234  } else if(var == 3) {
235  Solout.Resize(1);
236  int pos = Dimension() + 1;
237  Solout[0] = Sol[pos];//energy
238  return;
239  } else if(var == 4) {
240  Solout.Resize(1);
241  Solout[0] = Pressure(Sol);//pressure
242  return;
243  } else if(var == 5) {
244  int nstate = NStateVariables();
245  Solout.Resize(nstate);
246  for(int i=0;i<nstate;i++) Solout[i] = Sol[i];//(ro,ro*u,ro*v,ro*w,E)
247  return;
248  } else if(var == 6) {
249  int nstate = NStateVariables();
250  Solout.Resize(1);
251  REAL ro2 = Sol[0]*Sol[0];
252  REAL veloc = 0.0;
253  for(int i=1;i<nstate-1;i++) veloc += Sol[i]*Sol[i];//velocity vector
254  Solout[0] = sqrt(veloc/ro2);
255  return;
256  } else if(var == 7) {
257  // int nstate = NStateVariables();
258  Solout.Resize(1);
259  STATE cspeed;
260  STATE us;
261  TPZEulerConsLaw::cSpeed(Sol, fGamma, cspeed);
262  TPZEulerConsLaw::uRes(Sol, us);
263  Solout[0] = us / cspeed;
264  return;
265  } else {
266  //cout << "TPZEulerConsLaw::Solution variable in the base class\n";
267  TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
268  }
269 }
270 
271 
272 void TPZEulerConsLaw::SetDelta(REAL delta)
273 {
274  fArtDiff.SetDelta(delta);
275 }
276 
277 //------------------Differentiable variables setup
278 
279 #ifdef _AUTODIFF
280 
281 void TPZEulerConsLaw::PrepareFAD(
282  TPZVec<STATE> & sol, TPZFMatrix<STATE> & dsol,
283  TPZFMatrix<REAL> & phi, TPZFMatrix<REAL> & dphi,
284  TPZVec<FADREAL> & FADsol,
285  TPZVec<FADREAL> & FADdsol)
286 {
287  int nState = NStateVariables();
288  int nShape = phi.Rows();
289  int i_state, i_shape, k;
290  int nDer = nState * nShape;
291 
292  // initializing the differentiable variables
293  FADREAL defaultFAD(nDer, REAL(0.), REAL(0.));
294  if(defaultFAD.dx(0)==1.)PZError << "\nError: FAD doesn't have default constructor for parameters: (number of derivatives, default value, default derivative value) !";
295  FADsol.Resize(nState);
296  FADsol.Fill(defaultFAD);
297 
298  FADdsol.Resize(nState * fDim);
299  FADdsol.Fill(defaultFAD);
300 
301  // copying the solution and spatial derivative values
302  for(i_state = 0; i_state < nState; i_state++)
303  {
304  FADsol[i_state].val() = sol[i_state];
305  for(k = 0; k < fDim; k ++)
306  FADdsol[i_state * fDim + k].val() = dsol(k,i_state);
307  }
308 
309  // preparing the coefficient derivatives
310  for(i_shape = 0; i_shape < nShape; i_shape++)
311  for(i_state = 0; i_state < nState; i_state++)
312  {
313  int index = i_state + i_shape * nState;
314  FADsol[i_state].fastAccessDx(index)=phi(i_shape,0);
315  for(k = 0; k < fDim; k++)
316  FADdsol[i_state * fDim + k].fastAccessDx(index)=dphi(k,i_shape);
317  }
318 }
319 
320 void TPZEulerConsLaw::PrepareInterfaceFAD(
321  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
322  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
323  TPZVec<FADREAL> & FADsolL,
324  TPZVec<FADREAL> & FADsolR)
325 {
326  int nState = NStateVariables();
327  int nShapeL = phiL.Rows();
328  int nShapeR = phiR.Rows();
329  int i_state, i_shape;
330  int nDerL = nState * nShapeL;
331  int nDerR = nState * nShapeR;
332 
333  // initializing the differentiable variables
334  FADREAL defaultFAD(nDerL + nDerR, REAL(0.), REAL(0.));
335  if(defaultFAD.dx(0)==1.)PZError << "\nError: FAD doesn't have default constructor for parameters: (number of derivatives, default value, default derivative value) !";
336  FADsolL.Resize(nState);
337  FADsolL.Fill(defaultFAD);
338 
339  FADsolR.Resize(nState);
340  FADsolR.Fill(defaultFAD);
341 
342  // copying the solution and spatial derivatives values
343  for(i_state = 0; i_state < nState; i_state++)
344  {
345  FADsolL[i_state].val() = solL[i_state];
346  FADsolR[i_state].val() = solR[i_state];
347  }
348 
349  // preparing the coefficient derivatives
350  for(i_shape = 0; i_shape < nShapeL; i_shape++)
351  for(i_state = 0; i_state < nState; i_state++)
352  {
353  int index = i_state + i_shape * nState;
354  FADsolL[i_state].fastAccessDx(index)=phiL(i_shape,0);
355  //FADsolR[i_state].fastAccessDx(index)=phiL(i_shape,0);
356  }
357  for(i_shape = 0/*nShapeL*/; i_shape < /*nShapeL +*/ nShapeR; i_shape++)
358  for(i_state = 0; i_state < nState; i_state++)
359  {
360  int index = i_state + (i_shape + nShapeL) * nState;
361  //FADsolL[i_state].fastAccessDx(index)=phiR(i_shape,0);
362  FADsolR[i_state].fastAccessDx(index)=phiR(i_shape,0);
363  }
364 }
365 
366 template <class T>
367 void TPZEulerConsLaw::PrepareFastestInterfaceFAD(
368  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
369  TPZVec<T> & FADsolL,
370  TPZVec<T> & FADsolR)
371 {
372  int nState = solL.NElements();
373  int nVars = nState * 2;
374 
375  FADsolL.Resize(nState);
376  FADsolR.Resize(nState);
377 
378  for(int i = 0; i < nState; i++)
379  {
380  FADsolL[i] = solL[i];
381  FADsolL[i].diff(i, nVars);
382 
383  FADsolR[i] = solR[i];
384  FADsolR[i].diff(i + nState, nVars);
385  }
386 }
387 
388 #endif
389 
390 //----------------Contributions
391 
392 void TPZEulerConsLaw::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
393 {
394  int numbersol = data.sol.size();
395  if (numbersol != 1) {
396  DebugStop();
397  }
398 
399  // initial guesses for sol
400  // fForcingFunction is null at iterations > 0
401  if(fForcingFunction)
402  {
403  TPZVec<STATE> res;
404  int i, nState = NStateVariables();
405  fForcingFunction->Execute(data.x, res);
406  for(i = 0; i < nState; i++)
407  data.sol[0][i] = res[i];
408  }
409 
410  if(fContributionTime == Last_CT)
411  {
412  ContributeLast(data.x, data.jacinv,
413  data.sol[0], data.dsol[0],
414  weight,
415  data.phi, data.dphix,
416  ef);
417  return;
418  }
419 
420  if(fContributionTime == Advanced_CT)
421  {
422  ContributeAdv(data.x, data.jacinv,
423  data.sol[0], data.dsol[0],
424  weight,
425  data.phi, data.dphix,
426  ek, ef);
427  return;
428  }
429 
430  PZError << "TPZEulerConsLaw::Contribute> Unhandled Contribution Time";
431 }
432 
433 void TPZEulerConsLaw::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ef)
434 {
435  int numbersol = data.sol.size();
436  if (numbersol != 1) {
437  DebugStop();
438  }
439 
440  // initial guesses for sol
441  // fForcingFunction is null at iterations > 0
442  if(fForcingFunction)
443  {
444  TPZVec<STATE> res;
445  int i, nState = NStateVariables();
446  fForcingFunction->Execute(data.x, res);
447  for(i = 0; i < nState; i++)
448  data.sol[i] = res[i];
449  }
450 
451  if(fContributionTime == Last_CT)
452  {
453  ContributeLast(data.x, data.jacinv,
454  data.sol[0], data.dsol[0],
455  weight,
456  data.phi, data.dphix,
457  ef);
458  return;
459  }
460 
461  if(fContributionTime == Advanced_CT)
462  {
463  ContributeAdv(data.x, data.jacinv,
464  data.sol[0], data.dsol[0],
465  weight,
466  data.phi, data.dphix,
467  ef);
468  return;
469  }
470 
471  PZError << "TPZEulerConsLaw::Contribute> Unhandled Contribution Time";
472 }
473 
474 void TPZEulerConsLaw::ContributeLast(TPZVec<REAL> &x,TPZFMatrix<REAL> &jacinv,
475  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
476  REAL weight,
477  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
478  TPZFMatrix<STATE> &ef)
479 {
480  // contributing the explicit parcell of the residual to the
481  // rhs.
482 
483  // the parcell T2 is always explicit.
484  if(fResidualType == Residual_RT)
485  {
486  ContributeExplT2(x,sol,weight,phi,ef);
487  }
488 
489  // contributing volume-based quantities
490  // diffusive term
491  if (fDiff == Explicit_TD)
492  ContributeExplDiff(x, jacinv, sol,dsol,weight, phi, dphi, ef);
493 
494  // Volume Convective term
495  if (fConvVol == Explicit_TD)
496  ContributeExplConvVol(x, sol, weight, phi, dphi, ef);
497 
498 
499 }
500 
501 
502 void TPZEulerConsLaw::ContributeAdv(TPZVec<REAL> &x,TPZFMatrix<REAL> &jacinv,
503  TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,
504  REAL weight,
505  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
506  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
507 {
508  // contributing the implicit parcell of the residual to the
509  // rhs.
510 
511 
512  if(fResidualType == Residual_RT)
513  {
514  // the parcell T1 is always implicit.
515  ContributeImplT1(x,sol,dsol,weight, phi,dphi,ek,ef);
516  }
517 
518  // contributing volume-based quantities
519  // diffusive term
520  if (fDiff == Implicit_TD)
521  {
522  // if diffusive term is implicit
523  // then the FAD classes must be initialized
524 #ifdef _AUTODIFF
525 #ifdef FASTEST_IMPLICIT
526  ContributeFastestImplDiff(fDim, x, jacinv, sol, dsol,
527  phi, dphi, weight, ek, ef);
528 #else
529  TPZVec<FADREAL> FADsol, FADdsol;
530  PrepareFAD(sol, dsol, phi, dphi, FADsol, FADdsol);
531  ContributeImplDiff(x, jacinv, FADsol,FADdsol, weight, ek, ef);
532 #endif
533 #else
534  std::cout << "TPZEulerConsLaw::Contribute> Implicit diffusive contribution: _AUTODIFF directive not configured -> Using an approximation to the tgMatrix";
535  ContributeApproxImplDiff(x, jacinv, sol,dsol,weight,phi,dphi,ek,ef);
536 #endif
537  }else
538  {
539  if (fDiff == ApproxImplicit_TD)
540  ContributeApproxImplDiff(x, jacinv, sol,dsol,weight,phi,dphi,ek,ef);
541  }
542 
543  // Volume convective term
544  if (fConvVol == Implicit_TD)
545  ContributeImplConvVol(x,sol,dsol,weight,phi,dphi,ek,ef);
546  /*}else
547  {
548  // Flux_RT -> contribution only to the residual vector
549  if (fDiff == Implicit_TD)
550  ContributeExplDiff(x, jacinv, sol,dsol,weight, phi, dphi, ef);
551  if (fConvVol == Implicit_TD)
552  ContributeExplConvVol(x, sol, weight, phi, dphi, ef);
553  }*/
554 }
555 
556 void TPZEulerConsLaw::ContributeAdv(TPZVec<REAL> &x,TPZFMatrix<REAL> &jacinv,
557  TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,
558  REAL weight,
559  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
560  TPZFMatrix<STATE> &ef)
561 {
562  // contributing the implicit parcell of the residual to the
563  // rhs.
564 
565  if(fResidualType == Residual_RT)
566  {
567  ContributeExplT1(x,sol,dsol,weight, phi,dphi,ef);
568  }
569  // Flux_RT -> contribution only to the residual vector
570  if (fDiff == Implicit_TD || fDiff == ApproxImplicit_TD)
571  ContributeExplDiff(x, jacinv, sol,dsol,weight, phi, dphi, ef);
572  if (fConvVol == Implicit_TD)
573  ContributeExplConvVol(x, sol, weight, phi, dphi, ef);
574 
575 }
576 
577 void TPZEulerConsLaw::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright,
578  REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
579 {
580 #ifdef LOG4CXX
581  if(fluxroe->isDebugEnabled()){
582  std::stringstream sout;
583  sout << "solL " << dataleft.sol << std::endl << "solR " << dataright.sol << std::endl;
584  LOGPZ_DEBUG(fluxroe,sout.str().c_str());
585  LOGPZ_DEBUG(fluxappr,sout.str().c_str());
586  }
587 #endif
588  int numbersol = dataleft.sol.size();
589  if (numbersol != 1) {
590  DebugStop();
591  }
592 
593  // contributing face-based quantities
594  if (fConvFace == Implicit_TD && fContributionTime == Advanced_CT)
595  {
596  // if face contribution is implicit,
597  // then the FAD classes must be initialized
598 #ifdef _AUTODIFF
599 #ifdef FASTEST_IMPLICIT
600  ContributeFastestImplConvFace(fDim, data.x, dataleft.sol[0], dataright.sol[0],
601  weight, data.normal, dataleft.phi, dataright.phi, ek, ef);
602 #else
603  TPZVec<FADREAL> FADsolL, FADsolR;
604  PrepareInterfaceFAD(dataleft.sol[0], dataright.sol[0], dataleft.phi, dataright.phi, FADsolL, FADsolR);
605  ContributeImplConvFace(data.x,FADsolL,FADsolR, weight, data.normal, dataleft.phi, dataright.phi, ek, ef);
606 #ifdef LOG4CXX
607  if(fluxroe->isDebugEnabled()){
608  std::stringstream sout;
609  ek.Print("computed tangent matrix",sout);
610  ef.Print("computed rhs",sout);
611  LOGPZ_DEBUG(fluxroe,sout.str().c_str());
612  }
613 #endif
614 #endif
615 #else
616  // forcing explicit contribution and issueing an warning
617  std::cout << "TPZEulerConsLaw::ContributeInterface> Implicit face convective contribution: _AUTODIFF directive not configured";
618  // ContributeApproxImplConvFace(x,data.HSize,FADsolL,FADsolR, weight, normal, phiL, phiR, ek, ef);
619  ContributeExplConvFace(data.x,dataleft.sol[0],dataright.sol[0],weight,data.normal,dataleft.phi,dataright.phi,ef);
620 #endif
621  }
622  else if(fConvFace == ApproxImplicit_TD && fContributionTime == Advanced_CT)
623  {
624 #ifdef _AUTODIFF
625  TPZVec<FADREAL> FADsolL, FADsolR;
626  PrepareInterfaceFAD(dataleft.sol[0], dataright.sol[0], dataleft.phi, dataright.phi, FADsolL, FADsolR);
627  ContributeApproxImplConvFace(data.x,data.HSize,FADsolL,FADsolR, weight, data.normal, dataleft.phi, dataright.phi, ek, ef);
628 #endif
629 #ifdef LOG4CXX
630  if(fluxappr->isDebugEnabled()){
631  std::stringstream sout;
632  ek.Print("computed tangent matrix",sout);
633  ef.Print("computed rhs",sout);
634  LOGPZ_DEBUG(fluxappr,sout.str().c_str());
635  }
636 #endif
637  }
638 
639  if(fConvFace == Explicit_TD && fContributionTime == Last_CT)
640  {
641  ContributeExplConvFace(data.x,dataleft.sol[0],dataright.sol[0],weight,data.normal,dataleft.phi,dataright.phi,ef);
642  }
643 }
644 
645 
646 void TPZEulerConsLaw::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright,
647  REAL weight, TPZFMatrix<STATE> &ef)
648 {
649 
650  // contributing face-based quantities
651  if ((fConvFace == Implicit_TD && fContributionTime == Advanced_CT)
652  ||
653  (fConvFace == ApproxImplicit_TD && fContributionTime == Advanced_CT)
654  ||
655  (fConvFace == Explicit_TD && fContributionTime == Last_CT))
656  {
657  ContributeExplConvFace(data.x,dataleft.sol[0],dataright.sol[0],weight,data.normal,dataleft.phi,dataright.phi,ef);
658  }
659 
660 }
661 
662 void TPZEulerConsLaw::ContributeBC(TPZMaterialData &data,
663  REAL weight,
664  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef,
665  TPZBndCond &bc)
666 {
667  int phr = data.phi.Rows();
668  short in,jn,i,j;
669  int nstate = NStateVariables();
670  REAL v2[5];//m�imo nstate
671  for(i=0;i<nstate;i++) v2[i] = bc.Val2()(i,0);
672 
673  switch (bc.Type()) {
674  case 0 :// Dirichlet condition
675  for(in = 0 ; in < phr; in++) {
676  for(i = 0 ; i < nstate; i++)
677  ef(in*nstate+i,0) += gBigNumber * weight * v2[i] * data.phi(in,0);
678  for (jn = 0 ; jn < phr; jn++) {
679  for(i = 0 ; i < nstate; i++)
680  ek(in*nstate+i,jn*nstate+i) -= gBigNumber * weight * data.phi(in,0) * data.phi(jn,0);
681  }
682  }
683  break;
684  case 1 :// Neumann condition
685  for(in = 0 ; in < data.phi.Rows(); in++) {
686  for(i = 0 ; i < nstate; i++)
687  ef(in*nstate+i,0) += v2[i] * data.phi(in,0) * weight;
688  }
689  break;
690  case 2 :// condi�o mista
691  for(in = 0 ; in < data.phi.Rows(); in++) {
692  for(i = 0 ; i < nstate; i++)
693  ef(in*nstate+i, 0) += weight * v2[i] * data.phi(in, 0);
694  for (jn = 0 ; jn < data.phi.Rows(); jn++) {
695  for(i = 0 ; i < nstate; i++) for(j = 0 ; j < nstate; j++)
696  ek(in*nstate+i,jn*nstate+j) -= weight * bc.Val1()(i,j) * data.phi(in,0) * data.phi(jn,0);
697  }
698  }
699  }
700 }
701 
702 void TPZEulerConsLaw::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
703  REAL weight,
704  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,
705  TPZBndCond &bc)
706 {
707  int nstate = NStateVariables();
708  TPZVec<STATE> solR(nstate,0.);
709  int numbersol = dataleft.sol.size();
710  if (numbersol != 1) {
711  DebugStop();
712  }
713 
714  TPZFMatrix<STATE> dsolR(dataleft.dsol[0].Rows(), dataleft.dsol[0].Cols(),0.);
715  TPZFMatrix<REAL> phiR(0,0), dphiR(0,0);
716  int entropyFix;
717 
718  // contributing face-based quantities
719  if (fConvFace == Implicit_TD && fContributionTime == Advanced_CT)
720  {
721  // if face contribution is implicit,
722  // then the FAD classes must be initialized
723 #ifdef _AUTODIFF
724 
725  // if(bc.Type() ==5 && fDim == 2)
726  // {
727  // int entropyFix2;
728  // TPZManVector<REAL,5 > flux(nstate,0.);
729  // ComputeGhostState(solL, solR, normal, bc, entropyFix2);
730  // Roe_Flux<REAL>(solL, solR, normal, fGamma, flux, entropyFix2);
731  // REAL norflux = flux[1]*normal[1]-flux[2]*normal[0];
732  // REAL err = fabs(flux[0])+fabs(norflux)+fabs(flux[3]);
733  // if(err > 1.e-5)
734  // {
735  // cout << "fluxo de parede errado 1 err " << err << endl;
736  // }
737  // }
738 
739 #ifdef FASTEST_IMPLICIT
740  ContributeFastestBCInterface(fDim, data.x, dataleft.sol[0], dataleft.dsol[0],
741  weight, data.normal, dataleft.phi, data.dphixl, data.axesleft, ek, ef, bc);
742 #else
743  TPZVec<FADREAL> FADsolL, FADsolR;
744  PrepareInterfaceFAD(dataleft.sol[0], solR, dataleft.phi, phiR, FADsolL, FADsolR);
745  ComputeGhostState(FADsolL, FADsolR, data.normal, bc, entropyFix);
746  ContributeImplConvFace(data.x,FADsolL,FADsolR, weight, data.normal, dataleft.phi, phiR, ek, ef, entropyFix);
747 #ifdef LOG4CXX
748  if(fluxroe->isDebugEnabled()){
749  std::stringstream sout;
750  ek.Print("computed tangent matrix",sout);
751  ef.Print("computed rhs",sout);
752  LOGPZ_DEBUG(fluxroe,sout.str().c_str());
753  }
754 #endif
755 #endif
756 #else
757  // forcint explicit contribution and issueing an warning
758  std::cout << "TPZEulerConsLaw::ContributeInterface> Implicit face convective contribution: _AUTODIFF directive not configured";
759  // ComputeGhostState(FADsolL, FADsolR, normal, bc, entropyFix);
760  // ContributeApproxImplConvFace(x,solL,solR,weight,normal,phiL,phiR,ek,ef,entropyFix);
761 #endif
762  }
763  else if (fConvFace == ApproxImplicit_TD && fContributionTime == Advanced_CT)
764  {
765 #ifdef _AUTODIFF
766  TPZVec<FADREAL> FADsolL, FADsolR;
767  PrepareInterfaceFAD(dataleft.sol[0], solR, dataleft.phi, phiR, FADsolL, FADsolR);
768  ComputeGhostState(FADsolL, FADsolR, data.normal, bc, entropyFix);
769  ContributeApproxImplConvFace(data.x,data.HSize,FADsolL,FADsolR, weight, data.normal, dataleft.phi, phiR, ek, ef, entropyFix);
770 #ifdef LOG4CXX
771  if(fluxappr->isDebugEnabled()){
772  std::stringstream sout;
773  ek.Print("computed tangent matrix",sout);
774  ef.Print("computed rhs",sout);
775  LOGPZ_DEBUG(fluxappr,sout.str().c_str());
776  }
777 #endif
778 #else
779  // ComputeGhostState(solL, solR, normal, bc, entropyFix);
780  // ContributeApproxImplConvFace(x,data.HSize,FADsolL,FADsolR, weight, normal, phiL, phiR, ek, ef, entropyFix);
781 #endif
782  }
783 
784  if(fConvFace == Explicit_TD && fContributionTime == Last_CT)
785  {
786  ComputeGhostState(dataleft.sol[0], solR, data.normal, bc, entropyFix);
787  ContributeExplConvFace(data.x,dataleft.sol[0],solR,weight,data.normal,dataleft.phi,phiR,ef, entropyFix);
788  }
789 }
790 
791 void TPZEulerConsLaw::ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
792  REAL weight,
793  TPZFMatrix<STATE> &ef,
794  TPZBndCond &bc)
795 {
796  int numbersol = dataleft.sol.size();
797  if (numbersol != 1) {
798  DebugStop();
799  }
800 
801  int nstate = NStateVariables();
802  TPZVec<STATE> solR(nstate,0.);
803  TPZFMatrix<STATE> dsolR(dataleft.dsol[0].Rows(), dataleft.dsol[0].Cols(),0.);
804  TPZFMatrix<REAL> phiR(0,0), dphiR(0,0);
805  int entropyFix;
806 
807  if((fConvFace == Implicit_TD && fContributionTime == Advanced_CT)
808  ||
809  (fConvFace == ApproxImplicit_TD && fContributionTime == Advanced_CT)
810  ||
811  (fConvFace == Explicit_TD && fContributionTime == Last_CT))
812  {
813  ComputeGhostState(dataleft.sol[0], solR, data.normal, bc, entropyFix);
814 
815  if(fDim == 2)
816  {// flux tests
817  TPZManVector<REAL,3> normal2(2,0.);
818  TPZManVector<STATE,5> flux2(nstate,0.);
819  normal2[0] = -data.normal[0];
820  normal2[1] = -data.normal[1];
821  TPZManVector<STATE,5 > flux(nstate,0.);
822  Roe_Flux<STATE>(dataleft.sol[0], solR, data.normal, fGamma, flux);
823  Roe_Flux<STATE>(solR, dataleft.sol[0], normal2, fGamma, flux2);
824  STATE fluxs = fabs(flux[0]+flux2[0])+fabs(flux[1]+flux2[1])+fabs(flux[2]+flux2[2])+fabs(flux[3]+flux2[3]);
825  if(fluxs > 1.e-10)
826  {
827  std::cout << "Fluxo nao simetrico fluxs = " << fluxs << std::endl;
828  }
829 
830  if(bc.Type() ==5)
831  {
832  REAL norflux = flux[1]*data.normal[1]-flux[2]*data.normal[0];
833  REAL err = fabs(flux[0])+fabs(norflux)+fabs(flux[3]);
834  if(err > 1.e-5)
835  {
836  std::cout << "fluxo de parede errado 2 err " << err << std::endl;
837  Roe_Flux<STATE>(dataleft.sol[0],solR,data.normal,fGamma,flux);
838  } else
839  {
840  Roe_Flux<STATE>(dataleft.sol[0],solR,data.normal,fGamma,flux);
841  }
842  }
843  } // end of tests
844 
845  ContributeExplConvFace(data.x,dataleft.sol[0],solR,weight,data.normal,
846  dataleft.phi,phiR,ef,entropyFix);
847  }
848 }
849 
850 #ifdef _AUTODIFF
851 
852 void TPZEulerConsLaw::ContributeFastestBCInterface(int dim,
853  TPZVec<REAL> &x,
854  TPZVec<STATE> &solL, TPZFMatrix<STATE> &dsolL,
855  REAL weight, TPZVec<REAL> &normal,
856  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &dphiL, TPZFMatrix<REAL> &axesleft,
857  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
858 {
859 
860  switch(dim)
861  {
862  case(1):
863  ContributeFastestBCInterface_dim<1>(x, solL, dsolL,
864  weight, normal,
865  phiL, dphiL,axesleft,
866  ek, ef, bc);
867  break;
868  case(2):
869  ContributeFastestBCInterface_dim<2>(x, solL, dsolL,
870  weight, normal,
871  phiL, dphiL,axesleft,
872  ek, ef, bc);
873  break;
874  case(3):
875  ContributeFastestBCInterface_dim<3>(x, solL, dsolL,
876  weight, normal,
877  phiL, dphiL,axesleft,
878  ek, ef, bc);
879  break;
880  default:
881  PZError << "\nTPZEulerConsLaw::ContributeFastestBCInterface unhandled dimension\n";
882  exit(-1);
883  }
884 };
885 
886 
887 template <int dim>
888 void TPZEulerConsLaw::ContributeFastestBCInterface_dim(TPZVec<REAL> &x,
889  TPZVec<STATE> &solL, TPZFMatrix<STATE> &dsolL,
890  REAL weight, TPZVec<REAL> &normal,
891  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &dphiL,TPZFMatrix<REAL> &axesleft,
892  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
893 {
894 #ifdef _TFAD
895  typedef TFad<2*(dim+2), REAL> TFADREALInterface;
896 #endif
897 #ifdef _FAD
898  typedef Fad<REAL> TFADREALInterface;
899 #endif
900 #ifdef _TINYFAD
901  typedef TinyFad<2*(dim+2), REAL> TFADREALInterface;
902 #endif
903 
904  int entropyFix;
905 
906  int nstate = NStateVariables();
907  TPZVec<STATE> solR(nstate,0.);
908  TPZFMatrix<REAL> phiR(0,0);
909 
910  // initial guesses for left and right sol
911  // fForcingFunction is null at iterations > 0
912  if(fForcingFunction)
913  {
914  TPZVec<STATE> res;
915  fForcingFunction->Execute(x, res);
916  for(int i = 0; i < nstate; i++)
917  solL[i] = solR[i] = res[i];
918  }
919 
920  TPZVec<TFADREALInterface > FADsolL(nstate),
921  FADsolR(nstate);
922 
923  PrepareFastestInterfaceFAD(solL, solR, FADsolL, FADsolR);
924 
925  ComputeGhostState(FADsolL, FADsolR, normal, bc, entropyFix);
926 
927  ContributeFastestImplConvFace_T(x, FADsolL, FADsolR,
928  weight, normal,
929  phiL, phiR,
930  ek, ef,
931  entropyFix);
932 };
933 
934 #endif
935 
936 template <class T>
937 void TPZEulerConsLaw:: ComputeGhostState(TPZVec<T> &solL, TPZVec<T> &solR, TPZVec<REAL> &normal, TPZBndCond &bc, int & entropyFix)
938 {
939  entropyFix = 1;
940 
941  int nstate = NStateVariables();
942  T vpn=REAL(0.);
943  T us, un, c;
944  REAL Mach, temp;
945  int i;
946  //Riemann Invariants
947  T w1, w2, w5, uninf, usinf, cghost, usghost, unghost, p;
948  REAL cinf;
949 
950  switch (bc.Type()){
951  case 3://Dirichlet: nada a fazer a CC �a correta
952  for(i=0;i<nstate; i++) solR[i] = bc.Val2().operator()(i,0);
953  break;
954  case 4://recuperar valor da solu� MEF esquerda: saida livre
955  for(i=0;i<nstate; i++) solR[i] = solL[i];
956  break;
957  case 5://condi� de parede
958  for(i=1;i<nstate-1;i++) vpn += solL[i]*T(normal[i-1]);//v.n
959  for(i=1;i<nstate-1;i++) solR[i] = solL[i] - T(2.0*normal[i-1])*vpn;
960  solR[0] = solL[0];
961  solR[nstate-1] = solL[nstate-1];
962  entropyFix = 0;
963  break;
964  case 6://n� refletivas (campo distante)
965  for(i=0;i<nstate;i++) solR[i] = solL[i];
966  break;
967  case 7:// INLET (Dirichlet using Mach)
968  Mach = bc.Val2().operator()(1,0);
969  solR[0] = bc.Val2().operator()(0,0);//solL[0];
970  solR[nstate-1] = bc.Val2().operator()(nstate - 1,0);
971 
972  temp = Mach * Mach * fGamma * (fGamma - 1);
973  us = sqrt(2 * temp * solR[nstate-1] /
974  ( solR[0] * (2 + temp)) );
975 
976  for(i=1;i<nstate-1;i++) solR[i] = - us * normal[i-1];
977 
978  /* solR[nstate-1] = bc.Val2().operator()(nstate - 1,0)/(fGamma -1.) +
979  solR[0] * us * us / 2.;*/
980 
981  break;
982  case 8:// OUTLET (Dirichlet using Mach)
983  Mach = bc.Val2().operator()(1,0);
984  if(bc.Val2().operator()(0,0) == 0.)
985  {
986  solR[0] = solL[0];
987  }else
988  {
989  solR[0] = bc.Val2().operator()(0,0);//solL[0];
990  }
991 
992  temp = Mach * Mach * fGamma * (fGamma - 1);
993  us = sqrt(T(2.) * temp * solR[nstate-1] /
994  ( solR[0] * (T(2.) + temp)) );
995 
996  for(i=1;i<nstate-1;i++) solR[i] = us * normal[i-1];
997  break;
998  case 9:// INFLOW/OUTFLOW (dedpending on direction of internal
999  // velocity vector.
1000  // Inputs are in terms of primitive variables
1001  // rho, Mach, p
1002 
1003  // computing normal velocity and speed norm
1004  un = 0.;
1005  us = 0.;
1006  Mach = 0.;
1007  for(i = 1; i < nstate-1; i++)
1008  {
1009  un += solL[i]/solL[0]*normal[i-1];
1010  us += solL[i] * solL[i] / solL[0] / solL[0];
1011  Mach += bc.Val2()(i,0)*bc.Val2()(i,0);
1012  }
1013  us = sqrt(us);
1014  Mach = sqrt(Mach);
1015 
1016  // cinf = sqrt(gamma p / rho)
1017  cinf = sqrt(fGamma * bc.Val2()(nstate-1,0)/bc.Val2()(0,0));
1018 
1019  //computing the pressure
1020  // p = (gamma - 1.) * (rhoe - rho*vel^2 / 2)
1021  p = (fGamma - 1.) * (solL[nstate - 1] - solL[0] * us * us / T(2.));
1022  //c speed
1023  c = sqrt(fGamma * p / solL[0]);
1024 
1025  usinf = /*bc.Val2()(1,0)*/ Mach * cinf;
1026  uninf = un / us * usinf;
1027 
1028  if(un < REAL(0.))// Inflow
1029  {
1030  //if(normal[0]>0.) cout << "\ndirection error\n ";
1031 
1032  if(/*us>c*/ Mach >= 1.)
1033  {//supersonic
1034  // all Riemann invariants retain their imposed values
1035  solR[0] = bc.Val2()(0,0);
1036  // rho vel = rho * Mach * c * u_directioni / us
1037  for(i = 1; i < nstate-1; i++)
1038  solR[i] = (T)bc.Val2()(0,0) *
1039  usinf *
1040  solL[i]/solL[0] / us; // versor (direction)
1041  // rhoe = p / (gamma - 1) + rho * vel*vel/2
1042  solR[nstate-1] = T(bc.Val2()(nstate-1,0)) / T(fGamma - 1.) +
1043  T(bc.Val2()(0,0)) * usinf * usinf / T(2.);
1044  }
1045  else
1046  {//subsonic
1047  // Invariants w1 and w2 are imposed, w5 computed
1048  w1 = uninf - T(2.) * cinf/ T(fGamma - 1.);
1049  // Modified w2 invariant: w2 = p/rho^(gamma-1)
1050  // or w2 = c^2/(gamma * rho^(gamma-1))
1051  w2 = cinf * cinf / T(fGamma * pow((T)(bc.Val2()(0,0)), (T)(fGamma - 1.)));
1052  // w5 computed based on flow state
1053  w5 = un + T(2.) * c / T(fGamma - 1.);
1054 
1055  // computing ghost values
1056  cghost = (w5 - w1) * T((fGamma - 1.)/4.);
1057  solR[0] = pow(((T)cghost * cghost / (T(fGamma) * w2)), (T)(1./(fGamma - 1.)));
1058  unghost = (w1 + w5) / T(2.);
1059  usghost = us / un * unghost;
1060  for(i = 1; i < nstate - 1; i++)
1061  solR[i] = solR[0] // rho
1062  * usghost * // velocity
1063  bc.Val2()(i,0) / Mach; // element velocity component
1064 
1065 
1066  // rhoe = rho * (c^2 / (gamma(gamma -1)) + vel^2/2)
1067  solR[nstate - 1] = solR[0] * (
1068  cghost * cghost /T(fGamma * (fGamma - 1.)) +
1069  usghost * usghost / T(2.));
1070  /*
1071  // Invariants w1 and w2 are imposed, w5 computed
1072  w1 = uninf - T(2.) * cinf/ T(fGamma - 1.);
1073  // Modified w2 invariant: w2 = p/rho^(gamma-1)
1074  // or w2 = c^2/(gamma * rho^(gamma-1))
1075  w2 = cinf * cinf / T(fGamma * pow(bc.Val2()(0,0), fGamma - 1.));
1076  // w5 computed based on flow state
1077  w5 = un + T(2.) * c / T(fGamma - 1.);
1078 
1079  // computing ghost values
1080  cghost = (w5 - w1) * T((fGamma - 1.)/4.);
1081  solR[0] = pow(cghost * cghost / (T(fGamma) * w2), 1./(fGamma - 1.));
1082  unghost = (w1 + w5) / T(2.);
1083  for(i = 1; i < nstate - 1; i++)
1084  solR[i] = solR[0] // rho
1085  * unghost / un * // scale factor
1086  solL[i] / solL[0]; // element velocity component
1087  usghost = us / un * unghost;
1088 
1089  // rhoe = rho * (c^2 / (gamma(gamma -1)) + vel^2/2)
1090  solR[nstate - 1] = solR[0] * (
1091  cghost * cghost /T(fGamma * (fGamma - 1.)) +
1092  usghost * usghost / T(2.));*/
1093  }
1094  }else
1095  { // Outflow
1096  //if(normal[0]<0.) cout << "\ndirection error 2\n ";
1097  if(us>c)
1098  { // supersonic: no BC at all are required
1099  for(i = 0; i < nstate; i++)
1100  solR[i] = solL[i];
1101  }else
1102  { // subsonic outlet
1103  // only the condition w1 referring to the first Riemann invariant
1104  // is imposed. As a rule, the imposition of pressure is applied
1105  // instead.
1106 
1107  solR[0] = solL[0] * (T)pow(T(bc.Val2()(nstate-1,0))/p, T(1./fGamma));
1108 
1109  cghost = sqrt(fGamma * T(bc.Val2()(nstate-1,0))/ T(solR[0]));
1110 
1111  unghost = (c - cghost) * T(2./(fGamma - 1.));
1112  usghost = 0.;
1113  // ughost = u + 2.*(c - cghost)/(fGamma - 1.) * normal
1114  for(i = 1; i < nstate - 1; i++)
1115  {
1116  solR[i] = solR[0] * // rho
1117  (solL[i] / solL[0] + // element vel
1118  unghost * normal[i-1]); // ghost correction
1119  usghost += solR[i] * solR[i] / solR[0] / solR[0];
1120  }
1121  usghost = sqrt(usghost);
1122 
1123  // rhoe = p / (gamma - 1) + rho * vel*vel/2
1124  solR[nstate-1] = T(bc.Val2()(nstate-1,0)/(fGamma - 1.)) +
1125  solR[0] * usghost * usghost / T(2.);
1126 
1127  }
1128  }
1129  break;
1130 
1131 
1132  case 10:// Directional INFLOW
1133  // velocity vector.
1134  // Inputs are in terms of primitive variables
1135  // rho, Machx, Machy, Machz, p
1136 
1137  // computing normal velocity and speed norm
1138  un = 0.;
1139  us = 0.;
1140  Mach = 0.;
1141  for(i = 1; i < nstate-1; i++)
1142  {
1143  un += solL[i]/solL[0]*normal[i-1];
1144  us += solL[i] * solL[i] / solL[0] / solL[0];
1145  Mach += bc.Val2()(i,0)*bc.Val2()(i,0);
1146  }
1147  us = sqrt(us);
1148  Mach = sqrt(Mach);
1149 
1150  // cinf = sqrt(gamma p / rho)
1151  cinf = sqrt(fGamma * bc.Val2()(nstate-1,0)/bc.Val2()(0,0));
1152 
1153  //computing the pressure
1154  // p = (gamma - 1.) * (rhoe - rho*vel^2 / 2)
1155  p = (fGamma - 1.) * (solL[nstate - 1] - solL[0] * us * us / T(2.));
1156  //c speed
1157  c = sqrt(fGamma * p / solL[0]);
1158 
1159  usinf = /*bc.Val2()(1,0)*/ Mach * cinf;
1160  uninf = un / us * usinf;
1161 
1162  if(un < REAL(0.))// Inflow
1163  {
1164  if(/*us>c*/ Mach >= 1.)
1165  {//supersonic
1166  // all Riemann invariants retain their imposed values
1167  solR[0] = bc.Val2()(0,0);
1168  // rho vel = rho * Mach * c * u_directioni / us
1169  for(i = 1; i < nstate-1; i++)
1170  solR[i] = T(bc.Val2()(0,0)) *
1171  usinf *
1172  solL[i]/solL[0] / us; // versor (direction)
1173  // rhoe = p / (gamma - 1) + rho * vel*vel/2
1174  solR[nstate-1] = T(bc.Val2()(nstate-1,0)) / T(fGamma - 1.) +
1175  T(bc.Val2()(0,0)) * usinf * usinf / T(2.);
1176  }
1177  else
1178  {//subsonic
1179  // Invariants w1 and w2 are imposed, w5 computed
1180  w1 = uninf - T(2.) * cinf/ T(fGamma - 1.);
1181  // Modified w2 invariant: w2 = p/rho^(gamma-1)
1182  // or w2 = c^2/(gamma * rho^(gamma-1))
1183  w2 = cinf * cinf / T(fGamma * pow((T)(bc.Val2()(0,0)), (T)(fGamma - 1.))); // Unbelived - only is defined for pow(float, float) and pow(long double, long double) !!! Jorge
1184  // w5 computed based on flow state
1185  w5 = un + T(2.) * c / T(fGamma - 1.);
1186 
1187  // computing ghost values
1188  cghost = (w5 - w1) * T((fGamma - 1.)/4.);
1189  solR[0] = pow(cghost * cghost / (T(fGamma) * w2),(T)(1./(fGamma - 1.)));
1190  unghost = (w1 + w5) / T(2.);
1191  usghost = us / un * unghost;
1192  for(i = 1; i < nstate - 1; i++)
1193  solR[i] = solR[0] // rho
1194  * usghost * // velocity
1195  bc.Val2()(i,0) / Mach; // element velocity component
1196 
1197 
1198  // rhoe = rho * (c^2 / (gamma(gamma -1)) + vel^2/2)
1199  solR[nstate - 1] = solR[0] * (
1200  cghost * cghost /T(fGamma * (fGamma - 1.)) +
1201  usghost * usghost / T(2.));
1202  }
1203  }
1204  break;
1205 
1206  case 11:// SOLID WALL - No input needed
1207 
1208  entropyFix = 0;
1209  // Invariants w1 and w2 are imposed, w5 computed
1210  // This gives a set of closed BCs.
1211 
1212  // computing normal velocity and speed norm
1213  un = 0.;
1214  us = 0.;
1215  for(i = 1; i < nstate-1; i++)
1216  {
1217  un += solL[i]/solL[0]*normal[i-1];
1218  us += solL[i] * solL[i] / solL[0] / solL[0];
1219  }
1220  us = sqrt(us);
1221  //computing the pressure
1222  // p = (gamma - 1.) * (rhoe - rho*vel^2 / 2)
1223  p = (fGamma - 1.) * (solL[nstate - 1] - solL[0] * us * us / T(2.));
1224  //c speed
1225  c = sqrt(fGamma * p / solL[0]);
1226 
1227  // computing the ghost sound speed
1228  // cghost = c + (gamma - 1) * un / 2
1229  cghost = c + T((fGamma - 1.)/2.)*un;
1230  // ghost density
1231  // rhoghost = (cghost^2*rho^gamma/(gamma * p))^(1/gamma -1)
1232  solR[0] = pow(cghost * cghost * pow(T(solL[0]),T(fGamma))/(T(p * fGamma)), (T)(1./(fGamma - 1.)));
1233 
1234  // computing velocity vector
1235  usghost = 0.;
1236  for(i = 1; i < nstate-1; i++)
1237  {
1238  // vghost = u - (un * n)
1239  solR[i] = solR[0] * (solL[i]/ solL[0] - un * normal[i-1]);
1240  usghost += solR[i] * solR[i] / solR[0] / solR[0];
1241  }
1242  usghost = sqrt(usghost);
1243 
1244  // rhoe = rho * (c^2 / (gamma(gamma -1)) + vel^2/2)
1245  solR[nstate - 1] = solR[0] * (
1246  cghost * cghost /T(fGamma * (fGamma - 1.)) +
1247  usghost * usghost / T(2.));
1248  break;
1249  case 12: //Symmetry
1250  for(i=1;i<nstate-1;i++) vpn += solL[i]*T(normal[i-1]);//v.n
1251  for(i=1;i<nstate-1;i++) solR[i] = solL[i] - T(2.0*normal[i-1])*vpn;
1252  solR[0] = solL[0];
1253  solR[nstate-1] = solL[nstate-1];
1254  entropyFix = 1;
1255  break;
1256  default:
1257  for(i=0;i<nstate;i++) solR[i] = 0.;
1258  }
1259 }
1260 
1261 //------------------internal contributions
1262 
1263 void TPZEulerConsLaw::ContributeApproxImplDiff(TPZVec<REAL> &x,
1264  TPZFMatrix<REAL> &jacinv,
1265  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
1266  REAL weight,
1267  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
1268  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
1269 {
1270  // computing the determinant of the jacobian
1271  REAL deltaX = DeltaX(1./Det(jacinv));
1272 
1273  fArtDiff.ContributeApproxImplDiff(fDim, jacinv, sol, dsol, dphi,
1274  ek, ef, weight,
1275 #ifdef DIVTIMESTEP
1276  1.,
1277 #else
1278  TimeStep(),
1279 #endif
1280  deltaX
1281  );
1282 }
1283 
1284 void TPZEulerConsLaw::ContributeExplDiff(TPZVec<REAL> &x,
1285  TPZFMatrix<REAL> &jacinv,
1286  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
1287  REAL weight,
1288  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
1289  TPZFMatrix<STATE> &ef)
1290 {
1291  // computing the determinant of the jacobian
1292  REAL deltaX = DeltaX(1./Det(jacinv));
1293 
1294  fArtDiff.ContributeExplDiff(fDim, jacinv, sol, dsol, dphi,
1295  ef, weight,
1296 #ifdef DIVTIMESTEP
1297  1.,
1298 #else
1299  TimeStep(),
1300 #endif
1301  deltaX
1302  );
1303 }
1304 
1305 
1306 #ifdef _AUTODIFF
1307 void TPZEulerConsLaw::ContributeImplDiff(TPZVec<REAL> &x,
1308  TPZFMatrix<REAL> &jacinv,
1309  TPZVec<FADREAL> &sol,TPZVec<FADREAL> &dsol,
1310  REAL weight,
1311  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
1312 {
1313  // computing the determinant of the jacobian
1314  REAL deltaX = DeltaX(1./Det(jacinv));
1315 
1316  fArtDiff.ContributeImplDiff(fDim, jacinv, sol, dsol,
1317  ek, ef, weight,
1318 #ifdef DIVTIMESTEP
1319  1.,
1320 #else
1321  TimeStep(),
1322 #endif
1323  deltaX
1324  );
1325 }
1326 
1327 
1328 void TPZEulerConsLaw::ContributeFastestImplDiff(int dim, TPZVec<REAL> &x, TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol, TPZFMatrix<REAL> &phi,
1329  TPZFMatrix<REAL> &dphi, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
1330 {
1331  // computing the determinant of the jacobian
1332  REAL deltaX = DeltaX(1./Det(jacinv));
1333 
1334  fArtDiff.ContributeFastestImplDiff(dim, jacinv, sol, dsol, phi, dphi,
1335  ek, ef, weight,
1336 #ifdef DIVTIMESTEP
1337  1.,
1338 #else
1339  TimeStep(),
1340 #endif
1341  deltaX
1342  );
1343 }
1344 
1345 #endif
1346 
1347 void TPZEulerConsLaw::ContributeExplConvFace(TPZVec<REAL> &x,
1348  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
1349  REAL weight,TPZVec<REAL> &normal,
1350  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
1351  TPZFMatrix<STATE> &ef, int entropyFix)
1352 {
1353  int nState = NStateVariables();
1354  TPZVec<STATE> flux(nState,0.);
1355  Roe_Flux<STATE>(solL, solR, normal, fGamma, flux, entropyFix);
1356  int nShapeL = phiL.Rows(),
1357  nShapeR = phiR.Rows(),
1358  i_shape, i_state;
1359 
1360 #ifdef DIVTIMESTEP
1361  REAL constant = weight; // weight
1362 #else
1363  REAL constant = TimeStep() * weight; // deltaT * weight
1364 #endif
1365 
1366 
1367  // Contribution referring to the left element
1368  for(i_shape = 0; i_shape < nShapeL; i_shape ++)
1369  for(i_state = 0; i_state < nState; i_state++)
1370  ef(i_shape*nState + i_state,0) +=
1371  flux[i_state] * phiL(i_shape,0) * constant;
1372 
1373  // Contribution referring to the right element
1374  // REM: The contributions are negative in comparison
1375  // to the left elements: opposed normals
1376  for(i_shape = 0; i_shape < nShapeR; i_shape ++)
1377  for(i_state = 0; i_state < nState; i_state++)
1378  ef((nShapeL + i_shape)*nState + i_state,0) -=
1379  flux[i_state] * phiR(i_shape,0) * constant;
1380 }
1381 
1382 #ifdef _AUTODIFF
1383 
1384 void TPZEulerConsLaw::ContributeApproxImplConvFace(TPZVec<REAL> &x, REAL faceSize,
1385  TPZVec<FADREAL> &solL,TPZVec<FADREAL> &solR,
1386  REAL weight,TPZVec<REAL> &normal,
1387  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
1388  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, int entropyFix)
1389 {
1390  int nState = NStateVariables();
1391  TPZVec<FADREAL > flux(nState,REAL(0.));
1392  ApproxRoe_Flux(solL, solR, normal, fGamma, flux, entropyFix);
1393 
1394  // Testing whether Roe_Flux<REAL> gives the same result as Roe_Flux<FADREAL>
1395 
1396  /* TPZVec<REAL> solL2(nState,0.),solR2(nState,0.),flux2(nState,0.);
1397  int i;
1398  for(i=0; i<nState; i++)
1399  {
1400  solL2[i] = solL[i].val();
1401  solR2[i] = solR[i].val();
1402  }
1403  Roe_Flux(solL2, solR2, normal, fGamma, flux2,with_entropy_fix);
1404  REAL diff = fabs(flux[0].val()-flux2[0])+fabs(flux[1].val()-flux2[1])+fabs(flux[2].val()-flux2[2]);
1405  if(diff != 0.)
1406  {
1407  cout << "Roe<FADREAL> is different from Roe<REAL> diff " << diff << endl;
1408  }*/
1409  int nShapeL = phiL.Rows(),
1410  nShapeR = phiR.Rows(),
1411  i_shape, i_state, j,
1412  nDer = (nShapeL + nShapeR) * nState;
1413 
1414 
1415 #ifdef DIVTIMESTEP
1416  REAL constant = weight; // weight
1417 #else
1418  REAL constant = TimeStep() * weight; // deltaT * weight
1419 #endif
1420 
1421  // Contribution referring to the left element
1422  for(i_shape = 0; i_shape < nShapeL; i_shape ++)
1423  for(i_state = 0; i_state < nState; i_state++)
1424  {
1425  int index = i_shape*nState + i_state;
1426  ef(index,0) +=
1427  flux[i_state].val() * phiL(i_shape,0) * constant;
1428  for(j = 0; j < nDer; j++)
1429  ek(index, j) -= flux[i_state].dx/*fastAccessDx*/(j) *
1430  phiL(i_shape,0) * constant;
1431  }
1432 
1433  // Contribution referring to the right element
1434  // REM: The contributions are negative in comparison
1435  // to the left elements: opposed normals
1436  for(i_shape = 0; i_shape < nShapeR; i_shape ++)
1437  for(i_state = 0; i_state < nState; i_state++)
1438  {
1439  int index = (nShapeL + i_shape)*nState + i_state;
1440  ef(index,0) -=
1441  flux[i_state].val() * phiR(i_shape,0) * constant;
1442  for(j = 0; j < nDer; j++)
1443  ek(index, j) += flux[i_state].dx/*fastAccessDx*/(j) *
1444  phiR(i_shape,0) * constant;
1445  }
1446 
1447 }
1448 
1449 #endif
1450 
1451 void TPZEulerConsLaw::ContributeApproxImplConvFace(TPZVec<REAL> &x, REAL faceSize,
1452  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
1453  REAL weight,TPZVec<REAL> &normal,
1454  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
1455  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, int entropyFix
1456  )
1457 {
1458  int nState = NStateVariables();
1459  TPZManVector< TPZManVector<STATE,5> ,3> FL(3),FR(3);
1460  TPZManVector< STATE,5> FN(nState,0.);
1461  Flux(solL, FL[0], FL[1], FL[2]);
1462  Flux(solR, FR[0], FR[1], FR[2]);
1463  TPZFNMatrix<36> DFNL(nState,nState,0.),DFNR(nState,nState,0.);
1464 
1465  TPZManVector<STATE,7 > fluxroe(nState,0.);
1466  Roe_Flux(solL, solR, normal, fGamma, fluxroe,entropyFix);
1467  TPZManVector< TPZDiffMatrix<STATE> ,3> AL(3),AR(3);
1468  JacobFlux(fGamma, fDim, solL, AL);
1469  JacobFlux(fGamma, fDim, solR, AR);
1470 
1471  int nShapeL = phiL.Rows(),
1472  nShapeR = phiR.Rows(),
1473  i_shape, i_state, i, j_state, j_shape;
1474  // nDer = (nShapeL + nShapeR) * nState;
1475 
1476 
1477 #ifdef DIVTIMESTEP
1478  REAL constant = weight; // weight
1479 #else
1480  REAL constant = TimeStep() * weight; // deltaT * weight
1481 #endif
1482 
1483  for(i_state=0; i_state<nState; i_state++)
1484  {
1485  FN[i_state]= -(faceSize/constant)*(solR[i_state]-solL[i_state]);
1486  DFNL(i_state,i_state) = (faceSize/constant);
1487  DFNR(i_state,i_state) = -(faceSize/constant);
1488  for(i=0; i<fDim; i++)
1489  {
1490  FN[i_state]+=0.5*normal[i]*(FL[i][i_state]+FR[i][i_state]);
1491  for(j_state=0; j_state<nState; j_state++)
1492  {
1493  DFNL(i_state,j_state) +=0.5*normal[i]*(AL[i](i_state,j_state));
1494  DFNR(i_state,j_state) +=0.5*normal[i]*(AR[i](i_state,j_state));
1495  }
1496  }
1497  }
1498  // Contribution referring to the left element
1499  for(i_shape = 0; i_shape < nShapeL; i_shape ++)
1500  for(i_state = 0; i_state < nState; i_state++)
1501  {
1502  int index = i_shape*nState + i_state;
1503  ef(index,0) +=
1504  fluxroe[i_state] * phiL(i_shape,0) * constant;
1505  for(j_shape = 0; j_shape < nShapeL; j_shape++)
1506  {
1507 
1508  for(j_state = 0; j_state < nState; j_state++)
1509  {
1510  int jndex = j_shape*nState + j_state;
1511  ek(index,jndex) -= DFNL(i_state,j_state) *phiL(i_shape,0) * phiL(j_shape,0) * constant;
1512  }
1513  }
1514  for(j_shape = 0; j_shape < nShapeR; j_shape++)
1515  {
1516 
1517  for(j_state = 0; j_state < nState; j_state++)
1518  {
1519  int jndex = (nShapeL+j_shape)*nState + j_state;
1520  ek(index,jndex) -= DFNR(i_state,j_state) *phiL(i_shape,0) * phiR(j_shape,0) * constant;
1521  }
1522  }
1523  /* ek(index, j) -= flux[i_state].dx(j) *
1524  phiL(i_shape,0) * constant;*/
1525  }
1526 
1527  // Contribution referring to the right element
1528  // REM: The contributions are negative in comparison
1529  // to the left elements: opposed normals
1530  for(i_shape = 0; i_shape < nShapeR; i_shape ++)
1531  for(i_state = 0; i_state < nState; i_state++)
1532  {
1533  int index = (nShapeL + i_shape)*nState + i_state;
1534  ef(index,0) -=
1535  fluxroe[i_state] * phiR(i_shape,0) * constant;
1536  for(j_shape = 0; j_shape < nShapeL; j_shape++)
1537  {
1538 
1539  for(j_state = 0; j_state < nState; j_state++)
1540  {
1541  int jndex = j_shape*nState + j_state;
1542  ek(index,jndex) += DFNL(i_state,j_state) *phiR(i_shape,0) * phiL(j_shape,0) * constant;
1543  }
1544  }
1545  for(j_shape = 0; j_shape < nShapeR; j_shape++)
1546  {
1547 
1548  for(j_state = 0; j_state < nState; j_state++)
1549  {
1550  int jndex = (nShapeL+j_shape)*nState + j_state;
1551  ek(index,jndex) += DFNR(i_state,j_state) *phiR(i_shape,0) * phiR(j_shape,0) * constant;
1552  }
1553  }
1554  /* ek(index, j) += flux[i_state].dx(j) *
1555  phiR(i_shape,0) * constant;*/
1556  }
1557 }
1558 
1559 #ifdef _AUTODIFF
1560 
1561 void TPZEulerConsLaw::ContributeImplConvFace(TPZVec<REAL> &x,
1562  TPZVec<FADREAL> &solL,TPZVec<FADREAL> &solR,
1563  REAL weight,TPZVec<REAL> &normal,
1564  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
1565  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,
1566  int entropyFix)
1567 {
1568  int nState = NStateVariables();
1569  TPZVec<FADREAL > flux(nState,REAL(0.));
1570  Roe_Flux(solL, solR, normal, fGamma, flux, entropyFix);
1571 
1572  // Testing whether Roe_Flux<REAL> gives the same result as Roe_Flux<FADREAL>
1573 
1574  /* TPZVec<REAL> solL2(nState,0.),solR2(nState,0.),flux2(nState,0.);
1575  int i;
1576  for(i=0; i<nState; i++)
1577  {
1578  solL2[i] = solL[i].val();
1579  solR2[i] = solR[i].val();
1580  }
1581  Roe_Flux(solL2, solR2, normal, fGamma, flux2,with_entropy_fix);
1582  REAL diff = fabs(flux[0].val()-flux2[0])+fabs(flux[1].val()-flux2[1])+fabs(flux[2].val()-flux2[2]);
1583  if(diff != 0.)
1584  {
1585  cout << "Roe<FADREAL> is different from Roe<REAL> diff " << diff << endl;
1586  }*/
1587  int nShapeL = phiL.Rows(),
1588  nShapeR = phiR.Rows(),
1589  i_shape, i_state, j,
1590  nDer = (nShapeL + nShapeR) * nState;
1591 
1592 
1593 #ifdef DIVTIMESTEP
1594  REAL constant = weight; // weight
1595 #else
1596  REAL constant = TimeStep() * weight; // deltaT * weight
1597 #endif
1598 
1599  // Contribution referring to the left element
1600  for(i_shape = 0; i_shape < nShapeL; i_shape ++)
1601  for(i_state = 0; i_state < nState; i_state++)
1602  {
1603  int index = i_shape*nState + i_state;
1604  ef(index,0) +=
1605  flux[i_state].val() * phiL(i_shape,0) * constant;
1606  for(j = 0; j < nDer; j++)
1607  ek(index, j) -= flux[i_state].dx/*fastAccessDx*/(j) *
1608  phiL(i_shape,0) * constant;
1609  }
1610 
1611  // Contribution referring to the right element
1612  // REM: The contributions are negative in comparison
1613  // to the left elements: opposed normals
1614  for(i_shape = 0; i_shape < nShapeR; i_shape ++)
1615  for(i_state = 0; i_state < nState; i_state++)
1616  {
1617  int index = (nShapeL + i_shape)*nState + i_state;
1618  ef(index,0) -=
1619  flux[i_state].val() * phiR(i_shape,0) * constant;
1620  for(j = 0; j < nDer; j++)
1621  ek(index, j) += flux[i_state].dx/*fastAccessDx*/(j) *
1622  phiR(i_shape,0) * constant;
1623  }
1624 }
1625 
1626 void TPZEulerConsLaw::ContributeFastestImplConvFace(int dim,
1627  TPZVec<REAL> &x,
1628  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
1629  REAL weight,TPZVec<REAL> &normal,
1630  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
1631  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,
1632  int entropyFix)
1633 {
1634  switch(dim)
1635  {
1636  case(1):
1637  ContributeFastestImplConvFace_dim<1>(x, solL, solR, weight, normal,
1638  phiL, phiR, ek, ef, entropyFix);
1639  break;
1640  case(2):
1641  ContributeFastestImplConvFace_dim<2>(x, solL, solR, weight, normal,
1642  phiL, phiR, ek, ef, entropyFix);
1643  break;
1644  case(3):
1645  ContributeFastestImplConvFace_dim<3>(x, solL, solR, weight, normal,
1646  phiL, phiR, ek, ef, entropyFix);
1647  break;
1648  default:
1649  PZError << "\nTPZEulerConsLaw::ContributeFastestImplConvFace unhandled dimension\n";
1650  exit(-1);
1651  }
1652 }
1653 
1654 template <int dim>
1655 void TPZEulerConsLaw::ContributeFastestImplConvFace_dim(
1656  TPZVec<REAL> &x,
1657  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
1658  REAL weight,TPZVec<REAL> &normal,
1659  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
1660  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,
1661  int entropyFix)
1662 {
1663 #ifdef _TFAD
1664  typedef TFad<2*(dim+2), REAL> TFADREALInterface;
1665 #endif
1666 #ifdef _FAD
1667  typedef Fad<REAL> TFADREALInterface;
1668 #endif
1669 #ifdef _TINYFAD
1670  typedef TinyFad<2*(dim+2), REAL> TFADREALInterface;
1671 #endif
1672 
1673  int nstate = NStateVariables(dim);
1674  TPZVec< TFADREALInterface > FADsolL(nstate),
1675  FADsolR(nstate);
1676  PrepareFastestInterfaceFAD(solL, solR, FADsolL, FADsolR);
1677 
1678  ContributeFastestImplConvFace_T(x, FADsolL, FADsolR, weight, normal,
1679  phiL, phiR, ek, ef, entropyFix);
1680 }
1681 
1682 template <class T>
1683 void TPZEulerConsLaw::ContributeFastestImplConvFace_T(TPZVec<REAL> &x,
1684  TPZVec<T> &FADsolL,TPZVec<T> &FADsolR,
1685  REAL weight,TPZVec<REAL> &normal,
1686  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
1687  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,
1688  int entropyFix)
1689 {
1690  const int nState = NStateVariables();
1691 
1692  TPZVec<T> FADflux(nState);
1693 
1694  int nShapeL = phiL.Rows(),
1695  nShapeR = phiR.Rows(),
1696  i_shape, i_state, j, k,
1697  nDerL = nShapeL * nState;
1698 
1699 
1700 #ifdef DIVTIMESTEP
1701  REAL constant = weight; // weight
1702 #else
1703  REAL constant = TimeStep() * weight; // deltaT * weight
1704 #endif
1705 
1706 
1707  Roe_Flux(FADsolL, FADsolR, normal, fGamma, FADflux, entropyFix);
1708 
1709  // Contribution referring to the left element
1710  for(i_shape = 0; i_shape < nShapeL; i_shape ++)
1711  for(i_state = 0; i_state < nState; i_state++)
1712  {
1713  int index = i_shape*nState + i_state;
1714  ef(index,0) +=
1715  FADflux[i_state].val() * phiL(i_shape,0) * constant;
1716  for(k = 0; k < nState; k++)
1717  {
1718  for(j = 0; j < nShapeL; j++)
1719  ek(index, j * nState + k) -=
1720  FADflux[i_state].dx/*fastAccessDx*/(k) *
1721  phiL(j) * //df/dUl
1722  phiL(i_shape,0) * //test function
1723  constant;
1724  for(j = 0; j < nShapeR; j++)
1725  ek(index, j*nState + k + nDerL) -=
1726  FADflux[i_state]./*fastAccessDx*/dx(k + nState) *
1727  phiR(j) * //df/dUl
1728  phiL(i_shape,0) * //test function
1729  constant;
1730  }
1731  }
1732 
1733  // Contribution referring to the right element
1734  // REM: The contributions are negative in comparison
1735  // to the left elements: opposed normals
1736  for(i_shape = 0; i_shape < nShapeR; i_shape ++)
1737  for(i_state = 0; i_state < nState; i_state++)
1738  {
1739  int index = (nShapeL + i_shape) * nState + i_state;
1740  ef(index,0) -=
1741  FADflux[i_state].val() * phiR(i_shape,0) * constant;
1742  for(k = 0; k < nState; k++)
1743  {
1744  for(j = 0; j < nShapeL; j++)
1745  ek(index, j * nState + k) +=
1746  FADflux[i_state]./*fastAccessDx*/dx(k) *
1747  phiL(j) * //df/dUl
1748  phiR(i_shape,0) * //test function
1749  constant;
1750  for(j = 0; j < nShapeR; j++)
1751  ek(index, j * nState + k + nDerL) +=
1752  FADflux[i_state].dx/*fastAccessDx*/(k + nState) *
1753  phiR(j) * //df/dUl
1754  phiR(i_shape,0) * //test function
1755  constant;
1756  }
1757  }
1758 
1759 }
1760 
1761 #endif
1762 
1763 void TPZEulerConsLaw::ContributeExplConvVol(TPZVec<REAL> &x,
1764  TPZVec<STATE> &sol,
1765  REAL weight,
1766  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
1767  TPZFMatrix<STATE> &ef)
1768 {
1769  TPZVec< TPZVec<STATE> > F(3);
1770  Flux(sol, F[0], F[1], F[2]);
1771 
1772 
1773 #ifdef DIVTIMESTEP
1774  REAL constant = -weight; // -weight
1775 #else
1776  REAL constant = -TimeStep() * weight; // -deltaT * weight
1777 #endif
1778 
1779  int i_state, i_shape, nShape = phi.Rows(), k;
1780  int nState = NStateVariables();
1781 
1782  for(i_shape=0; i_shape<nShape; i_shape++)
1783  for(i_state = 0; i_state<nState; i_state++)
1784  {
1785  int index = i_state + i_shape * nState;
1786  for(k=0; k<fDim; k++)
1787  ef(index,0) += (F[k])[i_state] * dphi(k, i_shape) * constant;
1788  // ef(index) += F<scalar>GradPhi
1789  }
1790 
1791 }
1792 
1793 
1794 void TPZEulerConsLaw::ContributeImplConvVol(TPZVec<REAL> &x,
1795  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
1796  REAL weight,
1797  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
1798  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
1799 {
1800  TPZVec< TPZVec<STATE> > F(3);
1801  Flux(sol, F[0], F[1], F[2]);
1802 
1803 #ifdef DIVTIMESTEP
1804  REAL constant = -weight; // -weight
1805 #else
1806  REAL constant = -TimeStep() * weight; // -deltaT * weight
1807 #endif
1808 
1809  TPZVec< TPZDiffMatrix<STATE> > Ai(3);
1810  JacobFlux(fGamma, fDim, sol, Ai);
1811  int j_shape, j_state;
1812  int i_state, i_shape, nShape = phi.Rows(), k;
1813  int nState = NStateVariables();
1814 
1815  for(i_shape=0; i_shape<nShape; i_shape++)
1816  for(i_state = 0; i_state<nState; i_state++)
1817  {
1818  int index = i_state + i_shape * nState;
1819  for(k=0; k<fDim; k++)
1820  {
1821  // ef(index) += F<scalar>GradPhi
1822  ef(index,0) += (F[k])[i_state] * dphi(k, i_shape)
1823  * constant;
1824  for(j_shape = 0; j_shape < nShape; j_shape++)
1825  for(j_state = 0; j_state < nState; j_state++)
1826  ek(index, j_state + j_shape * nState) -=
1827  Ai[k](i_state, j_state) *
1828  dphi(k, i_shape) *
1829  phi(j_shape,0) *
1830  constant;
1831  // dConvVol/dU* =
1832  // dConvVol/dU * dU/dU*
1833  // dConvVol/dU = Ai<scalar>GradPhi
1834  // dU/dU* = phi(j)
1835  }
1836  }
1837 }
1838 
1839 
1840 void TPZEulerConsLaw::ContributeExplT1(TPZVec<REAL> &x,
1841  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
1842  REAL weight,
1843  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
1844  TPZFMatrix<STATE> &ef)
1845 {
1846  if(TimeStep()==0.)
1847  {
1848  std::cout << __PRETTY_FUNCTION__ << " Zero timestep bailing out\n";
1849  return;
1850  }
1851 
1852  int i_shape, ij_state;
1853  int nState = NStateVariables();
1854  int nShape = phi.Rows();
1855 
1856 #ifdef DIVTIMESTEP
1857  REAL constant = weight / TimeStep();
1858 #else
1859  REAL constant = weight;
1860 #endif
1861 
1862  for(i_shape = 0; i_shape < nShape; i_shape++)
1863  for(ij_state = 0; ij_state < nState; ij_state++)
1864  {
1865  int index = i_shape * nState + ij_state;
1866  // ef += sol*phi(i)
1867  ef(index, 0) +=
1868  sol[ij_state] * phi(i_shape,0) *
1869  constant;
1870  }
1871 }
1872 
1873 void TPZEulerConsLaw::ContributeImplT1(TPZVec<REAL> &x,
1874  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
1875  REAL weight,
1876  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
1877  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
1878 {
1879  if(TimeStep()==0.)
1880  {
1881  std::cout << __PRETTY_FUNCTION__ << " Zero timestep bailing out\n";
1882  return;
1883  }
1884 
1885  int i_shape, ij_state, j_shape;
1886  int nState = NStateVariables();
1887  int nShape = phi.Rows();
1888 
1889 #ifdef DIVTIMESTEP
1890  REAL constant = weight / TimeStep();
1891 #else
1892  REAL constant = weight;
1893 #endif
1894 
1895  for(i_shape = 0; i_shape < nShape; i_shape++)
1896  for(ij_state = 0; ij_state < nState; ij_state++)
1897  {
1898  int index = i_shape * nState + ij_state;
1899  // ef += sol*phi(i)
1900  ef(index, 0) +=
1901  sol[ij_state] * phi(i_shape,0) *
1902  constant;
1903  // ek += phi(i)*phi(j)
1904  for(j_shape = 0; j_shape < nShape; j_shape++)
1905  ek(index, j_shape * nState + ij_state) -=
1906  phi(i_shape,0) *
1907  phi(j_shape,0) *
1908  constant;
1909  }
1910 }
1911 
1912 void TPZEulerConsLaw::ContributeExplT2(TPZVec<REAL> &x,
1913  TPZVec<STATE> &sol,
1914  REAL weight,
1915  TPZFMatrix<REAL> &phi,
1916  TPZFMatrix<STATE> &ef)
1917 {
1918  if(TimeStep()==0.)
1919  {
1920  std::cout << __PRETTY_FUNCTION__ << " Zero timestep bailing out\n";
1921  return;
1922  }
1923 
1924  int i_shape, i_state;
1925  int nState = NStateVariables();
1926  int nShape = phi.Rows();
1927 
1928 #ifdef DIVTIMESTEP
1929  REAL constant = weight / TimeStep();
1930 #else
1931  REAL constant = weight;
1932 #endif
1933 
1934  for(i_shape = 0; i_shape < nShape; i_shape++)
1935  for(i_state = 0; i_state < nState; i_state++)
1936  {
1937  int index = i_shape * nState + i_state;
1938  // ef += sol*phi(i)
1939  ef(index, 0) -=
1940  sol[i_state] * phi(i_shape,0) *
1941  constant; // the T2 parcell is negative
1942  }
1943 }
1944 
1945 
1946 void TPZEulerConsLaw::Write(TPZStream &buf, int withclassid) const
1947 {
1948  TPZConservationLaw::Write(buf, 0);
1949  fArtDiff.Write(buf, 0);
1950  int tmp = static_cast < int > (fDiff);
1951  buf.Write(& tmp,1);
1952  tmp = static_cast<int>(fConvVol);
1953  buf.Write(& tmp,1);
1954  tmp = static_cast<int>(fConvFace);
1955  buf.Write(&tmp,1);
1956 }
1957 
1958 void TPZEulerConsLaw::Read(TPZStream &buf, void *context)
1959 {
1960  TPZConservationLaw::Read(buf, context);
1961  fArtDiff.Read(buf, context);
1962  int diff;
1963  buf.Read(&diff,1);
1964  fDiff = static_cast<TPZTimeDiscr>(diff);
1965  buf.Read(&diff,1);
1966  fConvVol = static_cast<TPZTimeDiscr>(diff);
1967  buf.Read(&diff,1);
1968  fConvFace = static_cast<TPZTimeDiscr>(diff);
1969 }
1970 
1983 template <class T>
1984 void TPZEulerConsLaw::ApproxRoe_Flux(TPZVec<T> &solL, TPZVec<T> &solR,
1985  TPZVec<REAL> & normal, REAL gamma,
1986  TPZVec<T> & flux, int entropyFix)
1987 {
1988  // Normals outgoing from the BC elements into the
1989  // mesh elements -> all the normals are opposited to
1990  // the common convention -> changing the left/right
1991  // elements and normals.
1992  int nState = solL.NElements();
1993  if(nState == 5)
1994  {
1995  ApproxRoe_Flux<T>(solL[0], solL[1], solL[2], solL[3], solL[4],
1996  solR[0], solR[1], solR[2], solR[3], solR[4],
1997  normal[0], normal[1], normal[2],
1998  gamma,
1999  flux[0], flux[1], flux[2], flux[3], flux[4], entropyFix);
2000 
2001  }else if(nState == 4)
2002  {
2003  ApproxRoe_Flux<T>(solL[0], solL[1], solL[2], solL[3],
2004  solR[0], solR[1], solR[2], solR[3],
2005  normal[0], normal[1],
2006  gamma,
2007  flux[0], flux[1], flux[2], flux[3], entropyFix);
2008  }else if(nState == 3)
2009  {
2010  //using the 2D expression for 1d problem
2011  T auxL = REAL(0.),
2012  auxR = REAL(0.),
2013  fluxaux = REAL(0.);
2014  auxL = flux[0];
2015  ApproxRoe_Flux<T>(solL[0], solL[1], auxL, solL[2],
2016  solR[0], solR[1], auxR, solR[2],
2017  normal[0], 0,
2018  gamma,
2019  flux[0], flux[1], fluxaux, flux[2], entropyFix);
2020  }else
2021  {
2022  PZError << "No flux on " << nState << " state variables.\n";
2023  }
2024 }
2025 
2026 #ifdef _AUTODIFF
2027 
2030 template <>
2031 void TPZEulerConsLaw::ApproxRoe_Flux<FADREAL>(const FADREAL & rho_f,
2032  const FADREAL & rhou_f,
2033  const FADREAL & rhov_f,
2034  const FADREAL & rhow_f,
2035  const FADREAL & rhoE_f,
2036  const FADREAL & rho_t,
2037  const FADREAL & rhou_t,
2038  const FADREAL & rhov_t,
2039  const FADREAL & rhow_t,
2040  const FADREAL & rhoE_t,
2041  const REAL nx,
2042  const REAL ny,
2043  const REAL nz,
2044  const REAL gam,
2045  FADREAL & flux_rho,
2046  FADREAL & flux_rhou,
2047  FADREAL & flux_rhov,
2048  FADREAL & flux_rhow,
2049  FADREAL & flux_rhoE, int entropyFix)
2050 {
2051 
2052  typedef FADREAL T;
2053  // REAL alpha1,alpha2,alpha3,alpha4,alpha5,alpha;
2054  REAL a1,a2,a3,a4,a5,b1,b2,b3,b4,b5;
2055  // REAL ep_t, ep_f, p_t, p_f;
2056  // REAL rhouv_t, rhouv_f, rhouw_t, rhouw_f, rhovw_t, rhovw_f;
2057  // REAL lambda_f, lambda_t;
2058  T delta_rho, delta_rhou, delta_rhov, delta_rhow, delta_rhoE;
2059  REAL hnx, hny, hnz;
2060  REAL tempo11, usc;
2061 
2062  flux_rho = 0;
2063  flux_rhou = 0;
2064  flux_rhov = 0;
2065  flux_rhow = 0;
2066  flux_rhoE = 0;
2067 
2068  REAL gam1 = gam - 1.0;
2069  REAL irho_f = REAL(1.0)/rho_f.val();
2070  REAL irho_t = REAL(1.0)/rho_t.val();
2071 
2072  //
2073  //.. Compute the ROE Averages
2074  //
2075  //.... some useful quantities
2076  REAL coef1 = sqrt(rho_f.val());
2077  REAL coef2 = sqrt(rho_t.val());
2078  REAL somme_coef = coef1 + coef2;
2079  REAL isomme_coef = REAL(1.0)/somme_coef;
2080  REAL u_f = rhou_f.val()*irho_f;
2081  REAL v_f = rhov_f.val()*irho_f;
2082  REAL w_f = rhow_f.val()*irho_f;
2083  REAL h_f = (gam * rhoE_f.val()*irho_f) - (.5*gam1) * (u_f * u_f + v_f * v_f + w_f * w_f);
2084  REAL u_t = rhou_t.val()*irho_t;
2085  REAL v_t = rhov_t.val()*irho_t;
2086  REAL w_t = rhow_t.val()*irho_t;
2087  REAL h_t = (gam * rhoE_t.val()*irho_t) - (.5*gam1) * (u_t * u_t + v_t * v_t + w_t * w_t);
2088 
2089  //.... averages
2090  //REAL rho_ave = coef1 * coef2;
2091  REAL u_ave = (coef1 * u_f + coef2 * u_t) * isomme_coef;
2092  REAL v_ave = (coef1 * v_f + coef2 * v_t) * isomme_coef;
2093  REAL w_ave = (coef1 * w_f + coef2 * w_t) * isomme_coef;
2094  REAL h_ave = (coef1 * h_f + coef2 * h_t) * isomme_coef;
2095  //
2096  //.. Compute Speed of sound
2097  REAL scal = u_ave * nx + v_ave * ny + w_ave * nz;
2098  REAL norme = sqrt(nx * nx + ny * ny + nz * nz);
2099  REAL inorme = REAL(1.0)/norme;
2100  REAL u2pv2pw2 = u_ave * u_ave + v_ave * v_ave + w_ave * w_ave;
2101  REAL c_speed = gam1 * (h_ave - REAL(0.5) * u2pv2pw2);
2102  if(c_speed < REAL(1e-6)) c_speed = 1e-6;// <!> zeroes the derivatives? // avoid division by 0 if critical
2103  c_speed = sqrt(c_speed);
2104  REAL c_speed2 = c_speed * norme;
2105  //
2106  //.. Compute the eigenvalues of the Jacobian matrix
2107  REAL eig_val1 = scal - c_speed2;
2108  REAL eig_val2 = scal;
2109  REAL eig_val3 = scal + c_speed2;
2110  //
2111  //.. Compute the ROE flux
2112  //.... In this part many tests upon the eigenvalues
2113  //.... are done to simplify calculations
2114  //.... Here we use the two formes of the ROE flux :
2115  //.... phi(Wl,Wr) = F(Wl) + A-(Wroe)(Wr - Wl)
2116  //.... phi(Wl,Wr) = F(Wr) - A+(Wroe)(Wr - Wl)
2117  //
2118  if(eig_val2 <= REAL(0.0)) {
2119  T irho_t = REAL(1.0)/rho_t;
2120  T u_t = rhou_t*irho_t;
2121  T v_t = rhov_t*irho_t;
2122  T w_t = rhow_t*irho_t;
2123  T h_t = (gam * rhoE_t*irho_t) - (.5*gam1) * (u_t * u_t + v_t * v_t + w_t * w_t);
2124 
2125  T p_t,ep_t,rhouv_t,rhouw_t,rhovw_t;
2126  REAL lambda_f,lambda_t;
2127  p_t = gam1 * (rhoE_t - REAL(0.5) * (rhou_t * rhou_t +
2128  rhov_t * rhov_t + rhow_t * rhow_t) * irho_t);
2129  ep_t = rhoE_t + p_t;
2130  rhouv_t = rhou_t * v_t;
2131  rhouw_t = rhou_t * w_t;
2132  rhovw_t = rhov_t * w_t;
2133  flux_rho = rhou_t * nx + rhov_t * ny + rhow_t * nz;
2134  flux_rhou = (rhou_t * u_t + p_t) * nx + rhouv_t * ny + rhouw_t * nz;
2135  flux_rhov = rhouv_t * nx + (rhov_t * v_t + p_t) * ny + rhovw_t * nz;
2136  flux_rhow = rhouw_t * nx + rhovw_t * ny + (rhow_t * w_t + p_t) * nz;
2137  flux_rhoE = ep_t * (u_t * nx + v_t * ny + w_t * nz);
2138  //
2139  //.... A Entropic modification
2140  //
2141  REAL p_f = gam1 * (rhoE_f.val() - REAL(0.5) * (rhou_f.val() * rhou_f.val() + rhov_f.val() * rhov_f.val()
2142  + rhow_f.val() * rhow_f.val()) * irho_f);
2143  lambda_f = u_f * nx + v_f * ny + w_f * nz + norme
2144  * sqrt(gam * p_f * irho_f);
2145  lambda_t = u_t.val() * nx + v_t.val() * ny + w_t.val() * nz + norme
2146  * sqrt(gam * p_t.val() * irho_t.val());
2147  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
2148  eig_val3 = lambda_t * (eig_val3 - lambda_f) / (lambda_t - lambda_f);
2149  }
2150  //
2151  if (eig_val3 > REAL(0.0)) {
2152  //.. In this case A+ is obtained by multiplying the last
2153  //.. colomne of T-1 with the last row of T with eig_val3 //Cedric
2154  T alpha1,alpha2,alpha3,alpha4,alpha5,alpha;
2155  delta_rho = rho_t - rho_f; //right - left
2156  delta_rhou = rhou_t - rhou_f; //**_t - **_f
2157  delta_rhov = rhov_t - rhov_f;
2158  delta_rhow = rhow_t - rhow_f;
2159  delta_rhoE = rhoE_t - rhoE_f;
2160  //
2161  scal = scal * inorme;
2162  hnx = nx * inorme;
2163  hny = ny * inorme;
2164  hnz = nz * inorme;
2165  usc = REAL(1.0)/c_speed;
2166  tempo11 = gam1 * usc;
2167  //.. Last columne of the matrix T-1
2168  a1 = usc;
2169  a2 = u_ave * usc + hnx;
2170  a3 = v_ave * usc + hny;
2171  a4 = w_ave * usc + hnz;
2172  a5 = REAL(0.5) * u2pv2pw2 * usc + REAL(2.5) * c_speed + scal;
2173  //.. Last row of the matrix T * eig_val3
2174  b1 = REAL(0.5) * (REAL(0.5) * tempo11 * u2pv2pw2 - scal);
2175  b2 = REAL(0.5) * (hnx - tempo11 * u_ave);
2176  b3 = REAL(0.5) * (hny - tempo11 * v_ave);
2177  b4 = REAL(0.5) * (hnz - tempo11 * w_ave);
2178  b5 = REAL(0.5) * tempo11;
2179  //
2180  alpha1 = b1 * delta_rho;
2181  alpha2 = b2 * delta_rhou;
2182  alpha3 = b3 * delta_rhov;
2183  alpha4 = b4 * delta_rhow;
2184  alpha5 = b5 * delta_rhoE;
2185  alpha = eig_val3 * (alpha1 + alpha2 + alpha3 + alpha4 + alpha5);
2186  //
2187  flux_rho -= a1 * alpha;
2188  flux_rhou -= a2 * alpha;
2189  flux_rhov -= a3 * alpha;
2190  flux_rhow -= a4 * alpha;
2191  flux_rhoE -= a5 * alpha;
2192  }
2193  }
2194  //
2195  if(eig_val2 > REAL(0.0)) {
2196  T p_f,ep_f,rhouv_f,rhovw_f,rhouw_f;
2197  T irho_f = REAL(1.0)/rho_f.val();
2198 
2199  T u_f = rhou_f.val()*irho_f;
2200  T v_f = rhov_f.val()*irho_f;
2201  T w_f = rhow_f.val()*irho_f;
2202  T h_f = (gam * rhoE_f.val()*irho_f) - (.5*gam1) * (u_f * u_f + v_f * v_f + w_f * w_f);
2203  p_f = gam1 * (rhoE_f - REAL(0.5) * (rhou_f * rhou_f +
2204  rhov_f * rhov_f + rhow_f * rhow_f) * irho_f);
2205  ep_f = rhoE_f + p_f;
2206  rhouv_f = rhou_f * v_f;
2207  rhouw_f = rhou_f * w_f;
2208  rhovw_f = rhov_f * w_f;
2209  flux_rho = rhou_f * nx + rhov_f * ny + rhow_f * nz;
2210  flux_rhou = (rhou_f * u_f + p_f) * nx + rhouv_f * ny + rhouw_f * nz;
2211  flux_rhov = rhouv_f * nx + (rhov_f * v_f + p_f) * ny + rhovw_f * nz;
2212  flux_rhow = rhouw_f * nx + rhovw_f * ny + (rhow_f * w_f + p_f) * nz;
2213  flux_rhoE = ep_f * (u_f * nx + v_f * ny + w_f * nz);
2214  //
2215  // A Entropic modification
2216  //
2217  REAL p_t = gam1 * (rhoE_t.val() - REAL(0.5) * (rhou_t.val() * rhou_t.val() +
2218  rhov_t.val() * rhov_t.val() + rhow_t.val() * rhow_t.val()) * irho_t);
2219  REAL lambda_f = u_f.val() * nx + v_f.val() * ny + w_f.val() * nz - norme
2220  * sqrt(gam * p_f.val() * irho_f.val());
2221  REAL lambda_t = u_t * nx + v_t * ny + w_t * nz - norme
2222  * sqrt(gam * p_t * irho_t);
2223  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
2224  eig_val1 = lambda_f * (lambda_t - eig_val1) / (lambda_t - lambda_f);
2225  }
2226  //
2227  if (eig_val1 < REAL(0.0)) {
2228  //.. In this case A+ is obtained by multiplying the first
2229  //.. columne of T-1 with the first row of T with eig_val1
2230  T alpha1,alpha2,alpha3,alpha4,alpha5,alpha;
2231  delta_rho = rho_t - rho_f;
2232  delta_rhou = rhou_t - rhou_f;
2233  delta_rhov = rhov_t - rhov_f;
2234  delta_rhow = rhow_t - rhow_f;
2235  delta_rhoE = rhoE_t - rhoE_f;
2236  //
2237  scal = scal * inorme;
2238  hnx = nx * inorme;
2239  hny = ny * inorme;
2240  hnz = nz * inorme;
2241  usc = REAL(1.0)/c_speed;
2242  tempo11 = gam1 * usc;
2243  //.. First colomne of the matrix T-1
2244  a1 = usc;
2245  a2 = u_ave * usc - hnx;
2246  a3 = v_ave * usc - hny;
2247  a4 = w_ave * usc - hnz;
2248  a5 = REAL(0.5) * u2pv2pw2 * usc + REAL(2.5) * c_speed - scal;
2249  //.. First row of the matrix T * eig_val1
2250  b1 = REAL(0.5) * (REAL(0.5) * tempo11 * u2pv2pw2 + scal);
2251  b2 = -REAL(0.5) * (hnx + tempo11 * u_ave);
2252  b3 = -REAL(0.5) * (hny + tempo11 * v_ave);
2253  b4 = -REAL(0.5) * (hnz + tempo11 * w_ave);
2254  b5 = REAL(0.5) * tempo11;
2255  //
2256  alpha1 = b1 * delta_rho;
2257  alpha2 = b2 * delta_rhou;
2258  alpha3 = b3 * delta_rhov;
2259  alpha4 = b4 * delta_rhow;
2260  alpha5 = b5 * delta_rhoE;
2261  alpha = eig_val1 * (alpha1 + alpha2 + alpha3 + alpha4 + alpha5);
2262  //
2263  flux_rho += a1 * alpha;
2264  flux_rhou += a2 * alpha;
2265  flux_rhov += a3 * alpha;
2266  flux_rhow += a4 * alpha;
2267  flux_rhoE += a5 * alpha;
2268  }
2269  }
2270 }
2271 /*
2272  REAL operator() (const FADREAL &a)
2273  {
2274  return a.val();
2275  }*/
2276 
2277 template <>
2278 void TPZEulerConsLaw::ApproxRoe_Flux<FADREAL>(const FADREAL & rho_f,
2279  const FADREAL & rhou_f,
2280  const FADREAL & rhov_f,
2281  const FADREAL & rhoE_f,
2282  const FADREAL & rho_t,
2283  const FADREAL & rhou_t,
2284  const FADREAL & rhov_t,
2285  const FADREAL & rhoE_t,
2286  const REAL nx,
2287  const REAL ny,
2288  const REAL gam,
2289  FADREAL &flux_rho,
2290  FADREAL &flux_rhou,
2291  FADREAL &flux_rhov,
2292  FADREAL &flux_rhoE, int entropyFix)
2293 {
2294  typedef FADREAL T;
2295  typedef REAL locREAL;
2296  locREAL rho_fv, rhou_fv, rhov_fv, rhoE_fv, rho_tv, rhou_tv, rhov_tv, rhoE_tv;
2297 
2298  rho_fv = rho_f.val();
2299  rhou_fv = rhou_f.val();
2300  rhov_fv = rhov_f.val();
2301  rhoE_fv = rhoE_f.val();
2302  rho_tv = rho_t.val();
2303  rhou_tv = rhou_t.val();
2304  rhov_tv = rhov_t.val();
2305  rhoE_tv = rhoE_t.val();
2306 
2307 
2308  // rho_fv = rho_f;
2309  // rhou_fv = rhou_f;
2310  // rhov_fv = rhov_f;
2311  // rhoE_fv = rhoE_f;
2312  // rho_tv = rho_t;
2313  // rhou_tv = rhou_t;
2314  // rhov_tv = rhov_t;
2315  // rhoE_tv = rhoE_t;
2316 
2317  locREAL c_speed2;
2318  //
2319  //.. Compute the eigenvalues of the Jacobian matrix
2320  locREAL eig_val1;
2321  locREAL eig_val2;
2322  locREAL eig_val3;
2323  locREAL u_fv = rhou_fv/rho_fv;
2324  locREAL v_fv = rhov_fv/rho_fv;
2325  locREAL u_tv = rhou_tv/rho_tv;
2326  locREAL v_tv = rhov_tv/rho_tv;
2327  REAL gam1 = gam - REAL(1.0);
2328  locREAL p_tv = gam1 * (rhoE_tv - REAL(0.5) * (rhou_tv * rhou_tv + rhov_tv * rhov_tv) / rho_tv);
2329  locREAL p_fv = gam1 * (rhoE_fv - REAL(0.5) * (rhou_fv * rhou_fv + rhov_fv * rhov_fv) / rho_fv);
2330  REAL norme = sqrt(nx * nx + ny * ny);
2331  locREAL scal,c_speed;
2332  locREAL u_ave,v_ave,u2pv2;
2333  {
2334 
2335  flux_rho = 0;
2336  flux_rhou = 0;
2337  flux_rhov = 0;
2338  flux_rhoE = 0;
2339 
2340  //REAL gam2 = gam * (gam - 1.0);
2341  //REAL igam = 1.0 / (gam - 1.0);
2342 
2343  //
2344  //.. Compute the ROE Averages
2345  //
2346  //.... some useful quantities
2347  locREAL coef1 = sqrt(rho_fv);
2348  locREAL coef2 = sqrt(rho_tv);
2349  locREAL somme_coef = coef1 + coef2;
2350  locREAL h_f = (gam * rhoE_fv/rho_fv) - (u_fv * u_fv + v_fv * v_fv) * (gam1 / REAL(2.0));
2351  locREAL h_t = (gam * rhoE_tv/rho_tv) - (u_tv * u_tv + v_tv * v_tv) * (gam1 / REAL(2.0));
2352 
2353  //.... averages
2354  //REAL rho_ave = coef1 * coef2;
2355  // cout << "Approx coef1 " << coef1 << "coef2 " << coef2 << "h_f " << h_f << "h_t " << h_t << endl;
2356  u_ave = (coef1 * u_fv + coef2 * u_tv) / somme_coef;
2357  v_ave = (coef1 * v_fv + coef2 * v_tv) / somme_coef;
2358  locREAL h_ave = (coef1 * h_f + coef2 * h_t) / somme_coef;
2359  //
2360  //.. Compute Speed of sound
2361  scal = u_ave * nx + v_ave * ny;
2362  u2pv2 = u_ave * u_ave + v_ave * v_ave;
2363  c_speed = gam1 * (h_ave - REAL(0.5) * u2pv2);
2364  // cout << "Approx c_speed " << c_speed << endl << "h_ave " << h_ave << endl << "u2pv2 " << u2pv2 << endl;
2365  if(c_speed < REAL(1e-6)) c_speed = REAL(1e-6); // avoid division by 0 if critical
2366  c_speed = sqrt(c_speed);
2367  c_speed2 = c_speed * norme;
2368  //
2369  //.. Compute the eigenvalues of the Jacobian matrix
2370  eig_val1 = scal - c_speed2;
2371  eig_val2 = scal;
2372  eig_val3 = scal + c_speed2;
2373  // cout << "Eigenvalues appox" << eig_val1 << endl << eig_val2 << endl << eig_val3 << endl;
2374  }
2375  //
2376  //.. Compute the ROE flux
2377  //.... In this part many tests upon the eigenvalues
2378  //.... are done to simplify calculations
2379  //.... Here we use the two formes of the ROE flux :
2380  //.... phi(Wl,Wr) = F(Wl) + A-(Wroe)(Wr - Wl)
2381  //.... phi(Wl,Wr) = F(Wr) - A+(Wroe)(Wr - Wl)
2382  //
2383  if(eig_val2 <= REAL(0.0)) {
2384  T p_t,ep_t;
2385  p_t = gam1 * (rhoE_t - REAL(0.5) * (rhou_t * rhou_t + rhov_t * rhov_t) / rho_t);
2386  ep_t = rhoE_t + p_t;
2387  T u_t = rhou_t/rho_t;
2388  T v_t = rhov_t/rho_t;
2389  T rhouv_t = rhou_t * v_t;
2390  flux_rho = rhou_t * nx + rhov_t * ny;
2391  flux_rhou = (rhou_t * u_t + p_t) * nx + rhouv_t * ny;
2392  flux_rhov = rhouv_t * nx + (rhov_t * v_t + p_t) * ny;
2393  flux_rhoE = ep_t * (u_t * nx + v_t * ny);
2394  //
2395  //.... A Entropic modification
2396  //
2397  locREAL lambda_f = u_fv * nx + v_fv * ny + norme * sqrt(gam * p_fv / rho_fv);
2398  locREAL lambda_t = u_tv * nx + v_tv * ny + norme
2399  * sqrt(gam * p_tv / rho_tv);
2400  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
2401  eig_val3 = lambda_t * (eig_val3 - lambda_f) / (lambda_t - lambda_f);
2402  }
2403  //
2404  if (eig_val3 > REAL(0.0)) {
2405  //.. In this case A+ is obtained by multiplying the last
2406  //.. colomne of T-1 with the last row of T with eig_val3
2407  T alpha1,alpha2,alpha3,alpha4,alpha;
2408  locREAL a1,a2,a3,a4,b1,b2,b3,b4;
2409  T delta_rho, delta_rhou,delta_rhov, delta_rhoE;
2410  delta_rho = rho_t - rho_f;
2411  delta_rhou = rhou_t - rhou_f;
2412  delta_rhov = rhov_t - rhov_f;
2413  delta_rhoE = rhoE_t - rhoE_f;
2414  //
2415  scal = scal / norme;
2416  REAL hnx = nx / norme;
2417  REAL hny = ny / norme;
2418  locREAL usc = REAL(1.0)/c_speed;
2419  locREAL tempo11 = gam1 * usc;
2420  //.. Last columne of the matrix T-1
2421  a1 = usc;
2422  a2 = u_ave * usc + hnx;
2423  a3 = v_ave * usc + hny;
2424  a4 = REAL(0.5) * u2pv2 * usc + REAL(2.5) * c_speed + scal;
2425  //.. Last row of the matrix T * eig_val3
2426  b1 = REAL(0.5) * eig_val3 * (REAL(0.5) * tempo11 * u2pv2 - scal);
2427  b2 = REAL(0.5) * eig_val3 * (hnx - tempo11 * u_ave);
2428  b3 = REAL(0.5) * eig_val3 * (hny - tempo11 * v_ave);
2429  b4 = REAL(0.5) * eig_val3 * tempo11;
2430  //
2431  alpha1 = a1 * b1 * delta_rho;
2432  alpha2 = a1 * b2 * delta_rhou;
2433  alpha3 = a1 * b3 * delta_rhov;
2434  alpha4 = a1 * b4 * delta_rhoE;
2435  alpha = alpha1 + alpha2 + alpha3 + alpha4;
2436  //
2437  flux_rho -= alpha;
2438  flux_rhou -= a2 * b1 * delta_rho + a2 * b2 * delta_rhou +
2439  a2 * b3 * delta_rhov + a2 * b4 * delta_rhoE;
2440  flux_rhov -= a3 * b1 * delta_rho + a3 * b2 * delta_rhou +
2441  a3 * b3 * delta_rhov + a3 * b4 * delta_rhoE;
2442  flux_rhoE -= a4 * b1 * delta_rho + a4 * b2 * delta_rhou +
2443  a4 * b3 * delta_rhov + a4 * b4 * delta_rhoE;
2444  }
2445  }
2446  //
2447  if(eig_val2 > REAL(0.0)) {
2448  T p_f,ep_f;
2449  T u_f = rhou_f/rho_f;
2450  T v_f = rhov_f/rho_f;
2451  p_f = gam1 * (rhoE_f - REAL(0.5) * (rhou_f * rhou_f +
2452  rhov_f * rhov_f) / rho_f);
2453  ep_f = rhoE_f + p_f;
2454  T rhouv_f = rhou_f * v_f;
2455  flux_rho = rhou_f * nx + rhov_f * ny;
2456  flux_rhou = (rhou_f * u_f + p_f) * nx + rhouv_f * ny;
2457  flux_rhov = rhouv_f * nx + (rhov_f * v_f + p_f) * ny;
2458  flux_rhoE = ep_f * (u_f * nx + v_f * ny);
2459  //
2460  // A Entropic modification
2461  //
2462  locREAL lambda_f = u_fv * nx + v_fv * ny - norme * sqrt(gam * p_fv / rho_fv);
2463  locREAL lambda_t = u_tv * nx + v_tv * ny - norme * sqrt(gam * p_tv / rho_tv);
2464  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
2465  eig_val1 = lambda_f * (lambda_t - eig_val1) / (lambda_t - lambda_f);
2466  }
2467  //
2468  if (eig_val1 < REAL(0.0)) {
2469  //.. In this case A+ is obtained by multiplying the first
2470  //.. columne of T-1 with the first row of T with eig_val1
2471  T alpha1,alpha2,alpha3,alpha4,alpha;
2472  locREAL a1,a2,a3,a4,b1,b2,b3,b4;
2473  T delta_rho, delta_rhou,delta_rhov, delta_rhoE;
2474  delta_rho = rho_t - rho_f;
2475  delta_rhou = rhou_t - rhou_f;
2476  delta_rhov = rhov_t - rhov_f;
2477  delta_rhoE = rhoE_t - rhoE_f;
2478  //
2479  scal = scal / norme;
2480  REAL hnx = nx / norme;
2481  REAL hny = ny / norme;
2482  locREAL usc = REAL(1.0)/c_speed;
2483  locREAL tempo11 = gam1 * usc;
2484  //.. First colomne of the matrix T-1
2485  a1 = usc;
2486  a2 = u_ave * usc - hnx;
2487  a3 = v_ave * usc - hny;
2488  a4 = REAL(0.5) * u2pv2 * usc + REAL(2.5) * c_speed - scal;
2489  //.. First row of the matrix T * eig_val1
2490  b1 = REAL(0.5) * eig_val1 * (REAL(0.5) * tempo11 * u2pv2 + scal);
2491  b2 = -REAL(0.5) * eig_val1 * (hnx + tempo11 * u_ave);
2492  b3 = -REAL(0.5) * eig_val1 * (hny + tempo11 * v_ave);
2493  b4 = REAL(0.5) * eig_val1 * tempo11;
2494  //
2495  alpha1 = a1 * b1 * delta_rho;
2496  alpha2 = a1 * b2 * delta_rhou;
2497  alpha3 = a1 * b3 * delta_rhov;
2498  alpha4 = a1 * b4 * delta_rhoE;
2499  alpha = alpha1 + alpha2 + alpha3 + alpha4;
2500  //
2501  flux_rho += alpha;
2502  flux_rhou += a2 * b1 * delta_rho + a2 * b2 * delta_rhou +
2503  a2 * b3 * delta_rhov + a2 * b4 * delta_rhoE;
2504  flux_rhov += a3 * b1 * delta_rho + a3 * b2 * delta_rhou +
2505  a3 * b3 * delta_rhov + a3 * b4 * delta_rhoE;
2506  flux_rhoE += a4 * b1 * delta_rho + a4 * b2 * delta_rhou +
2507  a4 * b3 * delta_rhov + a4 * b4 * delta_rhoE;
2508  }
2509  }
2510 }
2511 #endif
2512 
2515 int TPZEulerConsLaw::ClassId() const{
2516  return Hash("TPZEulerConsLaw") ^ TPZConservationLaw::ClassId() << 1;
2517 }
2518 
2519 
TPZEulerConsLaw::ContributeAdv
virtual void ContributeAdv(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:502
TPZMaterial::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout)
Returns the solution associated with the var index based on the finite element approximation.
Definition: TPZMaterial.cpp:200
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZEulerConsLaw
This material implements the weak statement of the compressible euler equations.
Definition: pzeulerconslaw.h:33
TPZConservationLaw::Read
void Read(TPZStream &buf, void *context) override
Read the element data from a stream.
Definition: pzconslaw.cpp:83
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
TPZConservationLaw::fResidualType
TPZResidualType fResidualType
Variable to indicate the type of residual to be computed by Assemble.
Definition: pzconslaw.h:268
TPZEulerConsLaw::JacobFlux
static void JacobFlux(REAL gamma, int dim, TPZVec< T > &U, TPZVec< TPZDiffMatrix< T > > &Ai)
Jacobian of the tensor flux of Euler.
Definition: pzeulerconslaw.h:615
TPZSavable.h
Contains declaration of the TPZSavable class which defines the interface to save and restore objects ...
TPZMaterialData::normal
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
Definition: pzmaterialdata.h:69
TPZEulerConsLaw::VariableIndex
virtual int VariableIndex(const std::string &name) override
Returns the relative index of a variable according to its name.
Definition: pzeulerconslaw.cpp:157
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZConservationLaw::SetTimeStep
void SetTimeStep(REAL timeStep)
Sets the time step used for time integration.
Definition: pzconslaw.h:302
TPZEulerConsLaw::ContributeBCInterface
virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to stiffness matrix and load vector at one BC integration point...
Definition: pzeulerconslaw.cpp:702
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
LeastSquares_AD
Use Least squares method to applied artificial diffusion term.
Definition: pzartdiff.h:46
TPZEulerConsLaw::ClassId
virtual int ClassId() const override
Class identificator.
Definition: pzeulerconslaw.cpp:2515
TPZEulerConsLaw::fConvFace
TPZTimeDiscr fConvFace
Definition: pzeulerconslaw.h:544
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZArtDiff::ContributeApproxImplDiff
void ContributeApproxImplDiff(int dim, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, TPZFMatrix< REAL > &dphix, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, REAL weight, REAL timeStep, REAL deltaX)
Contributes the diffusion term to the tangent matrix (ek-approximated) and residual vector (ef) ...
Definition: pzartdiff.cpp:675
TPZEulerConsLaw::ContributeExplConvFace
void ContributeExplConvFace(TPZVec< REAL > &x, TPZVec< STATE > &solL, TPZVec< STATE > &solR, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &phiR, TPZFMatrix< STATE > &ef, int entropyFix=1)
Definition: pzeulerconslaw.cpp:1347
ApproxImplicit_TD
Semi implicit time discretization.
Definition: pzconslaw.h:38
TPZEulerConsLaw::ComputeGhostState
void ComputeGhostState(TPZVec< T > &solL, TPZVec< T > &solR, TPZVec< REAL > &normal, TPZBndCond &bc, int &entropyFix)
Computes the ghost state variables bsed on the BC type.
Definition: pzeulerconslaw.cpp:937
TPZConservationLaw::Dimension
int Dimension() const override
Returns the dimension of the problem.
Definition: pzconslaw.h:276
TPZEulerConsLaw::SetDelta
void SetDelta(REAL delta)
Sets the delta parameter inside the artifficial diffusion term.
Definition: pzeulerconslaw.cpp:272
TPZArtDiff::DiffusionName
TPZString DiffusionName()
Returns the name of diffusive term.
Definition: pzartdiff.cpp:48
TPZMaterial::VariableIndex
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Definition: TPZMaterial.cpp:141
pzerror.h
Defines PZError.
Residual_RT
Definition: pzconslaw.h:61
TPZEulerConsLaw::fConvVol
TPZTimeDiscr fConvVol
Definition: pzeulerconslaw.h:544
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
TPZArtDiff::SetDelta
void SetDelta(REAL delta)
Sets the value for delta.
Definition: pzartdiff.cpp:110
pzeulerconslaw.h
Contains the TPZEulerConsLaw class which implements the weak statement of the compressible euler equa...
TPZEulerConsLaw::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Contributes to the residual vector the boundary conditions.
Definition: pzeulerconslaw.cpp:662
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZArtDiff::Write
void Write(TPZStream &buf, int withclassid) const override
Save the element data to a stream.
Definition: pzartdiff.cpp:829
TPZMaterialData::jacinv
TPZFNMatrix< 9, REAL > jacinv
value of the inverse of the jacobian at the integration point
Definition: pzmaterialdata.h:67
Fad
Definition: fad.h:54
TPZEulerConsLaw::Flux
void Flux(TPZVec< T > &U, TPZVec< T > &Fx, TPZVec< T > &Fy, TPZVec< T > &Fz)
tensor of the three-dimensional flux of Euler
Definition: pzeulerconslaw.h:556
TPZString::Str
const char * Str() const
Explicitly convertes a TPZString into a const null ended char string.
Definition: pzstring.cpp:70
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
tmp
AutoPointerMutexArrayInit tmp
Definition: tpzautopointer.cpp:42
TPZTimeDiscr
TPZTimeDiscr
Indicates the type of time discretization.
Definition: pzconslaw.h:34
TPZEulerConsLaw::uRes
static void uRes(TPZVec< T > &sol, T &us)
Definition: pzeulerconslaw.h:1340
TPZEulerConsLaw::ContributeApproxImplDiff
void ContributeApproxImplDiff(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1263
TPZEulerConsLaw::ContributeApproxImplConvFace
void ContributeApproxImplConvFace(TPZVec< REAL > &x, REAL faceSize, TPZVec< STATE > &solL, TPZVec< STATE > &solR, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &phiR, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, int entropyFix=1)
Definition: pzeulerconslaw.cpp:1451
TPZBndCond::Val2
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZMaterialData::phi
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
Definition: pzmaterialdata.h:55
TPZEulerConsLaw::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzeulerconslaw.cpp:1946
TPZEulerConsLaw::NFluxes
virtual int NFluxes() override
Returns the number of fluxes associated to this material.
Definition: pzeulerconslaw.cpp:211
TPZEulerConsLaw::ApproxRoe_Flux
static void ApproxRoe_Flux(TPZVec< T > &solL, TPZVec< T > &solR, TPZVec< REAL > &normal, REAL gamma, TPZVec< T > &flux, int entropyFix=1)
This flux encapsulates the two and three dimensional fluxes acquired from the Mouse program...
Definition: pzeulerconslaw.cpp:1984
TPZEulerConsLaw::Print
virtual void Print(std::ostream &out) override
Prints the state of internal variables.
Definition: pzeulerconslaw.cpp:102
TPZEulerConsLaw::Read
void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzeulerconslaw.cpp:1958
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZEulerConsLaw::Det
REAL Det(TPZFMatrix< REAL > &Mat)
Computes the determinant of a 2d or 3d matrix. Used by recompute the element size.
Definition: pzeulerconslaw.cpp:184
TPZEulerConsLaw::ContributeImplT1
void ContributeImplT1(TPZVec< REAL > &x, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1873
TPZMaterial::Print
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
Definition: TPZMaterial.cpp:103
pzelmat.h
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
TPZConservationLaw::TimeStep
REAL TimeStep()
Returns the value of the time step.
Definition: pzconslaw.h:307
TPZConservationLaw::fDim
int fDim
Dimension of the problem.
Definition: pzconslaw.h:243
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
pzstring.h
String implementation.
Advanced_CT
Definition: pzconslaw.h:51
TPZEulerConsLaw::ContributeExplT2
void ContributeExplT2(TPZVec< REAL > &x, TPZVec< STATE > &sol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1912
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
TPZEulerConsLaw::ContributeFastestBCInterface
void ContributeFastestBCInterface(int dim, TPZVec< REAL > &x, TPZVec< STATE > &solL, TPZFMatrix< STATE > &dsolL, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &dphiL, TPZFMatrix< REAL > &axesleft, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZEulerConsLaw::DeltaX
REAL DeltaX(REAL detJac)
Estimates the deltax (element diameter) based on the inverse of the jacobian.
Definition: pzeulerconslaw.cpp:179
TPZEulerConsLaw::fArtDiff
TPZArtDiff fArtDiff
diffusive term
Definition: pzeulerconslaw.h:541
TPZEulerConsLaw::SetTimeDiscr
void SetTimeDiscr(TPZTimeDiscr Diff, TPZTimeDiscr ConvVol, TPZTimeDiscr ConvFace)
Configures the time discretization of some contributions.
Definition: pzeulerconslaw.cpp:61
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TFad
Definition: tfad.h:64
TPZEulerConsLaw::ContributeLast
virtual void ContributeLast(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:474
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZArtDiff::Read
void Read(TPZStream &buf, void *context) override
Read the element data from a stream.
Definition: pzartdiff.cpp:838
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
test.res
string res
Definition: test.py:151
TPZEulerConsLaw::ContributeFastestBCInterface_dim
void ContributeFastestBCInterface_dim(TPZVec< REAL > &x, TPZVec< STATE > &solL, TPZFMatrix< STATE > &dsolL, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &dphiL, TPZFMatrix< REAL > &axesleft, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZConservationLaw::fCFL
REAL fCFL
CFL number.
Definition: pzconslaw.h:249
fastAccessDx
expr_ expr_ fastAccessDx(i) *cos(expr_.val())) FAD_FUNC_MACRO(TFadFuncTan
TPZMaterialData::HSize
REAL HSize
measure of the size of the element
Definition: pzmaterialdata.h:83
Last_CT
Definition: pzconslaw.h:50
TPZEulerConsLaw::ContributeInterface
virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Contributes to the residual vector and tangent matrix the face-based quantities.
Definition: pzeulerconslaw.cpp:577
TPZEulerConsLaw::ContributeExplT1
void ContributeExplT1(TPZVec< REAL > &x, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1840
TPZEulerConsLaw::NStateVariables
virtual int NStateVariables() const override
Object-based overload.
Definition: pzeulerconslaw.cpp:91
TPZMaterial::gBigNumber
static REAL gBigNumber
Big number to penalization method, used for Dirichlet conditions.
Definition: TPZMaterial.h:83
TPZEulerConsLaw::OptimalCFL
REAL OptimalCFL(int degree)
returns the best value for the CFL number based on the interpolation degree.
Definition: pzeulerconslaw.cpp:68
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
gamma
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
Definition: tpzgaussrule.cpp:567
TPZEulerConsLaw::ContributeExplDiff
void ContributeExplDiff(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1284
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZConservationLaw
Implements the interface for conservation laws, keeping track of the timestep as well.
Definition: pzconslaw.h:71
pzartdiff.h
Contains the TPZArtDiff class which implements a numerical diffusivity coefficient.
TPZEulerConsLaw::cSpeed
static void cSpeed(TPZVec< T > &sol, REAL gamma, T &c)
Evaluates the speed of sound in the fluid.
Definition: pzeulerconslaw.h:1318
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZConservationLaw::Write
void Write(TPZStream &buf, int withclassid) const override
Save the element data to a stream.
Definition: pzconslaw.cpp:74
FL
#define FL(A)
Function for dynamic cast of the material based on map A (second data)
Definition: pzflowcmesh.cpp:147
TPZConservationLaw::CFL
REAL CFL()
Returns the CFL number.
Definition: pzconslaw.h:281
TPZEulerConsLaw::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
Definition: pzeulerconslaw.cpp:218
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZEulerConsLaw::ContributeExplConvVol
void ContributeExplConvVol(TPZVec< REAL > &x, TPZVec< STATE > &sol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1763
TPZEulerConsLaw::fDiff
TPZTimeDiscr fDiff
variables indication whether the following terms are implicit
Definition: pzeulerconslaw.h:544
TPZEulerConsLaw::Roe_Flux
static void Roe_Flux(TPZVec< T > &solL, TPZVec< T > &solR, TPZVec< REAL > &normal, REAL gamma, TPZVec< T > &flux, int entropyFix=1)
This flux encapsulates the two and three dimensional fluxes acquired from the Mouse program...
Definition: pzeulerconslaw.h:872
TPZVec::Fill
void Fill(const T &copy, const int64_t from=0, const int64_t numelem=-1)
Will fill the elements of the vector with a copy object.
Definition: pzvec.h:460
TPZConservationLaw::fGamma
REAL fGamma
Ratio between specific heat is constant and the specific heat the constant volume of a polytropic gas...
Definition: pzconslaw.h:255
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
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.
Implicit_TD
Implicit time discretization.
Definition: pzconslaw.h:39
TPZConservationLaw::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzconslaw.h:230
TPZEulerConsLaw::ContributeImplConvVol
void ContributeImplConvVol(TPZVec< REAL > &x, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1794
TPZEulerConsLaw::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: pzeulerconslaw.cpp:169
TPZConservationLaw::fContributionTime
TPZContributeTime fContributionTime
Variable indicating the context of the solution.
Definition: pzconslaw.h:261
TPZEulerConsLaw::Pressure
static void Pressure(REAL gamma, int dim, T &press, TPZVec< T > &U)
Thermodynamic pressure determined by the law of an ideal gas.
Definition: pzeulerconslaw.h:817
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZArtDiffType
TPZArtDiffType
Enumerate to define the possible types of artificial diffusion term to stabilize the numerical scheme...
Definition: pzartdiff.h:43
TPZConservationLaw::Print
virtual void Print(std::ostream &out) override
Prints the state of internal variables.
Definition: pzconslaw.cpp:49
TPZArtDiff::ContributeExplDiff
void ContributeExplDiff(int dim, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, TPZFMatrix< REAL > &dphix, TPZFMatrix< STATE > &ef, REAL weight, REAL timeStep, REAL deltaX)
Contributes the diffusion term to the tangent matrix (ek-approximated) and residual vector (ef) ...
Definition: pzartdiff.cpp:707
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZEulerConsLaw::SetTimeStep
virtual REAL SetTimeStep(REAL maxveloc, REAL deltax, int degree) override
See declaration in base class.
Definition: pzeulerconslaw.cpp:73
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
Explicit_TD
Explicit time discretization. Can to be Euler method, Runge Kutta method, etc.
Definition: pzconslaw.h:37
TPZEulerConsLaw::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Contributes to the residual vector and tangent matrix the volume-based quantities.
Definition: pzeulerconslaw.cpp:392
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13