proteus  1.8.1
C/C++/Fortran libraries
subgridError.c
Go to the documentation of this file.
1 #include "subgridError.h"
2 
12 void calculateSubgridError_tauRes(int nElements_global,
13  int nQuadraturePoints_element,
14  int nDOF_trial_element,
15  double* tau,
16  double* pdeResidual,
17  double* dpdeResidual,
18  double* subgridError,
19  double* dsubgridError)
20 {
21  int eN,k,j;
22  for(eN=0;eN<nElements_global;eN++)
23  for (k=0;k<nQuadraturePoints_element;k++)
24  {
25  subgridError[eN*nQuadraturePoints_element+
26  k] =
27  -tau[eN*nQuadraturePoints_element+
28  k]
29  *
30  pdeResidual[eN*nQuadraturePoints_element+
31  k];
32  for (j=0;j<nDOF_trial_element;j++)
33  dsubgridError[eN*nQuadraturePoints_element*nDOF_trial_element+
34  k*nDOF_trial_element+
35  j]
36  =
37  -tau[eN*nQuadraturePoints_element+
38  k]
39  *
40  dpdeResidual[eN*nQuadraturePoints_element*nDOF_trial_element+
41  k*nDOF_trial_element+
42  j];
43  }
44 }
45 
49 void calculateSubgridError_ADR_tau_p(int nElements_global,
50  int nQuadraturePoints_element,
51  int nSpace,
52  double* elementDiameter,
53  double* dmt,
54  double* df,
55  double* a,
56  double* da,
57  double* grad_phi,
58  double* dphi,
59  double* dr,
60  double* pe,
61  double* cfl,
62  double* tau)
63 {
64  int eN,k,I,J,nSpace2=nSpace*nSpace;
65  double h,Vlin,Alin,A_II,num,den,cfl1,cfl2,vlin[3],alpha,beta;
66  for(eN=0;eN<nElements_global;eN++)
67  {
68  h = elementDiameter[eN];
69  for (k=0;k<nQuadraturePoints_element;k++)
70  {
71  Vlin = 0.0;
72  Alin = 0.0;
73  for(I=0;I<nSpace;I++)
74  {
75  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
76  k*nSpace +
77  I];
78  for(J=0;J<nSpace;J++)
79  vlin[I]
80  -=
81  da[eN*nQuadraturePoints_element*nSpace2 +
82  k*nSpace2 +
83  I*nSpace +
84  J]
85  *
86  grad_phi[eN*nQuadraturePoints_element*nSpace +
87  k*nSpace +
88  J];
89  }
90  for(I=0;I<nSpace;I++)
91  Vlin += vlin[I]*vlin[I];
92  Vlin = sqrt(Vlin);
93  /*max diagonal entry for A. This choice needs to be looked into*/
94  for(I=0;I<nSpace;I++)
95  {
96  A_II = a[eN*nQuadraturePoints_element*nSpace2 +
97  k*nSpace2 +
98  I*nSpace +
99  I];
100  Alin = (A_II > Alin) ? A_II : Alin;
101  }
102  Alin*=dphi[eN*nQuadraturePoints_element +
103  k];
104  cfl1 = Vlin/h;
105  cfl2 = Alin/(h*h);
106  if (cfl1 > cfl2)
107  cfl[eN*nQuadraturePoints_element +
108  k] = cfl1;
109  else
110  cfl[eN*nQuadraturePoints_element +
111  k] = cfl2;
112  num = 0.5*Vlin*h + 1.0e-8;
113  den = Alin + num*1.0e-8;
114  pe[eN*nQuadraturePoints_element +
115  k] = num/den;
116  tau[eN*nQuadraturePoints_element + k]=0.5*h*(
117  1.0/tanh(pe[eN*nQuadraturePoints_element+
118  k]) -
119  1.0/pe[eN*nQuadraturePoints_element+
120  k])/(Vlin+1.0e-8);
121 /* alpha = pe[eN*nQuadraturePoints_element + */
122 /* k]; */
123 /* beta = Vlin; */
124 /* num = h*h*h*(-3.0*exp(2.0*alpha)/alpha+exp(2.0*alpha)+3.0*exp(2.0*alpha)/(alpha*alpha) - 3.0/(alpha*alpha) - 1.0 - 1.0/alpha); */
125 /* den = 72.0*beta*(exp(2.0*alpha) - exp(2.0*alpha)/alpha + 1.0 + 1.0/alpha); */
126 /* tau[eN*nQuadraturePoints_element + k]=num/(den+1.0e-8); */
127 /* printf("tau %12.5e \n",tau[eN*nQuadraturePoints_element+k]); */
128  }
129  }
130 }
131 
132 void calculateSubgridError_ADR_tau_p_sd(int nElements_global,
133  int nQuadraturePoints_element,
134  int nSpace,
135  int* rowptr,
136  int* colind,
137  double* elementDiameter,
138  double* dmt,
139  double* df,
140  double* a,
141  double* da,
142  double* grad_phi,
143  double* dphi,
144  double* dr,
145  double* pe,
146  double* cfl,
147  double* tau)
148 {
149  int eN,k,I,J,m,nnz=rowptr[nSpace];
150  double h,Vlin,Alin,A_II,num,den,cfl1,cfl2,vlin[3],alpha,beta;
151  for(eN=0;eN<nElements_global;eN++)
152  {
153  h = elementDiameter[eN];
154  for (k=0;k<nQuadraturePoints_element;k++)
155  {
156  Vlin = 0.0;
157  Alin = 0.0;
158  for(I=0;I<nSpace;I++)
159  {
160  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
161  k*nSpace +
162  I];
163  for(J=0;J<nSpace;J++)
164  vlin[I]
165  -=
166  da[eN*nQuadraturePoints_element*nnz+
167  k*nnz+
168  m]
169  *
170  grad_phi[eN*nQuadraturePoints_element*nSpace +
171  k*nSpace +
172  colind[m]];
173  }
174  for(I=0;I<nSpace;I++)
175  Vlin += vlin[I]*vlin[I];
176  Vlin = sqrt(Vlin);
177  /*max diagonal entry for A. This choice needs to be looked into*/
178  for(I=0;I<nSpace;I++)
179  {
180  for(m=rowptr[I];m<rowptr[I+1];m++)
181  {
182  if (I==colind[m])
183  {
184  A_II = a[eN*nQuadraturePoints_element*nnz+
185  k*nnz+
186  m];
187  Alin = (A_II > Alin) ? A_II : Alin;
188  }
189  }
190  }
191  Alin*=dphi[eN*nQuadraturePoints_element +
192  k];
193  cfl1 = Vlin/h;
194  cfl2 = Alin/(h*h);
195  if (cfl1 > cfl2)
196  cfl[eN*nQuadraturePoints_element +
197  k] = cfl1;
198  else
199  cfl[eN*nQuadraturePoints_element +
200  k] = cfl2;
201  num = 0.5*Vlin*h + 1.0e-8;
202  den = Alin + num*1.0e-8;
203  pe[eN*nQuadraturePoints_element +
204  k] = num/den;
205  tau[eN*nQuadraturePoints_element + k]=0.5*h*(
206  1.0/tanh(pe[eN*nQuadraturePoints_element+
207  k]) -
208  1.0/pe[eN*nQuadraturePoints_element+
209  k])/(Vlin+1.0e-8);
210 /* alpha = pe[eN*nQuadraturePoints_element + */
211 /* k]; */
212 /* beta = Vlin; */
213 /* num = h*h*h*(-3.0*exp(2.0*alpha)/alpha+exp(2.0*alpha)+3.0*exp(2.0*alpha)/(alpha*alpha) - 3.0/(alpha*alpha) - 1.0 - 1.0/alpha); */
214 /* den = 72.0*beta*(exp(2.0*alpha) - exp(2.0*alpha)/alpha + 1.0 + 1.0/alpha); */
215 /* tau[eN*nQuadraturePoints_element + k]=num/(den+1.0e-8); */
216 /* printf("tau %12.5e \n",tau[eN*nQuadraturePoints_element+k]); */
217  }
218  }
219 }
220 
224 void calculateSubgridError_ADR_tau_1(int nElements_global,
225  int nQuadraturePoints_element,
226  int nSpace,
227  double* elementDiameter,
228  double* dmt,
229  double* df,
230  double* a,
231  double* da,
232  double* grad_phi,
233  double* dphi,
234  double* dr,
235  double* pe,
236  double* cfl,
237  double* tau)
238 {
239  int eN,k,I,J,nSpace2=nSpace*nSpace;
240  double h,Vlin,Alin,A_II,num,den,cfl1,cfl2,vlin[3],tauMax;
241  for(eN=0;eN<nElements_global;eN++)
242  {
243  h = elementDiameter[eN];
244  tauMax=0.0;
245  for (k=0;k<nQuadraturePoints_element;k++)
246  {
247  Vlin = 0.0;
248  Alin = 0.0;
249  for(I=0;I<nSpace;I++)
250  {
251  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
252  k*nSpace +
253  I];
254  for(J=0;J<nSpace;J++)
255  vlin[I]
256  +=
257  da[eN*nQuadraturePoints_element*nSpace2 +
258  k*nSpace2 +
259  I*nSpace +
260  J]
261  *
262  grad_phi[eN*nQuadraturePoints_element*nSpace +
263  k*nSpace +
264  J];
265  }
266  for(I=0;I<nSpace;I++)
267  Vlin += vlin[I]*vlin[I];
268  Vlin = sqrt(Vlin);
269  /*max diagonal entry for A. This choice needs to be looked into*/
270  for(I=0;I<nSpace;I++)
271  {
272  A_II = a[eN*nQuadraturePoints_element*nSpace2 +
273  k*nSpace2 +
274  I*nSpace +
275  I];
276  Alin = (A_II > Alin) ? A_II : Alin;
277  }
278  Alin*=dphi[eN*nQuadraturePoints_element +
279  k];
280  cfl1 = Vlin/h;
281  cfl2 = Alin/(h*h);
282  if (cfl1 > cfl2)
283  cfl[eN*nQuadraturePoints_element +
284  k] = cfl1;
285  else
286  cfl[eN*nQuadraturePoints_element +
287  k] = cfl2;
288  num = 0.5*Vlin*h + 1.0e-8;
289  den = Alin + num*1.0e-8;
290  pe[eN*nQuadraturePoints_element +
291  k] = num/den;
292  tau[eN*nQuadraturePoints_element + k]=1.0/((2.0*Vlin/h)+
293  (4.0*Alin/(h*h)) +
294  fabs(dr[eN*nQuadraturePoints_element +
295  k])+
296  fabs(dmt[eN*nQuadraturePoints_element +
297  k])+1.0e-8);
298 /* tauMax=fmax(tauMax,tau[eN*nQuadraturePoints_element + k]); */
299  }
300 /* for (k=0;k<nQuadraturePoints_element;k++) */
301 /* tau[eN*nQuadraturePoints_element + k]=tauMax; */
302  }
303 }
304 void calculateSubgridError_ADR_tau_1_sd(int nElements_global,
305  int nQuadraturePoints_element,
306  int nSpace,
307  int* rowptr,
308  int* colind,
309  double* elementDiameter,
310  double* dmt,
311  double* df,
312  double* a,
313  double* da,
314  double* grad_phi,
315  double* dphi,
316  double* dr,
317  double* pe,
318  double* cfl,
319  double* tau)
320 {
321  int eN,k,I,m,nnz=rowptr[nSpace];
322  double h,Vlin,Alin,A_II,num,den,cfl1,cfl2,vlin[3],tauMax;
323  for(eN=0;eN<nElements_global;eN++)
324  {
325  h = elementDiameter[eN];
326  tauMax=0.0;
327  for (k=0;k<nQuadraturePoints_element;k++)
328  {
329  Vlin = 0.0;
330  Alin = 0.0;
331  for(I=0;I<nSpace;I++)
332  {
333  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
334  k*nSpace +
335  I];
336  for(m=rowptr[I];m<rowptr[I+1];m++)
337  vlin[I]
338  +=
339  da[eN*nQuadraturePoints_element*nnz+
340  k*nnz+
341  m]
342  *
343  grad_phi[eN*nQuadraturePoints_element*nSpace +
344  k*nSpace +
345  colind[m]];
346  }
347  for(I=0;I<nSpace;I++)
348  Vlin += vlin[I]*vlin[I];
349  Vlin = sqrt(Vlin);
350  /*max diagonal entry for A. This choice needs to be looked into*/
351  for(I=0;I<nSpace;I++)
352  {
353  for(m=rowptr[I];m<rowptr[I+1];m++)
354  {
355  if(I==colind[m])
356  {
357  A_II = a[eN*nQuadraturePoints_element*nnz+
358  k*nnz+
359  m];
360  Alin = (A_II > Alin) ? A_II : Alin;
361  }
362  }
363  }
364  Alin*=dphi[eN*nQuadraturePoints_element +
365  k];
366  cfl1 = Vlin/h;
367  cfl2 = Alin/(h*h);
368  if (cfl1 > cfl2)
369  cfl[eN*nQuadraturePoints_element +
370  k] = cfl1;
371  else
372  cfl[eN*nQuadraturePoints_element +
373  k] = cfl2;
374  num = 0.5*Vlin*h + 1.0e-8;
375  den = Alin + num*1.0e-8;
376  pe[eN*nQuadraturePoints_element +
377  k] = num/den;
378  tau[eN*nQuadraturePoints_element + k]=1.0/((2.0*Vlin/h)+
379  (4.0*Alin/(h*h)) +
380  fabs(dr[eN*nQuadraturePoints_element +
381  k])+
382  fabs(dmt[eN*nQuadraturePoints_element +
383  k])+1.0e-8);
384 /* tauMax=fmax(tauMax,tau[eN*nQuadraturePoints_element + k]); */
385  }
386 /* for (k=0;k<nQuadraturePoints_element;k++) */
387 /* tau[eN*nQuadraturePoints_element + k]=tauMax; */
388  }
389 }
390 
394 void calculateSubgridError_ADR_tau_2(int nElements_global,
395  int nQuadraturePoints_element,
396  int nSpace,
397  double* elementDiameter,
398  double* dmt,
399  double* df,
400  double* a,
401  double* da,
402  double* grad_phi,
403  double* dphi,
404  double* dr,
405  double* pe,
406  double* cfl,
407  double* tau)
408 {
409  int eN,k,I,J,nSpace2=nSpace*nSpace;
410  double h,Vlin,Alin,A_II,num,den,cfl1,cfl2,vlin[3];
411  for(eN=0;eN<nElements_global;eN++)
412  {
413  h = elementDiameter[eN];
414  for (k=0;k<nQuadraturePoints_element;k++)
415  {
416  Vlin = 0.0;
417  Alin = 0.0;
418  for(I=0;I<nSpace;I++)
419  {
420  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
421  k*nSpace +
422  I];
423  for(J=0;J<nSpace;J++)
424  {
425  vlin[I]
426  -=
427  da[eN*nQuadraturePoints_element*nSpace2 +
428  k*nSpace2 +
429  I*nSpace +
430  J]
431  *
432  grad_phi[eN*nQuadraturePoints_element*nSpace +
433  k*nSpace +
434  J];
435  }
436  }
437  for(I=0;I<nSpace;I++)
438  Vlin += vlin[I]*vlin[I];
439  Vlin = sqrt(Vlin);
440  /*max diagonal entry for A. This choice needs to be looked into*/
441  for(I=0;I<nSpace;I++)
442  {
443  A_II = a[eN*nQuadraturePoints_element*nSpace2 +
444  k*nSpace2 +
445  I*nSpace +
446  I];
447  Alin = (A_II > Alin) ? A_II : Alin;
448  }
449  Alin*=dphi[eN*nQuadraturePoints_element +
450  k];
451  cfl1 = Vlin/h;
452  cfl2 = Alin/(h*h);
453  if (cfl1 > cfl2)
454  cfl[eN*nQuadraturePoints_element +
455  k] = cfl1;
456  else
457  cfl[eN*nQuadraturePoints_element +
458  k] = cfl2;
459  num = 0.5*Vlin*h + 1.0e-8;
460  den = Alin + num*1.0e-8;
461 /* num = 0.5*Vlin*h; */
462 /* den = Alin; */
463  tau[eN*nQuadraturePoints_element +
464  k]
465  =1.0/sqrt((2.0*Vlin/h)*(2.0*Vlin/h)+
466  9.0*(4.0*Alin/(h*h))*(4.0*Alin/(h*h)) +
467  dr[eN*nQuadraturePoints_element +
468  k]
469  *
470  dr[eN*nQuadraturePoints_element +
471  k]
472  +
473  dmt[eN*nQuadraturePoints_element +
474  k]
475  *
476  dmt[eN*nQuadraturePoints_element +
477  k]+1.0e-8);
478  pe[eN*nQuadraturePoints_element +
479  k] = num/den;
480  }
481  }
482 }
483 
484 void calculateSubgridError_ADR_tau_2_sd(int nElements_global,
485  int nQuadraturePoints_element,
486  int nSpace,
487  int* rowptr,
488  int* colind,
489  double* elementDiameter,
490  double* dmt,
491  double* df,
492  double* a,
493  double* da,
494  double* grad_phi,
495  double* dphi,
496  double* dr,
497  double* pe,
498  double* cfl,
499  double* tau)
500 {
501  int eN,k,I,m,nnz=rowptr[nSpace];
502  double h,Vlin,Alin,A_II,num,den,cfl1,cfl2,vlin[3];
503  for(eN=0;eN<nElements_global;eN++)
504  {
505  h = elementDiameter[eN];
506  for (k=0;k<nQuadraturePoints_element;k++)
507  {
508  Vlin = 0.0;
509  Alin = 0.0;
510  for(I=0;I<nSpace;I++)
511  {
512  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
513  k*nSpace +
514  I];
515  for(m=rowptr[I];m<rowptr[I+1];m++)
516  {
517  vlin[I]
518  -=
519  da[eN*nQuadraturePoints_element*nnz+
520  k*nnz+
521  m]
522  *
523  grad_phi[eN*nQuadraturePoints_element*nSpace +
524  k*nSpace +
525  colind[m]];
526  }
527  }
528  for(I=0;I<nSpace;I++)
529  Vlin += vlin[I]*vlin[I];
530  Vlin = sqrt(Vlin);
531  /*max diagonal entry for A. This choice needs to be looked into*/
532  for(I=0;I<nSpace;I++)
533  {
534  for(m=rowptr[I];m<rowptr[I+1];m++)
535  {
536  if(I==colind[m])
537  {
538  A_II = a[eN*nQuadraturePoints_element*nnz+
539  k*nnz+
540  m];
541  Alin = (A_II > Alin) ? A_II : Alin;
542  }
543  }
544  }
545  Alin*=dphi[eN*nQuadraturePoints_element +
546  k];
547  cfl1 = Vlin/h;
548  cfl2 = Alin/(h*h);
549  if (cfl1 > cfl2)
550  cfl[eN*nQuadraturePoints_element +
551  k] = cfl1;
552  else
553  cfl[eN*nQuadraturePoints_element +
554  k] = cfl2;
555  num = 0.5*Vlin*h + 1.0e-8;
556  den = Alin + num*1.0e-8;
557 /* num = 0.5*Vlin*h; */
558 /* den = Alin; */
559  tau[eN*nQuadraturePoints_element +
560  k]
561  =1.0/sqrt((2.0*Vlin/h)*(2.0*Vlin/h)+
562  9.0*(4.0*Alin/(h*h))*(4.0*Alin/(h*h)) +
563  dr[eN*nQuadraturePoints_element +
564  k]
565  *
566  dr[eN*nQuadraturePoints_element +
567  k]
568  +
569  dmt[eN*nQuadraturePoints_element +
570  k]
571  *
572  dmt[eN*nQuadraturePoints_element +
573  k]+1.0e-8);
574  pe[eN*nQuadraturePoints_element +
575  k] = num/den;
576  }
577  }
578 }
579 
580 void calculateSubgridError_ADR_generic_tau(int nElements_global,
581  int nQuadraturePoints_element,
582  int nSpace,
583  double* inverseJ,
584  double* dmt,
585  double* df,
586  double* a,
587  double* da,
588  double* grad_phi,
589  double* dphi,
590  double* dr,
591  double* pe,
592  double* cfl,
593  double* tau)
594 {
595  int eN,k,I,J,K,nSpace2=nSpace*nSpace;
596  double Vlin,Alin,A_II,cfl1,cfl2,vlin[3],G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0*36.0,g_I,g_dot_g;
597  for(eN=0;eN<nElements_global;eN++)
598  {
599  for (k=0;k<nQuadraturePoints_element;k++)
600  {
601  Vlin = 0.0;
602  Alin = 0.0;
603  for(I=0;I<nSpace;I++)
604  {
605  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
606  k*nSpace +
607  I];
608  for(J=0;J<nSpace;J++)
609  {
610  vlin[I]
611  -=
612  da[eN*nQuadraturePoints_element*nSpace2 +
613  k*nSpace2 +
614  I*nSpace +
615  J]
616  *
617  grad_phi[eN*nQuadraturePoints_element*nSpace +
618  k*nSpace +
619  J];
620  }
621  }
622  v_dot_Gv=0.0;
623  CI_nu2_G_ddot_G=0.0;
624  g_dot_g=0.0;
625  for (I=0;I<nSpace;I++)
626  {
627  Gv_I=0.0;
628  g_I=0.0;
629  for (J=0;J<nSpace;J++)
630  {
631  G_IJ=0.0;
632  for (K=0;K<nSpace;K++)
633  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
634  k*nSpace2+
635  K*nSpace+
636  I]
637  *
638  inverseJ[eN*nQuadraturePoints_element*nSpace2+
639  k*nSpace2+
640  K*nSpace+
641  J];
642  Gv_I += G_IJ
643  *
644  vlin[J];
645  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
646  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
647  k*nSpace2+
648  J*nSpace+
649  I];
650  }
651  g_dot_g += g_I*g_I;
652  v_dot_Gv += vlin[I]
653  *
654  Gv_I;
655 
656  }
657  /*max diagonal entry for A. This choice needs to be looked into*/
658  for(I=0;I<nSpace;I++)
659  {
660  A_II = a[eN*nQuadraturePoints_element*nSpace2 +
661  k*nSpace2 +
662  I*nSpace +
663  I];
664  Alin = (A_II > Alin) ? A_II : Alin;
665  }
666  Alin*=dphi[eN*nQuadraturePoints_element +
667  k];
668  CI_nu2_G_ddot_G*=CI*Alin*Alin;
669  cfl1 = 0.5*sqrt(v_dot_Gv);
670  cfl2 = 0.25*sqrt(CI_nu2_G_ddot_G);
671  if (cfl1 > cfl2)
672  cfl[eN*nQuadraturePoints_element +
673  k] = cfl1;
674  else
675  cfl[eN*nQuadraturePoints_element +
676  k] = cfl2;
677  tau[eN*nQuadraturePoints_element +
678  k]
679  =1.0/sqrt(v_dot_Gv +
680  CI_nu2_G_ddot_G +
681  4.0*dr[eN*nQuadraturePoints_element +
682  k]
683  *
684  dr[eN*nQuadraturePoints_element +
685  k]
686  +
687  4.0*dmt[eN*nQuadraturePoints_element +
688  k]
689  *
690  dmt[eN*nQuadraturePoints_element +
691  k]+1.0e-8);
692 /* printf("tau %12.5e \n",tau[eN*nQuadraturePoints_element+k]); */
693  pe[eN*nQuadraturePoints_element +
694  k] = sqrt(v_dot_Gv)/sqrt(CI_nu2_G_ddot_G);
695  }
696  }
697 }
698 
699 void calculateSubgridError_ADR_generic_tau_sd(int nElements_global,
700  int nQuadraturePoints_element,
701  int nSpace,
702  int* rowptr,
703  int* colind,
704  double* inverseJ,
705  double* dmt,
706  double* df,
707  double* a,
708  double* da,
709  double* grad_phi,
710  double* dphi,
711  double* dr,
712  double* pe,
713  double* cfl,
714  double* tau)
715 {
716  int eN,k,I,J,K,m,nnz=rowptr[nSpace],nSpace2=nSpace*nSpace;
717  double Vlin,Alin,A_II,cfl1,cfl2,vlin[3],G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0*36.0,g_I,g_dot_g;
718  for(eN=0;eN<nElements_global;eN++)
719  {
720  for (k=0;k<nQuadraturePoints_element;k++)
721  {
722  Vlin = 0.0;
723  Alin = 0.0;
724  for(I=0;I<nSpace;I++)
725  {
726  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
727  k*nSpace +
728  I];
729  for(m=rowptr[I];m<rowptr[I+1];m++)
730  {
731  vlin[I]
732  -=
733  da[eN*nQuadraturePoints_element*nnz+
734  k*nnz+
735  m]
736  *
737  grad_phi[eN*nQuadraturePoints_element*nSpace +
738  k*nSpace +
739  colind[m]];
740  }
741  }
742  v_dot_Gv=0.0;
743  CI_nu2_G_ddot_G=0.0;
744  g_dot_g=0.0;
745  for (I=0;I<nSpace;I++)
746  {
747  Gv_I=0.0;
748  g_I=0.0;
749  for (J=0;J<nSpace;J++)
750  {
751  G_IJ=0.0;
752  for (K=0;K<nSpace;K++)
753  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
754  k*nSpace2+
755  K*nSpace+
756  I]
757  *
758  inverseJ[eN*nQuadraturePoints_element*nSpace2+
759  k*nSpace2+
760  K*nSpace+
761  J];
762  Gv_I += G_IJ
763  *
764  vlin[J];
765  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
766  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
767  k*nSpace2+
768  J*nSpace+
769  I];
770  }
771  g_dot_g += g_I*g_I;
772  v_dot_Gv += vlin[I]
773  *
774  Gv_I;
775 
776  }
777  /*max diagonal entry for A. This choice needs to be looked into*/
778  for(I=0;I<nSpace;I++)
779  {
780  for(m=rowptr[I];m<rowptr[I+1];m++)
781  {
782  if (I==colind[m])
783  {
784  A_II = a[eN*nQuadraturePoints_element*nnz +
785  k*nnz + m];
786  Alin = (A_II > Alin) ? A_II : Alin;
787  }
788  }
789  }
790  Alin*=dphi[eN*nQuadraturePoints_element +
791  k];
792  CI_nu2_G_ddot_G*=CI*Alin*Alin;
793  cfl1 = 0.5*sqrt(v_dot_Gv);
794  cfl2 = 0.25*sqrt(CI_nu2_G_ddot_G);
795  if (cfl1 > cfl2)
796  cfl[eN*nQuadraturePoints_element +
797  k] = cfl1;
798  else
799  cfl[eN*nQuadraturePoints_element +
800  k] = cfl2;
801  tau[eN*nQuadraturePoints_element +
802  k]
803  =1.0/sqrt(v_dot_Gv +
804  CI_nu2_G_ddot_G +
805  4.0*dr[eN*nQuadraturePoints_element +
806  k]
807  *
808  dr[eN*nQuadraturePoints_element +
809  k]
810  +
811  4.0*dmt[eN*nQuadraturePoints_element +
812  k]
813  *
814  dmt[eN*nQuadraturePoints_element +
815  k]+1.0e-8);
816 /* printf("tau %12.5e \n",tau[eN*nQuadraturePoints_element+k]); */
817  pe[eN*nQuadraturePoints_element +
818  k] = sqrt(v_dot_Gv)/sqrt(CI_nu2_G_ddot_G);
819  }
820  }
821 }
822 
826 void calculateSubgridError_ADR_tau(int nElements_global,
827  int nQuadraturePoints_element,
828  int nSpace,
829  char stabilization,
830  double* elementDiameter,
831  double* dmt,
832  double* df,
833  double* a,
834  double* da,
835  double* grad_phi,
836  double* dphi,
837  double* dr,
838  double* pe,
839  double* cfl,
840  double* tau)
841 {
842  if(stabilization == '2')
843  calculateSubgridError_ADR_tau_2(nElements_global, nQuadraturePoints_element, nSpace, elementDiameter, dmt, df, a, da, grad_phi, dphi, dr, pe, cfl, tau);
844  else if(stabilization == '1')
845  calculateSubgridError_ADR_tau_1(nElements_global, nQuadraturePoints_element, nSpace, elementDiameter, dmt, df, a, da, grad_phi, dphi, dr, pe, cfl, tau);
846  else if(stabilization == 'p')
847  calculateSubgridError_ADR_tau_p(nElements_global, nQuadraturePoints_element, nSpace, elementDiameter, dmt, df, a, da, grad_phi, dphi, dr, pe, cfl, tau);
848 
849 }
850 
851 void calculateSubgridError_ADR_tau_sd(int nElements_global,
852  int nQuadraturePoints_element,
853  int nSpace,
854  int* rowptr,
855  int* colind,
856  char stabilization,
857  double* elementDiameter,
858  double* dmt,
859  double* df,
860  double* a,
861  double* da,
862  double* grad_phi,
863  double* dphi,
864  double* dr,
865  double* pe,
866  double* cfl,
867  double* tau)
868 {
869  if(stabilization == '2')
870  calculateSubgridError_ADR_tau_2_sd(nElements_global, nQuadraturePoints_element, nSpace, rowptr,colind,elementDiameter, dmt, df, a, da, grad_phi, dphi, dr, pe, cfl, tau);
871  else if(stabilization == '1')
872  calculateSubgridError_ADR_tau_1_sd(nElements_global, nQuadraturePoints_element, nSpace, rowptr,colind,elementDiameter, dmt, df, a, da, grad_phi, dphi, dr, pe, cfl, tau);
873  else if(stabilization == 'p')
874  calculateSubgridError_ADR_tau_p_sd(nElements_global, nQuadraturePoints_element, nSpace, rowptr,colind,elementDiameter, dmt, df, a, da, grad_phi, dphi, dr, pe, cfl, tau);
875 }
876 
880 void calculateSubgridError_A_tau_1(int nElements_global,
881  int nQuadraturePoints_element,
882  int nSpace,
883  double* elementDiameter,
884  double* dmt,
885  double* df,
886  double* cfl,
887  double* tau)
888 {
889  int eN,k,I;
890  double h,Vlin;
891  for(eN=0;eN<nElements_global;eN++)
892  {
893  h = elementDiameter[eN];
894  for (k=0;k<nQuadraturePoints_element;k++)
895  {
896  Vlin = 0.0;
897  for(I=0;I<nSpace;I++)
898  {
899  Vlin +=
900  df[eN*nQuadraturePoints_element*nSpace +
901  k*nSpace +
902  I]
903  *
904  df[eN*nQuadraturePoints_element*nSpace +
905  k*nSpace +
906  I];
907  }
908  Vlin = sqrt(Vlin);
909  cfl[eN*nQuadraturePoints_element +
910  k] = Vlin/h;
911  tau[eN*nQuadraturePoints_element +
912  k]
913  =1.0/((2.0*Vlin/h)+
914  fabs(dmt[eN*nQuadraturePoints_element +
915  k])+1.0e-8);
916  }
917  }
918 }
919 
923 void calculateSubgridError_A_tau_2(int nElements_global,
924  int nQuadraturePoints_element,
925  int nSpace,
926  double* elementDiameter,
927  double* dmt,
928  double* df,
929  double* cfl,
930  double* tau)
931 {
932  int eN,k,I;
933  double h,Vlin;
934  for(eN=0;eN<nElements_global;eN++)
935  {
936  h = elementDiameter[eN];
937  for (k=0;k<nQuadraturePoints_element;k++)
938  {
939  Vlin = 0.0;
940  for(I=0;I<nSpace;I++)
941  {
942  Vlin +=
943  df[eN*nQuadraturePoints_element*nSpace +
944  k*nSpace +
945  I]
946  *
947  df[eN*nQuadraturePoints_element*nSpace +
948  k*nSpace +
949  I];
950  }
951  Vlin = sqrt(Vlin);
952  cfl[eN*nQuadraturePoints_element +
953  k] = Vlin/h;
954  tau[eN*nQuadraturePoints_element +
955  k]
956  =1.0/sqrt((2.0*Vlin/h)*(2.0*Vlin/h)+
957  dmt[eN*nQuadraturePoints_element +
958  k]
959  *
960  dmt[eN*nQuadraturePoints_element +
961  k]+1.0e-8);
962 /* printf("tau[%i,%i] = %12.5e \n",eN,k,tau[eN*nQuadraturePoints_element + */
963 /* k]); */
964  }
965  }
966 }
967 
971 void calculateSubgridError_A_tau(int nElements_global,
972  int nQuadraturePoints_element,
973  int nSpace,
974  char stabilization,
975  double* elementDiameter,
976  double* dmt,
977  double* df,
978  double* cfl,
979  double* tau)
980 {
981  if(stabilization == '2')
982  calculateSubgridError_A_tau_2(nElements_global, nQuadraturePoints_element, nSpace, elementDiameter, dmt, df, cfl, tau);
983  else if(stabilization == '1')
984  calculateSubgridError_A_tau_1(nElements_global, nQuadraturePoints_element, nSpace, elementDiameter, dmt, df, cfl, tau);
985 }
986 
990 void calculateSubgridErrorStokes2D_GLS_velocity(int nElements_global,
991  int nQuadraturePoints_element,
992  int nDOF_trial_element,
993  int nSpace,
994  double* elementDiameter,
995  double* a,
996  double* pdeResidualU,
997  double* dpdeResidualU_dp,
998  double* dpdeResidualU_du,
999  double* pdeResidualV,
1000  double* dpdeResidualV_dp,
1001  double* dpdeResidualV_dv,
1002  double* subgridErrorU,
1003  double* dsubgridErrorU_dp,
1004  double* dsubgridErrorU_du,
1005  double* subgridErrorV,
1006  double* dsubgridErrorV_dp,
1007  double* dsubgridErrorV_dv)
1008 {
1009  int eN,k,nSpace2=nSpace*nSpace,j;
1010  double h,viscosity,tau1;
1011  for(eN=0;eN<nElements_global;eN++)
1012  {
1013  h = elementDiameter[eN];
1014  for (k=0;k<nQuadraturePoints_element;k++)
1015  {
1016  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
1017  k*nSpace2];
1018  tau1 = h*h/(12.0*viscosity);
1019  subgridErrorU[eN*nQuadraturePoints_element+
1020  k] =
1021  tau1
1022  *
1023  pdeResidualU[eN*nQuadraturePoints_element+
1024  k];
1025  subgridErrorV[eN*nQuadraturePoints_element+
1026  k] =
1027  tau1
1028  *
1029  pdeResidualV[eN*nQuadraturePoints_element+
1030  k];
1031  for (j=0;j<nDOF_trial_element;j++)
1032  {
1033  /* u */
1034  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1035  k*nDOF_trial_element+
1036  j] =
1037  tau1
1038  *
1039  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1040  k*nDOF_trial_element+
1041  j];
1042  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1043  k*nDOF_trial_element+
1044  j] =
1045  tau1
1046  *
1047  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1048  k*nDOF_trial_element+
1049  j];
1050  /* v */
1051  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1052  k*nDOF_trial_element+
1053  j] =
1054  tau1
1055  *
1056  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1057  k*nDOF_trial_element+
1058  j];
1059  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1060  k*nDOF_trial_element+
1061  j] =
1062  tau1
1063  *
1064  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1065  k*nDOF_trial_element+
1066  j];
1067  }
1068  }
1069  }
1070 }
1071 
1072 void calculateSubgridErrorStokes2D_GLS_velocity_sd(int nElements_global,
1073  int nQuadraturePoints_element,
1074  int nDOF_trial_element,
1075  int nSpace,
1076  double* elementDiameter,
1077  double* a,
1078  double* pdeResidualU,
1079  double* dpdeResidualU_dp,
1080  double* dpdeResidualU_du,
1081  double* pdeResidualV,
1082  double* dpdeResidualV_dp,
1083  double* dpdeResidualV_dv,
1084  double* subgridErrorU,
1085  double* dsubgridErrorU_dp,
1086  double* dsubgridErrorU_du,
1087  double* subgridErrorV,
1088  double* dsubgridErrorV_dp,
1089  double* dsubgridErrorV_dv)
1090 {
1091  int eN,k,nSpace2=nSpace*nSpace,j;
1092  double h,viscosity,tau1;
1093  for(eN=0;eN<nElements_global;eN++)
1094  {
1095  h = elementDiameter[eN];
1096  for (k=0;k<nQuadraturePoints_element;k++)
1097  {
1098  viscosity = a[eN*nQuadraturePoints_element*nSpace +
1099  k*nSpace];
1100  tau1 = h*h/(12.0*viscosity);
1101  subgridErrorU[eN*nQuadraturePoints_element+
1102  k] =
1103  tau1
1104  *
1105  pdeResidualU[eN*nQuadraturePoints_element+
1106  k];
1107  subgridErrorV[eN*nQuadraturePoints_element+
1108  k] =
1109  tau1
1110  *
1111  pdeResidualV[eN*nQuadraturePoints_element+
1112  k];
1113  for (j=0;j<nDOF_trial_element;j++)
1114  {
1115  /* u */
1116  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1117  k*nDOF_trial_element+
1118  j] =
1119  tau1
1120  *
1121  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1122  k*nDOF_trial_element+
1123  j];
1124  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1125  k*nDOF_trial_element+
1126  j] =
1127  tau1
1128  *
1129  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1130  k*nDOF_trial_element+
1131  j];
1132  /* v */
1133  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1134  k*nDOF_trial_element+
1135  j] =
1136  tau1
1137  *
1138  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1139  k*nDOF_trial_element+
1140  j];
1141  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1142  k*nDOF_trial_element+
1143  j] =
1144  tau1
1145  *
1146  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1147  k*nDOF_trial_element+
1148  j];
1149  }
1150  }
1151  }
1152 }
1153 
1157 void calculateSubgridErrorStokes3D_GLS_velocity(int nElements_global,
1158  int nQuadraturePoints_element,
1159  int nDOF_trial_element,
1160  int nSpace,
1161  double* elementDiameter,
1162  double* a,
1163  double* pdeResidualU,
1164  double* dpdeResidualU_dp,
1165  double* dpdeResidualU_du,
1166  double* pdeResidualV,
1167  double* dpdeResidualV_dp,
1168  double* dpdeResidualV_dv,
1169  double* pdeResidualW,
1170  double* dpdeResidualW_dp,
1171  double* dpdeResidualW_dw,
1172  double* subgridErrorU,
1173  double* dsubgridErrorU_dp,
1174  double* dsubgridErrorU_du,
1175  double* subgridErrorV,
1176  double* dsubgridErrorV_dp,
1177  double* dsubgridErrorV_dv,
1178  double* subgridErrorW,
1179  double* dsubgridErrorW_dp,
1180  double* dsubgridErrorW_dw)
1181 {
1182  int eN,k,nSpace2=nSpace*nSpace,j;
1183  double h,viscosity,tau1;
1184  for(eN=0;eN<nElements_global;eN++)
1185  {
1186  h = elementDiameter[eN];
1187  for (k=0;k<nQuadraturePoints_element;k++)
1188  {
1189  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
1190  k*nSpace2];
1191  tau1 = h*h/(12.0*viscosity);
1192  subgridErrorU[eN*nQuadraturePoints_element+
1193  k] =
1194  tau1
1195  *
1196  pdeResidualU[eN*nQuadraturePoints_element+
1197  k];
1198  subgridErrorV[eN*nQuadraturePoints_element+
1199  k] =
1200  tau1
1201  *
1202  pdeResidualV[eN*nQuadraturePoints_element+
1203  k];
1204  subgridErrorW[eN*nQuadraturePoints_element+
1205  k] =
1206  tau1
1207  *
1208  pdeResidualW[eN*nQuadraturePoints_element+
1209  k];
1210  for (j=0;j<nDOF_trial_element;j++)
1211  {
1212  /* u */
1213  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1214  k*nDOF_trial_element+
1215  j] =
1216  tau1
1217  *
1218  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1219  k*nDOF_trial_element+
1220  j];
1221  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1222  k*nDOF_trial_element+
1223  j] =
1224  tau1
1225  *
1226  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1227  k*nDOF_trial_element+
1228  j];
1229  /* v */
1230  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1231  k*nDOF_trial_element+
1232  j] =
1233  tau1
1234  *
1235  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1236  k*nDOF_trial_element+
1237  j];
1238  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1239  k*nDOF_trial_element+
1240  j] =
1241  tau1
1242  *
1243  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1244  k*nDOF_trial_element+
1245  j];
1246  /* w */
1247  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1248  k*nDOF_trial_element+
1249  j] =
1250  tau1
1251  *
1252  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1253  k*nDOF_trial_element+
1254  j];
1255  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1256  k*nDOF_trial_element+
1257  j] =
1258  tau1
1259  *
1260  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1261  k*nDOF_trial_element+
1262  j];
1263  }
1264  }
1265  }
1266 }
1267 
1268 void calculateSubgridErrorStokes3D_GLS_velocity_sd(int nElements_global,
1269  int nQuadraturePoints_element,
1270  int nDOF_trial_element,
1271  int nSpace,
1272  double* elementDiameter,
1273  double* a,
1274  double* pdeResidualU,
1275  double* dpdeResidualU_dp,
1276  double* dpdeResidualU_du,
1277  double* pdeResidualV,
1278  double* dpdeResidualV_dp,
1279  double* dpdeResidualV_dv,
1280  double* pdeResidualW,
1281  double* dpdeResidualW_dp,
1282  double* dpdeResidualW_dw,
1283  double* subgridErrorU,
1284  double* dsubgridErrorU_dp,
1285  double* dsubgridErrorU_du,
1286  double* subgridErrorV,
1287  double* dsubgridErrorV_dp,
1288  double* dsubgridErrorV_dv,
1289  double* subgridErrorW,
1290  double* dsubgridErrorW_dp,
1291  double* dsubgridErrorW_dw)
1292 {
1293  int eN,k,nSpace2=nSpace*nSpace,j;
1294  double h,viscosity,tau1;
1295  for(eN=0;eN<nElements_global;eN++)
1296  {
1297  h = elementDiameter[eN];
1298  for (k=0;k<nQuadraturePoints_element;k++)
1299  {
1300  viscosity = a[eN*nQuadraturePoints_element*nSpace +
1301  k*nSpace];
1302  tau1 = h*h/(12.0*viscosity);
1303  subgridErrorU[eN*nQuadraturePoints_element+
1304  k] =
1305  tau1
1306  *
1307  pdeResidualU[eN*nQuadraturePoints_element+
1308  k];
1309  subgridErrorV[eN*nQuadraturePoints_element+
1310  k] =
1311  tau1
1312  *
1313  pdeResidualV[eN*nQuadraturePoints_element+
1314  k];
1315  subgridErrorW[eN*nQuadraturePoints_element+
1316  k] =
1317  tau1
1318  *
1319  pdeResidualW[eN*nQuadraturePoints_element+
1320  k];
1321  for (j=0;j<nDOF_trial_element;j++)
1322  {
1323  /* u */
1324  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1325  k*nDOF_trial_element+
1326  j] =
1327  tau1
1328  *
1329  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1330  k*nDOF_trial_element+
1331  j];
1332  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1333  k*nDOF_trial_element+
1334  j] =
1335  tau1
1336  *
1337  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1338  k*nDOF_trial_element+
1339  j];
1340  /* v */
1341  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1342  k*nDOF_trial_element+
1343  j] =
1344  tau1
1345  *
1346  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1347  k*nDOF_trial_element+
1348  j];
1349  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1350  k*nDOF_trial_element+
1351  j] =
1352  tau1
1353  *
1354  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1355  k*nDOF_trial_element+
1356  j];
1357  /* w */
1358  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1359  k*nDOF_trial_element+
1360  j] =
1361  tau1
1362  *
1363  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1364  k*nDOF_trial_element+
1365  j];
1366  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1367  k*nDOF_trial_element+
1368  j] =
1369  tau1
1370  *
1371  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1372  k*nDOF_trial_element+
1373  j];
1374  }
1375  }
1376  }
1377 }
1378 
1382 void calculateSubgridErrorStokes2D_GLS_velocity_pressure(int nElements_global,
1383  int nQuadraturePoints_element,
1384  int nDOF_trial_element,
1385  int nSpace,
1386  double* elementDiameter,
1387  double* a,
1388  double* pdeResidualP,
1389  double* dpdeResidualP_du,
1390  double* dpdeResidualP_dv,
1391  double* pdeResidualU,
1392  double* dpdeResidualU_dp,
1393  double* dpdeResidualU_du,
1394  double* pdeResidualV,
1395  double* dpdeResidualV_dp,
1396  double* dpdeResidualV_dv,
1397  double* subgridErrorP,
1398  double* dsubgridErrorP_du,
1399  double* dsubgridErrorP_dv,
1400  double* subgridErrorU,
1401  double* dsubgridErrorU_dp,
1402  double* dsubgridErrorU_du,
1403  double* subgridErrorV,
1404  double* dsubgridErrorV_dp,
1405  double* dsubgridErrorV_dv)
1406 {
1407  int eN,k,nSpace2=nSpace*nSpace,j;
1408  double h,viscosity,tau1,tau2;
1409  for(eN=0;eN<nElements_global;eN++)
1410  {
1411  h = elementDiameter[eN];
1412  for (k=0;k<nQuadraturePoints_element;k++)
1413  {
1414  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
1415  k*nSpace2];
1416  /* GLS momentum */
1417  tau1 = h*h/(12.0*viscosity);
1418  subgridErrorU[eN*nQuadraturePoints_element+
1419  k] =
1420  tau1
1421  *
1422  pdeResidualU[eN*nQuadraturePoints_element+
1423  k];
1424  subgridErrorV[eN*nQuadraturePoints_element+
1425  k] =
1426  tau1
1427  *
1428  pdeResidualV[eN*nQuadraturePoints_element+
1429  k];
1430  /* GLS pressure */
1431  tau2 = 6.0*viscosity;
1432  subgridErrorP[eN*nQuadraturePoints_element+
1433  k] =
1434  tau2
1435  *pdeResidualP[eN*nQuadraturePoints_element+
1436  k];
1437  for (j=0;j<nDOF_trial_element;j++)
1438  {
1439  /* GLS momentum*/
1440  /* u */
1441  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1442  k*nDOF_trial_element+
1443  j] =
1444  tau1
1445  *
1446  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1447  k*nDOF_trial_element+
1448  j];
1449  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1450  k*nDOF_trial_element+
1451  j] =
1452  tau1
1453  *
1454  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1455  k*nDOF_trial_element+
1456  j];
1457  /* v */
1458  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1459  k*nDOF_trial_element+
1460  j] =
1461  tau1
1462  *
1463  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1464  k*nDOF_trial_element+
1465  j];
1466  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1467  k*nDOF_trial_element+
1468  j] =
1469  tau1
1470  *
1471  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1472  k*nDOF_trial_element+
1473  j];
1474  /* GLS pressure */
1475  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1476  k*nDOF_trial_element+
1477  j]
1478  =
1479  tau2
1480  *
1481  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1482  k*nDOF_trial_element+
1483  j];
1484  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1485  k*nDOF_trial_element+
1486  j]
1487  =
1488  tau2
1489  *
1490  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1491  k*nDOF_trial_element+
1492  j];
1493  }
1494  }
1495  }
1496 }
1497 
1498 void calculateSubgridErrorStokes2D_GLS_velocity_pressure_sd(int nElements_global,
1499  int nQuadraturePoints_element,
1500  int nDOF_trial_element,
1501  int nSpace,
1502  double* elementDiameter,
1503  double* a,
1504  double* pdeResidualP,
1505  double* dpdeResidualP_du,
1506  double* dpdeResidualP_dv,
1507  double* pdeResidualU,
1508  double* dpdeResidualU_dp,
1509  double* dpdeResidualU_du,
1510  double* pdeResidualV,
1511  double* dpdeResidualV_dp,
1512  double* dpdeResidualV_dv,
1513  double* subgridErrorP,
1514  double* dsubgridErrorP_du,
1515  double* dsubgridErrorP_dv,
1516  double* subgridErrorU,
1517  double* dsubgridErrorU_dp,
1518  double* dsubgridErrorU_du,
1519  double* subgridErrorV,
1520  double* dsubgridErrorV_dp,
1521  double* dsubgridErrorV_dv)
1522 {
1523  int eN,k,nSpace2=nSpace*nSpace,j;
1524  double h,viscosity,tau1,tau2;
1525  for(eN=0;eN<nElements_global;eN++)
1526  {
1527  h = elementDiameter[eN];
1528  for (k=0;k<nQuadraturePoints_element;k++)
1529  {
1530  viscosity = a[eN*nQuadraturePoints_element*nSpace +
1531  k*nSpace];
1532  /* GLS momentum */
1533  tau1 = h*h/(12.0*viscosity);
1534  subgridErrorU[eN*nQuadraturePoints_element+
1535  k] =
1536  tau1
1537  *
1538  pdeResidualU[eN*nQuadraturePoints_element+
1539  k];
1540  subgridErrorV[eN*nQuadraturePoints_element+
1541  k] =
1542  tau1
1543  *
1544  pdeResidualV[eN*nQuadraturePoints_element+
1545  k];
1546  /* GLS pressure */
1547  tau2 = 6.0*viscosity;
1548  subgridErrorP[eN*nQuadraturePoints_element+
1549  k] =
1550  tau2
1551  *pdeResidualP[eN*nQuadraturePoints_element+
1552  k];
1553  for (j=0;j<nDOF_trial_element;j++)
1554  {
1555  /* GLS momentum*/
1556  /* u */
1557  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1558  k*nDOF_trial_element+
1559  j] =
1560  tau1
1561  *
1562  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1563  k*nDOF_trial_element+
1564  j];
1565  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1566  k*nDOF_trial_element+
1567  j] =
1568  tau1
1569  *
1570  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1571  k*nDOF_trial_element+
1572  j];
1573  /* v */
1574  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1575  k*nDOF_trial_element+
1576  j] =
1577  tau1
1578  *
1579  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1580  k*nDOF_trial_element+
1581  j];
1582  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1583  k*nDOF_trial_element+
1584  j] =
1585  tau1
1586  *
1587  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1588  k*nDOF_trial_element+
1589  j];
1590  /* GLS pressure */
1591  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1592  k*nDOF_trial_element+
1593  j]
1594  =
1595  tau2
1596  *
1597  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1598  k*nDOF_trial_element+
1599  j];
1600  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1601  k*nDOF_trial_element+
1602  j]
1603  =
1604  tau2
1605  *
1606  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1607  k*nDOF_trial_element+
1608  j];
1609  }
1610  }
1611  }
1612 }
1613 
1617 void calculateSubgridErrorStokes3D_GLS_velocity_pressure(int nElements_global,
1618  int nQuadraturePoints_element,
1619  int nDOF_trial_element,
1620  int nSpace,
1621  double* elementDiameter,
1622  double* a,
1623  double* pdeResidualP,
1624  double* dpdeResidualP_du,
1625  double* dpdeResidualP_dv,
1626  double* dpdeResidualP_dw,
1627  double* pdeResidualU,
1628  double* dpdeResidualU_dp,
1629  double* dpdeResidualU_du,
1630  double* pdeResidualV,
1631  double* dpdeResidualV_dp,
1632  double* dpdeResidualV_dv,
1633  double* pdeResidualW,
1634  double* dpdeResidualW_dp,
1635  double* dpdeResidualW_dw,
1636  double* subgridErrorP,
1637  double* dsubgridErrorP_du,
1638  double* dsubgridErrorP_dv,
1639  double* dsubgridErrorP_dw,
1640  double* subgridErrorU,
1641  double* dsubgridErrorU_dp,
1642  double* dsubgridErrorU_du,
1643  double* subgridErrorV,
1644  double* dsubgridErrorV_dp,
1645  double* dsubgridErrorV_dv,
1646  double* subgridErrorW,
1647  double* dsubgridErrorW_dp,
1648  double* dsubgridErrorW_dw)
1649 {
1650  int eN,k,nSpace2=nSpace*nSpace,j;
1651  double h,viscosity,tau1,tau2;
1652  for(eN=0;eN<nElements_global;eN++)
1653  {
1654  h = elementDiameter[eN];
1655  for (k=0;k<nQuadraturePoints_element;k++)
1656  {
1657  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
1658  k*nSpace2];
1659  /* GLS momentum */
1660  tau1 = h*h/(12.0*viscosity);
1661  subgridErrorU[eN*nQuadraturePoints_element+
1662  k] =
1663  tau1
1664  *
1665  pdeResidualU[eN*nQuadraturePoints_element+
1666  k];
1667  subgridErrorV[eN*nQuadraturePoints_element+
1668  k] =
1669  tau1
1670  *
1671  pdeResidualV[eN*nQuadraturePoints_element+
1672  k];
1673  subgridErrorW[eN*nQuadraturePoints_element+
1674  k] =
1675  tau1
1676  *
1677  pdeResidualW[eN*nQuadraturePoints_element+
1678  k];
1679  /* GLS pressure */
1680  tau2 = 6.0*viscosity;
1681  subgridErrorP[eN*nQuadraturePoints_element+
1682  k] =
1683  tau2
1684  *pdeResidualP[eN*nQuadraturePoints_element+
1685  k];
1686  for (j=0;j<nDOF_trial_element;j++)
1687  {
1688  /* GLS momentum*/
1689  /* u */
1690  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1691  k*nDOF_trial_element+
1692  j] =
1693  tau1
1694  *
1695  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1696  k*nDOF_trial_element+
1697  j];
1698  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1699  k*nDOF_trial_element+
1700  j] =
1701  tau1
1702  *
1703  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1704  k*nDOF_trial_element+
1705  j];
1706  /* v */
1707  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1708  k*nDOF_trial_element+
1709  j] =
1710  tau1
1711  *
1712  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1713  k*nDOF_trial_element+
1714  j];
1715  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1716  k*nDOF_trial_element+
1717  j] =
1718  tau1
1719  *
1720  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1721  k*nDOF_trial_element+
1722  j];
1723  /* w */
1724  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1725  k*nDOF_trial_element+
1726  j] =
1727  tau1
1728  *
1729  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1730  k*nDOF_trial_element+
1731  j];
1732  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1733  k*nDOF_trial_element+
1734  j] =
1735  tau1
1736  *
1737  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1738  k*nDOF_trial_element+
1739  j];
1740  /* GLS pressure */
1741  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1742  k*nDOF_trial_element+
1743  j]
1744  =
1745  tau2
1746  *
1747  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1748  k*nDOF_trial_element+
1749  j];
1750  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1751  k*nDOF_trial_element+
1752  j]
1753  =
1754  tau2
1755  *
1756  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1757  k*nDOF_trial_element+
1758  j];
1759  }
1760  }
1761  }
1762 }
1763 void calculateSubgridErrorStokes3D_GLS_velocity_pressure_sd(int nElements_global,
1764  int nQuadraturePoints_element,
1765  int nDOF_trial_element,
1766  int nSpace,
1767  double* elementDiameter,
1768  double* a,
1769  double* pdeResidualP,
1770  double* dpdeResidualP_du,
1771  double* dpdeResidualP_dv,
1772  double* dpdeResidualP_dw,
1773  double* pdeResidualU,
1774  double* dpdeResidualU_dp,
1775  double* dpdeResidualU_du,
1776  double* pdeResidualV,
1777  double* dpdeResidualV_dp,
1778  double* dpdeResidualV_dv,
1779  double* pdeResidualW,
1780  double* dpdeResidualW_dp,
1781  double* dpdeResidualW_dw,
1782  double* subgridErrorP,
1783  double* dsubgridErrorP_du,
1784  double* dsubgridErrorP_dv,
1785  double* dsubgridErrorP_dw,
1786  double* subgridErrorU,
1787  double* dsubgridErrorU_dp,
1788  double* dsubgridErrorU_du,
1789  double* subgridErrorV,
1790  double* dsubgridErrorV_dp,
1791  double* dsubgridErrorV_dv,
1792  double* subgridErrorW,
1793  double* dsubgridErrorW_dp,
1794  double* dsubgridErrorW_dw)
1795 {
1796  int eN,k,nSpace2=nSpace*nSpace,j;
1797  double h,viscosity,tau1,tau2;
1798  for(eN=0;eN<nElements_global;eN++)
1799  {
1800  h = elementDiameter[eN];
1801  for (k=0;k<nQuadraturePoints_element;k++)
1802  {
1803  viscosity = a[eN*nQuadraturePoints_element*nSpace +
1804  k*nSpace];
1805  /* GLS momentum */
1806  tau1 = h*h/(12.0*viscosity);
1807  subgridErrorU[eN*nQuadraturePoints_element+
1808  k] =
1809  tau1
1810  *
1811  pdeResidualU[eN*nQuadraturePoints_element+
1812  k];
1813  subgridErrorV[eN*nQuadraturePoints_element+
1814  k] =
1815  tau1
1816  *
1817  pdeResidualV[eN*nQuadraturePoints_element+
1818  k];
1819  subgridErrorW[eN*nQuadraturePoints_element+
1820  k] =
1821  tau1
1822  *
1823  pdeResidualW[eN*nQuadraturePoints_element+
1824  k];
1825  /* GLS pressure */
1826  tau2 = 6.0*viscosity;
1827  subgridErrorP[eN*nQuadraturePoints_element+
1828  k] =
1829  tau2
1830  *pdeResidualP[eN*nQuadraturePoints_element+
1831  k];
1832  for (j=0;j<nDOF_trial_element;j++)
1833  {
1834  /* GLS momentum*/
1835  /* u */
1836  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1837  k*nDOF_trial_element+
1838  j] =
1839  tau1
1840  *
1841  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1842  k*nDOF_trial_element+
1843  j];
1844  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1845  k*nDOF_trial_element+
1846  j] =
1847  tau1
1848  *
1849  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1850  k*nDOF_trial_element+
1851  j];
1852  /* v */
1853  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1854  k*nDOF_trial_element+
1855  j] =
1856  tau1
1857  *
1858  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1859  k*nDOF_trial_element+
1860  j];
1861  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1862  k*nDOF_trial_element+
1863  j] =
1864  tau1
1865  *
1866  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1867  k*nDOF_trial_element+
1868  j];
1869  /* w */
1870  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1871  k*nDOF_trial_element+
1872  j] =
1873  tau1
1874  *
1875  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
1876  k*nDOF_trial_element+
1877  j];
1878  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1879  k*nDOF_trial_element+
1880  j] =
1881  tau1
1882  *
1883  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
1884  k*nDOF_trial_element+
1885  j];
1886  /* GLS pressure */
1887  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1888  k*nDOF_trial_element+
1889  j]
1890  =
1891  tau2
1892  *
1893  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
1894  k*nDOF_trial_element+
1895  j];
1896  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1897  k*nDOF_trial_element+
1898  j]
1899  =
1900  tau2
1901  *
1902  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
1903  k*nDOF_trial_element+
1904  j];
1905  }
1906  }
1907  }
1908 }
1909 
1910 
1914 void calculateSubgridErrorNavierStokes2D_GLS_tau(int nElements_global,
1915  int nQuadraturePoints_element,
1916  int nSpace,
1917  double hFactor,
1918  double* elementDiameter,
1919  double* dmt,
1920  double* dm,
1921  double* f,
1922  double* a,
1923  double* tau0,
1924  double* tau1,
1925  double* cfl)
1926 {
1927  int eN,k,nSpace2=nSpace*nSpace,I;
1928  double h,oneByAbsdt,density,viscosity,Re_max=0.0,CFL_max=0.0,nrm_v;
1929  for(eN=0;eN<nElements_global;eN++)
1930  {
1931  h = hFactor*elementDiameter[eN];
1932  //printf("h %12.5e \n",h);
1933  for (k=0;k<nQuadraturePoints_element;k++)
1934  {
1935  density = dm[eN*nQuadraturePoints_element+
1936  k];
1937  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
1938  k*nSpace2+
1939  nSpace+1];
1940  nrm_v=0.0;
1941  for(I=0;I<nSpace;I++)
1942  nrm_v+=f[eN*nQuadraturePoints_element*nSpace+
1943  k*nSpace+
1944  I]*
1945  f[eN*nQuadraturePoints_element*nSpace+
1946  k*nSpace+
1947  I];
1948  nrm_v = sqrt(nrm_v);
1949  Re_max = fmax(nrm_v*h/(viscosity/density),Re_max);
1950  CFL_max = fmax(nrm_v/h,CFL_max);
1951  cfl[eN*nQuadraturePoints_element+k] = nrm_v/h;
1952  if (0)//abs(dmt[eN*nQuadraturePoints_element+k]) > 10.0*(density*nrm_v/h))
1953  {
1954  oneByAbsdt = 10.0*(density*nrm_v/h);
1955  }
1956  else
1957  {
1958  oneByAbsdt = fabs(dmt[eN*nQuadraturePoints_element+k]);
1959  }
1960  tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) +
1961  2.0*density*nrm_v/h +
1962  oneByAbsdt);
1963  tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity +
1964  2.0*density*nrm_v*h+
1965  oneByAbsdt*h*h;
1966  printf("nrm_v %12.5e tau_v %12.5e tau_p %12.5e \n",nrm_v,tau0[eN*nQuadraturePoints_element+k],tau1[eN*nQuadraturePoints_element+k]);
1967 /* tau0[eN*nQuadraturePoints_element+k] = density/(4.0*viscosity/(h*h) + */
1968 /* 2.0*density*nrm_v/h + */
1969 /* oneByAbsdt); */
1970 /* tau1[eN*nQuadraturePoints_element+k] = (4.0*viscosity + */
1971 /* 2.0*density*nrm_v*h+ */
1972 /* oneByAbsdt*h*h)/density; */
1973  //printf("tau0 %12.5e \n",tau0[eN*nQuadraturePoints_element+k]);
1974  //printf("tau1 %12.5e \n",tau1[eN*nQuadraturePoints_element+k]);
1975 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
1976 /* 2.0*density*nrm_v/h); */
1977 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + 2.0*density*nrm_v*h; */
1978 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
1979 /* 2.0*density*nrm_v/h); */
1980 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
1981 /* 2.0*density*nrm_v*h; */
1982 /* if (fabs(dmt[eN*nQuadraturePoints_element+k]) < nrm_v/h) */
1983 /* { */
1984 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
1985 /* 2.0*density*nrm_v/h + */
1986 /* fabs(dmt[eN*nQuadraturePoints_element+k])); */
1987 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
1988 /* 2.0*density*nrm_v*h+ */
1989 /* fabs(dmt[eN*nQuadraturePoints_element+k])*h*h; */
1990 /* } */
1991 /* else */
1992 /* { */
1993 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) +2.0*density*nrm_v/h); */
1994 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + 2.0*density*nrm_v*h; */
1995 /* } */
1996 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
1997 /* 2.0*density*nrm_v/h); */
1998 /* if (fabs(dmt[eN*nQuadraturePoints_element+k]) > nrm_v/h) */
1999 /* { */
2000 /* tau1[eN*nQuadraturePoints_element+k] = (4.0*viscosity + 2.0*density*nrm_v*h); */
2001 /* } */
2002 /* else */
2003 /* tau1[eN*nQuadraturePoints_element+k] = h*h/tau0[eN*nQuadraturePoints_element+k]; */
2004 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2005 /* 2.0*density*nrm_v/h + */
2006 /* fabs(dmt[eN*nQuadraturePoints_element+k])); */
2007 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
2008 /* 2.0*density*nrm_v*h; */
2009 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2010 /* 2.0*density*nrm_v/h); */
2011 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
2012 /* 2.0*density*nrm_v*h+ */
2013 /* fabs(dmt[eN*nQuadraturePoints_element+k])*h*h; */
2014  }
2015  }
2016  /* printf("Element Re = %12.5e \n",Re_max);
2017  printf("Element CFL = %12.5e \n",CFL_max);*/
2018 }
2019 
2020 void calculateSubgridErrorNavierStokes2D_GLS_tau_sd(int nElements_global,
2021  int nQuadraturePoints_element,
2022  int nSpace,
2023  double hFactor,
2024  double* elementDiameter,
2025  double* dmt,
2026  double* dm,
2027  double* f,
2028  double* a,
2029  double* tau0,
2030  double* tau1,
2031  double* cfl)
2032 {
2033  int eN,k,nSpace2=nSpace*nSpace,I;
2034  double h,oneByAbsdt,density,viscosity,Re_max=0.0,CFL_max=0.0,nrm_v;
2035  for(eN=0;eN<nElements_global;eN++)
2036  {
2037  h = hFactor*elementDiameter[eN];
2038  //printf("h %12.5e \n",h);
2039  for (k=0;k<nQuadraturePoints_element;k++)
2040  {
2041  density = dm[eN*nQuadraturePoints_element+
2042  k];
2043  viscosity = a[eN*nQuadraturePoints_element*nSpace +
2044  k*nSpace+1];
2045  nrm_v=0.0;
2046  for(I=0;I<nSpace;I++)
2047  nrm_v+=f[eN*nQuadraturePoints_element*nSpace+
2048  k*nSpace+
2049  I]*
2050  f[eN*nQuadraturePoints_element*nSpace+
2051  k*nSpace+
2052  I];
2053  nrm_v = sqrt(nrm_v);
2054  Re_max = fmax(nrm_v*h/(viscosity/density),Re_max);
2055  CFL_max = fmax(nrm_v/h,CFL_max);
2056  cfl[eN*nQuadraturePoints_element+k] = nrm_v/h;
2057  if (0)//abs(dmt[eN*nQuadraturePoints_element+k]) > 10.0*(density*nrm_v/h))
2058  {
2059  oneByAbsdt = 10.0*(density*nrm_v/h);
2060  }
2061  else
2062  {
2063  oneByAbsdt = fabs(dmt[eN*nQuadraturePoints_element+k]);
2064  }
2065  tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) +
2066  2.0*density*nrm_v/h +
2067  oneByAbsdt);
2068  tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity +
2069  2.0*density*nrm_v*h+
2070  oneByAbsdt*h*h;
2071  /* printf("nrm_v %12.5e tau_v %12.5e tau_p %12.5e \n",nrm_v,tau0[eN*nQuadraturePoints_element+k],tau1[eN*nQuadraturePoints_element+k]); */
2072 /* tau0[eN*nQuadraturePoints_element+k] = density/(4.0*viscosity/(h*h) + */
2073 /* 2.0*density*nrm_v/h + */
2074 /* oneByAbsdt); */
2075 /* tau1[eN*nQuadraturePoints_element+k] = (4.0*viscosity + */
2076 /* 2.0*density*nrm_v*h+ */
2077 /* oneByAbsdt*h*h)/density; */
2078  //printf("tau0 %12.5e \n",tau0[eN*nQuadraturePoints_element+k]);
2079  //printf("tau1 %12.5e \n",tau1[eN*nQuadraturePoints_element+k]);
2080 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2081 /* 2.0*density*nrm_v/h); */
2082 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + 2.0*density*nrm_v*h; */
2083 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2084 /* 2.0*density*nrm_v/h); */
2085 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
2086 /* 2.0*density*nrm_v*h; */
2087 /* if (fabs(dmt[eN*nQuadraturePoints_element+k]) < nrm_v/h) */
2088 /* { */
2089 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2090 /* 2.0*density*nrm_v/h + */
2091 /* fabs(dmt[eN*nQuadraturePoints_element+k])); */
2092 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
2093 /* 2.0*density*nrm_v*h+ */
2094 /* fabs(dmt[eN*nQuadraturePoints_element+k])*h*h; */
2095 /* } */
2096 /* else */
2097 /* { */
2098 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) +2.0*density*nrm_v/h); */
2099 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + 2.0*density*nrm_v*h; */
2100 /* } */
2101 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2102 /* 2.0*density*nrm_v/h); */
2103 /* if (fabs(dmt[eN*nQuadraturePoints_element+k]) > nrm_v/h) */
2104 /* { */
2105 /* tau1[eN*nQuadraturePoints_element+k] = (4.0*viscosity + 2.0*density*nrm_v*h); */
2106 /* } */
2107 /* else */
2108 /* tau1[eN*nQuadraturePoints_element+k] = h*h/tau0[eN*nQuadraturePoints_element+k]; */
2109 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2110 /* 2.0*density*nrm_v/h + */
2111 /* fabs(dmt[eN*nQuadraturePoints_element+k])); */
2112 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
2113 /* 2.0*density*nrm_v*h; */
2114 /* tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h) + */
2115 /* 2.0*density*nrm_v/h); */
2116 /* tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity + */
2117 /* 2.0*density*nrm_v*h+ */
2118 /* fabs(dmt[eN*nQuadraturePoints_element+k])*h*h; */
2119  }
2120  }
2121  /* printf("Element Re = %12.5e \n",Re_max);
2122  printf("Element CFL = %12.5e \n",CFL_max);*/
2123 }
2124 
2125 void calculateSubgridErrorNavierStokes2D_generic_tau(int nElements_global,
2126  int nQuadraturePoints_element,
2127  int nSpace,
2128  double* inverseJ,
2129  double* dmt,
2130  double* dm,
2131  double* f,
2132  double* a,
2133  double* tau0,
2134  double* tau1,
2135  double* cfl)
2136 {
2137  int eN,k,nSpace2=nSpace*nSpace,I,J,K;
2138  double density,viscosity,Re_max=0.0,CFL_max=0.0,nrm_v,G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0*36.0,g_I,g_dot_g,den;
2139  for(eN=0;eN<nElements_global;eN++)
2140  {
2141  for (k=0;k<nQuadraturePoints_element;k++)
2142  {
2143  v_dot_Gv=0.0;
2144  CI_nu2_G_ddot_G=0.0;
2145  g_dot_g=0.0;
2146  for (I=0;I<nSpace;I++)
2147  {
2148  Gv_I=0.0;
2149  g_I=0.0;
2150  for (J=0;J<nSpace;J++)
2151  {
2152  G_IJ=0.0;
2153  for (K=0;K<nSpace;K++)
2154  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2155  k*nSpace2+
2156  K*nSpace+
2157  I]
2158  *
2159  inverseJ[eN*nQuadraturePoints_element*nSpace2+
2160  k*nSpace2+
2161  K*nSpace+
2162  J];
2163  Gv_I += G_IJ
2164  *
2165  f[eN*nQuadraturePoints_element*nSpace+
2166  k*nSpace+
2167  J];
2168  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
2169  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2170  k*nSpace2+
2171  J*nSpace+
2172  I];
2173  }
2174  g_dot_g += g_I*g_I;
2175  v_dot_Gv += f[eN*nQuadraturePoints_element*nSpace+
2176  k*nSpace+
2177  I]
2178  *
2179  Gv_I;
2180 
2181  }
2182  density = dm[eN*nQuadraturePoints_element+
2183  k];
2184  /*mwf add another version for Laplace form?*/
2185  /*cek this should now work for either form*/
2186  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
2187  k*nSpace2+
2188  nSpace+1];
2189  CI_nu2_G_ddot_G*=CI*viscosity*viscosity;
2190  cfl[eN*nQuadraturePoints_element+k] = sqrt(v_dot_Gv);
2191  den = sqrt(4.0*dmt[eN*nQuadraturePoints_element+k]*dmt[eN*nQuadraturePoints_element+k]+
2192  v_dot_Gv+
2193  CI_nu2_G_ddot_G+
2194  1.0e-8);
2195 /* den = sqrt(v_dot_Gv+ */
2196 /* CI_nu2_G_ddot_G+ */
2197 /* 1.0e-8); */
2198  if(den > 1.0e-8)
2199  {
2200  tau0[eN*nQuadraturePoints_element+k] = 1.0/den;
2201  tau1[eN*nQuadraturePoints_element+k] = sqrt(v_dot_Gv+CI_nu2_G_ddot_G)/g_dot_g;
2202  //tau1[eN*nQuadraturePoints_element+k] = 1.0/(tau0[eN*nQuadraturePoints_element+k]*g_dot_g);
2203  }
2204  else
2205  {
2206  tau0[eN*nQuadraturePoints_element+k] = 0.0;
2207  tau1[eN*nQuadraturePoints_element+k] = 0.0;
2208  }
2209 /* printf("%12.5e %12.5e \n",tau0[eN*nQuadraturePoints_element+k],tau1[eN*nQuadraturePoints_element+k]); */
2210  }
2211  }
2212 /* printf("Element Re = %12.5e \n",Re_max); */
2213 /* printf("Element CFL = %12.5e \n",CFL_max); */
2214 }
2215 void calculateSubgridErrorNavierStokes2D_generic_tau_sd(int nElements_global,
2216  int nQuadraturePoints_element,
2217  int nSpace,
2218  double* inverseJ,
2219  double* dmt,
2220  double* dm,
2221  double* f,
2222  double* a,
2223  double* tau0,
2224  double* tau1,
2225  double* cfl)
2226 {
2227  int eN,k,nSpace2=nSpace*nSpace,I,J,K;
2228  double density,viscosity,Re_max=0.0,CFL_max=0.0,nrm_v,G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0*36.0,g_I,g_dot_g,den;
2229  for(eN=0;eN<nElements_global;eN++)
2230  {
2231  for (k=0;k<nQuadraturePoints_element;k++)
2232  {
2233  v_dot_Gv=0.0;
2234  CI_nu2_G_ddot_G=0.0;
2235  g_dot_g=0.0;
2236  for (I=0;I<nSpace;I++)
2237  {
2238  Gv_I=0.0;
2239  g_I=0.0;
2240  for (J=0;J<nSpace;J++)
2241  {
2242  G_IJ=0.0;
2243  for (K=0;K<nSpace;K++)
2244  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2245  k*nSpace2+
2246  K*nSpace+
2247  I]
2248  *
2249  inverseJ[eN*nQuadraturePoints_element*nSpace2+
2250  k*nSpace2+
2251  K*nSpace+
2252  J];
2253  Gv_I += G_IJ
2254  *
2255  f[eN*nQuadraturePoints_element*nSpace+
2256  k*nSpace+
2257  J];
2258  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
2259  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2260  k*nSpace2+
2261  J*nSpace+
2262  I];
2263  }
2264  g_dot_g += g_I*g_I;
2265  v_dot_Gv += f[eN*nQuadraturePoints_element*nSpace+
2266  k*nSpace+
2267  I]
2268  *
2269  Gv_I;
2270 
2271  }
2272  density = dm[eN*nQuadraturePoints_element+
2273  k];
2274  /*mwf add another version for Laplace form?*/
2275  /*cek this should now work for either form*/
2276  viscosity = a[eN*nQuadraturePoints_element*nSpace +
2277  k*nSpace+1];
2278  CI_nu2_G_ddot_G*=CI*viscosity*viscosity;
2279  cfl[eN*nQuadraturePoints_element+k] = sqrt(v_dot_Gv);
2280  den = sqrt(4.0*dmt[eN*nQuadraturePoints_element+k]*dmt[eN*nQuadraturePoints_element+k]+
2281  v_dot_Gv+
2282  CI_nu2_G_ddot_G+
2283  1.0e-8);
2284 /* den = sqrt(v_dot_Gv+ */
2285 /* CI_nu2_G_ddot_G+ */
2286 /* 1.0e-8); */
2287  if(den > 1.0e-8)
2288  {
2289  tau0[eN*nQuadraturePoints_element+k] = 1.0/den;
2290  tau1[eN*nQuadraturePoints_element+k] = sqrt(v_dot_Gv+CI_nu2_G_ddot_G)/g_dot_g;
2291  //tau1[eN*nQuadraturePoints_element+k] = 1.0/(tau0[eN*nQuadraturePoints_element+k]*g_dot_g);
2292  }
2293  else
2294  {
2295  tau0[eN*nQuadraturePoints_element+k] = 0.0;
2296  tau1[eN*nQuadraturePoints_element+k] = 0.0;
2297  }
2298 /* printf("%12.5e %12.5e \n",tau0[eN*nQuadraturePoints_element+k],tau1[eN*nQuadraturePoints_element+k]); */
2299  }
2300  }
2301 /* printf("Element Re = %12.5e \n",Re_max); */
2302 /* printf("Element CFL = %12.5e \n",CFL_max); */
2303 }
2304 void calculateSubgridErrorNavierStokes2D_generic_withBodyForce_tau(int nElements_global,
2305  int nQuadraturePoints_element,
2306  int nSpace,
2307  double* inverseJ,
2308  double* dmt,
2309  double* dm,
2310  double* f,
2311  double* a,
2312  double* dr,
2313  double* tau0,
2314  double* tau1,
2315  double* cfl)
2316 {
2317  int eN,k,nSpace2=nSpace*nSpace,I,J,K;
2318  double h,density,viscosity,Re_max=0.0,CFL_max=0.0,nrm_v,G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0,g_I,g_dot_g,den;
2319  for(eN=0;eN<nElements_global;eN++)
2320  {
2321  for (k=0;k<nQuadraturePoints_element;k++)
2322  {
2323  v_dot_Gv=0.0;
2324  CI_nu2_G_ddot_G=0.0;
2325  g_dot_g=0.0;
2326  for (I=0;I<nSpace;I++)
2327  {
2328  Gv_I=0.0;
2329  g_I=0.0;
2330  for (J=0;J<nSpace;J++)
2331  {
2332  G_IJ=0.0;
2333  for (K=0;K<nSpace;K++)
2334  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2335  k*nSpace2+
2336  K*nSpace+
2337  I]
2338  *
2339  inverseJ[eN*nQuadraturePoints_element*nSpace2+
2340  k*nSpace2+
2341  K*nSpace+
2342  J];
2343  Gv_I += G_IJ
2344  *
2345  f[eN*nQuadraturePoints_element*nSpace+
2346  k*nSpace+
2347  J];
2348  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
2349  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2350  k*nSpace2+
2351  J*nSpace+
2352  I];
2353  }
2354  g_dot_g += g_I*g_I;
2355  v_dot_Gv += f[eN*nQuadraturePoints_element*nSpace+
2356  k*nSpace+
2357  I]
2358  *
2359  Gv_I;
2360 
2361  }
2362  density = dm[eN*nQuadraturePoints_element+
2363  k];
2364  viscosity = 0.5*a[eN*nQuadraturePoints_element*nSpace2 +
2365  k*nSpace2];
2366  CI_nu2_G_ddot_G*=CI*viscosity*viscosity;
2367  nrm_v = sqrt(v_dot_Gv);
2368  Re_max = fmax(nrm_v/(viscosity/density),Re_max);
2369  CFL_max = fmax(nrm_v,CFL_max);
2370  cfl[eN*nQuadraturePoints_element+k] = nrm_v/h;
2371  den = sqrt(4.0*dmt[eN*nQuadraturePoints_element+k]*dmt[eN*nQuadraturePoints_element+k]+
2372  4.0*dr[eN*nQuadraturePoints_element+k]*dr[eN*nQuadraturePoints_element+k]+
2373  v_dot_Gv+
2374  CI_nu2_G_ddot_G+
2375  1.0e-8);
2376  if(den > 1.0e-8)
2377  {
2378  tau0[eN*nQuadraturePoints_element+k] = 1.0/den;
2379  tau1[eN*nQuadraturePoints_element+k] = 1.0/(tau0[eN*nQuadraturePoints_element+k]*g_dot_g);
2380  }
2381  else
2382  {
2383  tau0[eN*nQuadraturePoints_element+k] = 0.0;
2384  tau1[eN*nQuadraturePoints_element+k] = 0.0;
2385  }
2386  /*mwf debug
2387  if (den < 1.0e-4 || 1)
2388  {
2389  printf("NSbodyForce tau eN=%d den = %12.5e dmt=%12.5e dr=%12.5e density=%12.5e viscosity= %12.5e\n",eN,den,
2390  dmt[eN*nQuadraturePoints_element+k],dr[eN*nQuadraturePoints_element+k],density,viscosity);
2391  }
2392  */
2393 /* printf("%12.5e %12.5e \n",tau0[eN*nQuadraturePoints_element+k],tau1[eN*nQuadraturePoints_element+k]); */
2394  }
2395  }
2396 /* printf("Element Re = %12.5e \n",Re_max); */
2397 /* printf("Element CFL = %12.5e \n",CFL_max); */
2398 }
2399 
2400 void calculateSubgridErrorNavierStokes2D_generic_withBodyForce_tau_sd(int nElements_global,
2401  int nQuadraturePoints_element,
2402  int nSpace,
2403  double* inverseJ,
2404  double* dmt,
2405  double* dm,
2406  double* f,
2407  double* a,
2408  double* dr,
2409  double* tau0,
2410  double* tau1,
2411  double* cfl)
2412 {
2413  int eN,k,nSpace2=nSpace*nSpace,I,J,K;
2414  double h,density,viscosity,Re_max=0.0,CFL_max=0.0,nrm_v,G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0,g_I,g_dot_g,den;
2415  for(eN=0;eN<nElements_global;eN++)
2416  {
2417  for (k=0;k<nQuadraturePoints_element;k++)
2418  {
2419  v_dot_Gv=0.0;
2420  CI_nu2_G_ddot_G=0.0;
2421  g_dot_g=0.0;
2422  for (I=0;I<nSpace;I++)
2423  {
2424  Gv_I=0.0;
2425  g_I=0.0;
2426  for (J=0;J<nSpace;J++)
2427  {
2428  G_IJ=0.0;
2429  for (K=0;K<nSpace;K++)
2430  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2431  k*nSpace2+
2432  K*nSpace+
2433  I]
2434  *
2435  inverseJ[eN*nQuadraturePoints_element*nSpace2+
2436  k*nSpace2+
2437  K*nSpace+
2438  J];
2439  Gv_I += G_IJ
2440  *
2441  f[eN*nQuadraturePoints_element*nSpace+
2442  k*nSpace+
2443  J];
2444  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
2445  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
2446  k*nSpace2+
2447  J*nSpace+
2448  I];
2449  }
2450  g_dot_g += g_I*g_I;
2451  v_dot_Gv += f[eN*nQuadraturePoints_element*nSpace+
2452  k*nSpace+
2453  I]
2454  *
2455  Gv_I;
2456 
2457  }
2458  density = dm[eN*nQuadraturePoints_element+
2459  k];
2460  viscosity = 0.5*a[eN*nQuadraturePoints_element*nSpace +
2461  k*nSpace];
2462  CI_nu2_G_ddot_G*=CI*viscosity*viscosity;
2463  nrm_v = sqrt(v_dot_Gv);
2464  Re_max = fmax(nrm_v/(viscosity/density),Re_max);
2465  CFL_max = fmax(nrm_v,CFL_max);
2466  cfl[eN*nQuadraturePoints_element+k] = nrm_v/h;
2467  den = sqrt(4.0*dmt[eN*nQuadraturePoints_element+k]*dmt[eN*nQuadraturePoints_element+k]+
2468  4.0*dr[eN*nQuadraturePoints_element+k]*dr[eN*nQuadraturePoints_element+k]+
2469  v_dot_Gv+
2470  CI_nu2_G_ddot_G+
2471  1.0e-8);
2472  if(den > 1.0e-8)
2473  {
2474  tau0[eN*nQuadraturePoints_element+k] = 1.0/den;
2475  tau1[eN*nQuadraturePoints_element+k] = 1.0/(tau0[eN*nQuadraturePoints_element+k]*g_dot_g);
2476  }
2477  else
2478  {
2479  tau0[eN*nQuadraturePoints_element+k] = 0.0;
2480  tau1[eN*nQuadraturePoints_element+k] = 0.0;
2481  }
2482  /*mwf debug
2483  if (den < 1.0e-4 || 1)
2484  {
2485  printf("NSbodyForce tau eN=%d den = %12.5e dmt=%12.5e dr=%12.5e density=%12.5e viscosity= %12.5e\n",eN,den,
2486  dmt[eN*nQuadraturePoints_element+k],dr[eN*nQuadraturePoints_element+k],density,viscosity);
2487  }
2488  */
2489 /* printf("%12.5e %12.5e \n",tau0[eN*nQuadraturePoints_element+k],tau1[eN*nQuadraturePoints_element+k]); */
2490  }
2491  }
2492 /* printf("Element Re = %12.5e \n",Re_max); */
2493 /* printf("Element CFL = %12.5e \n",CFL_max); */
2494 }
2495 
2499 void calculateSubgridErrorNavierStokes2D_GLS_tauRes(int nElements_global,
2500  int nQuadraturePoints_element,
2501  int nDOF_trial_element,
2502  int nSpace,
2503  double* tau0,
2504  double* tau1,
2505  double* pdeResidualP,
2506  double* dpdeResidualP_du,
2507  double* dpdeResidualP_dv,
2508  double* pdeResidualU,
2509  double* dpdeResidualU_dp,
2510  double* dpdeResidualU_du,
2511  double* dpdeResidualU_dv,
2512  double* pdeResidualV,
2513  double* dpdeResidualV_dp,
2514  double* dpdeResidualV_du,
2515  double* dpdeResidualV_dv,
2516  double* subgridErrorP,
2517  double* dsubgridErrorP_du,
2518  double* dsubgridErrorP_dv,
2519  double* subgridErrorU,
2520  double* dsubgridErrorU_dp,
2521  double* dsubgridErrorU_du,
2522  double* dsubgridErrorU_dv,
2523  double* subgridErrorV,
2524  double* dsubgridErrorV_dp,
2525  double* dsubgridErrorV_du,
2526  double* dsubgridErrorV_dv)
2527 {
2528  int eN,k,j;
2529  for(eN=0;eN<nElements_global;eN++)
2530  {
2531  for (k=0;k<nQuadraturePoints_element;k++)
2532  {
2533  /* GLS momentum */
2534  subgridErrorU[eN*nQuadraturePoints_element+
2535  k] =
2536  -tau0[eN*nQuadraturePoints_element+k]
2537  *
2538  pdeResidualU[eN*nQuadraturePoints_element+
2539  k];
2540  subgridErrorV[eN*nQuadraturePoints_element+
2541  k] =
2542  -tau0[eN*nQuadraturePoints_element+k]
2543  *
2544  pdeResidualV[eN*nQuadraturePoints_element+
2545  k];
2546  /* GLS pressure */
2547  subgridErrorP[eN*nQuadraturePoints_element+
2548  k] =
2549  -tau1[eN*nQuadraturePoints_element+k]
2550  *pdeResidualP[eN*nQuadraturePoints_element+
2551  k];
2552  for (j=0;j<nDOF_trial_element;j++)
2553  {
2554  /* GLS momentum*/
2555  /* u */
2556  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2557  k*nDOF_trial_element+
2558  j] =
2559  -tau0[eN*nQuadraturePoints_element+k]
2560  *
2561  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2562  k*nDOF_trial_element+
2563  j];
2564  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2565  k*nDOF_trial_element+
2566  j] =
2567  -tau0[eN*nQuadraturePoints_element+k]
2568  *
2569  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2570  k*nDOF_trial_element+
2571  j];
2572  dsubgridErrorU_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2573  k*nDOF_trial_element+
2574  j] =
2575  -tau0[eN*nQuadraturePoints_element+k]
2576  *
2577  dpdeResidualU_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2578  k*nDOF_trial_element+
2579  j];
2580  /* v */
2581  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2582  k*nDOF_trial_element+
2583  j] =
2584  -tau0[eN*nQuadraturePoints_element+k]
2585  *
2586  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2587  k*nDOF_trial_element+
2588  j];
2589  dsubgridErrorV_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2590  k*nDOF_trial_element+
2591  j] =
2592  -tau0[eN*nQuadraturePoints_element+k]
2593  *
2594  dpdeResidualV_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2595  k*nDOF_trial_element+
2596  j];
2597  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2598  k*nDOF_trial_element+
2599  j] =
2600  -tau0[eN*nQuadraturePoints_element+k]
2601  *
2602  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2603  k*nDOF_trial_element+
2604  j];
2605  /* GLS pressure */
2606  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2607  k*nDOF_trial_element+
2608  j]
2609  =
2610  -tau1[eN*nQuadraturePoints_element+k]
2611  *
2612  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2613  k*nDOF_trial_element+
2614  j];
2615  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2616  k*nDOF_trial_element+
2617  j]
2618  =
2619  -tau1[eN*nQuadraturePoints_element+k]
2620  *
2621  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2622  k*nDOF_trial_element+
2623  j];
2624  }
2625  }
2626  }
2627 }
2628 
2632 void calculateSubgridErrorStokes_GLS_tau(int nElements_global,
2633  int nQuadraturePoints_element,
2634  int nSpace,
2635  double* elementDiameter,
2636  double* pfac,
2637  double* a,
2638  double* tau0,
2639  double* tau1)
2640 {
2641  int eN,k,nSpace2=nSpace*nSpace;
2642  double h,viscosity,density;//caller should scale if kinematic viscosity
2643  for(eN=0;eN<nElements_global;eN++)
2644  {
2645  h = elementDiameter[eN];
2646  for (k=0;k<nQuadraturePoints_element;k++)
2647  {
2648  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
2649  k*nSpace2];
2650  density = pfac[eN*nQuadraturePoints_element*nSpace+k*nSpace+0];
2651  tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h));
2652  tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity/density;
2653  }
2654  }
2655 }
2656 
2657 void calculateSubgridErrorStokes_GLS_tau_sd(int nElements_global,
2658  int nQuadraturePoints_element,
2659  int nSpace,
2660  double* elementDiameter,
2661  double* pfac,
2662  double* a,
2663  double* tau0,
2664  double* tau1)
2665 {
2666  int eN,k;
2667  double h,viscosity, density;
2668  for(eN=0;eN<nElements_global;eN++)
2669  {
2670  h = elementDiameter[eN];
2671  for (k=0;k<nQuadraturePoints_element;k++)
2672  {
2673  viscosity = a[eN*nQuadraturePoints_element*nSpace*nSpace +
2674  k*nSpace*nSpace];
2675  density = pfac[eN*nQuadraturePoints_element*nSpace+k*nSpace+0];
2676  tau0[eN*nQuadraturePoints_element+k] = 1.0/(4.0*viscosity/(h*h));
2677  tau1[eN*nQuadraturePoints_element+k] = 4.0*viscosity/density;
2678  }
2679  }
2680 }
2681 
2685 void calculateSubgridErrorStokes2D_GLS_tauRes(int nElements_global,
2686  int nQuadraturePoints_element,
2687  int nDOF_trial_element,
2688  int nSpace,
2689  double* tau0,
2690  double* tau1,
2691  double* pdeResidualP,
2692  double* dpdeResidualP_du,
2693  double* dpdeResidualP_dv,
2694  double* pdeResidualU,
2695  double* dpdeResidualU_dp,
2696  double* dpdeResidualU_du,
2697  double* pdeResidualV,
2698  double* dpdeResidualV_dp,
2699  double* dpdeResidualV_dv,
2700  double* subgridErrorP,
2701  double* dsubgridErrorP_du,
2702  double* dsubgridErrorP_dv,
2703  double* subgridErrorU,
2704  double* dsubgridErrorU_dp,
2705  double* dsubgridErrorU_du,
2706  double* subgridErrorV,
2707  double* dsubgridErrorV_dp,
2708  double* dsubgridErrorV_dv)
2709 {
2710  int eN,k,j;
2711  for(eN=0;eN<nElements_global;eN++)
2712  {
2713  for (k=0;k<nQuadraturePoints_element;k++)
2714  {
2715  /* GLS momentum */
2716  subgridErrorU[eN*nQuadraturePoints_element+
2717  k] =
2718  -tau0[eN*nQuadraturePoints_element+k]
2719  *
2720  pdeResidualU[eN*nQuadraturePoints_element+
2721  k];
2722  subgridErrorV[eN*nQuadraturePoints_element+
2723  k] =
2724  -tau0[eN*nQuadraturePoints_element+k]
2725  *
2726  pdeResidualV[eN*nQuadraturePoints_element+
2727  k];
2728  /* GLS pressure */
2729  subgridErrorP[eN*nQuadraturePoints_element+
2730  k] =
2731  -tau1[eN*nQuadraturePoints_element+k]
2732  *pdeResidualP[eN*nQuadraturePoints_element+
2733  k];
2734  for (j=0;j<nDOF_trial_element;j++)
2735  {
2736  /* GLS momentum*/
2737  /* u */
2738  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2739  k*nDOF_trial_element+
2740  j] =
2741  -tau0[eN*nQuadraturePoints_element+k]
2742  *
2743  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2744  k*nDOF_trial_element+
2745  j];
2746  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2747  k*nDOF_trial_element+
2748  j] =
2749  -tau0[eN*nQuadraturePoints_element+k]
2750  *
2751  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2752  k*nDOF_trial_element+
2753  j];
2754  /* v */
2755  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2756  k*nDOF_trial_element+
2757  j] =
2758  -tau0[eN*nQuadraturePoints_element+k]
2759  *
2760  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2761  k*nDOF_trial_element+
2762  j];
2763  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2764  k*nDOF_trial_element+
2765  j] =
2766  -tau0[eN*nQuadraturePoints_element+k]
2767  *
2768  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2769  k*nDOF_trial_element+
2770  j];
2771  /* GLS pressure */
2772  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2773  k*nDOF_trial_element+
2774  j]
2775  =
2776  -tau1[eN*nQuadraturePoints_element+k]
2777  *
2778  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2779  k*nDOF_trial_element+
2780  j];
2781  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2782  k*nDOF_trial_element+
2783  j]
2784  =
2785  -tau1[eN*nQuadraturePoints_element+k]
2786  *
2787  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2788  k*nDOF_trial_element+
2789  j];
2790  }
2791  }
2792  }
2793 }
2797 void calculateSubgridErrorStokes3D_GLS_tauRes(int nElements_global,
2798  int nQuadraturePoints_element,
2799  int nDOF_trial_element,
2800  int nSpace,
2801  double* tau0,
2802  double* tau1,
2803  double* pdeResidualP,
2804  double* dpdeResidualP_du,
2805  double* dpdeResidualP_dv,
2806  double* dpdeResidualP_dw,
2807  double* pdeResidualU,
2808  double* dpdeResidualU_dp,
2809  double* dpdeResidualU_du,
2810  double* pdeResidualV,
2811  double* dpdeResidualV_dp,
2812  double* dpdeResidualV_dv,
2813  double* pdeResidualW,
2814  double* dpdeResidualW_dp,
2815  double* dpdeResidualW_dw,
2816  double* subgridErrorP,
2817  double* dsubgridErrorP_du,
2818  double* dsubgridErrorP_dv,
2819  double* dsubgridErrorP_dw,
2820  double* subgridErrorU,
2821  double* dsubgridErrorU_dp,
2822  double* dsubgridErrorU_du,
2823  double* subgridErrorV,
2824  double* dsubgridErrorV_dp,
2825  double* dsubgridErrorV_dv,
2826  double* subgridErrorW,
2827  double* dsubgridErrorW_dp,
2828  double* dsubgridErrorW_dw)
2829 {
2830  int eN,k,j;
2831  for(eN=0;eN<nElements_global;eN++)
2832  {
2833  for (k=0;k<nQuadraturePoints_element;k++)
2834  {
2835  /* GLS momentum */
2836  subgridErrorU[eN*nQuadraturePoints_element+
2837  k] =
2838  -tau0[eN*nQuadraturePoints_element+k]
2839  *
2840  pdeResidualU[eN*nQuadraturePoints_element+
2841  k];
2842  subgridErrorV[eN*nQuadraturePoints_element+
2843  k] =
2844  -tau0[eN*nQuadraturePoints_element+k]
2845  *
2846  pdeResidualV[eN*nQuadraturePoints_element+
2847  k];
2848  subgridErrorW[eN*nQuadraturePoints_element+
2849  k] =
2850  -tau0[eN*nQuadraturePoints_element+k]
2851  *
2852  pdeResidualW[eN*nQuadraturePoints_element+
2853  k];
2854  /* GLS pressure */
2855  subgridErrorP[eN*nQuadraturePoints_element+
2856  k] =
2857  -tau1[eN*nQuadraturePoints_element+k]
2858  *pdeResidualP[eN*nQuadraturePoints_element+
2859  k];
2860  for (j=0;j<nDOF_trial_element;j++)
2861  {
2862  /* GLS momentum*/
2863  /* u */
2864  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2865  k*nDOF_trial_element+
2866  j] =
2867  -tau0[eN*nQuadraturePoints_element+k]
2868  *
2869  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2870  k*nDOF_trial_element+
2871  j];
2872  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2873  k*nDOF_trial_element+
2874  j] =
2875  -tau0[eN*nQuadraturePoints_element+k]
2876  *
2877  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2878  k*nDOF_trial_element+
2879  j];
2880  /* v */
2881  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2882  k*nDOF_trial_element+
2883  j] =
2884  -tau0[eN*nQuadraturePoints_element+k]
2885  *
2886  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2887  k*nDOF_trial_element+
2888  j];
2889  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2890  k*nDOF_trial_element+
2891  j] =
2892  -tau0[eN*nQuadraturePoints_element+k]
2893  *
2894  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2895  k*nDOF_trial_element+
2896  j];
2897  /* w */
2898  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2899  k*nDOF_trial_element+
2900  j] =
2901  -tau0[eN*nQuadraturePoints_element+k]
2902  *
2903  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
2904  k*nDOF_trial_element+
2905  j];
2906  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
2907  k*nDOF_trial_element+
2908  j] =
2909  -tau0[eN*nQuadraturePoints_element+k]
2910  *
2911  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
2912  k*nDOF_trial_element+
2913  j];
2914  /* GLS pressure */
2915  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2916  k*nDOF_trial_element+
2917  j]
2918  =
2919  -tau1[eN*nQuadraturePoints_element+k]
2920  *
2921  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
2922  k*nDOF_trial_element+
2923  j];
2924  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2925  k*nDOF_trial_element+
2926  j]
2927  =
2928  -tau1[eN*nQuadraturePoints_element+k]
2929  *
2930  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
2931  k*nDOF_trial_element+
2932  j];
2933  dsubgridErrorP_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
2934  k*nDOF_trial_element+
2935  j]
2936  =
2937  -tau1[eN*nQuadraturePoints_element+k]
2938  *
2939  dpdeResidualP_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
2940  k*nDOF_trial_element+
2941  j];
2942  }
2943  }
2944  }
2945 }
2946 
2950 void calculateSubgridErrorNavierStokes3D_GLS_tauRes(int nElements_global,
2951  int nQuadraturePoints_element,
2952  int nDOF_trial_element,
2953  int nSpace,
2954  double* tau0,
2955  double* tau1,
2956  double* pdeResidualP,
2957  double* dpdeResidualP_du,
2958  double* dpdeResidualP_dv,
2959  double* dpdeResidualP_dw,
2960  double* pdeResidualU,
2961  double* dpdeResidualU_dp,
2962  double* dpdeResidualU_du,
2963  double* dpdeResidualU_dv,
2964  double* dpdeResidualU_dw,
2965  double* pdeResidualV,
2966  double* dpdeResidualV_dp,
2967  double* dpdeResidualV_du,
2968  double* dpdeResidualV_dv,
2969  double* dpdeResidualV_dw,
2970  double* pdeResidualW,
2971  double* dpdeResidualW_dp,
2972  double* dpdeResidualW_du,
2973  double* dpdeResidualW_dv,
2974  double* dpdeResidualW_dw,
2975  double* subgridErrorP,
2976  double* dsubgridErrorP_du,
2977  double* dsubgridErrorP_dv,
2978  double* dsubgridErrorP_dw,
2979  double* subgridErrorU,
2980  double* dsubgridErrorU_dp,
2981  double* dsubgridErrorU_du,
2982  double* dsubgridErrorU_dv,
2983  double* dsubgridErrorU_dw,
2984  double* subgridErrorV,
2985  double* dsubgridErrorV_dp,
2986  double* dsubgridErrorV_du,
2987  double* dsubgridErrorV_dv,
2988  double* dsubgridErrorV_dw,
2989  double* subgridErrorW,
2990  double* dsubgridErrorW_dp,
2991  double* dsubgridErrorW_du,
2992  double* dsubgridErrorW_dv,
2993  double* dsubgridErrorW_dw)
2994 {
2995  int eN,k,nSpace2=nSpace*nSpace,j;
2996  for(eN=0;eN<nElements_global;eN++)
2997  {
2998  for (k=0;k<nQuadraturePoints_element;k++)
2999  {
3000  /* GLS momentum */
3001  subgridErrorU[eN*nQuadraturePoints_element+
3002  k] =
3003  -tau0[eN*nQuadraturePoints_element+k]
3004  *
3005  pdeResidualU[eN*nQuadraturePoints_element+
3006  k];
3007  subgridErrorV[eN*nQuadraturePoints_element+
3008  k] =
3009  -tau0[eN*nQuadraturePoints_element+k]
3010  *
3011  pdeResidualV[eN*nQuadraturePoints_element+
3012  k];
3013  subgridErrorW[eN*nQuadraturePoints_element+
3014  k] =
3015  -tau0[eN*nQuadraturePoints_element+k]
3016  *
3017  pdeResidualW[eN*nQuadraturePoints_element+
3018  k];
3019  /* GLS pressure */
3020  subgridErrorP[eN*nQuadraturePoints_element+
3021  k] =
3022  -tau1[eN*nQuadraturePoints_element+k]
3023  *pdeResidualP[eN*nQuadraturePoints_element+
3024  k];
3025  for (j=0;j<nDOF_trial_element;j++)
3026  {
3027  /* GLS momentum*/
3028  /* u */
3029  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3030  k*nDOF_trial_element+
3031  j] =
3032  -tau0[eN*nQuadraturePoints_element+k]
3033  *
3034  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3035  k*nDOF_trial_element+
3036  j];
3037  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3038  k*nDOF_trial_element+
3039  j] =
3040  -tau0[eN*nQuadraturePoints_element+k]
3041  *
3042  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3043  k*nDOF_trial_element+
3044  j];
3045  dsubgridErrorU_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3046  k*nDOF_trial_element+
3047  j] =
3048  -tau0[eN*nQuadraturePoints_element+k]
3049  *
3050  dpdeResidualU_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3051  k*nDOF_trial_element+
3052  j];
3053  dsubgridErrorU_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3054  k*nDOF_trial_element+
3055  j] =
3056  -tau0[eN*nQuadraturePoints_element+k]
3057  *
3058  dpdeResidualU_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3059  k*nDOF_trial_element+
3060  j];
3061  /* v */
3062  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3063  k*nDOF_trial_element+
3064  j] =
3065  -tau0[eN*nQuadraturePoints_element+k]
3066  *
3067  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3068  k*nDOF_trial_element+
3069  j];
3070  dsubgridErrorV_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3071  k*nDOF_trial_element+
3072  j] =
3073  -tau0[eN*nQuadraturePoints_element+k]
3074  *
3075  dpdeResidualV_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3076  k*nDOF_trial_element+
3077  j];
3078  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3079  k*nDOF_trial_element+
3080  j] =
3081  -tau0[eN*nQuadraturePoints_element+k]
3082  *
3083  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3084  k*nDOF_trial_element+
3085  j];
3086  dsubgridErrorV_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3087  k*nDOF_trial_element+
3088  j] =
3089  -tau0[eN*nQuadraturePoints_element+k]
3090  *
3091  dpdeResidualV_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3092  k*nDOF_trial_element+
3093  j];
3094  /* w */
3095  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3096  k*nDOF_trial_element+
3097  j] =
3098  -tau0[eN*nQuadraturePoints_element+k]
3099  *
3100  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3101  k*nDOF_trial_element+
3102  j];
3103  dsubgridErrorW_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3104  k*nDOF_trial_element+
3105  j] =
3106  -tau0[eN*nQuadraturePoints_element+k]
3107  *
3108  dpdeResidualW_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3109  k*nDOF_trial_element+
3110  j];
3111  dsubgridErrorW_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3112  k*nDOF_trial_element+
3113  j] =
3114  -tau0[eN*nQuadraturePoints_element+k]
3115  *
3116  dpdeResidualW_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3117  k*nDOF_trial_element+
3118  j];
3119  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3120  k*nDOF_trial_element+
3121  j] =
3122  -tau0[eN*nQuadraturePoints_element+k]
3123  *
3124  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3125  k*nDOF_trial_element+
3126  j];
3127  /* GLS pressure */
3128  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3129  k*nDOF_trial_element+
3130  j]
3131  =
3132  -tau1[eN*nQuadraturePoints_element+k]
3133  *
3134  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3135  k*nDOF_trial_element+
3136  j];
3137  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3138  k*nDOF_trial_element+
3139  j]
3140  =
3141  -tau1[eN*nQuadraturePoints_element+k]
3142  *
3143  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3144  k*nDOF_trial_element+
3145  j];
3146  dsubgridErrorP_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3147  k*nDOF_trial_element+
3148  j]
3149  =
3150  -tau1[eN*nQuadraturePoints_element+k]
3151  *
3152  dpdeResidualP_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3153  k*nDOF_trial_element+
3154  j];
3155  }
3156  }
3157  }
3158 }
3159 
3163 void calculateSubgridErrorNavierStokes3D_GLS_velocity_pressure(int nElements_global,
3164  int nQuadraturePoints_element,
3165  int nDOF_trial_element,
3166  int nSpace,
3167  double* elementDiameter,
3168  double* dm,
3169  double* f,
3170  double* a,
3171  double* pdeResidualP,
3172  double* dpdeResidualP_du,
3173  double* dpdeResidualP_dv,
3174  double* dpdeResidualP_dw,
3175  double* pdeResidualU,
3176  double* dpdeResidualU_dp,
3177  double* dpdeResidualU_du,
3178  double* pdeResidualV,
3179  double* dpdeResidualV_dp,
3180  double* dpdeResidualV_dv,
3181  double* pdeResidualW,
3182  double* dpdeResidualW_dp,
3183  double* dpdeResidualW_dw,
3184  double* subgridErrorP,
3185  double* dsubgridErrorP_du,
3186  double* dsubgridErrorP_dv,
3187  double* dsubgridErrorP_dw,
3188  double* subgridErrorU,
3189  double* dsubgridErrorU_dp,
3190  double* dsubgridErrorU_du,
3191  double* subgridErrorV,
3192  double* dsubgridErrorV_dp,
3193  double* dsubgridErrorV_dv,
3194  double* subgridErrorW,
3195  double* dsubgridErrorW_dp,
3196  double* dsubgridErrorW_dw)
3197 {
3198  int eN,k,nSpace2=nSpace*nSpace,j;
3199  double h,viscosity,tau1,tau2;
3200  for(eN=0;eN<nElements_global;eN++)
3201  {
3202  h = elementDiameter[eN];
3203  for (k=0;k<nQuadraturePoints_element;k++)
3204  {
3205  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
3206  k*nSpace2];
3207  /* GLS momentum */
3208  tau1 = h*h/(12.0*viscosity);
3209  subgridErrorU[eN*nQuadraturePoints_element+
3210  k] =
3211  tau1
3212  *
3213  pdeResidualU[eN*nQuadraturePoints_element+
3214  k];
3215  subgridErrorV[eN*nQuadraturePoints_element+
3216  k] =
3217  tau1
3218  *
3219  pdeResidualV[eN*nQuadraturePoints_element+
3220  k];
3221  subgridErrorW[eN*nQuadraturePoints_element+
3222  k] =
3223  tau1
3224  *
3225  pdeResidualW[eN*nQuadraturePoints_element+
3226  k];
3227  /* GLS pressure */
3228  tau2 = 6.0*viscosity;
3229  subgridErrorP[eN*nQuadraturePoints_element+
3230  k] =
3231  tau2
3232  *pdeResidualP[eN*nQuadraturePoints_element+
3233  k];
3234  for (j=0;j<nDOF_trial_element;j++)
3235  {
3236  /* GLS momentum*/
3237  /* u */
3238  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3239  k*nDOF_trial_element+
3240  j] =
3241  tau1
3242  *
3243  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3244  k*nDOF_trial_element+
3245  j];
3246  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3247  k*nDOF_trial_element+
3248  j] =
3249  tau1
3250  *
3251  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3252  k*nDOF_trial_element+
3253  j];
3254  /* v */
3255  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3256  k*nDOF_trial_element+
3257  j] =
3258  tau1
3259  *
3260  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3261  k*nDOF_trial_element+
3262  j];
3263  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3264  k*nDOF_trial_element+
3265  j] =
3266  tau1
3267  *
3268  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3269  k*nDOF_trial_element+
3270  j];
3271  /* w */
3272  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3273  k*nDOF_trial_element+
3274  j] =
3275  tau1
3276  *
3277  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3278  k*nDOF_trial_element+
3279  j];
3280  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3281  k*nDOF_trial_element+
3282  j] =
3283  tau1
3284  *
3285  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3286  k*nDOF_trial_element+
3287  j];
3288  /* GLS pressure */
3289  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3290  k*nDOF_trial_element+
3291  j]
3292  =
3293  tau2
3294  *
3295  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3296  k*nDOF_trial_element+
3297  j];
3298  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3299  k*nDOF_trial_element+
3300  j]
3301  =
3302  tau2
3303  *
3304  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3305  k*nDOF_trial_element+
3306  j];
3307  }
3308  }
3309  }
3310 }
3311 
3312 void calculateSubgridErrorNavierStokes3D_GLS_velocity_pressure_sd(int nElements_global,
3313  int nQuadraturePoints_element,
3314  int nDOF_trial_element,
3315  int nSpace,
3316  double* elementDiameter,
3317  double* dm,
3318  double* f,
3319  double* a,
3320  double* pdeResidualP,
3321  double* dpdeResidualP_du,
3322  double* dpdeResidualP_dv,
3323  double* dpdeResidualP_dw,
3324  double* pdeResidualU,
3325  double* dpdeResidualU_dp,
3326  double* dpdeResidualU_du,
3327  double* pdeResidualV,
3328  double* dpdeResidualV_dp,
3329  double* dpdeResidualV_dv,
3330  double* pdeResidualW,
3331  double* dpdeResidualW_dp,
3332  double* dpdeResidualW_dw,
3333  double* subgridErrorP,
3334  double* dsubgridErrorP_du,
3335  double* dsubgridErrorP_dv,
3336  double* dsubgridErrorP_dw,
3337  double* subgridErrorU,
3338  double* dsubgridErrorU_dp,
3339  double* dsubgridErrorU_du,
3340  double* subgridErrorV,
3341  double* dsubgridErrorV_dp,
3342  double* dsubgridErrorV_dv,
3343  double* subgridErrorW,
3344  double* dsubgridErrorW_dp,
3345  double* dsubgridErrorW_dw)
3346 {
3347  int eN,k,j;
3348  double h,viscosity,tau1,tau2;
3349  for(eN=0;eN<nElements_global;eN++)
3350  {
3351  h = elementDiameter[eN];
3352  for (k=0;k<nQuadraturePoints_element;k++)
3353  {
3354  viscosity = a[eN*nQuadraturePoints_element*nSpace +
3355  k*nSpace];
3356  /* GLS momentum */
3357  tau1 = h*h/(12.0*viscosity);
3358  subgridErrorU[eN*nQuadraturePoints_element+
3359  k] =
3360  tau1
3361  *
3362  pdeResidualU[eN*nQuadraturePoints_element+
3363  k];
3364  subgridErrorV[eN*nQuadraturePoints_element+
3365  k] =
3366  tau1
3367  *
3368  pdeResidualV[eN*nQuadraturePoints_element+
3369  k];
3370  subgridErrorW[eN*nQuadraturePoints_element+
3371  k] =
3372  tau1
3373  *
3374  pdeResidualW[eN*nQuadraturePoints_element+
3375  k];
3376  /* GLS pressure */
3377  tau2 = 6.0*viscosity;
3378  subgridErrorP[eN*nQuadraturePoints_element+
3379  k] =
3380  tau2
3381  *pdeResidualP[eN*nQuadraturePoints_element+
3382  k];
3383  for (j=0;j<nDOF_trial_element;j++)
3384  {
3385  /* GLS momentum*/
3386  /* u */
3387  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3388  k*nDOF_trial_element+
3389  j] =
3390  tau1
3391  *
3392  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3393  k*nDOF_trial_element+
3394  j];
3395  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3396  k*nDOF_trial_element+
3397  j] =
3398  tau1
3399  *
3400  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3401  k*nDOF_trial_element+
3402  j];
3403  /* v */
3404  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3405  k*nDOF_trial_element+
3406  j] =
3407  tau1
3408  *
3409  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3410  k*nDOF_trial_element+
3411  j];
3412  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3413  k*nDOF_trial_element+
3414  j] =
3415  tau1
3416  *
3417  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3418  k*nDOF_trial_element+
3419  j];
3420  /* w */
3421  dsubgridErrorW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3422  k*nDOF_trial_element+
3423  j] =
3424  tau1
3425  *
3426  dpdeResidualW_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3427  k*nDOF_trial_element+
3428  j];
3429  dsubgridErrorW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3430  k*nDOF_trial_element+
3431  j] =
3432  tau1
3433  *
3434  dpdeResidualW_dw[eN*nQuadraturePoints_element*nDOF_trial_element+
3435  k*nDOF_trial_element+
3436  j];
3437  /* GLS pressure */
3438  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3439  k*nDOF_trial_element+
3440  j]
3441  =
3442  tau2
3443  *
3444  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3445  k*nDOF_trial_element+
3446  j];
3447  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3448  k*nDOF_trial_element+
3449  j]
3450  =
3451  tau2
3452  *
3453  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3454  k*nDOF_trial_element+
3455  j];
3456  }
3457  }
3458  }
3459 }
3460 
3464 void calculateSubgridErrorStokes2D_1(int nElements_global,
3465  int nQuadraturePoints_element,
3466  int nDOF_trial_element,
3467  int nSpace,
3468  double* elementDiameter,
3469  double* u,
3470  double* v,
3471  double* a,
3472  double* pdeResidualP,
3473  double* dpdeResidualP_du,
3474  double* dpdeResidualP_dv,
3475  double* pdeResidualU,
3476  double* dpdeResidualU_dp,
3477  double* dpdeResidualU_du,
3478  double* pdeResidualV,
3479  double* dpdeResidualV_dp,
3480  double* dpdeResidualV_dv,
3481  double* subgridErrorP,
3482  double* dsubgridErrorP_dp,
3483  double* dsubgridErrorP_du,
3484  double* dsubgridErrorP_dv,
3485  double* subgridErrorU,
3486  double* dsubgridErrorU_dp,
3487  double* dsubgridErrorU_du,
3488  double* dsubgridErrorU_dv,
3489  double* subgridErrorV,
3490  double* dsubgridErrorV_dp,
3491  double* dsubgridErrorV_du,
3492  double* dsubgridErrorV_dv)
3493 {
3494  int eN,k,nSpace2=nSpace*nSpace,j;
3495  double h,velocity_norm,viscosity,tau1,tau2,U,V;
3496  for(eN=0;eN<nElements_global;eN++)
3497  {
3498  h = elementDiameter[eN];
3499  for (k=0;k<nQuadraturePoints_element;k++)
3500  {
3501  U = u[eN*nQuadraturePoints_element+
3502  k];
3503  V = v[eN*nQuadraturePoints_element+
3504  k];
3505  velocity_norm = sqrt(U*U + V*V);
3506  viscosity = a[eN*nQuadraturePoints_element*nSpace2 +
3507  k*nSpace2];
3508  tau1 = 1.0/(4.0*viscosity/(h*h)+2.0*velocity_norm/h+1.0e-8);
3509  tau2 = 4.0*viscosity+2.0*velocity_norm*h;
3510 /* subgridErrorP[eN*nQuadraturePoints_element+ */
3511 /* k] = */
3512 /* tau1 */
3513 /* * */
3514 /* (pdeResidualU[eN*nQuadraturePoints_element+ */
3515 /* k] + */
3516 /* pdeResidualV[eN*nQuadraturePoints_element+ */
3517 /* k]); */
3518 /* subgridErrorU[eN*nQuadraturePoints_element+ */
3519 /* k] = */
3520 /* tau1 */
3521 /* * */
3522 /* pdeResidualU[eN*nQuadraturePoints_element+ */
3523 /* k] */
3524 /* + */
3525 /* tau2 */
3526 /* * */
3527 /* pdeResidualP[eN*nQuadraturePoints_element+ */
3528 /* k]; */
3529 /* subgridErrorV[eN*nQuadraturePoints_element+ */
3530 /* k] = */
3531 /* tau1 */
3532 /* * */
3533 /* pdeResidualV[eN*nQuadraturePoints_element+ */
3534 /* k] */
3535 /* + */
3536 /* tau2 */
3537 /* * */
3538 /* pdeResidualP[eN*nQuadraturePoints_element+ */
3539 /* k]; */
3540  /* GLS momentum */
3541  tau1 = h*h/(12.0*viscosity);
3542  subgridErrorU[eN*nQuadraturePoints_element+
3543  k] =
3544  tau1
3545  *
3546  pdeResidualU[eN*nQuadraturePoints_element+
3547  k];
3548  subgridErrorV[eN*nQuadraturePoints_element+
3549  k] =
3550  tau1
3551  *
3552  pdeResidualV[eN*nQuadraturePoints_element+
3553  k];
3554  /* GLS pressure */
3555  tau2 = 6.0*viscosity;
3556  subgridErrorP[eN*nQuadraturePoints_element+
3557  k] =
3558  tau2
3559  *pdeResidualP[eN*nQuadraturePoints_element+
3560  k];
3561  for (j=0;j<nDOF_trial_element;j++)
3562  {
3563  /* /\* p *\/ */
3564 /* dsubgridErrorP_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3565 /* k*nDOF_trial_element+ */
3566 /* j] */
3567 /* = */
3568 /* tau1 */
3569 /* * */
3570 /* (dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3571 /* k*nDOF_trial_element+ */
3572 /* j] + */
3573 /* dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3574 /* k*nDOF_trial_element+ */
3575 /* j]); */
3576 /* dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3577 /* k*nDOF_trial_element+ */
3578 /* j] */
3579 /* = */
3580 /* tau1 */
3581 /* * */
3582 /* dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3583 /* k*nDOF_trial_element+ */
3584 /* j]; */
3585 /* dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3586 /* k*nDOF_trial_element+ */
3587 /* j] */
3588 /* = */
3589 /* tau1 */
3590 /* * */
3591 /* dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3592 /* k*nDOF_trial_element+ */
3593 /* j]; */
3594 /* /\* u *\/ */
3595 /* dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3596 /* k*nDOF_trial_element+ */
3597 /* j] = */
3598 /* tau1 */
3599 /* * */
3600 /* dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3601 /* k*nDOF_trial_element+ */
3602 /* j]; */
3603 /* dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3604 /* k*nDOF_trial_element+ */
3605 /* j] = */
3606 /* tau1 */
3607 /* * */
3608 /* dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3609 /* k*nDOF_trial_element+ */
3610 /* j] */
3611 /* + */
3612 /* tau2 */
3613 /* * */
3614 /* dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3615 /* k*nDOF_trial_element+ */
3616 /* j]; */
3617 /* dsubgridErrorU_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3618 /* k*nDOF_trial_element+ */
3619 /* j] = */
3620 /* tau2 */
3621 /* * */
3622 /* dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3623 /* k*nDOF_trial_element+ */
3624 /* j]; */
3625 /* /\* v *\/ */
3626 /* dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3627 /* k*nDOF_trial_element+ */
3628 /* j] = */
3629 /* tau1 */
3630 /* * */
3631 /* dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3632 /* k*nDOF_trial_element+ */
3633 /* j]; */
3634 /* dsubgridErrorV_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3635 /* k*nDOF_trial_element+ */
3636 /* j] = */
3637 /* tau2 */
3638 /* * */
3639 /* dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3640 /* k*nDOF_trial_element+ */
3641 /* j]; */
3642 /* dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3643 /* k*nDOF_trial_element+ */
3644 /* j] = */
3645 /* tau1 */
3646 /* * */
3647 /* dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3648 /* k*nDOF_trial_element+ */
3649 /* j] */
3650 /* + */
3651 /* tau2 */
3652 /* * */
3653 /* dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3654 /* k*nDOF_trial_element+ */
3655 /* j]; */
3656  /* GLS momentum*/
3657  /* u */
3658  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3659  k*nDOF_trial_element+
3660  j] =
3661  tau1
3662  *
3663  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3664  k*nDOF_trial_element+
3665  j];
3666  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3667  k*nDOF_trial_element+
3668  j] =
3669  tau1
3670  *
3671  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3672  k*nDOF_trial_element+
3673  j];
3674  /* v */
3675  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3676  k*nDOF_trial_element+
3677  j] =
3678  tau1
3679  *
3680  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3681  k*nDOF_trial_element+
3682  j];
3683  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3684  k*nDOF_trial_element+
3685  j] =
3686  tau1
3687  *
3688  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3689  k*nDOF_trial_element+
3690  j];
3691  /* GLS pressure */
3692  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3693  k*nDOF_trial_element+
3694  j]
3695  =
3696  tau2
3697  *
3698  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3699  k*nDOF_trial_element+
3700  j];
3701  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3702  k*nDOF_trial_element+
3703  j]
3704  =
3705  tau2
3706  *
3707  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3708  k*nDOF_trial_element+
3709  j];
3710  }
3711  }
3712  }
3713 }
3714 
3715 void calculateSubgridErrorStokes2D_1_sd(int nElements_global,
3716  int nQuadraturePoints_element,
3717  int nDOF_trial_element,
3718  int nSpace,
3719  double* elementDiameter,
3720  double* u,
3721  double* v,
3722  double* a,
3723  double* pdeResidualP,
3724  double* dpdeResidualP_du,
3725  double* dpdeResidualP_dv,
3726  double* pdeResidualU,
3727  double* dpdeResidualU_dp,
3728  double* dpdeResidualU_du,
3729  double* pdeResidualV,
3730  double* dpdeResidualV_dp,
3731  double* dpdeResidualV_dv,
3732  double* subgridErrorP,
3733  double* dsubgridErrorP_dp,
3734  double* dsubgridErrorP_du,
3735  double* dsubgridErrorP_dv,
3736  double* subgridErrorU,
3737  double* dsubgridErrorU_dp,
3738  double* dsubgridErrorU_du,
3739  double* dsubgridErrorU_dv,
3740  double* subgridErrorV,
3741  double* dsubgridErrorV_dp,
3742  double* dsubgridErrorV_du,
3743  double* dsubgridErrorV_dv)
3744 {
3745  int eN,k,j;
3746  double h,velocity_norm,viscosity,tau1,tau2,U,V;
3747  for(eN=0;eN<nElements_global;eN++)
3748  {
3749  h = elementDiameter[eN];
3750  for (k=0;k<nQuadraturePoints_element;k++)
3751  {
3752  U = u[eN*nQuadraturePoints_element+
3753  k];
3754  V = v[eN*nQuadraturePoints_element+
3755  k];
3756  velocity_norm = sqrt(U*U + V*V);
3757  viscosity = a[eN*nQuadraturePoints_element*nSpace +
3758  k*nSpace];
3759  tau1 = 1.0/(4.0*viscosity/(h*h)+2.0*velocity_norm/h+1.0e-8);
3760  tau2 = 4.0*viscosity+2.0*velocity_norm*h;
3761 /* subgridErrorP[eN*nQuadraturePoints_element+ */
3762 /* k] = */
3763 /* tau1 */
3764 /* * */
3765 /* (pdeResidualU[eN*nQuadraturePoints_element+ */
3766 /* k] + */
3767 /* pdeResidualV[eN*nQuadraturePoints_element+ */
3768 /* k]); */
3769 /* subgridErrorU[eN*nQuadraturePoints_element+ */
3770 /* k] = */
3771 /* tau1 */
3772 /* * */
3773 /* pdeResidualU[eN*nQuadraturePoints_element+ */
3774 /* k] */
3775 /* + */
3776 /* tau2 */
3777 /* * */
3778 /* pdeResidualP[eN*nQuadraturePoints_element+ */
3779 /* k]; */
3780 /* subgridErrorV[eN*nQuadraturePoints_element+ */
3781 /* k] = */
3782 /* tau1 */
3783 /* * */
3784 /* pdeResidualV[eN*nQuadraturePoints_element+ */
3785 /* k] */
3786 /* + */
3787 /* tau2 */
3788 /* * */
3789 /* pdeResidualP[eN*nQuadraturePoints_element+ */
3790 /* k]; */
3791  /* GLS momentum */
3792  tau1 = h*h/(12.0*viscosity);
3793  subgridErrorU[eN*nQuadraturePoints_element+
3794  k] =
3795  tau1
3796  *
3797  pdeResidualU[eN*nQuadraturePoints_element+
3798  k];
3799  subgridErrorV[eN*nQuadraturePoints_element+
3800  k] =
3801  tau1
3802  *
3803  pdeResidualV[eN*nQuadraturePoints_element+
3804  k];
3805  /* GLS pressure */
3806  tau2 = 6.0*viscosity;
3807  subgridErrorP[eN*nQuadraturePoints_element+
3808  k] =
3809  tau2
3810  *pdeResidualP[eN*nQuadraturePoints_element+
3811  k];
3812  for (j=0;j<nDOF_trial_element;j++)
3813  {
3814  /* /\* p *\/ */
3815 /* dsubgridErrorP_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3816 /* k*nDOF_trial_element+ */
3817 /* j] */
3818 /* = */
3819 /* tau1 */
3820 /* * */
3821 /* (dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3822 /* k*nDOF_trial_element+ */
3823 /* j] + */
3824 /* dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3825 /* k*nDOF_trial_element+ */
3826 /* j]); */
3827 /* dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3828 /* k*nDOF_trial_element+ */
3829 /* j] */
3830 /* = */
3831 /* tau1 */
3832 /* * */
3833 /* dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3834 /* k*nDOF_trial_element+ */
3835 /* j]; */
3836 /* dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3837 /* k*nDOF_trial_element+ */
3838 /* j] */
3839 /* = */
3840 /* tau1 */
3841 /* * */
3842 /* dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3843 /* k*nDOF_trial_element+ */
3844 /* j]; */
3845 /* /\* u *\/ */
3846 /* dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3847 /* k*nDOF_trial_element+ */
3848 /* j] = */
3849 /* tau1 */
3850 /* * */
3851 /* dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3852 /* k*nDOF_trial_element+ */
3853 /* j]; */
3854 /* dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3855 /* k*nDOF_trial_element+ */
3856 /* j] = */
3857 /* tau1 */
3858 /* * */
3859 /* dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3860 /* k*nDOF_trial_element+ */
3861 /* j] */
3862 /* + */
3863 /* tau2 */
3864 /* * */
3865 /* dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3866 /* k*nDOF_trial_element+ */
3867 /* j]; */
3868 /* dsubgridErrorU_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3869 /* k*nDOF_trial_element+ */
3870 /* j] = */
3871 /* tau2 */
3872 /* * */
3873 /* dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3874 /* k*nDOF_trial_element+ */
3875 /* j]; */
3876 /* /\* v *\/ */
3877 /* dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3878 /* k*nDOF_trial_element+ */
3879 /* j] = */
3880 /* tau1 */
3881 /* * */
3882 /* dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3883 /* k*nDOF_trial_element+ */
3884 /* j]; */
3885 /* dsubgridErrorV_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3886 /* k*nDOF_trial_element+ */
3887 /* j] = */
3888 /* tau2 */
3889 /* * */
3890 /* dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3891 /* k*nDOF_trial_element+ */
3892 /* j]; */
3893 /* dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3894 /* k*nDOF_trial_element+ */
3895 /* j] = */
3896 /* tau1 */
3897 /* * */
3898 /* dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3899 /* k*nDOF_trial_element+ */
3900 /* j] */
3901 /* + */
3902 /* tau2 */
3903 /* * */
3904 /* dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+ */
3905 /* k*nDOF_trial_element+ */
3906 /* j]; */
3907  /* GLS momentum*/
3908  /* u */
3909  dsubgridErrorU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3910  k*nDOF_trial_element+
3911  j] =
3912  tau1
3913  *
3914  dpdeResidualU_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3915  k*nDOF_trial_element+
3916  j];
3917  dsubgridErrorU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3918  k*nDOF_trial_element+
3919  j] =
3920  tau1
3921  *
3922  dpdeResidualU_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3923  k*nDOF_trial_element+
3924  j];
3925  /* v */
3926  dsubgridErrorV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3927  k*nDOF_trial_element+
3928  j] =
3929  tau1
3930  *
3931  dpdeResidualV_dp[eN*nQuadraturePoints_element*nDOF_trial_element+
3932  k*nDOF_trial_element+
3933  j];
3934  dsubgridErrorV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3935  k*nDOF_trial_element+
3936  j] =
3937  tau1
3938  *
3939  dpdeResidualV_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3940  k*nDOF_trial_element+
3941  j];
3942  /* GLS pressure */
3943  dsubgridErrorP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3944  k*nDOF_trial_element+
3945  j]
3946  =
3947  tau2
3948  *
3949  dpdeResidualP_du[eN*nQuadraturePoints_element*nDOF_trial_element+
3950  k*nDOF_trial_element+
3951  j];
3952  dsubgridErrorP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3953  k*nDOF_trial_element+
3954  j]
3955  =
3956  tau2
3957  *
3958  dpdeResidualP_dv[eN*nQuadraturePoints_element*nDOF_trial_element+
3959  k*nDOF_trial_element+
3960  j];
3961  }
3962  }
3963  }
3964 }
3965 
3966 void calculateSubgridErrorShallowWater1D(int nElements_global,
3967  int nQuadraturePoints_element,
3968  double g,
3969  double* elementDiameter,
3970  double* h,
3971  double* hu,
3972  double* cfl_1,
3973  double* cfl_2)
3974 {
3975  int eN,k;
3976  double u,c;
3977  for(eN=0;eN<nElements_global;eN++)
3978  {
3979  for (k=0;k<nQuadraturePoints_element;k++)
3980  {
3981  u = hu[eN*nQuadraturePoints_element+k]/fmax(h[eN*nQuadraturePoints_element+k],1.0e-8);
3982  c = sqrt(fabs(g*h[eN*nQuadraturePoints_element+k]));
3983  cfl_1[eN*nQuadraturePoints_element+k] = (u-c)/elementDiameter[eN];
3984  cfl_2[eN*nQuadraturePoints_element+k] = (u+c)/elementDiameter[eN];
3985  }
3986  }
3987 }
3988 void calculateSubgridErrorShallowWater2D(int nElements_global,
3989  int nQuadraturePoints_element,
3990  double g,
3991  double* elementDiameter,
3992  double* h,
3993  double* hu,
3994  double* hv,
3995  double* cfl_1,
3996  double* cfl_2,
3997  double* cfl_3)
3998 {
3999  int eN,k;
4000  double u,v,speed,c;
4001  for(eN=0;eN<nElements_global;eN++)
4002  {
4003  for (k=0;k<nQuadraturePoints_element;k++)
4004  {
4005  u = hu[eN*nQuadraturePoints_element+k]/fmax(h[eN*nQuadraturePoints_element+k],1.0e-8);
4006  v = hv[eN*nQuadraturePoints_element+k]/fmax(h[eN*nQuadraturePoints_element+k],1.0e-8);
4007  speed = sqrt(u*u + v*v);
4008  c = sqrt(fabs(g*h[eN*nQuadraturePoints_element+k]));
4009  cfl_1[eN*nQuadraturePoints_element+k] = (speed-c)/elementDiameter[eN];
4010  cfl_2[eN*nQuadraturePoints_element+k] = speed/elementDiameter[eN];
4011  cfl_3[eN*nQuadraturePoints_element+k] = (speed+c)/elementDiameter[eN];
4012  }
4013  }
4014 }
4015 
4016 /*this should reduce to harari's tau * dt in the linear case*/
4017 void calculateSubgridError_Harari_tau_sd(int nElements_global,
4018  int nQuadraturePoints_element,
4019  int nSpace,
4020  double dt,
4021  int* rowptr,
4022  int* colind,
4023  double* elementDiameter,
4024  double* a,
4025  double* tau)
4026 {
4027  int eN,k,I,m,nnz=rowptr[nSpace];
4028  double h,A_max,A_II,tauMax,heps;
4029  for(eN=0;eN<nElements_global;eN++)
4030  {
4031  h = elementDiameter[eN];
4032  tauMax=0.0;
4033  for (k=0;k<nQuadraturePoints_element;k++)
4034  {
4035  A_max = 0.0;
4036  /*max diagonal entry for A. This choice needs to be looked into*/
4037  for(I=0;I<nSpace;I++)
4038  {
4039  for(m=rowptr[I];m<rowptr[I+1];m++)
4040  {
4041  if(I==colind[m])
4042  {
4043  A_II = a[eN*nQuadraturePoints_element*nnz+
4044  k*nnz+
4045  m];
4046  A_max = (A_II > A_max) ? A_II : A_max;
4047  }
4048  }
4049  }
4050  if (A_max > 0.0)
4051  {
4052  heps = h/(A_max*dt);
4053  heps = (heps > 20.0) ? 20.0 : heps;
4054  }
4055  tau[eN*nQuadraturePoints_element + k]=dt + 6.0*A_max*dt*dt/(h*h)*(1.0-cosh(heps))/(2.0+cosh(heps));
4056  }
4057  }
4058 }
4059 
4060 
4064 void calculateSubgridErrorGradient_tauRes(int nElements_global,
4065  int nQuadraturePoints_element,
4066  int nSpace,
4067  double* tau_gradient,
4068  double* grad_pdeResidual,
4069  double* grad_subgridError)
4070 {
4071  int eN,k,j,I;
4072  for(eN=0;eN<nElements_global;eN++)
4073  for (k=0;k<nQuadraturePoints_element;k++)
4074  {
4075  for (I=0; I < nSpace; I++)
4076  {
4077  grad_subgridError[eN*nQuadraturePoints_element*nSpace+
4078  k*nSpace + I] =
4079  -tau_gradient[eN*nQuadraturePoints_element+
4080  k]
4081  *
4082  grad_pdeResidual[eN*nQuadraturePoints_element*nSpace+
4083  k*nSpace + I];
4084  }
4085 
4086  }
4087 }
4088 
4089 
4090 
4091 void calculateSubgridError_ADR_Sangalli_tau(int nElements_global,
4092  int nQuadraturePoints_element,
4093  int nSpace,
4094  double* inverseJ,
4095  double* dmt,
4096  double* df,
4097  double* a,
4098  double* da,
4099  double* grad_phi,
4100  double* dphi,
4101  double* dr,
4102  double* pe,
4103  double* cfl,
4104  double* tau,
4105  double* tau_gradient)
4106 {
4107  int eN,k,I,J,K,nSpace2=nSpace*nSpace;
4108  double Vlin,Alin,A_II,cfl1,cfl2,Pe,Da,vlin[3],G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0*36.0,g_I,g_dot_g;
4109  double h_e,tau_00,tau_11,t_00,t_11,coshGamma;
4110  for(eN=0;eN<nElements_global;eN++)
4111  {
4112  for (k=0;k<nQuadraturePoints_element;k++)
4113  {
4114  Vlin = 0.0;
4115  Alin = 0.0;
4116  for(I=0;I<nSpace;I++)
4117  {
4118  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
4119  k*nSpace +
4120  I];
4121  for(J=0;J<nSpace;J++)
4122  {
4123  vlin[I]
4124  -=
4125  da[eN*nQuadraturePoints_element*nSpace2 +
4126  k*nSpace2 +
4127  I*nSpace +
4128  J]
4129  *
4130  grad_phi[eN*nQuadraturePoints_element*nSpace +
4131  k*nSpace +
4132  J];
4133  }
4134  }
4135  v_dot_Gv=0.0;
4136  CI_nu2_G_ddot_G=0.0;
4137  g_dot_g=0.0;
4138  for (I=0;I<nSpace;I++)
4139  {
4140  Gv_I=0.0;
4141  g_I=0.0;
4142  for (J=0;J<nSpace;J++)
4143  {
4144  G_IJ=0.0;
4145  for (K=0;K<nSpace;K++)
4146  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
4147  k*nSpace2+
4148  K*nSpace+
4149  I]
4150  *
4151  inverseJ[eN*nQuadraturePoints_element*nSpace2+
4152  k*nSpace2+
4153  K*nSpace+
4154  J];
4155  Gv_I += G_IJ
4156  *
4157  vlin[J];
4158  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
4159  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
4160  k*nSpace2+
4161  J*nSpace+
4162  I];
4163  }
4164  g_dot_g += g_I*g_I;
4165  v_dot_Gv += vlin[I]
4166  *
4167  Gv_I;
4168 
4169  }
4170  /*max diagonal entry for A. This choice needs to be looked into*/
4171  for(I=0;I<nSpace;I++)
4172  {
4173  A_II = a[eN*nQuadraturePoints_element*nSpace2 +
4174  k*nSpace2 +
4175  I*nSpace +
4176  I];
4177  Alin = (A_II > Alin) ? A_II : Alin;
4178  }
4179  Alin*=dphi[eN*nQuadraturePoints_element +
4180  k];
4181  CI_nu2_G_ddot_G*=CI*Alin*Alin;
4182  /*need to check this*/
4183  h_e = sqrt(g_dot_g);
4184  cfl1 = 0.5*sqrt(v_dot_Gv);
4185  cfl2 = 0.25*sqrt(CI_nu2_G_ddot_G);
4186  Pe = sqrt(v_dot_Gv)/sqrt(CI_nu2_G_ddot_G);
4187  Da = (dmt[eN*nQuadraturePoints_element + k] + dr[eN*nQuadraturePoints_element + k])/(cfl1 + 1.0e-12);
4188 
4189 
4190  /*handle different cases here?*/
4191  if (fabs(Da) < 1.0e-8)
4192  {
4193  tau_11 = 0.0;
4194  tau_00 = 1.0/sqrt(v_dot_Gv +
4195  CI_nu2_G_ddot_G +
4196  +1.0e-8);
4197  }
4198  else if (Da > Pe*10.0)
4199  {
4200  tau_00 = 0.5/dmt[eN*nQuadraturePoints_element + k];
4201  tau_11 = h_e*h_e/(12.0*dmt[eN*nQuadraturePoints_element + k]);
4202  }
4203  else if (cfl1 < 1.0e-8)
4204  {
4205  if (Pe*(2.0*Da + Pe) < 0.0)
4206  coshGamma = cos(sqrt(Pe*(-2.0*Da - Pe)));
4207  else if (Pe*(2.0*Da + Pe) > 1000.0)
4208  coshGamma = cosh(100.0);
4209  else
4210  coshGamma = cosh(sqrt(Pe*(2.0*Da + Pe)));
4211 
4212  tau_00 = -1.0/dmt[eN*nQuadraturePoints_element + k]*1.0/(-2.0 - Pe*Da/(coshGamma - 1.0 - 2.0*Pe*Da) + 1.0e-8);
4213  tau_11 = -h_e*h_e/(dmt[eN*nQuadraturePoints_element + k]*6.0)*(-1.0 - 6.0/(2.0*Pe*Da) + 3.0*coshGamma -Pe*Da/(-2.0 + 2.0*coshGamma - Pe*Da));
4214  }
4215  else
4216  {
4217  t_00 = fmin(Pe/6.0,fmin(0.5/Da,0.5));
4218  t_11 = fmin(Pe/60.0,fmin(-1.0/(12.0*Da),1.0/24.0));
4219  t_11*= -1.0;
4220  tau_00 = t_00/cfl1;
4221  tau_11 = t_11/cfl1/(h_e*h_e);
4222  }
4223  assert(tau_11 <= 0.0);
4224  assert(tau_00 >= 0.0);
4225 
4226  tau[eN*nQuadraturePoints_element + k] = tau_00;
4227  tau_gradient[eN*nQuadraturePoints_element + k] = tau_11;
4228 
4229  if (cfl1 > cfl2)
4230  cfl[eN*nQuadraturePoints_element +
4231  k] = cfl1;
4232  else
4233  cfl[eN*nQuadraturePoints_element +
4234  k] = cfl2;
4235 
4236  pe[eN*nQuadraturePoints_element +
4237  k] = sqrt(v_dot_Gv)/sqrt(CI_nu2_G_ddot_G);
4238 
4239 
4240  }
4241  }
4242 }
4243 /*note computes cfl as always advection cfl*/
4244 void calculateSubgridError_ADR_Sangalli_tau_sd(int nElements_global,
4245  int nQuadraturePoints_element,
4246  int nSpace,
4247  int* rowptr,
4248  int* colind,
4249  double* inverseJ,
4250  double* dmt,
4251  double* df,
4252  double* a,
4253  double* da,
4254  double* grad_phi,
4255  double* dphi,
4256  double* dr,
4257  double* pe,
4258  double* cfl,
4259  double* tau,
4260  double* tau_gradient)
4261 {
4262  int eN,k,I,J,K,m,nnz=rowptr[nSpace],nSpace2=nSpace*nSpace;
4263  double Vlin,Alin,A_II,cfl1,cfl2,vlin[3],G_IJ,v_dot_Gv,CI_nu2_G_ddot_G,Gv_I,CI=36.0*36.0,g_I,g_dot_g;
4264  double Da,h_e,tau_00,tau_11,t_00,t_11,coshGamma,Pe,sinhAlpha,coshAlpha,PeEval,gamma2,gam2Eval,DaEval;
4265  for(eN=0;eN<nElements_global;eN++)
4266  {
4267  for (k=0;k<nQuadraturePoints_element;k++)
4268  {
4269  Vlin = 0.0;
4270  Alin = 0.0;
4271  for(I=0;I<nSpace;I++)
4272  {
4273  vlin[I] = df[eN*nQuadraturePoints_element*nSpace +
4274  k*nSpace +
4275  I];
4276  for(m=rowptr[I];m<rowptr[I+1];m++)
4277  {
4278  vlin[I]
4279  -=
4280  da[eN*nQuadraturePoints_element*nnz+
4281  k*nnz+
4282  m]
4283  *
4284  grad_phi[eN*nQuadraturePoints_element*nSpace +
4285  k*nSpace +
4286  colind[m]];
4287  }
4288  }
4289  v_dot_Gv=0.0;
4290  CI_nu2_G_ddot_G=0.0;
4291  g_dot_g=0.0;
4292  for (I=0;I<nSpace;I++)
4293  {
4294  Gv_I=0.0;
4295  g_I=0.0;
4296  for (J=0;J<nSpace;J++)
4297  {
4298  G_IJ=0.0;
4299  for (K=0;K<nSpace;K++)
4300  G_IJ +=inverseJ[eN*nQuadraturePoints_element*nSpace2+
4301  k*nSpace2+
4302  K*nSpace+
4303  I]
4304  *
4305  inverseJ[eN*nQuadraturePoints_element*nSpace2+
4306  k*nSpace2+
4307  K*nSpace+
4308  J];
4309  Gv_I += G_IJ
4310  *
4311  vlin[J];
4312  CI_nu2_G_ddot_G+=G_IJ*G_IJ;
4313  g_I+=inverseJ[eN*nQuadraturePoints_element*nSpace2+
4314  k*nSpace2+
4315  J*nSpace+
4316  I];
4317  }
4318  g_dot_g += g_I*g_I;
4319  v_dot_Gv += vlin[I]
4320  *
4321  Gv_I;
4322 
4323  }
4324  /*max diagonal entry for A. This choice needs to be looked into*/
4325  for(I=0;I<nSpace;I++)
4326  {
4327  for(m=rowptr[I];m<rowptr[I+1];m++)
4328  {
4329  if (I==colind[m])
4330  {
4331  A_II = a[eN*nQuadraturePoints_element*nnz +
4332  k*nnz + m];
4333  Alin = (A_II > Alin) ? A_II : Alin;
4334  }
4335  }
4336  }
4337  Alin*=dphi[eN*nQuadraturePoints_element +
4338  k];
4339  CI_nu2_G_ddot_G*=CI*Alin*Alin;
4340  /*need to check this*/
4341  h_e = 1.0/sqrt(g_dot_g);
4342  cfl1 = sqrt(v_dot_Gv);/*mwf not sure about 0.5 here 0.5*sqrt(v_dot_Gv);*/
4343  cfl2 = 0.25*sqrt(CI_nu2_G_ddot_G);
4344  Pe = sqrt(v_dot_Gv)/(sqrt(CI_nu2_G_ddot_G)+1.0e-12);
4345  Da = (dmt[eN*nQuadraturePoints_element + k] + dr[eN*nQuadraturePoints_element + k])/(cfl1 + 1.0e-12);
4346 
4347 
4348  /*handle different cases here?*/
4349  if (fabs(Da) < 1.0e-8)
4350  {
4351  tau_11 = 0.0;
4352  tau_00 = 1.0/sqrt(v_dot_Gv +
4353  CI_nu2_G_ddot_G +
4354  +1.0e-8);
4355  /*mwf debug
4356  printf("Sangalli Da ~ 0 case Da= %12.5e Pe=%12.5e cfl=%12.5e tau_11=%12.5e tau_00=%12.5e \n",Da,Pe,cfl1,tau_11,tau_00);
4357  */
4358  }
4359  else if (Da > Pe*10.0)
4360  {
4361  tau_00 = 0.5/(dmt[eN*nQuadraturePoints_element + k]+dr[eN*nQuadraturePoints_element + k]);
4362  tau_11 = -h_e*h_e/(12.0*(dmt[eN*nQuadraturePoints_element + k] + dr[eN*nQuadraturePoints_element + k]));
4363  /*mwf debug
4364  printf("Sangalli Da >> 1 case h_e= %12.5e dmt=%12.5e dr=%12.5e Da= %12.5e Pe=%12.5e cfl=%12.5e tau_11=%12.5e tau_00=%12.5e \n",h_e,
4365  dmt[eN*nQuadraturePoints_element + k],dr[eN*nQuadraturePoints_element + k],
4366  Da,Pe,cfl1,tau_11,tau_00);
4367  */
4368  }
4369  else if (cfl1 < 1.0e-8)
4370  {
4371  if (Pe*(2.0*Da + Pe) < 0.0)
4372  coshGamma = cos(sqrt(Pe*(-2.0*Da - Pe)));
4373  else if (Pe*(2.0*Da + Pe) > 1000.0)
4374  coshGamma = cosh(100.0);
4375  else
4376  coshGamma = cosh(sqrt(Pe*(2.0*Da + Pe)));
4377 
4378  tau_00 = -1.0/dmt[eN*nQuadraturePoints_element + k]*1.0/(-2.0 - Pe*Da/(coshGamma - 1.0 - 2.0*Pe*Da) + 1.0e-8);
4379  tau_11 = -h_e*h_e/(dmt[eN*nQuadraturePoints_element + k]*6.0)*(-1.0 - 6.0/(2.0*Pe*Da) + 3.0*coshGamma -Pe*Da/(-2.0 + 2.0*coshGamma - Pe*Da));
4380  /*mwf debug
4381  printf("Sangalli cfl ~ 0 case Da= %12.5e Pe=%12.5e cfl=%12.5e tau_11=%12.5e tau_00=%12.5e \n",Da,Pe,cfl1,tau_11,tau_00);
4382  */
4383  }
4384  else if (Da >= 0.0 && 0)
4385  {
4386  t_00 = fmin(Pe/6.0,fmin(0.5/Da,0.5));
4387  t_11 = fmin(Pe/60.0,fmin(1.0/(12.0*Da),1.0/24.0));
4388  t_11*= -1.0;
4389  tau_00 = t_00/cfl1;
4390  tau_11 = h_e*h_e*t_11/cfl1;
4391  /*mwf debug */
4392  printf("Sangalli general reduced formula Da= %12.5e Pe=%12.5e cfl=%12.5e t_11=%12.5e tau_11=%12.5e t_00=%12.5e tau_00=%12.5e \n",Da,Pe,cfl1,t_11,tau_11,t_00,tau_00);
4393 
4394  }
4395  else
4396  {
4397  PeEval = fmin(Pe,10.0);
4398  DaEval = -Da;
4399  if (DaEval < -10.0)
4400  DaEval = -10.0;
4401  if (DaEval > 10.0)
4402  DaEval = 10.0;
4403  gamma2 = PeEval*(-2.0*DaEval + PeEval);
4404 
4405  if (gamma2 < 0.0)
4406  coshGamma = cos(sqrt(PeEval*(2.0*DaEval - PeEval)));
4407  else
4408  coshGamma = cosh(sqrt(gamma2));
4409  sinhAlpha = sinh(PeEval);
4410  coshAlpha = cosh(PeEval);
4411 
4412  t_00 = 1.0/(-2.0*DaEval + DaEval*DaEval*sinhAlpha/(-coshAlpha + coshGamma + DaEval*sinhAlpha));
4413  t_11 = (-3.0 -DaEval*DaEval + 3.0*DaEval/(PeEval+1.0e-12) + DaEval*(3.0*DaEval * coshGamma + sinhAlpha*(-3.0 + DaEval*DaEval))/
4414  (-2.0*coshAlpha + 2.0*coshGamma + DaEval*sinhAlpha))/(6.0*DaEval*DaEval*DaEval + 1.0e-12);
4415  /*go ahead and switch sign to match Sangalli convention?*/
4416 /* if (Da > 0.0) */
4417 /* DaEval = fmin(Da,10.0); */
4418 /* if (Da < 0.0) */
4419 /* DaEval = fmax(Da,-10.0); */
4420 /* gamma2 = PeEval*(2.0*DaEval + PeEval); */
4421 /* gam2Eval= gamma2; */
4422 /* if (gamma2 < 0.0) */
4423 /* coshGamma = cos(sqrt(PeEval*(-2.0*DaEval - PeEval))); */
4424 /* else */
4425 /* coshGamma = cosh(sqrt(gam2Eval)); */
4426 /* sinhAlpha = sinh(PeEval); */
4427 /* coshAlpha = cosh(PeEval); */
4428 
4429 /* t_00 = 1.0/(2.0*DaEval + DaEval*DaEval*sinhAlpha/(-coshAlpha + coshGamma - DaEval*sinhAlpha + 1.0e-12) + 1.0e-12); */
4430 /* t_11 = -1.0/(6.0*DaEval*DaEval*DaEval+1.0e-12)* */
4431 /* (-3.0 - DaEval*DaEval -3.0*DaEval/(PeEval+1.0e-12) */
4432 /* -DaEval*(-3.0*DaEval*coshGamma +sinhAlpha*(-3.0+DaEval*DaEval))/(-2.0*coshAlpha + 2.0*coshGamma -DaEval*sinhAlpha)); */
4433 
4434 /* t_11 = (DaEval*(sinhAlpha*(DaEval*DaEval - 3.) - 3.*DaEval*coshGamma)/(2.*coshGamma - sinhAlpha*DaEval - 2*coshAlpha) + DaEval*DaEval + 3.*DaEval/PeEval + 3.)/(6.*DaEval*DaEval*DaEval); */
4435 
4436 
4437 
4438 
4439  tau_00 = t_00/cfl1;
4440  tau_11 = h_e*h_e*t_11/cfl1;
4441  /*mwf debug*/
4442  printf("\t Sangalli general h_e= %12.5e Da= %12.5e DaEval= %12.5e Pe=%12.5e PeEval= %12.5e cfl=%12.5e \n\t gamma2= %12.5e coshGamma=%12.5e sinhAlpha=%12.5e coshAlpha=%12.5e \n\t t_11=%12.5e tau_11=%12.5e t_00= %12.5e tau_00=%12.5e \n",h_e,Da,DaEval,Pe,PeEval,cfl1,gamma2,coshGamma,sinhAlpha,coshAlpha,t_11,tau_11,t_00,tau_00);
4443 
4444  }
4445 
4446  /*mwf try to catch eval errors*/
4447  if (tau_11 > 0.0)
4448  {
4449  /*mwf debug*/
4450  printf("Sangalli tau_11 = %12.5e > 0 switching ",tau_11);
4451  t_11 = fmin(Pe/60.0,fmin(1.0/(12.0*Da),1.0/24.0));
4452  t_11*= -1.0;
4453  tau_11 = h_e*h_e*t_11/cfl1;
4454  /*mwf debug*/
4455  printf("tau_11 = %12.5e \n",tau_11);
4456 
4457 
4458  }
4459  assert(tau_11 <= 0.0);
4460  assert(tau_00 >= 0.0);
4461 
4462  tau[eN*nQuadraturePoints_element + k] = tau_00;
4463  tau_gradient[eN*nQuadraturePoints_element + k] = tau_11;
4464 
4465  if (cfl1 > cfl2)
4466  cfl[eN*nQuadraturePoints_element +
4467  k] = cfl1;
4468  else
4469  cfl[eN*nQuadraturePoints_element +
4470  k] = cfl2;
4471 
4472  pe[eN*nQuadraturePoints_element +
4473  k] = sqrt(v_dot_Gv)/sqrt(CI_nu2_G_ddot_G);
4474 
4475  }
4476  }
4477 }
4478 
calculateSubgridError_A_tau
void calculateSubgridError_A_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, char stabilization, double *elementDiameter, double *dmt, double *df, double *cfl, double *tau)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation.
Definition: subgridError.c:971
calculateSubgridError_ADR_Sangalli_tau
void calculateSubgridError_ADR_Sangalli_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, double *inverseJ, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau, double *tau_gradient)
Definition: subgridError.c:4091
calculateSubgridError_ADR_tau_1
void calculateSubgridError_ADR_tau_1(int nElements_global, int nQuadraturePoints_element, int nSpace, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:224
calculateSubgridErrorStokes2D_GLS_tauRes
void calculateSubgridErrorStokes2D_GLS_tauRes(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *tau0, double *tau1, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:2685
calculateSubgridError_tauRes
void calculateSubgridError_tauRes(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, double *tau, double *pdeResidual, double *dpdeResidual, double *subgridError, double *dsubgridError)
Calculate the ASGS subgrid error given tau and the strong residual.
Definition: subgridError.c:12
calculateSubgridErrorStokes2D_GLS_velocity_pressure_sd
void calculateSubgridErrorStokes2D_GLS_velocity_pressure_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv)
Definition: subgridError.c:1498
calculateSubgridErrorNavierStokes3D_GLS_tauRes
void calculateSubgridErrorNavierStokes3D_GLS_tauRes(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *tau0, double *tau1, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *dpdeResidualP_dw, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *dpdeResidualU_dv, double *dpdeResidualU_dw, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_du, double *dpdeResidualV_dv, double *dpdeResidualV_dw, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_du, double *dpdeResidualW_dv, double *dpdeResidualW_dw, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *dsubgridErrorP_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *dsubgridErrorU_dv, double *dsubgridErrorU_dw, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_du, double *dsubgridErrorV_dv, double *dsubgridErrorV_dw, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_du, double *dsubgridErrorW_dv, double *dsubgridErrorW_dw)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:2950
f
Double f
Definition: Headers.h:64
calculateSubgridError_ADR_tau_1_sd
void calculateSubgridError_ADR_tau_1_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, int *rowptr, int *colind, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Definition: subgridError.c:304
calculateSubgridErrorNavierStokes2D_generic_withBodyForce_tau
void calculateSubgridErrorNavierStokes2D_generic_withBodyForce_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, double *inverseJ, double *dmt, double *dm, double *f, double *a, double *dr, double *tau0, double *tau1, double *cfl)
Definition: subgridError.c:2304
calculateSubgridErrorStokes_GLS_tau
void calculateSubgridErrorStokes_GLS_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, double *elementDiameter, double *pfac, double *a, double *tau0, double *tau1)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:2632
calculateSubgridErrorNavierStokes3D_GLS_velocity_pressure
void calculateSubgridErrorNavierStokes3D_GLS_velocity_pressure(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *dm, double *f, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *dpdeResidualP_dw, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_dw, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *dsubgridErrorP_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_dw)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:3163
calculateSubgridError_ADR_tau_p_sd
void calculateSubgridError_ADR_tau_p_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, int *rowptr, int *colind, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Definition: subgridError.c:132
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
subgridError.h
C implementations of subgrid error calculations.
calculateSubgridErrorStokes3D_GLS_velocity_pressure
void calculateSubgridErrorStokes3D_GLS_velocity_pressure(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *dpdeResidualP_dw, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_dw, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *dsubgridErrorP_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_dw)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:1617
calculateSubgridErrorStokes2D_GLS_velocity_sd
void calculateSubgridErrorStokes2D_GLS_velocity_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv)
Definition: subgridError.c:1072
calculateSubgridErrorGradient_tauRes
void calculateSubgridErrorGradient_tauRes(int nElements_global, int nQuadraturePoints_element, int nSpace, double *tau_gradient, double *grad_pdeResidual, double *grad_subgridError)
Calculate the ASGS subgrid error given tau and the strong residual.
Definition: subgridError.c:4064
calculateSubgridErrorStokes3D_GLS_tauRes
void calculateSubgridErrorStokes3D_GLS_tauRes(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *tau0, double *tau1, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *dpdeResidualP_dw, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_dw, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *dsubgridErrorP_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_dw)
3D version of Stokes GLS tau
Definition: subgridError.c:2797
num
Int num
Definition: Headers.h:32
calculateSubgridErrorStokes3D_GLS_velocity_pressure_sd
void calculateSubgridErrorStokes3D_GLS_velocity_pressure_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *dpdeResidualP_dw, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_dw, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *dsubgridErrorP_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_dw)
Definition: subgridError.c:1763
calculateSubgridErrorStokes_GLS_tau_sd
void calculateSubgridErrorStokes_GLS_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, double *elementDiameter, double *pfac, double *a, double *tau0, double *tau1)
Definition: subgridError.c:2657
calculateSubgridError_ADR_tau_2_sd
void calculateSubgridError_ADR_tau_2_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, int *rowptr, int *colind, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Definition: subgridError.c:484
calculateSubgridErrorStokes3D_GLS_velocity
void calculateSubgridErrorStokes3D_GLS_velocity(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_dw)
Calculate the subgrid error for velocity in 3D Stokes equation with a GLS-like formula.
Definition: subgridError.c:1157
U
Double U
Definition: Headers.h:88
v
Double v
Definition: Headers.h:95
calculateSubgridErrorShallowWater1D
void calculateSubgridErrorShallowWater1D(int nElements_global, int nQuadraturePoints_element, double g, double *elementDiameter, double *h, double *hu, double *cfl_1, double *cfl_2)
Definition: subgridError.c:3966
calculateSubgridError_ADR_tau_2
void calculateSubgridError_ADR_tau_2(int nElements_global, int nQuadraturePoints_element, int nSpace, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:394
calculateSubgridErrorStokes2D_1_sd
void calculateSubgridErrorStokes2D_1_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *u, double *v, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *subgridErrorP, double *dsubgridErrorP_dp, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *dsubgridErrorU_dv, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_du, double *dsubgridErrorV_dv)
Definition: subgridError.c:3715
pe
Double pe
Definition: Headers.h:75
calculateSubgridErrorNavierStokes2D_generic_withBodyForce_tau_sd
void calculateSubgridErrorNavierStokes2D_generic_withBodyForce_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, double *inverseJ, double *dmt, double *dm, double *f, double *a, double *dr, double *tau0, double *tau1, double *cfl)
Definition: subgridError.c:2400
u
Double u
Definition: Headers.h:89
c
Double c
Definition: Headers.h:54
calculateSubgridErrorStokes2D_1
void calculateSubgridErrorStokes2D_1(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *u, double *v, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *subgridErrorP, double *dsubgridErrorP_dp, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *dsubgridErrorU_dv, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_du, double *dsubgridErrorV_dv)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:3464
calculateSubgridErrorShallowWater2D
void calculateSubgridErrorShallowWater2D(int nElements_global, int nQuadraturePoints_element, double g, double *elementDiameter, double *h, double *hu, double *hv, double *cfl_1, double *cfl_2, double *cfl_3)
Definition: subgridError.c:3988
calculateSubgridError_ADR_tau_p
void calculateSubgridError_ADR_tau_p(int nElements_global, int nQuadraturePoints_element, int nSpace, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:49
calculateSubgridErrorStokes2D_GLS_velocity
void calculateSubgridErrorStokes2D_GLS_velocity(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv)
Calculate the subgrid error for velocity in 2D Stokes equation with a GLS-like formula.
Definition: subgridError.c:990
calculateSubgridErrorNavierStokes2D_generic_tau_sd
void calculateSubgridErrorNavierStokes2D_generic_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, double *inverseJ, double *dmt, double *dm, double *f, double *a, double *tau0, double *tau1, double *cfl)
Definition: subgridError.c:2215
calculateSubgridErrorNavierStokes3D_GLS_velocity_pressure_sd
void calculateSubgridErrorNavierStokes3D_GLS_velocity_pressure_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *dm, double *f, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *dpdeResidualP_dw, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_dw, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *dsubgridErrorP_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_dw)
Definition: subgridError.c:3312
calculateSubgridError_ADR_generic_tau
void calculateSubgridError_ADR_generic_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, double *inverseJ, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Definition: subgridError.c:580
calculateSubgridErrorNavierStokes2D_GLS_tau_sd
void calculateSubgridErrorNavierStokes2D_GLS_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, double hFactor, double *elementDiameter, double *dmt, double *dm, double *f, double *a, double *tau0, double *tau1, double *cfl)
Definition: subgridError.c:2020
calculateSubgridErrorNavierStokes2D_GLS_tau
void calculateSubgridErrorNavierStokes2D_GLS_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, double hFactor, double *elementDiameter, double *dmt, double *dm, double *f, double *a, double *tau0, double *tau1, double *cfl)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:1914
calculateSubgridError_ADR_tau_sd
void calculateSubgridError_ADR_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, int *rowptr, int *colind, char stabilization, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Definition: subgridError.c:851
calculateSubgridErrorStokes3D_GLS_velocity_sd
void calculateSubgridErrorStokes3D_GLS_velocity_sd(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *pdeResidualW, double *dpdeResidualW_dp, double *dpdeResidualW_dw, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv, double *subgridErrorW, double *dsubgridErrorW_dp, double *dsubgridErrorW_dw)
Definition: subgridError.c:1268
calculateSubgridError_A_tau_2
void calculateSubgridError_A_tau_2(int nElements_global, int nQuadraturePoints_element, int nSpace, double *elementDiameter, double *dmt, double *df, double *cfl, double *tau)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:923
calculateSubgridError_ADR_generic_tau_sd
void calculateSubgridError_ADR_generic_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, int *rowptr, int *colind, double *inverseJ, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Definition: subgridError.c:699
calculateSubgridErrorStokes2D_GLS_velocity_pressure
void calculateSubgridErrorStokes2D_GLS_velocity_pressure(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *elementDiameter, double *a, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_dv, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_dv)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:1382
calculateSubgridError_Harari_tau_sd
void calculateSubgridError_Harari_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, double dt, int *rowptr, int *colind, double *elementDiameter, double *a, double *tau)
Definition: subgridError.c:4017
calculateSubgridErrorNavierStokes2D_GLS_tauRes
void calculateSubgridErrorNavierStokes2D_GLS_tauRes(int nElements_global, int nQuadraturePoints_element, int nDOF_trial_element, int nSpace, double *tau0, double *tau1, double *pdeResidualP, double *dpdeResidualP_du, double *dpdeResidualP_dv, double *pdeResidualU, double *dpdeResidualU_dp, double *dpdeResidualU_du, double *dpdeResidualU_dv, double *pdeResidualV, double *dpdeResidualV_dp, double *dpdeResidualV_du, double *dpdeResidualV_dv, double *subgridErrorP, double *dsubgridErrorP_du, double *dsubgridErrorP_dv, double *subgridErrorU, double *dsubgridErrorU_dp, double *dsubgridErrorU_du, double *dsubgridErrorU_dv, double *subgridErrorV, double *dsubgridErrorV_dp, double *dsubgridErrorV_du, double *dsubgridErrorV_dv)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:2499
calculateSubgridError_ADR_tau
void calculateSubgridError_ADR_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, char stabilization, double *elementDiameter, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation.
Definition: subgridError.c:826
calculateSubgridError_ADR_Sangalli_tau_sd
void calculateSubgridError_ADR_Sangalli_tau_sd(int nElements_global, int nQuadraturePoints_element, int nSpace, int *rowptr, int *colind, double *inverseJ, double *dmt, double *df, double *a, double *da, double *grad_phi, double *dphi, double *dr, double *pe, double *cfl, double *tau, double *tau_gradient)
Definition: subgridError.c:4244
calculateSubgridErrorNavierStokes2D_generic_tau
void calculateSubgridErrorNavierStokes2D_generic_tau(int nElements_global, int nQuadraturePoints_element, int nSpace, double *inverseJ, double *dmt, double *dm, double *f, double *a, double *tau0, double *tau1, double *cfl)
Definition: subgridError.c:2125
calculateSubgridError_A_tau_1
void calculateSubgridError_A_tau_1(int nElements_global, int nQuadraturePoints_element, int nSpace, double *elementDiameter, double *dmt, double *df, double *cfl, double *tau)
Calculate the stabilization parameter for the scalar advection-diffusion-reaction equation using the ...
Definition: subgridError.c:880
nnz
#define nnz
Definition: Richards.h:9
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20