NeoPZ
pzartdiff.h
Go to the documentation of this file.
1 
6 #ifndef TPZARTDIFF_HH
7 #define TPZARTDIFF_HH
8 
9 
10 #include "pzvec.h"
11 #include "pzfmatrix.h"
12 #include "pzdiffmatrix.h"
13 #include "pzstring.h"
14 #include "TPZCompElDisc.h"
15 #include "TPZSavable.h"
16 
17 #ifdef _AUTODIFF
18 #include "fadType.h"
19 #endif
20 
21 
23 inline REAL val(STATE & number)
24 {
25  return number;
26 }
27 
43 enum TPZArtDiffType
44 {
45  None_AD = -1,
46  LeastSquares_AD = 1,
47  Bornhaus_AD = 2,
48  SUPG_AD = 3,
49  TrnLeastSquares_AD = 4
50 };
51 
62 class TPZArtDiff : public TPZSavable
63 {
64 public:
66  TPZArtDiff();
68  TPZArtDiff(TPZArtDiffType type, REAL gamma, REAL CFL = 0., REAL delta = 0.);
69 
70  virtual ~TPZArtDiff();
71 
77  TPZArtDiffType ArtDiffType();
78 
83  void SetArtDiffType(TPZArtDiffType type);
84 
86  TPZString DiffusionName();
87 
89  REAL OptimalDelta();
90 
106  REAL Delta(REAL deltaX, TPZVec<STATE> & sol);
107 
109  void SetDelta(REAL delta);
110 
112  REAL OptimalCFL(int degree = TPZCompEl::GetgOrder());
113 
115  void Write(TPZStream &buf, int withclassid) const override;
116 
118  void Read(TPZStream &buf, void *context) override;
119 
121 public:
122 int ClassId() const override;
123 
124 
126  template< class T >
127  static void Pressure(REAL gamma, int dim, T& press, TPZVec<T> &U);
128 
141  template <class T>
142  void ODotOperator(TPZVec<REAL> &dphi, TPZVec<TPZDiffMatrix<T> > &M, TPZDiffMatrix<T> &Result);
143 
150  template <class T>
151  void ODotOperator(TPZVec<REAL> &dphi, TPZVec<TPZVec<T> > &TauDiv, TPZVec<T> &Result);
152 
153 
163  void Divergent(TPZFMatrix<STATE> &dsol,
164  TPZFMatrix<REAL> & dphi,
165  TPZVec<TPZDiffMatrix<STATE> > & Ai,
166  TPZVec<STATE> & Div,
167  TPZDiffMatrix<STATE> * dDiv);
168 
169 #ifdef _AUTODIFF
170 
183  template <class T>
184  void Divergent(TPZFMatrix<STATE> &dsol,
185  TPZFMatrix<REAL> & phi,
186  TPZFMatrix<REAL> & dphi,
187  TPZVec<TPZDiffMatrix<T> > & Ai,
188  TPZVec<STATE> & Div,
189  TPZDiffMatrix<STATE> * dDiv);
190 
191 
192 
199  void Divergent(TPZVec<FADREAL> &dsol,
200  TPZVec<TPZDiffMatrix<FADREAL> > & Ai,
201  TPZVec<FADREAL> & Div);
202 #endif
203 
214  template <class T>
215  void ComputeTau(int dim,
216  TPZFMatrix<REAL> &jacinv,
217  TPZVec<T> & Sol,
218  TPZVec<TPZDiffMatrix<T> > &Ai,
219  TPZVec<TPZDiffMatrix<T> > &Tau);
220 
225  template <class T>
226  inline static void RotMatrix(TPZVec<T> & sol, T & us, TPZDiffMatrix<T> &Rot, TPZDiffMatrix<T> &RotT);
227 
228  template <class T>
229  inline static void MMatrix(TPZVec<T> & sol, T & us, REAL gamma, TPZDiffMatrix<T> &M, TPZDiffMatrix<T> &Mi);
230 
231  template <class T>
232  inline static void RMMatrix(TPZVec<T> & sol, T & us, REAL gamma, TPZDiffMatrix<T> &RTM, TPZDiffMatrix<T> &RMi);
233 
234  template <class T>
235  static void EigenSystemSUPG(TPZVec<T> & sol, T & us, T & c, REAL gamma, TPZDiffMatrix<T> &X, TPZDiffMatrix<T> &Xi, TPZDiffMatrix<T> &Lambda);
236 
237  template <class T>
238  static void EigenSystemBornhaus(TPZVec<T> & sol, T & us, T & c, REAL gamma, TPZVec<REAL> & aaS, TPZDiffMatrix<T> &Y, TPZDiffMatrix<T> &Yi, TPZDiffMatrix<T> &Lambda);
239 
240  template <class T>
241  static void ContributeBornhaus(TPZVec<T> & sol, T & us, T & c, REAL gamma, TPZVec<REAL> & aaS, TPZDiffMatrix<T> &Mat);
242 
245 private:
246 
247  template <class T>
248  void SUPG(int dim, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau);
249 
250  template <class T>
251  void LS(int dim, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau);
252 
253  template <class T>
254  void LST(int dim, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau);
255 
256  template <class T>
257  void Bornhaus(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau);
258 
259 public:
260 
273  template <class T>
274  void PrepareDiff(int dim,
275  TPZFMatrix<REAL> &jacinv, TPZVec<T> &U,
276  TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau);
277 
290  void PrepareFastDiff(int dim,
291  TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol,
292  TPZFMatrix<STATE> &dsol, TPZFMatrix<REAL> & dphi,
293  TPZVec<TPZVec<STATE> > & TauDiv,
294  TPZVec<TPZDiffMatrix<STATE> > * pTaudDiv = NULL);
295 
296 #ifdef _AUTODIFF
297 
305  void PrepareFastDiff(int dim,
306  TPZFMatrix<REAL> &jacinv, TPZVec<FADREAL> &sol,
307  TPZVec<FADREAL> &dsol, TPZVec<TPZVec<FADREAL> > & TauDiv);
308 
309  template <int dim>
310  void PrepareFastestDiff(TPZFMatrix<REAL> &jacinv,
311  TPZVec<STATE> &sol,
312  TPZFMatrix<STATE> &dsol,
313  TPZFMatrix<REAL> &phi,
314  TPZFMatrix<REAL> &dphi,
315  TPZVec<TPZVec<STATE> > & TauDiv,
316  TPZVec<TPZDiffMatrix<STATE> > & dTauDiv);
319 #endif
320 
338  void ContributeApproxImplDiff(int dim,
339  TPZFMatrix<REAL> &jacinv,
340  TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,
341  TPZFMatrix<REAL> &dphix,
342  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef,
343  REAL weight, REAL timeStep,
344  REAL deltaX);
345 
358  void ContributeExplDiff(int dim,
359  TPZFMatrix<REAL> &jacinv,
360  TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,
361  TPZFMatrix<REAL> &dphix,
362  TPZFMatrix<STATE> &ef,
363  REAL weight, REAL timeStep,
364  REAL deltaX);
365 
366 #ifdef _AUTODIFF
367 
380  void ContributeImplDiff(int dim,
381  TPZFMatrix<REAL> &jacinv,
382  TPZVec<FADREAL> &sol, TPZVec<FADREAL> &dsol,
383  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef,
384  REAL weight, REAL timeStep,
385  REAL deltaX);
386 
401  void ContributeFastestImplDiff(int dim, TPZFMatrix<REAL> &jacinv,
402  TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,
403  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
404  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef,
405  REAL weight, REAL timeStep,
406  REAL deltaX);
407 
408  template <int dim>
409  void ContributeFastestImplDiff_dim(
410  TPZFMatrix<REAL> &jacinv,
411  TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,
412  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
413  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef,
414  REAL weight, REAL timeStep,
415  REAL deltaX);
416 
417 #endif
418 
420 private:
421 
423  TPZArtDiffType fArtDiffType;
424 
426  REAL fGamma;
427 
429  REAL fDelta;
430 
432  REAL fCFL;
433 
434 
435 };
436 
437 
438 template <class T>
439 void TPZArtDiff::RotMatrix(TPZVec<T> & sol, T & us, TPZDiffMatrix<T> &Rot, TPZDiffMatrix<T> &RotT)
440 {
441  int nstate = sol.NElements();
442  int dim = nstate - 2;
443 
444  T u, v, w, u2, v2, w2, uspuInv, usInv;
445 
446  switch (dim)
447  {
448  case (2):
449  u = sol[1]/sol[0];
450  v = sol[2]/sol[0];
451  u2 = u * u;
452  v2 = v * v;
453  usInv= T(1.)/us;
454 
455  Rot.Redim(nstate,nstate);
456  Rot(0,0) = 1.;
457  Rot(3,3) = 1.;
458  Rot(1,1) = u * usInv;
459  Rot(1,2) = v * usInv;
460  Rot(2,1) = - Rot(1,2);
461  Rot(2,2) = Rot(1,1);
462 
463  Rot.Transpose(RotT);
464  break;
465  case (3):
466  u = sol[1]/sol[0];
467  v = sol[2]/sol[0];
468  w = sol[3]/sol[0];
469  u2 = u * u;
470  v2 = v * v;
471  w2 = w * w;
472  uspuInv = T(1.)/(us + u);
473  usInv = T(1.)/us;
474  uspuInv *= usInv;
475 
476  Rot.Redim(nstate,nstate);
477  Rot(0,0) = 1.;
478  Rot(4,4) = 1.;
479  Rot(1,1) = u * usInv;
480  Rot(1,2) = v * usInv;
481  Rot(1,3) = w * usInv;
482  Rot(2,1) = -Rot(1,2);
483  Rot(2,2) = Rot(1,1) + w2 * uspuInv;
484  Rot(2,3) = - w*v * uspuInv;
485  Rot(3,1) = -Rot(1,3);
486  Rot(3,2) = Rot(2,3);
487  Rot(3,3) = Rot(1,1) + v2 * uspuInv;
488 
489  Rot.Transpose(RotT);
490  break;
491 
492  default:
493  PZError << "TPZArtDiff::RotMatrix Error: Invalid Dimension\n";
494  }
495 
496 }
497 
498 
499 template <class T>
500 void TPZArtDiff::MMatrix(TPZVec<T> & sol, T & us, REAL gamma, TPZDiffMatrix<T> &M, TPZDiffMatrix<T> &Mi)
501 {
502  int nstate = sol.NElements();
503  int dim = nstate - 2;
504 
505  T rhoInv, usInv;
506 
507  switch (dim)
508  {
509  case (2):
510  usInv= T(1.)/us;
511  rhoInv = T(1.)/sol[0];
512 
513  M.Redim(nstate,nstate);
514  M(0,0) = 1.;
515  M(1,0) = us;
516  M(1,1) = sol[0];
517  M(2,2) = sol[0];
518  M(3,0) = us * us /2.;
519  M(3,1) = sol[0] * us;
520  M(3,3) = T(1.)/(gamma-1.);
521 
522  Mi.Redim(nstate,nstate);
523  Mi(0,0) = 1.;
524  Mi(1,0) = - rhoInv * us;
525  Mi(1,1) = rhoInv;
526  Mi(2,2) = rhoInv;
527  Mi(3,0) = T((gamma-1.)/2.) * us * us;
528  Mi(3,1) = T(1.-gamma)*us;
529  Mi(3,3) = gamma-1.;
530 
531  break;
532  case(3):
533  usInv= T(1.)/us;
534  rhoInv = T(1.)/sol[0];
535 
536  M.Redim(nstate,nstate);
537  M(0,0) = 1.;
538  M(1,0) = us;
539  M(1,1) = sol[0];
540  M(2,2) = sol[0];
541  M(3,3) = sol[0];
542  M(4,0) = us * us /2.;
543  M(4,1) = sol[0] * us;
544  M(4,4) = T(1.)/(gamma-1.);
545 
546  Mi.Redim(nstate,nstate);
547  Mi(0,0) = 1.;
548  Mi(1,0) = - rhoInv * us;
549  Mi(1,1) = rhoInv;
550  Mi(2,2) = rhoInv;
551  Mi(3,3) = rhoInv;
552  Mi(4,0) = T((gamma-1.)/2.) * us * us;
553  Mi(4,1) = T(1.-gamma)*us;
554  Mi(4,4) = gamma-1.;
555  break;
556 
557  default:
558  PZError << "TPZArtDiff::MMatrix Error: Invalid Dimension\n";
559  }
560 }
561 
562 
563 template <class T>
564 void TPZArtDiff::RMMatrix(TPZVec<T> & sol, T & us, REAL gamma, TPZDiffMatrix<T> &RTM, TPZDiffMatrix<T> &RMi)
565 {
566  int nstate = sol.NElements();
567  int dim = nstate - 2;
568 
569  T u, v, w, u2, v2, w2,
570  rhoInv, usInv, uspu, uspuInv;
571  T kt, st;
572 
573  switch (dim)
574  {
575  case (2):
576  u = sol[1]/sol[0];
577  usInv= T(1.)/us;
578  rhoInv = T(1.)/sol[0];
579  // kt = u/us, st = v/us
580  kt = sol[1]/sol[0]/us;
581  st = sol[2]/sol[0]/us;
582 
583  RTM.Redim(nstate,nstate);
584  RTM(0,0) = 1.;
585  RTM(1,0) = us * kt;
586  RTM(1,1) = sol[0] * kt;
587  RTM(1,2) = -sol[0] * st;
588  RTM(2,0) = us * st;
589  RTM(2,1) = -RTM(1,2);
590  RTM(2,2) = RTM(1,1);
591  RTM(3,0) = us * us /2.;
592  RTM(3,1) = sol[0] * us;
593  RTM(3,3) = 1./(gamma-1.);
594 
595  RMi.Redim(nstate,nstate);
596  RMi(0,0) = 1.;
597  RMi(1,0) = - rhoInv * us;
598  RMi(1,1) = rhoInv * kt;
599  RMi(1,2) = rhoInv * st;
600  RMi(2,1) = - RMi(1,2);
601  RMi(2,2) = RMi(1,1);
602  RMi(3,0) = T((gamma-1.)/2.) * us * us;
603  RMi(3,1) = T(1.-gamma)*us*kt;
604  RMi(3,2) = T(1.-gamma)*us*st;
605  RMi(3,3) = gamma-1.;
606 
607  break;
608  case(3):
609  u = sol[1]/sol[0];
610  v = sol[2]/sol[0];
611  w = sol[3]/sol[0];
612  u2 = u * u;
613  v2 = v * v;
614  w2 = w * w;
615  uspu = us + u;
616  uspuInv = T(1.)/uspu;
617  usInv = T(1.)/us;
618  kt = T(1.)/(sol[0] * us * (uspu * u + v2 + w2));
619 #ifdef NOTDEFINED
620  RTM.Redim(nstate,nstate);
621  RTM(0,0) = 1.;
622  RTM(1,0) = u2 * usInv;
623  RTM(1,1) = sol[1]/*rhou*/*usInv;
624  RTM(1,2) = sol[2]/*rhov*/*usInv;
625  RTM(1,3) = sol[3]/*rhow*/*usInv;
626  RTM(2,0) = - u * v * usInv;
627  RTM(2,1) = - RTM(1,2);
628  RTM(2,2) = RTM(1,1) + sol[3] * w * uspuInv * usInv;
629  RTM(2,3) = - RTM(1,2) * w * uspuInv;
630  RTM(3,0) = - u * w * usInv;
631  RTM(3,1) = -RTM(1,3);
632  RTM(3,2) = RTM(2,3);
633  RTM(3,3) = RTM(1,1) + sol[2] * v * uspuInv * usInv;
634  RTM(4,0) = u2/2.;
635  RTM(4,1) = sol[1];
636  RTM(4,4) = 1./(gamma-1.);
637 
638  RMi.Redim(nstate,nstate);
639  RMi(0,0) = 1.;
640  RMi(1,0) = - u / sol[0];
641  RMi(1,1) = - RMi(1,0) * usInv;
642  RMi(1,2) = - v * usInv / sol[0];
643  RMi(1,3) = - w * usInv / sol[0];
644  RMi(2,1) = - RMi(1,2);
645  RMi(2,2) = (u * v2 + uspu * (u2+w2)) * kt;
646  RMi(2,3) = (u - uspu) * v * w * kt;
647  RMi(3,1) = - RMi(1,3);
648  RMi(3,2) = RMi(2,3);
649  RMi(3,3) = (uspu * (u2+v2) + u* w2) * kt;
650  RMi(4,0) = u2 * T((gamma-1.)/2.);
651  RMi(4,1) = - u2 * usInv * T(gamma - 1.);
652  RMi(4,2) = u * v * usInv * T(gamma -1.);
653  RMi(4,3) = u * w * usInv * T(gamma -1.);
654  RMi(4,4) = gamma -1.;
655 #endif
656 
657  RTM.Redim(nstate,nstate);
658  RTM(0,0) = 1.;
659  RTM(1,0) = u;
660  RTM(1,1) = sol[1]/*rhou*/*usInv;
661  RTM(1,2) = - sol[2]/*rhov*/*usInv;
662  RTM(1,3) = - sol[3]/*rhow*/*usInv;
663  RTM(2,0) = v;
664  RTM(2,1) = - RTM(1,2);
665  RTM(2,2) = RTM(1,1) + sol[3] * w * uspuInv * usInv;
666  RTM(2,3) = RTM(1,2) * w * uspuInv;
667  RTM(3,0) = w;
668  RTM(3,1) = -RTM(1,3);
669  RTM(3,2) = RTM(2,3);
670  RTM(3,3) = RTM(1,1) + sol[2] * v * uspuInv * usInv;
671  RTM(4,0) = us * us /2.;
672  RTM(4,1) = sol[0] * us;
673  RTM(4,4) = 1./(gamma-1.);
674 
675  RMi.Redim(nstate,nstate);
676  RMi(0,0) = 1.;
677  RMi(1,0) = - us / sol[0];
678  RMi(1,1) = u * usInv / sol[0];
679  RMi(1,2) = v * usInv / sol[0];
680  RMi(1,3) = w * usInv / sol[0];
681  RMi(2,1) = - RMi(1,2);
682  RMi(2,2) = (u * v2 + uspu * (u2+w2)) * kt;
683  RMi(2,3) = (u - uspu) * v * w * kt;
684  RMi(3,1) = - RMi(1,3);
685  RMi(3,2) = RMi(2,3);
686  RMi(3,3) = (uspu * (u2+v2) + u* w2) * kt;
687  RMi(4,0) = us * us * T((gamma-1.)/2.);
688  RMi(4,1) = u * T(1.-gamma);
689  RMi(4,2) = v * T(1.-gamma);
690  RMi(4,3) = w * T(1.-gamma);
691  RMi(4,4) = gamma -1.;
692 
693  break;
694 
695  default:
696  PZError << "TPZArtDiff::MMatrix Error: Invalid Dimension\n";
697  }
698 }
699 
700 
701 
702 template <class T>
703 void TPZArtDiff::EigenSystemSUPG(TPZVec<T> & sol, T & us, T & c, REAL gamma, TPZDiffMatrix<T> &X, TPZDiffMatrix<T> &Xi, TPZDiffMatrix<T> &Lambda)
704 {
705  int nstate = sol.NElements();
706  int dim = nstate - 2;
707 
708  T us2, c2, s, rho_c_us;
709 
710  c2 = c * c;
711  us2 = us * us;
712  rho_c_us = sol[0] * c * us;
713 
714  switch (dim)
715  {
716  case (2):
717  s = sqrt(c2 + us2 * REAL(16.));
718 
719  X.Redim(nstate,nstate);
720  X(0,0) = 1.;
721  X(0,2) = REAL(4.) * sol[0] * us / c2;
722  X(0,3) = X(0,2);
723  X(1,2) = - s/c -1.;
724  X(1,3) = s/c -1.;
725  X(2,1) = 1.;
726  X(3,2) = REAL(4.) * sol[0] * us;
727  X(3,3) = X(3,2);
728 
729  Xi.Redim(nstate,nstate);
730  Xi(0,0) = 1.;
731  Xi(0,3) = REAL(- 1.) / c2;
732  Xi(1,2) = 1.;
733  Xi(2,1) = - c / (REAL(2.) * s);
734  Xi(2,3) = (s-c)/(REAL(8.) * sol[0] * s * us);
735  Xi(3,1) = -Xi(2,1);
736  Xi(3,3) = (s+c)/(REAL(8.) * sol[0] * s * us);
737 
738  Lambda.Redim(nstate,nstate);
739  Lambda(0,0) = REAL(1.)/us;
740  Lambda(1,1) = REAL(1.)/sqrt(us2 + c2);
741  Lambda(2,2) = REAL(1.)/sqrt(us2 + REAL(1.5) * c2 - REAL(.5) * c * s);
742  Lambda(3,3) = 1./sqrt(us2 + REAL(1.5) * c2 + REAL(.5) * c * s);
743 
744  break;
745  case(3):
746 
747  s = sqrt(c2 + us2 * REAL(4.));
748 
749  X.Redim(nstate,nstate);
750  X(0,0) = 1.;
751  X(0,3) = REAL(1.) / c2;
752  X(0,4) = X(0,3);
753  X(1,3) = (-s-c)/(REAL(2.) * rho_c_us);
754  X(1,4) = ( s-c)/(REAL(2.) * rho_c_us);
755  X(2,2) = 1.;
756  X(3,1) = 1.;
757  X(4,3) = 1.;
758  X(4,4) = 1.;
759 
760 
761  Xi.Redim(nstate,nstate);
762  Xi(0,0) = 1.;
763  Xi(0,4) = REAL(- 1.) / c2;
764  Xi(1,3) = 1.;
765  Xi(2,2) = 1.;
766  Xi(3,1) = - rho_c_us / s;
767  Xi(3,4) = REAL(.5) - c / (REAL(2.) * s);
768  Xi(4,1) = -Xi(3,1);
769  Xi(4,4) = REAL(.5) + c / (REAL(2.) * s);
770 
771  Lambda.Redim(nstate,nstate);
772  Lambda(0,0) = REAL(1.)/us;
773  Lambda(1,1) = REAL(1.)/sqrt(us2 + c2);
774  Lambda(2,2) = Lambda(1,1);
775  Lambda(3,3) = REAL(1.)/sqrt(us2 + REAL(2.) * c2 - c * s);
776  Lambda(4,4) = REAL(1.)/sqrt(us2 + REAL(2.) * c2 + c * s);
777  break;
778 
779  default:
780  PZError << "TPZArtDiff::EigenSystemSUPG Error: Invalid Dimension\n";
781  }
782 }
783 
784 
785 
786 
787 template <class T>
788 void TPZArtDiff::EigenSystemBornhaus(TPZVec<T> & sol, T & us, T & c, REAL gamma, TPZVec<REAL> & aaS, TPZDiffMatrix<T> &Y, TPZDiffMatrix<T> &Yi, TPZDiffMatrix<T> &Lambda)
789 {
790  int nstate = sol.NElements();
791  int dim = nstate - 2;
792 
793  T k, us2, c2, rho_c, temp1, temp2, temp3, k2;
794 
795  c2 = c * c;
796  us2 = us * us;
797  rho_c = sol[0] * c;
798 
799  switch (dim)
800  {
801  case (2):
802  k2 = aaS[0]*aaS[0] + aaS[1]*aaS[1];
803  k = sqrt(k2);
804 
805 
806  Y.Redim(nstate,nstate);
807  Y(0,1) = 1.;
808  Y(0,2) = REAL(1.) / c2;
809  Y(0,3) = Y(0,2);
810  Y(1,0) = -aaS[1] ;
811  Y(1,2) = -aaS[0] / (k * rho_c);
812  Y(1,3) = -Y(1,2);
813  Y(2,0) = aaS[0];//1.;
814  Y(2,2) = -aaS[1] / (k * rho_c);
815  Y(2,3) = -Y(2,2);
816  Y(3,2) = 1.;
817  Y(3,3) = 1.;
818 
819  Yi.Redim(nstate,nstate);
820  Yi(0,1) = - /*aaS[0] * */aaS[1] / k2;
821  Yi(0,2) = /*aaS[0] * */aaS[0] / k2;
822  Yi(1,0) = 1.;
823  Yi(1,3) = REAL(-1.) / c2;
824  Yi(2,1) = - aaS[0] * rho_c / (REAL(2.) * k);
825  Yi(2,2) = - aaS[1] * rho_c / (REAL(2.) * k);
826  Yi(2,3) = .5;
827  Yi(3,1) = -Yi(2,1);
828  Yi(3,2) = -Yi(2,2);
829  Yi(3,3) = .5;
830 
831  temp1 = aaS[0] * us;
832  if(val(temp1) < 0)temp1 = -temp1;
833  temp2 = aaS[0] * us - k * c;
834  if(val(temp2) < 0)temp2 = -temp2;
835  temp3 = aaS[0] * us + k * c;
836  if(val(temp3) < 0)temp3 = -temp3;
837 
838  Lambda.Redim(nstate,nstate);
839  Lambda(0,0) = temp1;
840  Lambda(1,1) = temp1;
841  Lambda(2,2) = temp2;
842  Lambda(3,3) = temp3;
843 
844  break;
845  case (3):
846  k2 = aaS[0]*aaS[0] + aaS[1]*aaS[1] + aaS[2]*aaS[2];
847  k = sqrt(k2);
848 
849  Y.Redim(nstate,nstate);
850  Y(0,2) = 1.;
851  Y(0,3) = REAL(1.) / c2;
852  Y(0,4) = Y(0,3);
853  Y(1,0) = -aaS[2]/aaS[0];
854  Y(1,1) = -aaS[1]/aaS[0];
855  Y(1,3) = -aaS[0]/(k * rho_c);
856  Y(1,4) = -Y(1,3);
857  Y(2,1) = 1.;
858  Y(2,3) = -aaS[1]/(k * rho_c);
859  Y(2,4) = -Y(2,3);
860  Y(3,0) = 1.;
861  Y(3,3) = -aaS[2]/(k * rho_c);
862  Y(3,4) = -Y(3,3);
863  Y(4,3) = 1.;
864  Y(4,4) = 1.;
865 
866  Yi.Redim(nstate,nstate);
867  Yi(0,1) = - aaS[0] * aaS[2] / k2;
868  Yi(0,2) = - aaS[1] * aaS[2] / k2;
869  Yi(0,3) = (aaS[0] * aaS[0] + aaS[1] * aaS[1]) / k2;
870  Yi(1,1) = - aaS[0] * aaS[1] / k2;
871  Yi(1,2) = (aaS[0] * aaS[0] + aaS[2] * aaS[2]) / k2;
872  Yi(1,3) = Yi(0,2);
873  Yi(2,0) = 1.;
874  Yi(2,4) = REAL(-1.) / c2;
875  Yi(3,1) = - aaS[0] * rho_c / (REAL(2.) * k);
876  Yi(3,2) = - aaS[1] * rho_c / (REAL(2.) * k);
877  Yi(3,3) = - aaS[2] * rho_c / (REAL(2.) * k);
878  Yi(3,4) = .5;
879  Yi(4,1) = - Yi(3,1);
880  Yi(4,2) = - Yi(3,2);
881  Yi(4,3) = - Yi(3,3);
882  Yi(4,4) = .5;
883 
884  temp1 = aaS[0] * us;
885  if(val(temp1) < 0)temp1 = -temp1;
886  temp2 = aaS[0] * us - k * c;
887  if(val(temp2) < 0)temp2 = -temp2;
888  temp3 = aaS[0] * us + k * c;
889  if(val(temp3) < 0)temp3 = -temp3;
890 
891  Lambda.Redim(nstate,nstate);
892  Lambda(0,0) = temp1;
893  Lambda(1,1) = temp1;
894  Lambda(2,2) = temp1;
895  Lambda(3,3) = temp2;
896  Lambda(4,4) = temp3;
897 
898  break;
899  default:
900  PZError << "TPZArtDiff::EigenSystemBornhaus Error: Invalid Dimension\n";
901  }
902 }
903 
904 
905 
906 template <class T>
907 void TPZArtDiff::ContributeBornhaus(TPZVec<T> & sol, T & us, T & c, REAL gamma, TPZVec<REAL> & aaS, TPZDiffMatrix<T> &Mat)
908 {
909  int nstate = sol.NElements();
910  int dim = nstate - 2;
911 
912  T k, l1, l3, l4, l5, lstar,
913  lsum, ldiff, twoCK, c2, temp1, temp2/*, temp3*/, k2;
914  //T us2;
915 
916  c2 = c * c;
917  //us2 = us * us;
918  //rho_c = sol[0] * c;
919 
920  switch (dim)
921  {
922  case (2):
923  k2 = aaS[0]*aaS[0] + aaS[1]*aaS[1];
924  k = sqrt(k2);
925  twoCK = c * T(k * REAL(2.));
926 
927  l1 = aaS[0] * us;
928  if(val(l1) < 0.)l1 = -l1;
929  l3 = aaS[0] * us - k * c;
930  if(val(l3) < 0.)l3 = -l3;
931  l4 = aaS[0] * us + k * c;
932  if(val(l4) < 0.)l4 = -l4;
933 
934  temp1 = (l4 - l3) * sol[0] / twoCK;
935  temp2 = l3 + l4 - REAL(2.)*l1;
936  Mat(0,0) += l1;
937  Mat(0,1) += aaS[0] * temp1;
938  Mat(0,2) += aaS[1] * temp1;
939  Mat(0,3) += temp2 / (REAL(2.) * c2);
940  Mat(1,1) += (aaS[1]*aaS[1] * l1 + aaS[0]*aaS[0] * (l3+l4) *REAL(.5))/k2;
941  Mat(1,2) += aaS[0]*aaS[1] * temp2/(REAL(2.) * k2);
942  Mat(1,3) += aaS[0]*(l4-l3)/(c * k * sol[0] * REAL(2.));
943  Mat(2,1) += aaS[0]*aaS[1] * temp2/(REAL(2.) * k2);
944  Mat(2,2) += (aaS[0]*aaS[0] * l1 + aaS[1]*aaS[1] * (l3+l4) *REAL(.5))/k2;
945  Mat(2,3) += aaS[1]*(l4-l3)/(c * k * sol[0] * REAL(2.));
946  Mat(3,1) += aaS[0] * c2 * temp1;
947  Mat(3,2) += aaS[1] * c2 * temp1;
948  Mat(3,3) += (l3 + l4)/REAL(2.);
949 
950 
951  break;
952  case (3):
953  k2 = aaS[0]*aaS[0] + aaS[1]*aaS[1] + aaS[2]*aaS[2];
954  k = sqrt(k2);
955 
956  l1 = aaS[0] * us;
957  if(val(l1) < 0.)l1 = -l1;
958  l4 = aaS[0] * us - k * c;
959  if(val(l3) < 0.)l3 = -l3;
960  l5 = aaS[0] * us + k * c;
961  if(val(l4) < 0.)l4 = -l4;
962 
963  lsum = l5 + l4;
964  ldiff= l5 - l4;
965  lstar = lsum - REAL(2.)*l1;
966 
967  twoCK = c * k * REAL(2.);
968 
969  Mat.Redim(nstate,nstate);
970  Mat(0,0) += l1;
971  Mat(0,1) += aaS[0] * sol[0] * ldiff / twoCK;
972  Mat(0,2) += aaS[1] * sol[0] * ldiff / twoCK;
973  Mat(0,3) += aaS[2] * sol[0] * ldiff / twoCK;
974  Mat(0,4) += lstar / ( c2 * REAL(2.) );
975  Mat(1,1) += aaS[0] * aaS[0] * lstar / (k2 * REAL(2.)) + l1;
976  Mat(1,2) += aaS[0] * aaS[1] * lstar / (k2 * REAL(2.));
977  Mat(1,3) += aaS[0] * aaS[2] * lstar / (k2 * REAL(2.));
978  Mat(1,4) += aaS[0] / twoCK / sol[0] * ldiff;
979  Mat(2,1) += aaS[0] * aaS[1] * lstar / (k2 * REAL(2.));
980  Mat(2,2) += aaS[1] * aaS[1] * lstar / (k2 * REAL(2.)) + l1;
981  Mat(2,3) += aaS[1] * aaS[2] * lstar / (k2 * REAL(2.));
982  Mat(2,4) += aaS[1] / twoCK / sol[0] * ldiff;
983  Mat(3,1) += aaS[0] * aaS[2] * lstar / (k2 * REAL(2.));
984  Mat(3,2) += aaS[1] * aaS[2] * lstar / (k2 * REAL(2.));
985  Mat(3,3) += aaS[2] * aaS[2] * lstar / (k2 * REAL(2.)) + l1;
986  Mat(3,4) += aaS[2] / twoCK / sol[0] * ldiff;
987  Mat(4,1) += aaS[0] * c * sol[0] * ldiff / (k * REAL(2.));
988  Mat(4,2) += aaS[1] * c * sol[0] * ldiff / (k * REAL(2.));
989  Mat(4,3) += aaS[2] * c * sol[0] * ldiff / (k * REAL(2.));
990  Mat(4,4) += lstar/REAL(2.) + l1;
991 
992  break;
993  default:
994  PZError << "TPZArtDiff::EigenSystemBornhaus Error: Invalid Dimension\n";
995  }
996 }
997 
998 #endif
999 
TPZSavable.h
Contains declaration of the TPZSavable class which defines the interface to save and restore objects ...
TPZArtDiff::LST
void LST(int dim, TPZVec< T > &sol, TPZVec< TPZDiffMatrix< T > > &Ai, TPZVec< TPZDiffMatrix< T > > &Tau)
Definition: pzartdiff.cpp:391
TPZArtDiff::RMMatrix
static void RMMatrix(TPZVec< T > &sol, T &us, REAL gamma, TPZDiffMatrix< T > &RTM, TPZDiffMatrix< T > &RMi)
Definition: pzartdiff.h:564
pzdiffmatrix.h
Contains TPZDiffMatrix class which to hold the flux derivatives A B C and diffusive matrix coefficien...
None_AD
No Artificial diffusion term is considered.
Definition: pzartdiff.h:45
LeastSquares_AD
Use Least squares method to applied artificial diffusion term.
Definition: pzartdiff.h:46
pzvec.h
Templated vector implementation.
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
TPZArtDiff::DiffusionName
TPZString DiffusionName()
Returns the name of diffusive term.
Definition: pzartdiff.cpp:48
TPZDiffMatrix::Transpose
TPZDiffMatrix< T > & Transpose(TPZDiffMatrix< T > &matrix)
Transposes the matrix onto the parameter object.
Definition: pzdiffmatrix.h:207
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
TPZArtDiff::Write
void Write(TPZStream &buf, int withclassid) const override
Save the element data to a stream.
Definition: pzartdiff.cpp:829
TPZArtDiff::ODotOperator
void ODotOperator(TPZVec< REAL > &dphi, TPZVec< TPZDiffMatrix< T > > &M, TPZDiffMatrix< T > &Result)
Operation product point in the diffusion term.
Definition: pzartdiff.cpp:122
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZArtDiff::fDelta
REAL fDelta
Scalar coefficient of the element in the diffusion term.
Definition: pzartdiff.h:429
TPZArtDiff::SUPG
void SUPG(int dim, TPZVec< T > &sol, TPZVec< TPZDiffMatrix< T > > &Ai, TPZVec< TPZDiffMatrix< T > > &Tau)
Definition: pzartdiff.cpp:323
TPZArtDiff
This class adds to the term of diffusion to the variacional formulation of the differential equation...
Definition: pzartdiff.h:62
TPZArtDiff::Delta
REAL Delta(REAL deltaX, TPZVec< STATE > &sol)
Returns the stored value for delta.
Definition: pzartdiff.cpp:82
SUPG_AD
Use Supg method to applied artificial diffusion term.
Definition: pzartdiff.h:48
TPZArtDiff::Divergent
void Divergent(TPZFMatrix< STATE > &dsol, TPZFMatrix< REAL > &dphi, TPZVec< TPZDiffMatrix< STATE > > &Ai, TPZVec< STATE > &Div, TPZDiffMatrix< STATE > *dDiv)
Evaluates the divergent of F.
Definition: pzartdiff.cpp:224
TPZArtDiff::fGamma
REAL fGamma
Ratio between specific heat is constant and the specific heat the constant volume of a polytropic gas...
Definition: pzartdiff.h:426
pzstring.h
String implementation.
TPZDiffMatrix::Redim
void Redim(const int64_t rows, const int64_t cols)
Resizes and zeroes the matrix.
Definition: pzdiffmatrix.h:191
TPZArtDiff::EigenSystemBornhaus
static void EigenSystemBornhaus(TPZVec< T > &sol, T &us, T &c, REAL gamma, TPZVec< REAL > &aaS, TPZDiffMatrix< T > &Y, TPZDiffMatrix< T > &Yi, TPZDiffMatrix< T > &Lambda)
Definition: pzartdiff.h:788
TPZArtDiff::EigenSystemSUPG
static void EigenSystemSUPG(TPZVec< T > &sol, T &us, T &c, REAL gamma, TPZDiffMatrix< T > &X, TPZDiffMatrix< T > &Xi, TPZDiffMatrix< T > &Lambda)
Definition: pzartdiff.h:703
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
Bornhaus_AD
Use Bornhaus method to applied artificial diffusion term.
Definition: pzartdiff.h:47
TPZArtDiff::LS
void LS(int dim, TPZVec< T > &sol, TPZVec< TPZDiffMatrix< T > > &Ai, TPZVec< TPZDiffMatrix< T > > &Tau)
Definition: pzartdiff.cpp:382
TPZArtDiff::Read
void Read(TPZStream &buf, void *context) override
Read the element data from a stream.
Definition: pzartdiff.cpp:838
TPZArtDiff::~TPZArtDiff
virtual ~TPZArtDiff()
Definition: pzartdiff.cpp:34
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
TPZArtDiff::fArtDiffType
TPZArtDiffType fArtDiffType
Kind of artificial diffusion term to apply.
Definition: pzartdiff.h:423
TPZCompEl::GetgOrder
static int GetgOrder()
Set the default value of the interpolation order.
Definition: pzcompel.h:830
TPZCompElDisc.h
Contains declaration of TPZCompelDisc class which implements a computational element for discontinuou...
TrnLeastSquares_AD
Use Transpose of the Least Squares method to applied artificial diffusion term.
Definition: pzartdiff.h:49
TPZString
Implements strings as stack. Utility.
Definition: pzstring.h:18
gamma
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
Definition: tpzgaussrule.cpp:567
TPZArtDiff::PrepareFastDiff
void PrepareFastDiff(int dim, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, TPZFMatrix< REAL > &dphi, TPZVec< TPZVec< STATE > > &TauDiv, TPZVec< TPZDiffMatrix< STATE > > *pTaudDiv=NULL)
Prepares the data to compute the diffusive term as fast as possible, sparing operations.
Definition: pzartdiff.cpp:494
TPZArtDiff::Bornhaus
void Bornhaus(int dim, TPZFMatrix< REAL > &jacinv, TPZVec< T > &sol, TPZVec< TPZDiffMatrix< T > > &Ai, TPZVec< TPZDiffMatrix< T > > &Tau)
Definition: pzartdiff.cpp:400
TPZArtDiff::OptimalCFL
REAL OptimalCFL(int degree=TPZCompEl::GetgOrder())
Returns the best value for based on the interpolation degree.
Definition: pzartdiff.cpp:115
TPZArtDiff::fCFL
REAL fCFL
number
Definition: pzartdiff.h:432
TPZArtDiff::OptimalDelta
REAL OptimalDelta()
Returns the best value for delta based on the interpolation degree.
Definition: pzartdiff.cpp:75
TPZArtDiff::MMatrix
static void MMatrix(TPZVec< T > &sol, T &us, REAL gamma, TPZDiffMatrix< T > &M, TPZDiffMatrix< T > &Mi)
Definition: pzartdiff.h:500
TPZArtDiff::SetArtDiffType
void SetArtDiffType(TPZArtDiffType type)
Configures the type of artificial diffusion.
Definition: pzartdiff.cpp:43
TPZArtDiff::Pressure
static void Pressure(REAL gamma, int dim, T &press, TPZVec< T > &U)
pressure
TPZArtDiff::TPZArtDiff
TPZArtDiff()
Simple constructor.
Definition: pzartdiff.cpp:25
TPZArtDiff::PrepareDiff
void PrepareDiff(int dim, TPZFMatrix< REAL > &jacinv, TPZVec< T > &U, TPZVec< TPZDiffMatrix< T > > &Ai, TPZVec< TPZDiffMatrix< T > > &Tau)
Computes the common values A B C and Tau vector of matrixes for contributions.
Definition: pzartdiff.cpp:487
TPZDiffMatrix
Matrix class to hold the flux derivatives A B C and diffusive matrix coefficients. Matrix.
Definition: pzdiffmatrix.h:27
TPZArtDiff::ClassId
int ClassId() const override
Class identificator.
Definition: pzartdiff.cpp:848
TPZArtDiff::ComputeTau
void ComputeTau(int dim, TPZFMatrix< REAL > &jacinv, TPZVec< T > &Sol, TPZVec< TPZDiffMatrix< T > > &Ai, TPZVec< TPZDiffMatrix< T > > &Tau)
Computes the diffusive term according to the name.
Definition: pzartdiff.cpp:297
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
TPZArtDiff::ContributeBornhaus
static void ContributeBornhaus(TPZVec< T > &sol, T &us, T &c, REAL gamma, TPZVec< REAL > &aaS, TPZDiffMatrix< T > &Mat)
Definition: pzartdiff.h:907
TPZSavable
This class defines the interface to save and restore objects from TPZStream objects. Persistency.
Definition: TPZSavable.h:67
TPZArtDiffType
TPZArtDiffType
Enumerate to define the possible types of artificial diffusion term to stabilize the numerical scheme...
Definition: pzartdiff.h:43
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
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZArtDiff::RotMatrix
static void RotMatrix(TPZVec< T > &sol, T &us, TPZDiffMatrix< T > &Rot, TPZDiffMatrix< T > &RotT)
Definition: pzartdiff.h:439
TPZArtDiff::ArtDiffType
TPZArtDiffType ArtDiffType()
Returns the type of artifical diffusion.
Definition: pzartdiff.cpp:38
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13