NeoPZ
pzartdiff.cpp
Go to the documentation of this file.
1 
6 #include "pzartdiff.h"
7 #include "TPZDiffusionConsLaw.h"
8 #include "pzfmatrix.h"
9 #include "pzvec.h"
10 #include "pzreal.h"
11 #include "pzeulerconslaw.h"
12 
13 #define FASTESTDIFF
14 
15 TPZArtDiff::TPZArtDiff(TPZArtDiffType type, REAL gamma, REAL CFL, REAL delta):
16 TPZRegisterClassId(&TPZArtDiff::ClassId),
17 fArtDiffType(type),
18 fGamma(gamma),
19 fDelta(delta),
20 fCFL(CFL)
21 {
22 
23 }
24 
25 TPZArtDiff::TPZArtDiff():
26 TPZRegisterClassId(&TPZArtDiff::ClassId),
27 fArtDiffType(None_AD),
28 fGamma(1.4),
29 fDelta(0.),
30 fCFL(1.)
31 {
32 }
33 
34 TPZArtDiff::~TPZArtDiff()
35 {
36 }
37 
38 TPZArtDiffType TPZArtDiff::ArtDiffType()
39 {
40  return fArtDiffType;
41 }
42 
43 void TPZArtDiff::SetArtDiffType(TPZArtDiffType type)
44 {
45  fArtDiffType = type;
46 }
47 
48 TPZString TPZArtDiff::DiffusionName()
49 {
50  TPZString rtstr;
51  switch(fArtDiffType)
52  {
53  case SUPG_AD:
54  rtstr = "SUPG";
55  return rtstr;
56  // break;
57  case LeastSquares_AD:
58  rtstr = "LeastSquares";
59  return rtstr;
60  // break;
61  case Bornhaus_AD:
62  rtstr = "Bornhaus";
63  return rtstr;
64  // break;
65  case TrnLeastSquares_AD:
66  rtstr = "TrnLeastSquares";
67  return rtstr;
68  // break;
69  default:
70  PZError << "Unknown artificial diffusion term (" << fArtDiffType << ")";
71  return rtstr;
72  }
73 }
74 
75 REAL TPZArtDiff::OptimalDelta()
76 {
77  REAL cfl = OptimalCFL(/*degree*/);
78  REAL delta = ( (10./3.)*cfl*cfl - (2./3.)*cfl + 1./10. );
79  return delta;
80 }
81 
82 REAL TPZArtDiff::Delta(REAL deltax, TPZVec<STATE> & sol)
83 {
84  REAL dX = 0.;
85 
86  if(fDelta > 0.)dX = deltax * fDelta;
87  if(fDelta < 0.)dX =-fDelta;
88  if(fDelta == 0.)dX = OptimalDelta();
89 
90  STATE us, c, lambdaMax;
91 
92  switch(fArtDiffType)
93  {
94  case SUPG_AD:
95  case Bornhaus_AD:
96  return dX / 2.;
97  // break;
98  case LeastSquares_AD:
99  case TrnLeastSquares_AD:
100  TPZEulerConsLaw::uRes(sol, us);
101  TPZEulerConsLaw::cSpeed(sol, fGamma, c);
102  lambdaMax = us+c;
103  return dX /2. / lambdaMax;
104  default:
105  PZError << "Unknown artificial diffusion term (" << fArtDiffType << ")";
106  return 0.;
107  }
108 }
109 
110 void TPZArtDiff::SetDelta(REAL delta)
111 {
112  fDelta=delta;
113 }
114 
115 REAL TPZArtDiff::OptimalCFL(int degree)
116 {
117  if(fCFL != 0.) return fCFL;
118  return (1.0/(2.0*(REAL)degree+1.0));
119 }
120 
121 template <class T>
122 void TPZArtDiff::ODotOperator(TPZVec<REAL> &dphi, TPZVec<TPZDiffMatrix<T> > &M, TPZDiffMatrix<T> &Result){
123 
124  int dim = M.NElements();
125  int size = dphi.NElements();
126  if(size<1 || size>3){
127  PZError << "TPZArtDiff::PointOperator: error data size";
128  }
129 
130  Result.Redim(M[0].fRows, M[0].fCols());
131 
132  int i;
133  for (i=0;i<dim;i++)Result.Add(M[i], dphi[i]);
134 }
135 
136 template <class T>
137 void TPZArtDiff::ODotOperator(TPZVec<REAL> &dphi, TPZVec<TPZVec<T> > &TauDiv, TPZVec<T> &Result){
138 
139  int dim = TauDiv.NElements();
140  int size = dphi.NElements();
141  int neq = TauDiv[0].NElements();
142  if(size<1 || size>3){
143  PZError << "TPZArtDiff::PointOperator: error data size";
144  }
145 
146  Result.Resize(neq);
147  Result.Fill(REAL(0.));
148 
149  int i, k;
150  for(k=0;k<dim;k++)
151  for(i=0;i<neq;i++)Result[i] += TauDiv[k][i] * dphi[k];
152 }
153 
154 #ifdef _AUTODIFF
155 
156 template <class T>
157 void TPZArtDiff::Divergent(TPZFMatrix<STATE> &dsol,
158  TPZFMatrix<REAL> & phi,
159  TPZFMatrix<REAL> & dphi,
160  TPZVec<TPZDiffMatrix<T> > & Ai,
161  TPZVec<STATE> & Div,
162  TPZDiffMatrix<STATE> * dDiv)
163 {
164  int nstate = Ai[0].Cols();
165  int dim = nstate - 2;
166  int nshape = dphi.Cols();
167  Div.Resize(nstate);
168  Div = (STATE(0.));
169 
170  int i, j, k;
171 
172  // computing the divergent:
173  // A.du/dx + B.du/dy + C.du/dz
174  for(k=0;k<dim;k++)
175  for(i=0;i<nstate; i++)
176  for(j=0;j<nstate;j++)
177  {
178  Div[i] += Ai[k](i,j).val() * dsol(k,j);
179  }
180 
181  if(!dDiv)return;
182 
183  // computing an approximation to the divergent derivative:
184 
185  // dDiv/dUj ~= A.d2U/dUidx + B.d2U/dUidy + C.d2U/dUidz
186 
187  dDiv->Redim(nstate, nstate * nshape);
188  int l;
189  REAL buff;
190  for(l=0;l<nshape;l++)
191  for(j=0;j<nstate;j++)
192  for(i=0;i<nstate; i++)
193  {
194  buff =0.;
195  for(k=0;k<dim;k++)
196  {
197  buff+=Ai[k](i,j).val()*dphi(k,l);
198  }
199  dDiv->operator()(i,j+l*nstate)=buff;
200  }
201 
202  TPZVec<T> ADiv(nstate);
203  T temp;
204  for( k = 0; k < dim; k++)
205  {
206  //Ai[k].Multiply(Div, ADiv);
207  for(i = 0; i < nstate; i++)
208  {
209  temp = T(0.);
210  for(j = 0; j < nstate; j++)
211  temp += Ai[k](i,j) * (T)dsol(k,j);//[j];
212  ADiv[i] = temp;
213  }
214  for(l=0;l<nshape;l++)
215  for(j=0;j<nstate;j++)
216  for(i=0;i<nstate; i++)
217  dDiv->operator()(i,j+l*nstate) +=
218  ADiv[i]./*fastAccessDx*/dx(j) * phi(l,0);
219  }
220 }
221 
222 #endif // _AUTODIFF
223 
224 void TPZArtDiff::Divergent(TPZFMatrix<STATE> &dsol,
225  TPZFMatrix<REAL> & dphi,
226  TPZVec<TPZDiffMatrix<STATE> > & Ai,
227  TPZVec<STATE> & Div,
228  TPZDiffMatrix<STATE> * dDiv)
229 {
230  int nstate = Ai[0].Cols();
231  int dim = nstate - 2;
232  int nshape = dphi.Cols();
233  Div.Resize(nstate);
234  Div.Fill(0.);
235 
236  int i, j, k;
237 
238  // computing the divergent:
239  // A.du/dx + B.du/dy + C.du/dz
240  for(k=0;k<dim;k++)
241  for(i=0;i<nstate; i++)
242  for(j=0;j<nstate;j++)
243  {
244  Div[i]+=Ai[k](i,j)*dsol(k,j);
245  }
246 
247  if(!dDiv)return;
248  // computing an approximation to the divergent derivative:
249 
250  // dDiv/dUj ~= A.d2U/dUidx + B.d2U/dUidy + C.d2U/dUidz
251  dDiv->Redim(nstate, nstate * nshape);
252  int l;
253  REAL buff;
254  for(l=0;l<nshape;l++)
255  for(j=0;j<nstate;j++)
256  for(i=0;i<nstate; i++)
257  {
258  buff =0.;
259  for(k=0;k<dim;k++)
260  {
261  buff+=Ai[k](i,j)*dphi(k,l);
262  }
263  dDiv->operator()(i,j+l*nstate)=buff;
264  }
265 
266 }
267 
268 
269 #ifdef _AUTODIFF
270 void TPZArtDiff::Divergent(TPZVec<FADREAL> &dsol,
271  TPZVec<TPZDiffMatrix<FADREAL> > & Ai,
272  TPZVec<FADREAL> & Div)
273 {
274  int nstate = Ai[0].Cols();
275  int dim = nstate - 2;
276  Div.Resize(nstate);
277 
278  int i, j, k;
279 
280  // computing the divergent:
281  // A.du/dx + B.du/dy + C.du/dz
282 
283  for(i=0;i<nstate; i++)
284  {
285  Div[i]=0.;
286  for(j=0;j<nstate;j++)
287  for(k=0;k<dim;k++)
288  Div[i]+=Ai[k](i,j)*dsol[j*dim+k];
289  }
290 }
291 
292 #endif
293 
294 //----------------------Tau tensor
295 
296 template <class T>
297 void TPZArtDiff::ComputeTau(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > &Ai,
298  TPZVec<TPZDiffMatrix<T> > &Tau)
299 {
300  Tau.Resize(dim);
301 
302  switch(fArtDiffType)
303  {
304  case SUPG_AD:
305  SUPG(dim, sol, Ai, Tau);
306  break;
307  case LeastSquares_AD:
308  LS(dim, sol, Ai, Tau);
309  break;
310  case Bornhaus_AD:
311  Bornhaus(dim, jacinv, sol, Ai, Tau);
312  break;
313  case TrnLeastSquares_AD:
314  LST(dim, sol, Ai, Tau);
315  break;
316  default:
317  PZError << "Unknown artificial diffision term (" << fArtDiffType << ")";
318  }
319 }
320 
321 
322 template <class T>
323 void TPZArtDiff::SUPG(int dim, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau){
324 
325 #ifdef FASTESTDIFF
326 
327  TPZDiffMatrix<T> RTM, RMi,
328  X, Xi,
329  Temp, INVA2B2,
330  LambdaSUPG;
331  T us, c;
332 
333  TPZEulerConsLaw::uRes(sol, us);
334  TPZEulerConsLaw::cSpeed(sol, 1.4, c);
335 
336  RMMatrix(sol, us, fGamma, RTM, RMi);
337 
338  EigenSystemSUPG(sol, us, c, fGamma, X, Xi, LambdaSUPG);
339 
340 
341  RTM. Multiply(X, Temp);
342  Temp. Multiply(LambdaSUPG, INVA2B2);
343  INVA2B2.Multiply(Xi, Temp);
344  Temp. Multiply(RMi, INVA2B2);
345 
346  for(int i = 0; i < Ai.NElements();i++)
347  {
348  Ai[i].Multiply(INVA2B2, Tau[i]);
349  }
350 #else
351 
352  TPZDiffMatrix<T> Rot, RotT,
353  X, Xi,
354  M, Mi,
355  Temp, INVA2B2,
356  LambdaSUPG;
357  T us, c;
358 
359  TPZEulerConsLaw::uRes(sol, us);
360  TPZEulerConsLaw::cSpeed(sol, 1.4, c);
361 
362  RotMatrix(sol, us, Rot, RotT);
363  MMatrix(sol, us, fGamma, M, Mi);
364  EigenSystemSUPG(sol, us, c, fGamma, X, Xi, LambdaSUPG);
365 
366  RotT. Multiply(M, Temp);
367  Temp. Multiply(X, INVA2B2);
368  INVA2B2.Multiply(LambdaSUPG, Temp);
369  Temp. Multiply(Xi, INVA2B2);
370  INVA2B2.Multiply(Mi, Temp);
371  Temp. Multiply(Rot, INVA2B2);
372 
373  for(int i = 0; i < Ai.NElements();i++)
374  {
375  Ai[i].Multiply(INVA2B2, Tau[i]);
376  }
377 
378 #endif
379 }
380 
381 template <class T>
382 void TPZArtDiff::LS(int dim, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau){
383  int i;
384  for(i = 0; i < dim; i++)
385  {
386  Ai[i].Transpose(Tau[i]);
387  }
388 }
389 
390 template <class T>
391 void TPZArtDiff::LST(int dim, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau){
392  int i;
393  for(i = 0; i < dim; i++)
394  {
395  Tau[i] = Ai[i];
396  }
397 }
398 
399 template <class T>
400 void TPZArtDiff::Bornhaus(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau){
401 
402 #ifdef FASTESTDIFF
403 
404  int i, j;
405  int nstate = Ai[0].Rows();
406 
407  TPZDiffMatrix<T> RTM, RMi,
408  Y, Yi,
409  Temp, Temp2,
410  BornhausTau(nstate, nstate),
411  LambdaBornhaus;
412  T us, c;
413  TPZVec<REAL> alphas(dim,0.);
414 
415  TPZEulerConsLaw::uRes(sol, us);
416  TPZEulerConsLaw::cSpeed(sol, 1.4, c);
417 
418  RMMatrix(sol, us, fGamma, RTM, RMi);
419 
420  BornhausTau.Redim(nstate, nstate);
421 
422  for( i = 0; i < dim; i++)
423  {
424  for( j = 0; j < dim; j++)
425  alphas[j] = jacinv(i,j);
426  ContributeBornhaus(sol, us, c, fGamma, alphas, BornhausTau);
427  }
428 
429  RTM. Multiply(BornhausTau, Temp);
430  Temp.Multiply(RMi, BornhausTau);
431 
432  for( i = 0; i < dim;i++)
433  {
434  Ai[i].Multiply(BornhausTau, Tau[i]);
435  }
436 
437 #else
438 
439  int i, j;
440  int nstate = Ai[0].Rows();
441 
442  TPZDiffMatrix<T> Rot, RotT,
443  Y, Yi,
444  M, Mi,
445  Temp, Temp2,
446  BornhausTau(nstate, nstate),
447  LambdaBornhaus;
448  T us, c;
449  TPZVec<REAL> alphas(dim,0.);
450 
451  TPZEulerConsLaw::uRes(sol, us);
452  TPZEulerConsLaw::cSpeed(sol, 1.4, c);
453 
454  RotMatrix(sol, us, Rot, RotT);
455  MMatrix(sol, us, fGamma, M, Mi);
456 
457  for( i = 0; i < dim; i++)
458  {
459  for( j = 0; j < dim; j++)
460  alphas[j] = jacinv(i,j);
461 
462  EigenSystemBornhaus(sol, us, c, fGamma, alphas,
463  Y, Yi, LambdaBornhaus);
464 
465  Y. Multiply(LambdaBornhaus, Temp);
466  Temp.Multiply(Yi, Temp2);
467 
468  BornhausTau.Add(Temp2);
469  }
470 
471  RotT. Multiply(M, Temp);
472  Temp. Multiply(BornhausTau, Temp2);
473  Temp2. Multiply(Mi, Temp);
474  Temp. Multiply(Rot, BornhausTau);
475 
476  BornhausTau.Inverse();
477  for( i = 0; i < dim;i++)
478  {
479  Ai[i].Multiply(BornhausTau, Tau[i]);
480  }
481 #endif
482 
483 }
484 
485 //------------------ Diff setup
486 template <class T>
487 void TPZArtDiff::PrepareDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<T> &U,
488  TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau)
489 {
490  TPZEulerConsLaw::JacobFlux(fGamma, dim, U, Ai);
491  ComputeTau(dim, jacinv, U, Ai, Tau);
492 }
493 
494 void TPZArtDiff::PrepareFastDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol,
495  TPZFMatrix<STATE> &dsol, TPZFMatrix<REAL> & dphi,
496  TPZVec<TPZVec<STATE> > & TauDiv, TPZVec<TPZDiffMatrix<STATE> > * pTaudDiv)
497 {
498  TPZVec<TPZDiffMatrix<STATE> > Ai;
499  TPZVec<TPZDiffMatrix<STATE> > Tau;
500 
501  TPZEulerConsLaw::JacobFlux(fGamma, dim, sol, Ai);
502  ComputeTau(dim, jacinv, sol, Ai, Tau);
503 
504  TPZVec<STATE> Div;
505 
506  TPZDiffMatrix<STATE> * pdDiv = NULL;
507  TPZDiffMatrix<STATE> dDiv;
508  if(pTaudDiv) pdDiv = & dDiv;
509 
510  //Computing the divergent
511  Divergent(dsol, dphi, Ai, Div, pdDiv);
512  TauDiv.Resize(dim);
513  if(pTaudDiv)pTaudDiv->Resize(dim);
514 
515  // computing Tau.Div = {Tx.Div, Ty.Div, Tz.Div}
516  // and Tau.dDiv = {Tx.dDiv, Ty.dDiv, Tz.dDiv}, if requested
517  int k;
518  for(k=0;k<dim;k++)
519  {
520  Tau[k].Multiply(Div, TauDiv[k]);
521  if(pTaudDiv)Tau[k].Multiply(dDiv, pTaudDiv->operator[](k));
522 
523  }
524 }
525 
526 #ifdef _AUTODIFF
527 void TPZArtDiff::PrepareFastDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<FADREAL> &sol,
528  TPZVec<FADREAL> &dsol, TPZVec<TPZVec<FADREAL> > & TauDiv)
529 {
530  TPZVec<TPZDiffMatrix<FADREAL> > Ai;
531  TPZVec<TPZDiffMatrix<FADREAL> > Tau;
532 
533  TPZEulerConsLaw::JacobFlux(fGamma, dim, sol, Ai);
534  ComputeTau(dim, jacinv, sol, Ai, Tau);
535 
536  // #define TEST_PARTIAL_DIFF
537  // Uncomment line above to zero the derivatives of
538  // Tensors Ai and Tau, so that partial and complete
539  // derivatives should be the same.
540 #ifdef TEST_PARTIAL_DIFF
541  //Zeroeing the derivatives... comparison tests for FAD and Approximate derivative methods
542  int i, j, t, l;
543  for(t=0;t<3;t++)
544  for(i=0;i<5;i++)
545  for(j=0;j<5;j++)
546  { Ai[t](i,j).diff(0,30);
547  Tau[t](i,j).diff(0,30);
548  for(l=0;l<30;l++)
549  {
550  Ai[t](i,j).dx/*fastAccessDx*/(l)=0.;
551  Tau[t](i,j).dx/*fastAccessDx*/(l)=0.;
552  }
553  }
554 #endif
555  TPZVec<FADREAL> Div;
556 
557  //Computing the divergent
558  Divergent(dsol, Ai, Div);
559  TauDiv.Resize(dim);
560 
561  // computing Tau.Div = {Tx.Div, Ty.Div, Tz.Div}
562  int k;
563  for(k=0;k<dim;k++)
564  Tau[k].Multiply(Div, TauDiv[k]);
565 }
566 
567 template <int dim>
568 void TPZArtDiff::PrepareFastestDiff(TPZFMatrix<REAL> &jacinv,
569  TPZVec<STATE> &sol,
570  TPZFMatrix<STATE> &dsol,
571  TPZFMatrix<REAL> &phi,
572  TPZFMatrix<REAL> &dphi,
573  TPZVec<TPZVec<STATE> > & TauDiv,
574  TPZVec<TPZDiffMatrix<STATE> > & dTauDiv)
575 {
576 #ifdef _TFAD
577  typedef TFad<dim+2, REAL> TFADREALdim;
578 #endif
579 #ifdef _FAD
580  typedef Fad<REAL> TFADREALdim;
581 #endif
582 #ifdef _TINYFAD
583  typedef TinyFad<dim+2, REAL> TFADREALdim;
584 #endif
585 
586  const int nstate = sol.NElements();
587  const int nshape = phi.Rows();
588 
589  TPZVec<TFADREALdim > FADsol(nstate);
590  TPZVec<TPZDiffMatrix<TFADREALdim> > FADAi(dim);
591  TPZVec<TPZDiffMatrix<REAL> > Ai(dim);
592  TPZVec<TPZDiffMatrix<TFADREALdim> > FADTau(dim);
593  TPZVec<TPZDiffMatrix<STATE> > Tau(dim);
594  TPZVec<TPZVec<TFADREALdim> > FADTauDiv(dim);
595  TFADREALdim temp;
596 
597  int i, j, k, l;
598  for(i = 0; i < nstate; i++)
599  {
600  FADsol[i] = sol[i];
601  FADsol[i].diff(i, nstate);
602  }
603 
604  TPZEulerConsLaw::JacobFlux(fGamma, dim, FADsol, FADAi);
605  ComputeTau(dim, jacinv, FADsol, FADAi, FADTau);
606 
607  for( k = 0; k < dim; k++)
608  {
609  Tau[k].Redim(nstate, nstate);
610  Ai [k].Redim(nstate, nstate);
611  for(i = 0; i < nstate; i++)
612  for( j = 0; j < nstate; j++)
613  {
614  Tau[k](i,j) = FADTau[k](i,j).val();
615  Ai [k](i,j) = FADAi [k](i,j).val();
616  }
617  }
618 
619  TPZVec<STATE> Div;
620  TPZDiffMatrix<STATE> dDiv;
621 
622  //Computing the divergent with derivatives
623  Divergent(dsol, phi, dphi, FADAi, Div, &dDiv);
624 
625  TauDiv. Resize(dim);
626  dTauDiv.Resize(dim);
627 
628  //Computing Tau * Div and DTau * Div
629  for(k=0;k<dim;k++)
630  {
631  TauDiv [k].Resize(nstate);
632  dTauDiv[k].Redim(nstate, nstate);
633  //FADTau[k].Multiply(Div, FADTauDiv[k]);
634  FADTauDiv[k].Resize(nstate);
635  for(i = 0; i < nstate; i++)
636  {
637  temp = 0.;
638  for(j = 0; j < nstate; j++)
639  temp += FADTau[k](i,j) * ((REAL)Div[j]);
640  FADTauDiv[k][i] = temp;
641  }//
642 
643  // copying data using REAL
644  dTauDiv[k].Redim(nstate, nstate * nshape);
645  for(i = 0; i < nstate; i++)
646  {
647  TauDiv[k][i] = FADTauDiv[k][i].val();
648  for(j = 0; j < nstate; j++)
649  {
650  for(l = 0; l < nshape; l++)
651  dTauDiv[k](i,l * nstate + j) =
652  FADTauDiv[k][i].dx/*fastAccessDx*/(j) * phi(l,0);
653  }
654  }
655 
656  }
657 
658  TPZDiffMatrix<STATE> TaudDiv_k; // temporary storage
659 
660  for(k = 0; k < dim; k++)
661  {
662  Tau[k].Multiply(dDiv, TaudDiv_k);
663  dTauDiv[k].Add(TaudDiv_k, 1.);
664  }
665 
666  //Computing DTauDiv = DTau * Div + Tau * DDiv
667 }
668 
669 
670 #endif
671 
672 
673 //-----------------Contribute
674 
675 void TPZArtDiff::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)
676 {
677  REAL delta = Delta(deltaX, sol);
678  REAL constant = /*-*/ weight * delta * timeStep;
679 
680  TPZVec<TPZVec<STATE> > TauDiv;
681  TPZVec<TPZDiffMatrix<STATE> > * pTaudDiv = NULL;
682  TPZVec<TPZDiffMatrix<STATE> > TaudDiv;
683 
684  pTaudDiv = & TaudDiv;
685 
686  PrepareFastDiff(dim, jacinv, sol, dsol, dphix, TauDiv, pTaudDiv);
687 
688  int i, j, k, l;
689  int nshape = dphix.Cols();
690  int nstate = dim + 2;
691  int neq = nstate*nshape;
692 
693  REAL buff;
694 
695  // ODotProduct speeded up
696  for(l=0;l<nshape;l++)
697  for(i=0;i<nstate;i++)
698  for(k=0;k<dim;k++)
699  {
700  buff = dphix(k,l) * constant;
701  ef(i+l*nstate,0) += buff * TauDiv[k][i];
702  for(j=0;j<neq;j++)
703  ek(i+l*nstate,j) -= buff * TaudDiv[k](i,j);
704  }
705 }
706 
707 void TPZArtDiff::ContributeExplDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol, TPZFMatrix<REAL> &dphix, TPZFMatrix<STATE> &ef, REAL weight, REAL timeStep, REAL deltaX)
708 {
709  REAL delta = Delta(deltaX, sol);
710  REAL constant = /*-*/ weight * delta * timeStep;
711 
712  TPZVec<TPZVec<STATE> > TauDiv;
713 
714  PrepareFastDiff(dim, jacinv, sol, dsol, dphix, TauDiv, NULL);
715 
716  int i, k, l;
717  int nshape = dphix.Cols();
718  int nstate = dim + 2;
719 
720  // ODotProduct speeded up
721  for(l=0;l<nshape;l++)
722  for(i=0;i<nstate;i++)
723  for(k=0;k<dim;k++)
724  ef(i+l*nstate,0) += dphix(k,l) * TauDiv[k][i] * constant;
725 }
726 
727 
728 #ifdef _AUTODIFF
729 
730 void TPZArtDiff::ContributeImplDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<FADREAL> &sol, TPZVec<FADREAL> &dsol, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, REAL weight, REAL timeStep, REAL deltaX)
731 {
732  TPZVec<STATE> solReal(sol.NElements());
733  int i;
734  for(i = 0; i < sol.NElements(); i++)
735  solReal[i] = sol[i].val();
736 
737  REAL delta = Delta(deltaX, solReal);
738  REAL constant = /*-*/ delta * weight * timeStep;
739 
740  TPZVec<TPZVec<FADREAL> > TauDiv;
741 
742  PrepareFastDiff(dim, jacinv, sol, dsol, TauDiv);
743 
744  TPZVec<FADREAL> Diff;
745  TPZVec<REAL> gradv(dim);
746 
747  int j, k, l;
748  int nstate = dim + 2;
749  int neq = sol[0].size();
750  int nshape = neq/nstate;
751 
752  for(l=0;l<nshape;l++)
753  {
754  for(k=0;k<dim;k++)
755  gradv[k] = dsol[k].dx/*fastAccessDx*/(/*k+*/l*nstate);// always retrieving this information from the first state variable...
756  ODotOperator(gradv, TauDiv, Diff);
757  for(i=0;i<nstate;i++)
758  {
759  ef(i+l*nstate,0) += constant * Diff[i].val();
760  for(j=0;j<neq;j++)
761  ek(i+l*nstate, j) -= constant * Diff[i].dx/*fastAccessDx*/(j);
762  }
763  }
764 }
765 
766 template <int dim>
767 void TPZArtDiff::ContributeFastestImplDiff_dim(TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, REAL weight, REAL timeStep, REAL deltaX)
768 {
769  REAL delta = Delta(deltaX, sol);
770  REAL constant = /*-*/ weight * delta * timeStep;
771  REAL buff;
772 
773  TPZVec<TPZVec<STATE> > TauDiv;
774  TPZVec<TPZDiffMatrix<STATE> > dTauDiv;
775 
776  PrepareFastestDiff<dim>( jacinv, sol, dsol, phi, dphi, TauDiv, dTauDiv);
777 
778  int i, j, k, l;
779  int nshape = dphi.Cols();
780  int nstate = dim + 2;
781  int neq = nstate * nshape;
782 
783  // ODotProduct speeded up
784 
785  for(l=0;l<nshape;l++)
786  for(i=0;i<nstate;i++)
787  for(k=0;k<dim;k++)
788  {
789  buff = dphi(k,l) * constant;
790  ef(i+l*nstate,0) += buff * TauDiv[k][i];
791  for(j=0;j<neq;j++)
792  ek(i+l*nstate,j) -= buff * dTauDiv[k](i,j);
793  }
794 }
795 
796 void TPZArtDiff::ContributeFastestImplDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, REAL weight, REAL timeStep, REAL deltaX)
797 {
798  switch(dim)
799  {
800  case(1):
801  ContributeFastestImplDiff_dim<1>(jacinv, sol, dsol,
802  phi, dphi, ek, ef,
803  weight, timeStep, deltaX);
804  break;
805  case(2):
806  ContributeFastestImplDiff_dim<2>(jacinv, sol, dsol,
807  phi, dphi, ek, ef,
808  weight, timeStep, deltaX);
809  break;
810  case(3):
811  ContributeFastestImplDiff_dim<3>(jacinv, sol, dsol,
812  phi, dphi, ek, ef,
813  weight, timeStep, deltaX);
814  break;
815  default:
816  PZError << "\nTPZArtDiff::ContributeFastestImplDiff unhandled dimension\n";
817  exit(-1);
818  }
819 }
820 
821 template< class T >
822 void TPZArtDiff::Pressure(REAL gamma, int dim, T& press, TPZVec<T> &U)
823 {
824  TPZEulerConsLaw::Pressure(gamma, dim, press, U);
825 }
826 
827 #endif
828 
829 void TPZArtDiff::Write(TPZStream &buf, int withclassid) const
830 {
831  int tmp = static_cast<int>(fArtDiffType);
832  buf.Write(&tmp,1);
833  buf.Write(&fGamma, 1);
834  buf.Write(&fDelta, 1);
835  buf.Write(&fCFL, 1);
836 }
837 
838 void TPZArtDiff::Read(TPZStream &buf, void *context)
839 {
840  int ArtDiffType;
841  buf.Read(&ArtDiffType, 1);
842  fArtDiffType = static_cast<TPZArtDiffType>(ArtDiffType);
843  buf.Read(&fGamma, 1);
844  buf.Read(&fDelta, 1);
845  buf.Read(&fCFL, 1);
846 }
847 
848 int TPZArtDiff::ClassId() const{
849  return Hash("TPZArtDiff");
850 }
851 
TPZDiffMatrix::Add
void Add(TPZDiffMatrix< T > &matrix, const T &scale=T(1.))
Adds element by element.
Definition: pzdiffmatrix.h:318
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
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
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
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...
TPZArtDiff::Write
void Write(TPZStream &buf, int withclassid) const override
Save the element data to a stream.
Definition: pzartdiff.cpp:829
Fad
Definition: fad.h:54
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
tmp
AutoPointerMutexArrayInit tmp
Definition: tpzautopointer.cpp:42
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
TPZEulerConsLaw::uRes
static void uRes(TPZVec< T > &sol, T &us)
Definition: pzeulerconslaw.h:1340
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
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
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
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
dx
expr_ dx(i) *cos(expr_.val())
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
TFad
Definition: tfad.h:64
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
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZDiffusionConsLaw.h
Contains the TPZDiffusionConsLaw class which implements a Euler equation where is introduced a diffus...
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
TPZArtDiff::fArtDiffType
TPZArtDiffType fArtDiffType
Kind of artificial diffusion term to apply.
Definition: pzartdiff.h:423
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
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
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
pzartdiff.h
Contains the TPZArtDiff class which implements a numerical diffusivity coefficient.
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
TPZEulerConsLaw::cSpeed
static void cSpeed(TPZVec< T > &sol, REAL gamma, T &c)
Evaluates the speed of sound in the fluid.
Definition: pzeulerconslaw.h:1318
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
TPZDiffMatrix::Multiply
void Multiply(TPZVec< T > &In, TPZVec< T > &Out, const T &scale=T(1.))
Multiplies the matrix by a correspondent TPZVec vector. Dimensions are checked.
Definition: pzdiffmatrix.h:256
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
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
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
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
pzreal.h
Contains the declaration of TPZFlopCounter class and TPZCounter struct.
TPZArtDiff::ContributeBornhaus
static void ContributeBornhaus(TPZVec< T > &sol, T &us, T &c, REAL gamma, TPZVec< REAL > &aaS, TPZDiffMatrix< T > &Mat)
Definition: pzartdiff.h:907
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
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
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13