proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
subgridError.c
Go to the documentation of this file.
1#include "subgridError.h"
2
8
12void 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
49void 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
132void 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
224void 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}
304void 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
394void 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
484void 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
580void 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
699void 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
826void 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
851void 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
880void 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
923void 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
971void 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
990void 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
1072void 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
1157void 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
1268void 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
1382void 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
1498void 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
1617void 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}
1763void 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
1914void 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
2020void 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
2125void 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}
2215void 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}
2304void 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
2400void 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
2499void 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
2632void 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
2657void 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
2685void 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}
2794
2797void 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
2950void 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
3163void 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
3312void 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
3464void 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
3715void 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
3966void 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}
3988void 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*/
4017void 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
4064void 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
4091void 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*/
4244void 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
4480
U
Double U
Definition Headers.h:88
f
Double f
Definition Headers.h:64
u
Double u
Definition Headers.h:89
num
Int num
Definition Headers.h:32
v
Double v
Definition Headers.h:95
pe
Double pe
Definition Headers.h:75
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
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
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
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
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_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
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
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_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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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_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
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
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
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
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_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
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
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
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
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
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
c
#define c(i)
Definition jf.h:21
nnz
#define nnz
Definition m_comp_co2.h:19
subgridError.h
C implementations of subgrid error calculations.
Generated on for proteus by doxygen 1.18.0