proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
SW2D.h
Go to the documentation of this file.
1#ifndef SW2D_H
2#define SW2D_H
3#include <cmath>
4#include <iostream>
5#include "CompKernel.h"
6#include "ModelFactory.h"
7#include "ArgumentsDict.h"
8#include "xtensor-python/pyarray.hpp"
9
10namespace py = pybind11;
11
12//cek todo
13//2. Get stabilization right
14//3. Add Riemann solvers for external flux
15//4. Add Riemann solvers for internal flux and DG terms
16//5. Try other choices of variables h,hu,hv, Bova-Carey symmetrization?
17namespace proteus
18{
19 class SW2D_base
20 {
21 public:
22 virtual ~SW2D_base(){}
23 virtual void calculateResidual(arguments_dict& args)=0;
24 virtual void calculateJacobian(arguments_dict& args)=0;
25 virtual void calculateResidual_supg(arguments_dict& args)=0;
26 virtual void calculateJacobian_supg(arguments_dict& args)=0;
27 };
28
29 template<class CompKernelType,
30 int nSpace,
31 int nQuadraturePoints_element,
32 int nDOF_mesh_trial_element,
33 int nDOF_trial_element,
34 int nDOF_test_element,
35 int nQuadraturePoints_elementBoundary>
36 class SW2D : public SW2D_base
37 {
38 public:
39 const int nDOF_test_X_trial_element;
40 CompKernelType ck;
41 SW2D():
42 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
43 ck()
44 {
45 std::cout<<"Constructing SW2D<CompKernelTemplate<"
46 <<nSpace<<","
47 <<nQuadraturePoints_element<<","
48 <<nDOF_mesh_trial_element<<","
49 <<nDOF_trial_element<<","
50 <<nDOF_test_element<<","
51 <<nQuadraturePoints_elementBoundary<<">());"
52 <<std::endl<<std::flush;
53 }
54
55 inline
56 void evaluateCoefficients(const double nu,
57 const double g,
58 const double grad_b[nSpace],
59 const double& h,
60 const double& u,
61 const double& v,
62 double& mass_acc,
63 double& dmass_acc_h,
64 double& mom_u_acc,
65 double& dmom_u_acc_h,
66 double& dmom_u_acc_u,
67 double& mom_v_acc,
68 double& dmom_v_acc_h,
69 double& dmom_v_acc_v,
70 double mass_adv[nSpace],
71 double dmass_adv_h[nSpace],
72 double dmass_adv_u[nSpace],
73 double dmass_adv_v[nSpace],
74 double mom_u_adv[nSpace],
75 double dmom_u_adv_h[nSpace],
76 double dmom_u_adv_u[nSpace],
77 double dmom_u_adv_v[nSpace],
78 double mom_v_adv[nSpace],
79 double dmom_v_adv_h[nSpace],
80 double dmom_v_adv_u[nSpace],
81 double dmom_v_adv_v[nSpace],
82 double mom_u_diff_ten[nSpace],
83 double mom_v_diff_ten[nSpace],
84 double mom_uv_diff_ten[1],
85 double mom_vu_diff_ten[1],
86 double& mom_u_source,
87 double& dmom_u_source_h,
88 double& mom_v_source,
89 double& dmom_v_source_h)
90 {
91 //mass accumulation
92 mass_acc = h;
93 dmass_acc_h = 1.0;
94
95 //u momentum accumulation
96 mom_u_acc=h*u;
97 dmom_u_acc_h=u;
98 dmom_u_acc_u=h;
99
100 //v momentum accumulation
101 mom_v_acc=h*v;
102 dmom_v_acc_h=v;
103 dmom_v_acc_v=h;
104
105 //mass advective flux
106 mass_adv[0]=h*u;
107 mass_adv[1]=h*v;
108
109 dmass_adv_h[0]=u;
110 dmass_adv_h[1]=v;
111
112 dmass_adv_u[0]=h;
113 dmass_adv_u[1]=0.0;
114
115 dmass_adv_v[0]=0.0;
116 dmass_adv_v[1]=h;
117
118 //u momentum advective flux
119 mom_u_adv[0]=h*u*u + 0.5*g*h*h;
120 mom_u_adv[1]=h*u*v;
121
122 dmom_u_adv_h[0]=u*u + g*h;
123 dmom_u_adv_h[1]=u*v;
124
125 dmom_u_adv_u[0]=h*2.0*u;
126 dmom_u_adv_u[1]=h*v;
127
128 dmom_u_adv_v[0]=0.0;
129 dmom_u_adv_v[1]=h*u;
130
131 //v momentum advective_flux
132 mom_v_adv[0]=h*v*u;
133 mom_v_adv[1]=h*v*v + 0.5*g*h*h;
134
135 dmom_v_adv_h[0]=v*u;
136 dmom_v_adv_h[1]=v*v + g*h;
137
138 dmom_v_adv_u[0]=h*v;
139 dmom_v_adv_u[1]=0.0;
140
141 dmom_v_adv_v[0]=h*u;
142 dmom_v_adv_v[1]=h*2.0*v;
143
144 //u momentum diffusion tensor
145 mom_u_diff_ten[0] = 2.0*nu;
146 mom_u_diff_ten[1] = nu;
147
148 mom_uv_diff_ten[0]=nu;
149
150 //v momentum diffusion tensor
151 mom_v_diff_ten[0] = nu;
152 mom_v_diff_ten[1] = 2.0*nu;
153
154 mom_vu_diff_ten[0]=nu;
155
156 //momentum sources
157 mom_u_source = g*h*grad_b[0];
158 dmom_u_source_h = g*grad_b[0];
159
160 mom_v_source = g*h*grad_b[1];
161 dmom_v_source_h = g*grad_b[1];
162 }
163
164 inline
165 void calculateSubgridError_tau(const double& elementDiameter,
166 const double& nu,
167 const double& g,
168 const double& h,
169 const double& u,
170 const double& v,
171 double tau[9],
172 double& cfl)
173 {
174 //todo go back through rx,ry, etc. and the formula for tau in physical variables once this is known.
175 double rx[9],ry[9],rxInv[9],ryInv[9],c=sqrt(fmax(1.0e-8,g*h)),lambdax[3],lambday[3],tauxHat[3],tauyHat[3],taux[9],tauy[9],tauc[9],cflx,cfly,L[9],hStar=fmax(1.0e-8,h);
176 //eigenvalues and eigenvectors for conservation variables h,hu,hv
177 lambdax[0] = u - c;
178 lambdax[1] = u;
179 lambdax[2] = u + c;
180 if (u > 0.0)
181 cflx = (u+c)/elementDiameter;
182 else
183 cflx = fabs(u-c)/elementDiameter;
184
185 double rn=sqrt(1.0+(u-c)*(u-c) + v*v);
186 rx[0*3+0] = 1.0/rn;
187 rx[1*3+0] = (u - c)/rn;
188 rx[2*3+0] = v/rn;
189
190 rx[0*3+1] = 0.0;
191 rx[1*3+1] = 0.0;
192 rx[2*3+1] = 1.0;
193
194 rn = sqrt(1.0 + (u+c)*(u+c) + v*v);
195 rx[0*3+2] = 1.0/rn;
196 rx[1*3+2] = (u + c)/rn;
197 rx[2*3+2] = v/rn;
198
199 rxInv[0*3+0] = 1.0 - (c - u)/(2.0*c);
200 rxInv[1*3+0] = -v;
201 rxInv[2*3+0] = (c-u)/(2.0*c);
202
203 rxInv[0*3+1] = -1.0/(2.0*c);
204 rxInv[1*3+1] = 0.0;
205 rxInv[2*3+1] = 1.0/(2.0*c);
206
207 rxInv[0*3+2] = 0.0;
208 rxInv[1*3+2] = 1.0;
209 rxInv[2*3+2] = 0.0;
210
211 /* //cek hack assume no viscosity for now */
212 /* num = 0.5*fabs(lambdax[0])*elementDiameter + 1.0e-8; */
213 /* den = nu + num*1.0e-8; */
214 /* pe = num/den; */
215 /* tauxHat[0] = 0.5*h*(1.0/tanh(pe) - 1.0/pe)/(Vlin+1.0e-8); */
216
217 tauxHat[0] = 0.5*elementDiameter/(fabs(lambdax[0])+1.0e-8);
218 tauxHat[1] = 0.5*elementDiameter/(fabs(lambdax[1])+1.0e-8);
219 tauxHat[2] = 0.5*elementDiameter/(fabs(lambdax[2])+1.0e-8);
220
221 lambday[0] = v - c;
222 lambday[1] = v;
223 lambday[2] = v + c;
224
225 rn=sqrt(1.0 + u*u+(v-c)*(v-c));
226 ry[0*3+0] = 1.0/rn;
227 ry[1*3+0] = u/rn;
228 ry[2*3+0] = (v - c)/rn;
229
230 ry[0*3+1] = 0.0;
231 ry[1*3+1] = -1.0;
232 ry[2*3+1] = 0.0;
233
234 rn = sqrt(1.0 + u*u + (v+c)*(v+c));
235 ry[0*3+2] = 1.0/rn;
236 ry[1*3+2] = u/rn;
237 ry[2*3+2] = (v + c)/rn;
238
239 ryInv[0*3+0] = 1.0 - (c - v)/(2*c);
240 ryInv[1*3+0] = u;
241 ryInv[2*3+0] = (c-v)/(2*c);
242
243 ryInv[0*3+1] = 0.0;
244 ryInv[1*3+1] = -1.0;
245 ryInv[2*3+1] = 0.0;
246
247 ryInv[0*3+2] = -1.0/(2*c);
248 ryInv[1*3+2] = 0.0;
249 ryInv[2*3+2] = 1.0/(2*c);
250
251 //cek hack assume no viscosity for now
252 tauyHat[0] = 0.5*elementDiameter/(fabs(lambday[0])+1.0e-8);
253 tauyHat[1] = 0.5*elementDiameter/(fabs(lambday[1])+1.0e-8);
254 tauyHat[2] = 0.5*elementDiameter/(fabs(lambday[2])+1.0e-8);
255 if (v > 0.0)
256 cfly = (v+c)/elementDiameter;
257 else
258 cfly = fabs(v-c)/elementDiameter;
259 cfl = sqrt(cflx*cflx+cfly*cfly);//hack, conservative estimate
260
261 //project back to solution variables
262 //initialize, hack do with memset
263 double tmpx[9],tmpy[9];
264 for (int i=0;i<3;i++)
265 for (int j=0;j<3;j++)
266 {
267 taux[i*3+j] = 0.0;
268 tauy[i*3+j] = 0.0;
269 tauc[i*3+j] = 0.0;
270 tau[i*3+j] = 0.0;
271 tmpx[i*3+j] = 0.0;
272 tmpy[i*3+j] = 0.0;
273 L[i*3+j] = 0.0;
274 /* double Ix=0,Iy=0.0; */
275 /* for (int k=0;k<3;k++) */
276 /* { */
277 /* Ix += rx[i*3+k]*rxInv[k*3+j]; */
278 /* Iy += ry[i*3+k]*ryInv[k*3+j]; */
279 /* } */
280 /* std::cout<<i<<'\t'<<j<<'\t'<<Ix<<'\t'<<Iy<<std::endl; */
281 }
282 //transform from characteristic variables to conservation variables
283 for (int i=0;i<3;i++)
284 for (int j=0;j<3;j++)
285 for (int m=0;m<3;m++)
286 {
287 //taux[i*3+j] += rx[i*3+m]*tauxHat[m]*rxInv[m*3+j];
288 //tauy[i*3+j] += ry[i*3+m]*tauyHat[m]*ryInv[m*3+j];
289 if (m==j)
290 {
291 tmpx[i*3+m] += rx[i*3+m]*tauxHat[m];
292 tmpy[i*3+m] += ry[i*3+m]*tauyHat[m];
293 }
294 }
295 for (int i=0;i<3;i++)
296 for (int j=0;j<3;j++)
297 for (int m=0;m<3;m++)
298 {
299 //taux[i*3+j] += tmpx[i*3+m]*rxInv[m*3 + j];
300 //tauy[i*3+j] += tmpy[i*3+m]*ryInv[m*3 + j];
301 taux[i*3+j] += tmpx[i*3+m]*rx[j*3 + m];
302 tauy[i*3+j] += tmpy[i*3+m]*ry[j*3 + m];
303 }
304 //matrix norm
305 for (int i=0;i<3;i++)
306 for (int j=0;j<3;j++)
307 tauc[i*3+j] = sqrt(taux[i*3+j]*taux[i*3+j] + tauy[i*3+j]*tauy[i*3+j]);
308 //transform from conservation variables to solution variables h,u,v
309 L[0*3+0] = 1.0;
310 L[1*3+0] = -u/hStar;
311 L[1*3+1] = 1.0/hStar;
312 L[2*3+0] = -v/hStar;
313 L[2*3+2] = 1.0/hStar;
314 for (int i=0;i<3;i++)
315 for (int j=0;j<3;j++)
316 {
317 for (int k=0;k<3;k++)
318 tau[i*3+j] += L[i*3+k]*tauc[k*3+j];
319 //cek hack, turn it off until I get the ASGS stuff straigtened out
320 tau[i*3+j] = 0.0;
321 }
322 /* std::cout<<"tau"<<std::endl; */
323 /* for (int i=0;i<3;i++) */
324 /* { */
325 /* for (int j=0;j<3;j++) */
326 /* { */
327 /* std::cout<<tau[i*3+j]<<'\t'; */
328 /* } */
329 /* std::cout<<std::endl; */
330 /* } */
331 }
332 inline
333 void calculateSubgridError_tau_supg(const double& elementDiameter,
334 const double& nu,
335 const double& g,
336 const double& h,
337 const double& u,
338 const double& v,
339 const double& alpha,
340 const double& area,
341 double tau_x[9],
342 double tau_y[9],
343 double& cfl)
344 {
345 double c=sqrt(fmax(1.0e-8,g*h)),lambdax[3],lambday[3],cflx,cfly,hStar=fmax(1.0e-8,h),dx,dy,a,ainv;
346 int i=0;
347 //eigenvalues for x and y, using Berger and Stockstill notation where possible
348 //tau_x = \alpha \Delta x \hat{A}, tau_y = \alpha \Delta y \hat{B}
349 lambdax[0] = u + c;
350 lambdax[1] = u - c;
351 lambdax[2] = u;
352 if (u > 0.0)
353 cflx = (u+c)/elementDiameter;
354 else
355 cflx = fabs(u-c)/elementDiameter;
356 lambday[0] = v + c;
357 lambday[1] = v - c;
358 lambday[2] = v;
359 if (v > 0.0)
360 cfly = (v+c)/elementDiameter;
361 else
362 cfly = fabs(v-c)/elementDiameter;
363
364 cfl = sqrt(cflx*cflx+cfly*cfly);
365
366 dx = sqrt(fabs(area)); dy = sqrt(fabs(area));
367 a = sqrt(u*u + v*v + c*c);
368 ainv = 1.0/(a+1.0e-8);
369 //\hat{A}\alpha\Delta x
370 tau_x[0*3+0]= u; tau_x[0*3+1]=h; tau_x[0*3+2]=0.;
371 tau_x[1*3+0]= h*g; tau_x[1*3+1]=h*u; tau_x[1*3+2]=0.;
372 tau_x[2*3+0]= 0.; tau_x[2*3+1]=0.; tau_x[2*3+2]=u*h;
373 for (i=0; i < 9; i++)
374 tau_x[i] *= alpha*dx*ainv;
375 //\hat{B}\alpha\Delta y
376 tau_y[0*3+0]= v; tau_y[0*3+1]=0; tau_y[0*3+2]=h;
377 tau_y[1*3+0]= 0.; tau_y[1*3+1]=h*v; tau_y[1*3+2]=0.;
378 tau_y[2*3+0]= h*g; tau_y[2*3+1]=0.; tau_y[2*3+2]=v*h;
379 for (i=0; i < 9; i++)
380 tau_y[i] *= alpha*dy*ainv;
381
382 }
383
384 inline
385 void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_h,
386 const int& isDOFBoundary_u,
387 const int& isDOFBoundary_v,
388 const int& isFluxBoundary_h,
389 const int& isFluxBoundary_u,
390 const int& isFluxBoundary_v,
391 const double n[nSpace],
392 const double& bc_h,
393 const double bc_f_mass[nSpace],
394 const double bc_f_umom[nSpace],
395 const double bc_f_vmom[nSpace],
396 const double& bc_flux_mass,
397 const double& bc_flux_umom,
398 const double& bc_flux_vmom,
399 const double& h,
400 const double f_mass[nSpace],
401 const double f_umom[nSpace],
402 const double f_vmom[nSpace],
403 const double df_mass_dh[nSpace],
404 const double df_mass_du[nSpace],
405 const double df_mass_dv[nSpace],
406 const double df_umom_dh[nSpace],
407 const double df_umom_du[nSpace],
408 const double df_umom_dv[nSpace],
409 const double df_vmom_dh[nSpace],
410 const double df_vmom_du[nSpace],
411 const double df_vmom_dv[nSpace],
412 double& flux_mass,
413 double& flux_umom,
414 double& flux_vmom,
415 double* velocity)
416 {
417 //cek todo, need to do the Riemann solve
418 /* double flowDirection; */
419 flux_mass = 0.0;
420 flux_umom = 0.0;
421 flux_vmom = 0.0;
422 /* flowDirection=n[0]*f_mass[0]+n[1]*f_mass[1]; */
423 /* if (isDOFBoundary_u != 1) */
424 /* { */
425 /* flux_mass += n[0]*f_mass[0]; */
426 /* velocity[0] = f_mass[0]; */
427 /* if (flowDirection >= 0.0) */
428 /* { */
429 /* flux_umom += n[0]*f_umom[0]; */
430 /* flux_vmom += n[0]*f_vmom[0]; */
431 /* } */
432 /* } */
433 /* else */
434 /* { */
435 /* flux_mass += n[0]*bc_f_mass[0]; */
436 /* velocity[0] = bc_f_mass[0]; */
437 /* //cek still upwind the advection for Dirichlet? */
438 /* if (flowDirection >= 0.0) */
439 /* { */
440 /* flux_umom += n[0]*f_umom[0]; */
441 /* flux_vmom += n[0]*f_vmom[0]; */
442 /* } */
443 /* else */
444 /* { */
445 /* flux_umom+=n[0]*bc_f_umom[0]; */
446 /* flux_vmom+=n[0]*bc_f_vmom[0]; */
447 /* } */
448 /* } */
449 /* if (isDOFBoundary_v != 1) */
450 /* { */
451 /* flux_mass+=n[1]*f_mass[1]; */
452 /* velocity[1] = f_mass[1]; */
453 /* if (flowDirection >= 0.0) */
454 /* { */
455 /* flux_umom+=n[1]*f_umom[1]; */
456 /* flux_vmom+=n[1]*f_vmom[1]; */
457 /* } */
458 /* } */
459 /* else */
460 /* { */
461 /* flux_mass+=n[1]*bc_f_mass[1]; */
462 /* velocity[1] = bc_f_mass[1]; */
463 /* //cek still upwind the advection for Dirichlet? */
464 /* if (flowDirection >= 0.0) */
465 /* { */
466 /* flux_umom+=n[1]*f_umom[1]; */
467 /* flux_vmom+=n[1]*f_vmom[1]; */
468 /* } */
469 /* else */
470 /* { */
471 /* flux_umom+=n[1]*bc_f_umom[1]; */
472 /* flux_vmom+=n[1]*bc_f_vmom[1]; */
473 /* } */
474 /* } */
475 /* else */
476 /* { */
477 /* flux_mass +=n[2]*bc_f_mass[2]; */
478 /* velocity[2] = bc_f_mass[2]; */
479 /* //cek still upwind the advection for Dirichlet? */
480 /* if (flowDirection >= 0.0) */
481 /* { */
482 /* flux_umom+=n[2]*f_umom[2]; */
483 /* flux_vmom+=n[2]*f_vmom[2]; */
484 /* } */
485 /* else */
486 /* { */
487 /* flux_umom+=n[2]*bc_f_umom[2]; */
488 /* flux_vmom+=n[2]*bc_f_vmom[2]; */
489 /* } */
490 /* } */
491 /* if (isDOFBoundary_h == 1) */
492 /* { */
493 /* flux_umom+= n[0]*(bc_h-p)*oneByRho; */
494 /* flux_vmom+= n[1]*(bc_h-p)*oneByRho; */
495 /* } */
496 if (isFluxBoundary_h == 1)
497 {
498 //cek todo, not sure if we'll need this for SW2
499 //velocity[0] += (bc_flux_mass - flux_mass)*n[0];
500 //velocity[1] += (bc_flux_mass - flux_mass)*n[1];
501 //velocity[2] += (bc_flux_mass - flux_mass)*n[2];
502 flux_mass = bc_flux_mass;
503 }
504 if (isFluxBoundary_u == 1)
505 {
506 flux_umom = bc_flux_umom;
507 }
508 if (isFluxBoundary_v == 1)
509 {
510 flux_vmom = bc_flux_vmom;
511 }
512 }
513
514 inline
515 void exteriorNumericalAdvectiveFluxDerivatives(const int& isDOFBoundary_h,
516 const int& isDOFBoundary_u,
517 const int& isDOFBoundary_v,
518 const int& isFluxBoundary_h,
519 const int& isFluxBoundary_u,
520 const int& isFluxBoundary_v,
521 const double n[nSpace],
522 const double& bc_h,
523 const double bc_f_mass[nSpace],
524 const double bc_f_umom[nSpace],
525 const double bc_f_vmom[nSpace],
526 const double& bc_flux_mass,
527 const double& bc_flux_umom,
528 const double& bc_flux_vmom,
529 const double& h,
530 const double f_mass[nSpace],
531 const double f_umom[nSpace],
532 const double f_vmom[nSpace],
533 const double df_mass_du[nSpace],
534 const double df_mass_dv[nSpace],
535 const double df_umom_dh[nSpace],
536 const double df_umom_du[nSpace],
537 const double df_umom_dv[nSpace],
538 const double df_vmom_dh[nSpace],
539 const double df_vmom_du[nSpace],
540 const double df_vmom_dv[nSpace],
541 double& dflux_mass_dh,
542 double& dflux_mass_du,
543 double& dflux_mass_dv,
544 double& dflux_umom_dh,
545 double& dflux_umom_du,
546 double& dflux_umom_dv,
547 double& dflux_vmom_dh,
548 double& dflux_vmom_du,
549 double& dflux_vmom_dv)
550 {
551 double flowDirection;
552 dflux_mass_dh = 0.0;
553 dflux_mass_du = 0.0;
554 dflux_mass_dv = 0.0;
555
556 dflux_umom_dh = 0.0;
557 dflux_umom_du = 0.0;
558 dflux_umom_dv = 0.0;
559
560 dflux_vmom_dh = 0.0;
561 dflux_vmom_du = 0.0;
562 dflux_vmom_dv = 0.0;
563
564 flowDirection=n[0]*f_mass[0]+n[1]*f_mass[1];
565 /* if (isDOFBoundary_u != 1) */
566 /* { */
567 /* dflux_mass_du += n[0]*df_mass_du[0]; */
568 /* if (flowDirection >= 0.0) */
569 /* { */
570 /* dflux_umom_du += n[0]*df_umom_du[0]; */
571 /* dflux_vmom_du += n[0]*df_vmom_du[0]; */
572 /* dflux_vmom_dv += n[0]*df_vmom_dv[0]; */
573 /* } */
574 /* } */
575 /* else */
576 /* { */
577 /* //cek still upwind the advection for Dirichlet? */
578 /* if (flowDirection >= 0.0) */
579 /* { */
580 /* dflux_umom_du += n[0]*df_umom_du[0]; */
581 /* dflux_vmom_du += n[0]*df_vmom_du[0]; */
582 /* dflux_vmom_dv += n[0]*df_vmom_dv[0]; */
583 /* } */
584 /* else */
585 /* { */
586 /* if (isDOFBoundary_v != 1) */
587 /* dflux_vmom_dv += n[0]*df_vmom_dv[0]; */
588 /* } */
589 /* } */
590 /* if (isDOFBoundary_v != 1) */
591 /* { */
592 /* dflux_mass_dv += n[1]*df_mass_dv[1]; */
593 /* if (flowDirection >= 0.0) */
594 /* { */
595 /* dflux_umom_du += n[1]*df_umom_du[1]; */
596 /* dflux_umom_dv += n[1]*df_umom_dv[1]; */
597 /* dflux_vmom_dv += n[1]*df_vmom_dv[1]; */
598 /* } */
599 /* } */
600 /* else */
601 /* { */
602 /* //cek still upwind the advection for Dirichlet? */
603 /* if (flowDirection >= 0.0) */
604 /* { */
605 /* dflux_umom_du += n[1]*df_umom_du[1]; */
606 /* dflux_umom_dv += n[1]*df_umom_dv[1]; */
607 /* dflux_vmom_dv += n[1]*df_vmom_dv[1]; */
608 /* } */
609 /* else */
610 /* { */
611 /* if (isDOFBoundary_u != 1) */
612 /* dflux_umom_du += n[1]*df_umom_du[1]; */
613 /* } */
614 /* } */
615 /* else */
616 /* { */
617 /* //cek still upwind the advection for Dirichlet? */
618 /* if (flowDirection >= 0.0) */
619 /* { */
620 /* dflux_umom_du += n[2]*df_umom_du[2]; */
621 /* dflux_umom_dw += n[2]*df_umom_dw[2]; */
622 /* dflux_vmom_dv += n[2]*df_vmom_dv[2]; */
623 /* } */
624 /* else */
625 /* { */
626 /* if (isDOFBoundary_u != 1) */
627 /* dflux_umom_du += n[2]*df_umom_du[2]; */
628 /* if (isDOFBoundary_v != 1) */
629 /* dflux_vmom_dv += n[2]*df_vmom_dv[2]; */
630 /* } */
631 /* } */
632 /* if (isDOFBoundary_h == 1) */
633 /* { */
634 /* dflux_umom_dp= -n[0]*oneByRho; */
635 /* dflux_vmom_dp= -n[1]*oneByRho; */
636 /* } */
637 if (isFluxBoundary_h == 1)
638 {
639 dflux_mass_dh = 0.0;
640 dflux_mass_du = 0.0;
641 dflux_mass_dv = 0.0;
642 }
643 if (isFluxBoundary_u == 1)
644 {
645 dflux_umom_dh = 0.0;
646 dflux_umom_du = 0.0;
647 dflux_umom_dv = 0.0;
648 }
649 if (isFluxBoundary_v == 1)
650 {
651 dflux_vmom_dh = 0.0;
652 dflux_vmom_du = 0.0;
653 dflux_vmom_dv = 0.0;
654 }
655 }
656
657 /* inline */
658 /* void exteriorNumericalDiffusiveFlux(const double& eps, */
659 /* int* rowptr, */
660 /* int* colind, */
661 /* const int& isDOFBoundary, */
662 /* const int& isFluxBoundary, */
663 /* const double n[nSpace], */
664 /* double* bc_a, */
665 /* const double& bc_u, */
666 /* const double& bc_flux, */
667 /* double* a, */
668 /* const double grad_phi[nSpace], */
669 /* const double& u, */
670 /* const double& penalty, */
671 /* double& flux) */
672 /* { */
673 /* double diffusiveVelocityComponent_I,penaltyFlux,max_a; */
674 /* if(isDOFBoundary == 1) */
675 /* { */
676 /* flux = 0.0; */
677 /* max_a=0.0; */
678 /* for(int I=0;I<nSpace;I++) */
679 /* { */
680 /* diffusiveVelocityComponent_I=0.0; */
681 /* for(int m=rowptr[I];m<rowptr[I+1];m++) */
682 /* { */
683 /* diffusiveVelocityComponent_I -= a[m]*grad_phi[colind[m]]; */
684 /* max_a = fmax(max_a,a[m]); */
685 /* } */
686 /* flux+= diffusiveVelocityComponent_I*n[I]; */
687 /* } */
688 /* penaltyFlux = max_a*penalty*(u-bc_u); */
689 /* flux += penaltyFlux; */
690 /* } */
691 /* else if(isFluxBoundary == 1) */
692 /* { */
693 /* flux = bc_flux; */
694 /* } */
695 /* else */
696 /* { */
697 /* std::cerr<<"warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0"<<std::endl; */
698 /* flux = 0.0; */
699 /* } */
700 /* } */
701
702 /* inline */
703 /* double ExteriorNumericalDiffusiveFluxJacobian(const double& eps, */
704 /* int* rowptr, */
705 /* int* colind, */
706 /* const int& isDOFBoundary, */
707 /* const double n[nSpace], */
708 /* double* a, */
709 /* const double& v, */
710 /* const double grad_v[nSpace], */
711 /* const double& penalty) */
712 /* { */
713 /* double dvel_I,tmp=0.0,max_a=0.0; */
714 /* if(isDOFBoundary >= 1) */
715 /* { */
716 /* for(int I=0;I<nSpace;I++) */
717 /* { */
718 /* dvel_I=0.0; */
719 /* for(int m=rowptr[I];m<rowptr[I+1];m++) */
720 /* { */
721 /* dvel_I -= a[m]*grad_v[colind[m]]; */
722 /* max_a = fmax(max_a,a[m]); */
723 /* } */
724 /* tmp += dvel_I*n[I]; */
725 /* } */
726 /* tmp +=max_a*penalty*v; */
727 /* } */
728 /* return tmp; */
729 /* } */
730
731 void calculateResidual(arguments_dict& args)
732 {
733 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
734 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
735 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
736 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
737 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
738 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
739 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
740 xt::pyarray<double>& h_trial_ref = args.array<double>("h_trial_ref");
741 xt::pyarray<double>& h_grad_trial_ref = args.array<double>("h_grad_trial_ref");
742 xt::pyarray<double>& h_test_ref = args.array<double>("h_test_ref");
743 xt::pyarray<double>& h_grad_test_ref = args.array<double>("h_grad_test_ref");
744 xt::pyarray<double>& vel_trial_ref = args.array<double>("vel_trial_ref");
745 xt::pyarray<double>& vel_grad_trial_ref = args.array<double>("vel_grad_trial_ref");
746 xt::pyarray<double>& vel_test_ref = args.array<double>("vel_test_ref");
747 xt::pyarray<double>& vel_grad_test_ref = args.array<double>("vel_grad_test_ref");
748 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
749 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
750 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
751 xt::pyarray<double>& h_trial_trace_ref = args.array<double>("h_trial_trace_ref");
752 xt::pyarray<double>& h_grad_trial_trace_ref = args.array<double>("h_grad_trial_trace_ref");
753 xt::pyarray<double>& h_test_trace_ref = args.array<double>("h_test_trace_ref");
754 xt::pyarray<double>& h_grad_test_trace_ref = args.array<double>("h_grad_test_trace_ref");
755 xt::pyarray<double>& vel_trial_trace_ref = args.array<double>("vel_trial_trace_ref");
756 xt::pyarray<double>& vel_grad_trial_trace_ref = args.array<double>("vel_grad_trial_trace_ref");
757 xt::pyarray<double>& vel_test_trace_ref = args.array<double>("vel_test_trace_ref");
758 xt::pyarray<double>& vel_grad_test_trace_ref = args.array<double>("vel_grad_test_trace_ref");
759 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
760 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
761 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
762 int nElements_global = args.scalar<int>("nElements_global");
763 double useRBLES = args.scalar<double>("useRBLES");
764 double useMetrics = args.scalar<double>("useMetrics");
765 double alphaBDF = args.scalar<double>("alphaBDF");
766 double nu = args.scalar<double>("nu");
767 double g = args.scalar<double>("g");
768 double shockCapturingCoefficient = args.scalar<double>("shockCapturingCoefficient");
769 xt::pyarray<int>& h_l2g = args.array<int>("h_l2g");
770 xt::pyarray<int>& vel_l2g = args.array<int>("vel_l2g");
771 xt::pyarray<double>& b_dof = args.array<double>("b_dof");
772 xt::pyarray<double>& h_dof = args.array<double>("h_dof");
773 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
774 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
775 xt::pyarray<double>& h_dof_sge = args.array<double>("h_dof_sge");
776 xt::pyarray<double>& u_dof_sge = args.array<double>("u_dof_sge");
777 xt::pyarray<double>& v_dof_sge = args.array<double>("v_dof_sge");
778 xt::pyarray<double>& q_mass_acc = args.array<double>("q_mass_acc");
779 xt::pyarray<double>& q_mom_u_acc = args.array<double>("q_mom_u_acc");
780 xt::pyarray<double>& q_mom_v_acc = args.array<double>("q_mom_v_acc");
781 xt::pyarray<double>& q_mass_adv = args.array<double>("q_mass_adv");
782 xt::pyarray<double>& q_mass_acc_beta_bdf = args.array<double>("q_mass_acc_beta_bdf");
783 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.array<double>("q_mom_u_acc_beta_bdf");
784 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.array<double>("q_mom_v_acc_beta_bdf");
785 xt::pyarray<double>& q_velocity_sge = args.array<double>("q_velocity_sge");
786 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
787 xt::pyarray<double>& q_numDiff_h = args.array<double>("q_numDiff_h");
788 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
789 xt::pyarray<double>& q_numDiff_v = args.array<double>("q_numDiff_v");
790 xt::pyarray<double>& q_numDiff_h_last = args.array<double>("q_numDiff_h_last");
791 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
792 xt::pyarray<double>& q_numDiff_v_last = args.array<double>("q_numDiff_v_last");
793 xt::pyarray<int>& sdInfo_u_u_rowptr = args.array<int>("sdInfo_u_u_rowptr");
794 xt::pyarray<int>& sdInfo_u_u_colind = args.array<int>("sdInfo_u_u_colind");
795 xt::pyarray<int>& sdInfo_u_v_rowptr = args.array<int>("sdInfo_u_v_rowptr");
796 xt::pyarray<int>& sdInfo_u_v_colind = args.array<int>("sdInfo_u_v_colind");
797 xt::pyarray<int>& sdInfo_v_v_rowptr = args.array<int>("sdInfo_v_v_rowptr");
798 xt::pyarray<int>& sdInfo_v_v_colind = args.array<int>("sdInfo_v_v_colind");
799 xt::pyarray<int>& sdInfo_v_u_rowptr = args.array<int>("sdInfo_v_u_rowptr");
800 xt::pyarray<int>& sdInfo_v_u_colind = args.array<int>("sdInfo_v_u_colind");
801 int offset_h = args.scalar<int>("offset_h");
802 int offset_u = args.scalar<int>("offset_u");
803 int offset_v = args.scalar<int>("offset_v");
804 int stride_h = args.scalar<int>("stride_h");
805 int stride_u = args.scalar<int>("stride_u");
806 int stride_v = args.scalar<int>("stride_v");
807 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
808 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
809 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
810 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
811 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
812 xt::pyarray<int>& isDOFBoundary_h = args.array<int>("isDOFBoundary_h");
813 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
814 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
815 xt::pyarray<int>& isAdvectiveFluxBoundary_h = args.array<int>("isAdvectiveFluxBoundary_h");
816 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
817 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.array<int>("isAdvectiveFluxBoundary_v");
818 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
819 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.array<int>("isDiffusiveFluxBoundary_v");
820 xt::pyarray<double>& ebqe_bc_h_ext = args.array<double>("ebqe_bc_h_ext");
821 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.array<double>("ebqe_bc_flux_mass_ext");
822 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.array<double>("ebqe_bc_flux_mom_u_adv_ext");
823 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.array<double>("ebqe_bc_flux_mom_v_adv_ext");
824 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
825 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.array<double>("ebqe_bc_flux_u_diff_ext");
826 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
827 xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
828 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.array<double>("ebqe_bc_flux_v_diff_ext");
829 xt::pyarray<double>& q_velocity = args.array<double>("q_velocity");
830 xt::pyarray<double>& ebqe_velocity = args.array<double>("ebqe_velocity");
831 xt::pyarray<double>& flux = args.array<double>("flux");
832 xt::pyarray<double>& elementResidual_h_save = args.array<double>("elementResidual_h");
833 //
834 //loop over elements to compute volume integrals and load them into element and global residual
835 //
836 double globalConservationError=0.0,tauSum=0.0;
837 for(int eN=0;eN<nElements_global;eN++)
838 {
839 //declare local storage for element residual and initialize
840 double elementResidual_h[nDOF_test_element],
841 elementResidual_u[nDOF_test_element],
842 elementResidual_v[nDOF_test_element];
843 for (int i=0;i<nDOF_test_element;i++)
844 {
845 int eN_i = eN*nDOF_test_element+i;
846 elementResidual_h_save.data()[eN_i]=0.0;
847 elementResidual_h[i]=0.0;
848 elementResidual_u[i]=0.0;
849 elementResidual_v[i]=0.0;
850 }//i
851 //
852 //loop over quadrature points and compute integrands
853 //
854 for(int k=0;k<nQuadraturePoints_element;k++)
855 {
856 //compute indices and declare local storage
857 int eN_k = eN*nQuadraturePoints_element+k,
858 eN_k_nSpace = eN_k*nSpace,
859 eN_nDOF_trial_element = eN*nDOF_trial_element;
860 double b=0.0,h=0.0,u=0.0,v=0.0,h_sge=0.0,u_sge=0.0,v_sge=0.0,
861 grad_b[nSpace],grad_h[nSpace],grad_u[nSpace],grad_v[nSpace],
862 mass_acc=0.0,
863 dmass_acc_h=0.0,
864 mom_u_acc=0.0,
865 dmom_u_acc_h=0.0,
866 dmom_u_acc_u=0.0,
867 mom_v_acc=0.0,
868 dmom_v_acc_h=0.0,
869 dmom_v_acc_v=0.0,
870 mass_adv[nSpace],
871 dmass_adv_h[nSpace],
872 dmass_adv_u[nSpace],
873 dmass_adv_v[nSpace],
874 dmass_adv_h_sge[nSpace],
875 dmass_adv_u_sge[nSpace],
876 dmass_adv_v_sge[nSpace],
877 mom_u_adv[nSpace],
878 dmom_u_adv_h[nSpace],
879 dmom_u_adv_u[nSpace],
880 dmom_u_adv_v[nSpace],
881 dmom_u_adv_h_sge[nSpace],
882 dmom_u_adv_u_sge[nSpace],
883 dmom_u_adv_v_sge[nSpace],
884 mom_v_adv[nSpace],
885 dmom_v_adv_h[nSpace],
886 dmom_v_adv_u[nSpace],
887 dmom_v_adv_v[nSpace],
888 dmom_v_adv_h_sge[nSpace],
889 dmom_v_adv_u_sge[nSpace],
890 dmom_v_adv_v_sge[nSpace],
891 mom_u_diff_ten[nSpace],
892 mom_v_diff_ten[nSpace],
893 mom_uv_diff_ten[1],
894 mom_vu_diff_ten[1],
895 mom_u_source=0.0,
896 dmom_u_source_h=0.0,
897 mom_v_source=0.0,
898 dmom_v_source_h=0.0,
899 mass_acc_t=0.0,
900 dmass_acc_h_t=0.0,
901 mom_u_acc_t=0.0,
902 dmom_u_acc_h_t=0.0,
903 dmom_u_acc_u_t=0.0,
904 mom_v_acc_t=0.0,
905 dmom_v_acc_h_t=0.0,
906 dmom_v_acc_v_t=0.0,
907 tau[9],
908 pdeResidual_h=0.0,
909 pdeResidual_u=0.0,
910 pdeResidual_v=0.0,
911 Lstar_h_h[nDOF_test_element],
912 Lstar_u_h[nDOF_test_element],
913 Lstar_v_h[nDOF_test_element],
914 Lstar_h_u[nDOF_test_element],
915 Lstar_u_u[nDOF_test_element],
916 Lstar_v_u[nDOF_test_element],
917 Lstar_h_v[nDOF_test_element],
918 Lstar_u_v[nDOF_test_element],
919 Lstar_v_v[nDOF_test_element],
920 subgridError_h=0.0,
921 subgridError_u=0.0,
922 subgridError_v=0.0,
923 jac[nSpace*nSpace],
924 jacDet,
925 jacInv[nSpace*nSpace],
926 h_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
927 h_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
928 h_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
929 dV,x,y,xt,yt,
930 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
931 //get jacobian, etc for mapping reference element
932 ck.calculateMapping_element(eN,
933 k,
934 mesh_dof.data(),
935 mesh_l2g.data(),
936 mesh_trial_ref.data(),
937 mesh_grad_trial_ref.data(),
938 jac,
939 jacDet,
940 jacInv,
941 x,y);
942 ck.calculateMappingVelocity_element(eN,
943 k,
944 mesh_velocity_dof.data(),
945 mesh_l2g.data(),
946 mesh_trial_ref.data(),
947 xt,yt);
948 //xt=0.0;yt=0.0;
949 //std::cout<<"xt "<<xt<<'\t'<<yt<<'\t'<<std::endl;
950 //get the physical integration weight
951 dV = fabs(jacDet)*dV_ref.data()[k];
952 //get the trial function gradients
953 ck.gradTrialFromRef(&h_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,h_grad_trial);
954 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
955 //get the solution
956 ck.valFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],b);
957 ck.valFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h);
958 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u);
959 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v);
960 ck.valFromDOF(h_dof_sge.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h_sge);
961 ck.valFromDOF(u_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u_sge);
962 ck.valFromDOF(v_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v_sge);
963 //get the solution gradients
964 ck.gradFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_b);
965 ck.gradFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_h);
966 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
967 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
968 //precalculate test function products with integration weights
969 for (int j=0;j<nDOF_trial_element;j++)
970 {
971 h_test_dV[j] = h_test_ref.data()[k*nDOF_trial_element+j]*dV;
972 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
973 for (int I=0;I<nSpace;I++)
974 {
975 h_grad_test_dV[j*nSpace+I] = h_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
976 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
977 }
978 }
979 //save velocity at quadrature points for other models to use
980 q_velocity.data()[eN_k_nSpace+0]=u;
981 q_velocity.data()[eN_k_nSpace+1]=v;
982 //
983 //calculate pde coefficients at quadrature points
984 //
985 evaluateCoefficients(nu,
986 g,
987 grad_b,
988 h,
989 u,
990 v,
991 mass_acc,
992 dmass_acc_h,
993 mom_u_acc,
994 dmom_u_acc_h,
995 dmom_u_acc_u,
996 mom_v_acc,
997 dmom_v_acc_h,
998 dmom_v_acc_v,
999 mass_adv,
1000 dmass_adv_h,
1001 dmass_adv_u,
1002 dmass_adv_v,
1003 mom_u_adv,
1004 dmom_u_adv_h,
1005 dmom_u_adv_u,
1006 dmom_u_adv_v,
1007 mom_v_adv,
1008 dmom_v_adv_h,
1009 dmom_v_adv_u,
1010 dmom_v_adv_v,
1011 mom_u_diff_ten,
1012 mom_v_diff_ten,
1013 mom_uv_diff_ten,
1014 mom_vu_diff_ten,
1015 mom_u_source,
1016 dmom_u_source_h,
1017 mom_v_source,
1018 dmom_v_source_h);
1019 //
1020 //save momentum for time history and velocity for subgrid error
1021 //
1022 q_mass_acc.data()[eN_k] = mass_acc;
1023 q_mom_u_acc.data()[eN_k] = mom_u_acc;
1024 q_mom_v_acc.data()[eN_k] = mom_v_acc;
1025 //subgrid error uses grid scale discharge
1026 q_mass_adv.data()[eN_k_nSpace+0] = h*u;
1027 q_mass_adv.data()[eN_k_nSpace+1] = h*v;
1028 //
1029 //moving mesh
1030 //
1031 //transform the continuity equation as if the accumulation term was d(1)/dt
1032 /* mass_adv[0] -= MOVING_DOMAIN*mass_acc*xt; */
1033 /* mass_adv[1] -= MOVING_DOMAIN*mass_acc*yt; */
1034 /* dmass_adv_h[0] -= MOVING_DOMAIN*dmass_acc_h*xt; */
1035 /* dmass_adv_h[1] -= MOVING_DOMAIN*dmass_acc_h*yt; */
1036
1037 /* mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*xt; */
1038 /* mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt; */
1039 /* dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*xt; */
1040 /* dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt; */
1041
1042 /* mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*xt; */
1043 /* mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt; */
1044 /* dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*xt; */
1045 /* dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt; */
1046
1047 //
1048 //calculate time derivative at quadrature points
1049 //
1050 ck.bdf(alphaBDF,
1051 q_mass_acc_beta_bdf.data()[eN_k],
1052 mass_acc,
1053 dmass_acc_h,
1054 mass_acc_t,
1055 dmass_acc_h_t);
1056 ck.bdfC2(alphaBDF,
1057 q_mom_u_acc_beta_bdf.data()[eN_k],
1058 mom_u_acc,
1059 dmom_u_acc_h,
1060 dmom_u_acc_u,
1061 mom_u_acc_t,
1062 dmom_u_acc_h_t,
1063 dmom_u_acc_u_t);
1064 ck.bdfC2(alphaBDF,
1065 q_mom_v_acc_beta_bdf.data()[eN_k],
1066 mom_v_acc,
1067 dmom_v_acc_h,
1068 dmom_v_acc_v,
1069 mom_v_acc_t,
1070 dmom_v_acc_h_t,
1071 dmom_v_acc_v_t);
1072 //
1073 //calculate subgrid error (strong residual and adjoint)
1074 //
1075 //calculate strong residual
1076 dmass_adv_h_sge[0] = dmass_adv_h[0];
1077 dmass_adv_h_sge[1] = dmass_adv_h[1];
1078 dmass_adv_u_sge[0] = dmass_adv_u[0];
1079 dmass_adv_u_sge[1] = dmass_adv_u[1];
1080 dmass_adv_v_sge[0] = dmass_adv_v[0];
1081 dmass_adv_v_sge[1] = dmass_adv_v[1];
1082 dmom_u_adv_h_sge[0] = dmom_u_adv_h[0];
1083 dmom_u_adv_h_sge[1] = dmom_u_adv_h[1];
1084 dmom_u_adv_u_sge[0] = dmom_u_adv_u[0];
1085 dmom_u_adv_u_sge[1] = dmom_u_adv_u[1];
1086 dmom_u_adv_v_sge[0] = dmom_u_adv_v[0];
1087 dmom_u_adv_v_sge[1] = dmom_u_adv_v[1];
1088 dmom_v_adv_h_sge[0] = dmom_v_adv_h[0];
1089 dmom_v_adv_h_sge[1] = dmom_v_adv_h[1];
1090 dmom_v_adv_u_sge[0] = dmom_v_adv_u[0];
1091 dmom_v_adv_u_sge[1] = dmom_v_adv_u[1];
1092 dmom_v_adv_v_sge[0] = dmom_v_adv_v[0];
1093 dmom_v_adv_v_sge[1] = dmom_v_adv_v[1];
1094
1095 /* //mass advective flux */
1096 /* dmass_adv_h_sge[0]=u_sge; */
1097 /* dmass_adv_h_sge[1]=v_sge; */
1098
1099 /* dmass_adv_u_sge[0]=0.0;//h_sge; */
1100 /* dmass_adv_u_sge[1]=0.0; */
1101
1102 /* dmass_adv_v_sge[0]=0.0; */
1103 /* dmass_adv_v_sge[1]=0.0;//h_sge; */
1104
1105 /* //u momentum advective flux */
1106 /* dmom_u_adv_h_sge[0]=0.0;//u_sge*u_sge + g*h_sge; */
1107 /* dmom_u_adv_h_sge[1]=0.0;//u_sge*v_sge; */
1108
1109 /* dmom_u_adv_u_sge[0]=h_sge*u_sge;//h_sge*2.0*u_sge; */
1110 /* dmom_u_adv_u_sge[1]=h_sge*v_sge; */
1111
1112 /* dmom_u_adv_v_sge[0]=0.0; */
1113 /* dmom_u_adv_v_sge[1]=0.0;//h_sge*u_sge; */
1114
1115 /* //v momentum advective_flux */
1116 /* dmom_v_adv_h_sge[0]=0.0;//v_sge*u_sge; */
1117 /* dmom_v_adv_h_sge[1]=0.0;//v_sge*v_sge + g*h_sge; */
1118
1119 /* dmom_v_adv_u_sge[0]=0.0;//h_sge*v_sge; */
1120 /* dmom_v_adv_u_sge[1]=0.0; */
1121
1122 /* dmom_v_adv_v_sge[0]=h_sge*u_sge; */
1123 /* dmom_v_adv_v_sge[1]=h_sge*v_sge;//h_sge*2.0*v_sge; */
1124
1125
1126 //full linearization, lagged
1127
1128 //mass advective flux
1129 dmass_adv_h_sge[0]=u_sge;
1130 dmass_adv_h_sge[1]=v_sge;
1131
1132 dmass_adv_u_sge[0]=h_sge;
1133 dmass_adv_u_sge[1]=0.0;
1134
1135 dmass_adv_v_sge[0]=0.0;
1136 dmass_adv_v_sge[1]=h_sge;
1137
1138 //u momentum advective flux
1139 dmom_u_adv_h_sge[0]=u_sge*u_sge + g*h_sge;
1140 dmom_u_adv_h_sge[1]=u_sge*v_sge;
1141
1142 dmom_u_adv_u_sge[0]=h_sge*2.0*u_sge;
1143 dmom_u_adv_u_sge[1]=h_sge*v_sge;
1144
1145 dmom_u_adv_v_sge[0]=0.0;
1146 dmom_u_adv_v_sge[1]=h_sge*u_sge;
1147
1148 //v momentum advective_flux
1149 dmom_v_adv_h_sge[0]=v_sge*u_sge;
1150 dmom_v_adv_h_sge[1]=v_sge*v_sge + g*h_sge;
1151
1152 dmom_v_adv_u_sge[0]=h_sge*v_sge;
1153 dmom_v_adv_u_sge[1]=0.0;
1154
1155 dmom_v_adv_v_sge[0]=h_sge*u_sge;
1156 dmom_v_adv_v_sge[1]=h_sge*2.0*v_sge;
1157
1158 pdeResidual_h = ck.Mass_strong(mass_acc_t) +
1159 ck.Advection_strong(dmass_adv_h_sge,grad_h) +
1160 ck.Advection_strong(dmass_adv_u_sge,grad_u) +
1161 ck.Advection_strong(dmass_adv_v_sge,grad_v);
1162
1163 pdeResidual_u = ck.Mass_strong(mom_u_acc_t) +
1164 ck.Advection_strong(dmom_u_adv_h_sge,grad_h) +
1165 ck.Advection_strong(dmom_u_adv_u_sge,grad_u) +
1166 ck.Advection_strong(dmom_u_adv_v_sge,grad_v) +
1167 ck.Reaction_strong(mom_u_source);
1168
1169 pdeResidual_v = ck.Mass_strong(mom_v_acc_t) +
1170 ck.Advection_strong(dmom_v_adv_h_sge,grad_h) +
1171 ck.Advection_strong(dmom_v_adv_u_sge,grad_u) +
1172 ck.Advection_strong(dmom_v_adv_v_sge,grad_v) +
1173 ck.Reaction_strong(mom_v_source);
1174
1175 calculateSubgridError_tau(elementDiameter.data()[eN],
1176 nu,
1177 g,
1178 h_sge,
1179 u_sge,
1180 v_sge,
1181 tau,
1182 q_cfl.data()[eN_k]);
1183 for (int i=0;i<9;i++)
1184 tauSum += tau[i];
1185
1186 subgridError_h = - tau[0*3+0]*pdeResidual_h - tau[0*3+1]*pdeResidual_u - tau[0*3+2]*pdeResidual_v;
1187 subgridError_u = - tau[1*3+0]*pdeResidual_h - tau[1*3+1]*pdeResidual_u - tau[1*3+2]*pdeResidual_v;
1188 subgridError_v = - tau[2*3+0]*pdeResidual_h - tau[2*3+1]*pdeResidual_u - tau[2*3+2]*pdeResidual_v;
1189
1190 //adjoint times the test functions
1191 for (int i=0;i<nDOF_test_element;i++)
1192 {
1193 int i_nSpace = i*nSpace;
1194 Lstar_h_h[i]=ck.Advection_adjoint(dmass_adv_h_sge,&h_grad_test_dV[i_nSpace]);
1195 Lstar_u_h[i]=ck.Advection_adjoint(dmass_adv_u_sge,&h_grad_test_dV[i_nSpace]);
1196 Lstar_v_h[i]=ck.Advection_adjoint(dmass_adv_v_sge,&h_grad_test_dV[i_nSpace]);
1197
1198 Lstar_h_u[i]=ck.Advection_adjoint(dmom_u_adv_h_sge,&vel_grad_test_dV[i_nSpace]) +
1199 ck.Reaction_adjoint(dmom_u_source_h,vel_test_dV[i]);
1200 Lstar_u_u[i]=ck.Advection_adjoint(dmom_u_adv_u_sge,&vel_grad_test_dV[i_nSpace]);
1201 Lstar_v_u[i]=ck.Advection_adjoint(dmom_u_adv_v_sge,&vel_grad_test_dV[i_nSpace]);
1202
1203 Lstar_h_v[i]=ck.Advection_adjoint(dmom_v_adv_h_sge,&vel_grad_test_dV[i_nSpace]) +
1204 ck.Reaction_adjoint(dmom_v_source_h,vel_test_dV[i]);
1205 Lstar_u_v[i]=ck.Advection_adjoint(dmom_v_adv_u_sge,&vel_grad_test_dV[i_nSpace]);
1206 Lstar_v_v[i]=ck.Advection_adjoint(dmom_v_adv_v_sge,&vel_grad_test_dV[i_nSpace]);
1207
1208 }
1209
1210 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v);
1211 /* double */
1212 double norm_grad = 1.0;
1213 q_numDiff_u.data()[eN_k] = 0.5*elementDiameter.data()[eN]*norm_Rv/(norm_grad+1.0e-8);
1214 q_numDiff_v.data()[eN_k] = q_numDiff_u.data()[eN_k];
1215
1216 /* ck.calculateNumericalDiffusion(1.0, */
1217 /* elementDiameter.data()[eN], */
1218 /* pdeResidual_h, */
1219 /* grad_h, */
1220 /* q_numDiff_h.data()[eN_k]); */
1221
1222 //update element residual
1223
1224 for(int i=0;i<nDOF_test_element;i++)
1225 {
1226 int i_nSpace=i*nSpace;
1227
1228 elementResidual_h[i] += ck.Mass_weak(mass_acc_t,h_test_dV[i]) +
1229 ck.Advection_weak(mass_adv,&h_grad_test_dV[i_nSpace]) +
1230 ck.SubgridError(subgridError_h,Lstar_h_h[i]) +
1231 ck.SubgridError(subgridError_u,Lstar_u_h[i]) +
1232 ck.SubgridError(subgridError_v,Lstar_v_h[i]) +
1233 ck.NumericalDiffusion(q_numDiff_h_last.data()[eN_k],grad_h,&h_grad_test_dV[i_nSpace]);
1234
1235 elementResidual_u[i] += ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
1236 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
1237 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_u_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
1238 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
1239 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
1240 ck.SubgridError(subgridError_h,Lstar_h_u[i]) +
1241 ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
1242 ck.SubgridError(subgridError_v,Lstar_v_u[i]) +
1243 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]);
1244
1245 elementResidual_v[i] += ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
1246 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
1247 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_v_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
1248 ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
1249 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
1250 ck.SubgridError(subgridError_h,Lstar_h_v[i]) +
1251 ck.SubgridError(subgridError_u,Lstar_u_v[i]) +
1252 ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
1253 ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]);
1254 }
1255 }
1256
1257 //load element into global residual and save element residual
1258
1259 for(int i=0;i<nDOF_test_element;i++)
1260 {
1261 int eN_i=eN*nDOF_test_element+i;
1262
1263 elementResidual_h_save.data()[eN_i] += elementResidual_h[i];
1264
1265 globalResidual.data()[offset_h+stride_h*h_l2g.data()[eN_i]]+=elementResidual_h[i];
1266 globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
1267 globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
1268 }
1269 }
1270 std::cout<<"tauSum = "<<tauSum<<std::endl;
1271 /* */
1272 /* loop over exterior element boundaries to calculate surface integrals and load into element and global residuals */
1273 /* */
1274 /* ebNE is the Exterior element boundary INdex */
1275 /* ebN is the element boundary INdex */
1276 /* eN is the element index */
1277 /* for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) */
1278 /* { */
1279 /* int ebN = exteriorElementBoundariesArray.data()[ebNE], */
1280 /* eN = elementBoundaryElementsArray.data()[ebN*2+0], */
1281 /* ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0], */
1282 /* eN_nDOF_trial_element = eN*nDOF_trial_element; */
1283 /* double elementResidual_h[nDOF_test_element], */
1284 /* elementResidual_u[nDOF_test_element], */
1285 /* elementResidual_v[nDOF_test_element], */
1286 /* eps_rho,eps_mu; */
1287 /* for (int i=0;i<nDOF_test_element;i++) */
1288 /* { */
1289 /* elementResidual_h[i]=0.0; */
1290 /* elementResidual_u[i]=0.0; */
1291 /* elementResidual_v[i]=0.0; */
1292 /* } */
1293 /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
1294 /* { */
1295 /* int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb, */
1296 /* ebNE_kb_nSpace = ebNE_kb*nSpace, */
1297 /* ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb, */
1298 /* ebN_local_kb_nSpace = ebN_local_kb*nSpace; */
1299 /* double h_ext=0.0, */
1300 /* u_ext=0.0, */
1301 /* v_ext=0.0, */
1302 /* grad_h_ext[nSpace], */
1303 /* grad_u_ext[nSpace], */
1304 /* grad_v_ext[nSpace], */
1305 /* mom_u_acc_ext=0.0, */
1306 /* dmom_u_acc_u_ext=0.0, */
1307 /* mom_v_acc_ext=0.0, */
1308 /* dmom_v_acc_v_ext=0.0, */
1309 /* mass_adv_ext[nSpace], */
1310 /* dmass_adv_u_ext[nSpace], */
1311 /* dmass_adv_v_ext[nSpace], */
1312 /* mom_u_adv_ext[nSpace], */
1313 /* dmom_u_adv_u_ext[nSpace], */
1314 /* dmom_u_adv_v_ext[nSpace], */
1315 /* mom_v_adv_ext[nSpace], */
1316 /* dmom_v_adv_u_ext[nSpace], */
1317 /* dmom_v_adv_v_ext[nSpace], */
1318 /* mom_u_diff_ten_ext[nSpace], */
1319 /* mom_v_diff_ten_ext[nSpace], */
1320 /* mom_uv_diff_ten_ext[1], */
1321 /* mom_vu_diff_ten_ext[1], */
1322 /* mom_u_source_ext=0.0, */
1323 /* mom_v_source_ext=0.0, */
1324 /* mom_u_ham_ext=0.0, */
1325 /* dmom_u_ham_grad_h_ext[nSpace], */
1326 /* mom_v_ham_ext=0.0, */
1327 /* dmom_v_ham_grad_h_ext[nSpace], */
1328 /* dmom_u_adv_h_ext[nSpace], */
1329 /* dmom_v_adv_h_ext[nSpace], */
1330 /* flux_mass_ext=0.0, */
1331 /* flux_mom_u_adv_ext=0.0, */
1332 /* flux_mom_v_adv_ext=0.0, */
1333 /* flux_mom_u_diff_ext=0.0, */
1334 /* flux_mom_v_diff_ext=0.0, */
1335 /* bc_h_ext=0.0, */
1336 /* bc_u_ext=0.0, */
1337 /* bc_v_ext=0.0, */
1338 /* bc_mom_u_acc_ext=0.0, */
1339 /* bc_dmom_u_acc_u_ext=0.0, */
1340 /* bc_mom_v_acc_ext=0.0, */
1341 /* bc_dmom_v_acc_v_ext=0.0, */
1342 /* bc_mass_adv_ext[nSpace], */
1343 /* bc_dmass_adv_u_ext[nSpace], */
1344 /* bc_dmass_adv_v_ext[nSpace], */
1345 /* bc_mom_u_adv_ext[nSpace], */
1346 /* bc_dmom_u_adv_u_ext[nSpace], */
1347 /* bc_dmom_u_adv_v_ext[nSpace], */
1348 /* bc_mom_v_adv_ext[nSpace], */
1349 /* bc_dmom_v_adv_u_ext[nSpace], */
1350 /* bc_dmom_v_adv_v_ext[nSpace], */
1351 /* bc_mom_u_diff_ten_ext[nSpace], */
1352 /* bc_mom_v_diff_ten_ext[nSpace], */
1353 /* bc_mom_uv_diff_ten_ext[1], */
1354 /* bc_mom_vu_diff_ten_ext[1], */
1355 /* bc_mom_u_source_ext=0.0, */
1356 /* bc_mom_v_source_ext=0.0, */
1357 /* bc_mom_u_ham_ext=0.0, */
1358 /* bc_dmom_u_ham_grad_h_ext[nSpace], */
1359 /* bc_mom_v_ham_ext=0.0, */
1360 /* bc_dmom_v_ham_grad_h_ext[nSpace], */
1361 /* jac_ext[nSpace*nSpace], */
1362 /* jacDet_ext, */
1363 /* jacInv_ext[nSpace*nSpace], */
1364 /* boundaryJac[nSpace*(nSpace-1)], */
1365 /* metricTensor[(nSpace-1)*(nSpace-1)], */
1366 /* metricTensorDetSqrt, */
1367 /* dS,h_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element], */
1368 /* h_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace], */
1369 /* normal[3],x_ext,y_ext,xt_ext,yt_ext,integralScaling, */
1370 /* G[nSpace*nSpace],G_dd_G,tr_G,h_penalty; */
1371 /* compute information about mapping from reference element to physical element */
1372 /* ck.calculateMapping_elementBoundary(eN, */
1373 /* ebN_local, */
1374 /* kb, */
1375 /* ebN_local_kb, */
1376 /* mesh_dof.data(), */
1377 /* mesh_l2g.data(), */
1378 /* mesh_trial_trace_ref.data(), */
1379 /* mesh_grad_trial_trace_ref.data(), */
1380 /* boundaryJac_ref.data(), */
1381 /* jac_ext, */
1382 /* jacDet_ext, */
1383 /* jacInv_ext, */
1384 /* boundaryJac, */
1385 /* metricTensor, */
1386 /* metricTensorDetSqrt, */
1387 /* normal_ref.data(), */
1388 /* normal, */
1389 /* x_ext,y_ext); */
1390 /* ck.calculateMappingVelocity_elementBoundary(eN, */
1391 /* ebN_local, */
1392 /* kb, */
1393 /* ebN_local_kb, */
1394 /* mesh_velocity_dof.data(), */
1395 /* mesh_l2g.data(), */
1396 /* mesh_trial_trace_ref.data(), */
1397 /* xt_ext,yt_ext, */
1398 /* normal, */
1399 /* boundaryJac, */
1400 /* metricTensor, */
1401 /* integralScaling); */
1402 /* xt_ext=0.0;yt_ext=0.0; */
1403 /* std::cout<<"xt_ext "<<xt_ext<<'\t'<<yt_ext<<'\t'<<std::endl; */
1404 /* std::cout<<"integralScaling - metricTensorDetSrt ==============================="<<integralScaling-metricTensorDetSqrt<<std::endl; */
1405 /* dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb]; */
1406 /* compute shape and solution information */
1407 /* shape */
1408 /* ck.gradTrialFromRef(&h_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,h_grad_trial_trace); */
1409 /* ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace); */
1410 /* solution and gradients */
1411 /* ck.valFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],h_ext); */
1412 /* ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext); */
1413 /* ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext); */
1414 /* ck.gradFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial_trace,grad_h_ext); */
1415 /* ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext); */
1416 /* ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext); */
1417 /* precalculate test function products with integration weights */
1418 /* for (int j=0;j<nDOF_trial_element;j++) */
1419 /* { */
1420 /* h_test_dS[j] = h_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
1421 /* vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
1422 /* } */
1423 /* // */
1424 /* //debugging section for finite element calculations on exterior */
1425 /* // */
1426 /* std::cout<<"ebNE = "<<ebNE<<" kb = "<<kb<<std::endl; */
1427 /* for (int j=0;j<nDOF_trial_element;j++) */
1428 /* { */
1429 /* std::cout<<"h_trial_trace["<<j<<"] "<<h_trial_trace_ref.data()[ebN_local_kb*nDOF_trial_element+j]<<std::endl; */
1430 /* std::cout<<"vel_trial_trace["<<j<<"] "<<vel_trial_trace_ref.data()[ebN_local_kb*nDOF_trial_element+j]<<std::endl; */
1431 /* std::cout<<"h_test_dS["<<j<<"] "<<h_test_dS[j]<<std::endl; */
1432 /* std::cout<<"vel_test_dS["<<j<<"] "<<vel_test_dS[j]<<std::endl; */
1433 /* for (int I=0;I<nSpace;I++) */
1434 /* { */
1435 /* std::cout<<"h_grad_trial_trace["<<j<<","<<I<<"] "<<h_grad_trial_trace[j*nSpace+I]<<std::endl; */
1436 /* std::cout<<"vel_grad_trial_trace["<<j<<","<<I<<"] "<<vel_grad_trial_trace[j*nSpace+I]<<std::endl; */
1437 /* } */
1438 /* } */
1439 /* std::cout<<"h_ext "<<h_ext<<std::endl; */
1440 /* std::cout<<"u_ext "<<u_ext<<std::endl; */
1441 /* std::cout<<"v_ext "<<v_ext<<std::endl; */
1442 /* for(int I=0;I<nSpace;I++) */
1443 /* { */
1444 /* std::cout<<"grad_h_ext["<<I<<"] "<<grad_h_ext[I]<<std::endl; */
1445 /* std::cout<<"grad_u_ext["<<I<<"] "<<grad_u_ext[I]<<std::endl; */
1446 /* std::cout<<"grad_v_ext["<<I<<"] "<<grad_v_ext[I]<<std::endl; */
1447 /* } */
1448 /* */
1449 /* load the boundary values */
1450 /* */
1451 /* bc_h_ext = isDOFBoundary_h.data()[ebNE_kb]*ebqe_bc_h_ext.data()[ebNE_kb]+(1-isDOFBoundary_h.data()[ebNE_kb])*h_ext; */
1452 /* bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext; */
1453 /* bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext; */
1454 /* */
1455 /* calculate the pde coefficients using the solution and the boundary values for the solution */
1456 /* */
1457 /* cek debug */
1458 /* eps_rho=0.1; */
1459 /* eps_mu=0.1; */
1460 /* evaluateCoefficients(eps_rho, */
1461 /* eps_mu, */
1462 /* sigma, */
1463 /* rho_0, */
1464 /* nu_0, */
1465 /* rho_1, */
1466 /* nu_1, */
1467 /* g, */
1468 /* h_ext, */
1469 /* grad_h_ext, */
1470 /* u_ext, */
1471 /* v_ext, */
1472 /* mom_u_acc_ext, */
1473 /* dmom_u_acc_u_ext, */
1474 /* mom_v_acc_ext, */
1475 /* dmom_v_acc_v_ext, */
1476 /* mass_adv_ext, */
1477 /* dmass_adv_u_ext, */
1478 /* dmass_adv_v_ext, */
1479 /* mom_u_adv_ext, */
1480 /* dmom_u_adv_u_ext, */
1481 /* dmom_u_adv_v_ext, */
1482 /* mom_v_adv_ext, */
1483 /* dmom_v_adv_u_ext, */
1484 /* dmom_v_adv_v_ext, */
1485 /* mom_u_diff_ten_ext, */
1486 /* mom_v_diff_ten_ext, */
1487 /* mom_uv_diff_ten_ext, */
1488 /* mom_vu_diff_ten_ext, */
1489 /* mom_u_source_ext, */
1490 /* mom_v_source_ext, */
1491 /* mom_u_ham_ext, */
1492 /* dmom_u_ham_grad_h_ext, */
1493 /* mom_v_ham_ext, */
1494 /* dmom_v_ham_grad_h_ext); */
1495 /* evaluateCoefficients(eps_rho, */
1496 /* eps_mu, */
1497 /* sigma, */
1498 /* rho_0, */
1499 /* nu_0, */
1500 /* rho_1, */
1501 /* nu_1, */
1502 /* g, */
1503 /* bc_h_ext, */
1504 /* grad_h_ext,cek should't be used */
1505 /* bc_u_ext, */
1506 /* bc_v_ext, */
1507 /* bc_mom_u_acc_ext, */
1508 /* bc_dmom_u_acc_u_ext, */
1509 /* bc_mom_v_acc_ext, */
1510 /* bc_dmom_v_acc_v_ext, */
1511 /* bc_mass_adv_ext, */
1512 /* bc_dmass_adv_u_ext, */
1513 /* bc_dmass_adv_v_ext, */
1514 /* bc_mom_u_adv_ext, */
1515 /* bc_dmom_u_adv_u_ext, */
1516 /* bc_dmom_u_adv_v_ext, */
1517 /* bc_mom_v_adv_ext, */
1518 /* bc_dmom_v_adv_u_ext, */
1519 /* bc_dmom_v_adv_v_ext, */
1520 /* bc_mom_u_diff_ten_ext, */
1521 /* bc_mom_v_diff_ten_ext, */
1522 /* bc_mom_uv_diff_ten_ext, */
1523 /* bc_mom_vu_diff_ten_ext, */
1524 /* bc_mom_u_source_ext, */
1525 /* bc_mom_v_source_ext, */
1526 /* bc_mom_u_ham_ext, */
1527 /* bc_dmom_u_ham_grad_h_ext, */
1528 /* bc_mom_v_ham_ext, */
1529 /* bc_dmom_v_ham_grad_h_ext); */
1530 /* */
1531 /* moving domain */
1532 /* */
1533 /* mass_adv_ext[0] -= MOVING_DOMAIN*xt_ext; */
1534 /* mass_adv_ext[1] -= MOVING_DOMAIN*yt_ext; */
1535
1536 /* mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext; */
1537 /* mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext; */
1538 /* dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext; */
1539 /* dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext; */
1540
1541 /* mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext; */
1542 /* mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext; */
1543 /* dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext; */
1544 /* dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext; */
1545
1546 /* bc's */
1547 /* bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext; */
1548 /* bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext; */
1549
1550 /* bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext; */
1551 /* bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext; */
1552
1553 /* */
1554 /* calculate the numerical fluxes */
1555 /* */
1556 /* cek debug */
1557 /* ebqe_penalty_ext.data()[ebNE_kb] = 10.0; */
1558 /* */
1559 /* ck.calculateGScale(G,normal,h_penalty); */
1560 /* h_penalty = 10.0/h_penalty; */
1561 /* cek debug, do it the old way */
1562 /* h_penalty = 100.0/elementDiameter.data()[eN]; */
1563 /* exteriorNumericalAdvectiveFlux(isDOFBoundary_h.data()[ebNE_kb], */
1564 /* isDOFBoundary_u.data()[ebNE_kb], */
1565 /* isDOFBoundary_v.data()[ebNE_kb], */
1566 /* isAdvectiveFluxBoundary_h.data()[ebNE_kb], */
1567 /* isAdvectiveFluxBoundary_u.data()[ebNE_kb], */
1568 /* isAdvectiveFluxBoundary_v.data()[ebNE_kb], */
1569 /* dmom_u_ham_grad_h_ext[0],=1/rho, */
1570 /* normal, */
1571 /* bc_h_ext, */
1572 /* bc_mass_adv_ext, */
1573 /* bc_mom_u_adv_ext, */
1574 /* bc_mom_v_adv_ext, */
1575 /* ebqe_bc_flux_mass_ext.data()[ebNE_kb], */
1576 /* ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb], */
1577 /* ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb], */
1578 /* h_ext, */
1579 /* mass_adv_ext, */
1580 /* mom_u_adv_ext, */
1581 /* mom_v_adv_ext, */
1582 /* dmass_adv_u_ext, */
1583 /* dmass_adv_v_ext, */
1584 /* dmom_u_adv_h_ext, */
1585 /* dmom_u_adv_u_ext, */
1586 /* dmom_u_adv_v_ext, */
1587 /* dmom_v_adv_h_ext, */
1588 /* dmom_v_adv_u_ext, */
1589 /* dmom_v_adv_v_ext, */
1590 /* flux_mass_ext, */
1591 /* flux_mom_u_adv_ext, */
1592 /* flux_mom_v_adv_ext, */
1593 /* &ebqe_velocity.data()[ebNE_kb_nSpace]); */
1594 /* cek todo need to switch to full stress and add adjoint consistency */
1595 /* exteriorNumericalDiffusiveFlux(eps_rho, */
1596 /* sdInfo_u_u_rowptr.data(), */
1597 /* sdInfo_u_u_colind.data(), */
1598 /* isDOFBoundary_u.data()[ebNE_kb], */
1599 /* isDiffusiveFluxBoundary_u.data()[ebNE_kb], */
1600 /* normal, */
1601 /* bc_mom_u_diff_ten_ext, */
1602 /* bc_u_ext, */
1603 /* ebqe_bc_flux_u_diff_ext.data()[ebNE_kb], */
1604 /* mom_u_diff_ten_ext, */
1605 /* grad_u_ext, */
1606 /* u_ext, */
1607 /* h_penalty,ebqe_penalty_ext.data()[ebNE_kb], */
1608 /* flux_mom_u_diff_ext); */
1609 /* exteriorNumericalDiffusiveFlux(eps_rho, */
1610 /* sdInfo_v_v_rowptr.data(), */
1611 /* sdInfo_v_v_colind.data(), */
1612 /* isDOFBoundary_v.data()[ebNE_kb], */
1613 /* isDiffusiveFluxBoundary_v.data()[ebNE_kb], */
1614 /* normal, */
1615 /* bc_mom_v_diff_ten_ext, */
1616 /* bc_v_ext, */
1617 /* ebqe_bc_flux_v_diff_ext.data()[ebNE_kb], */
1618 /* mom_v_diff_ten_ext, */
1619 /* grad_v_ext, */
1620 /* v_ext, */
1621 /* h_penalty,ebqe_penalty_ext.data()[ebNE_kb], */
1622 /* flux_mom_v_diff_ext); */
1623 /* flux.data()[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext; */
1624 /* flux.data()[ebN*nQuadraturePoints_elementBoundary+kb] = 0.0;//cek debug */
1625 /* */
1626 /* update residuals */
1627 /* */
1628 /* for (int i=0;i<nDOF_test_element;i++) */
1629 /* { */
1630 /* elementResidual_h[i] += ck.ExteriorElementBoundaryFlux(flux_mass_ext,h_test_dS[i]); */
1631 /* globalConservationError += ck.ExteriorElementBoundaryFlux(flux_mass_ext,h_test_dS[i]); */
1632
1633 /* elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+ */
1634 /* ck.ExteriorElementBoundaryFlux(flux_mom_u_diff_ext,vel_test_dS[i]); */
1635
1636 /* elementResidual_v[i] += ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) + */
1637 /* ck.ExteriorElementBoundaryFlux(flux_mom_v_diff_ext,vel_test_dS[i]); */
1638
1639 /* }i */
1640 /* }kb */
1641 /* */
1642 /* update the element and global residual storage */
1643 /* */
1644 /* for (int i=0;i<nDOF_test_element;i++) */
1645 /* { */
1646 /* int eN_i = eN*nDOF_test_element+i; */
1647
1648 /* elementResidual_h_save.data()[eN_i] += elementResidual_h[i]; */
1649
1650 /* globalResidual.data()[offset_h+stride_h*h_l2g.data()[eN_i]]+=elementResidual_h[i]; */
1651 /* globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i]; */
1652 /* globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i]; */
1653 /* }i */
1654 /* // */
1655 /* //debug */
1656 /* // */
1657 /* for(int i=0;i<nDOF_test_element;i++) */
1658 /* { */
1659 /* std::cout<<"ebNE "<<ebNE<<" i "<<i<<std::endl; */
1660 /* std::cout<<"r_h"<<elementResidual_h[i]<<std::endl; */
1661 /* std::cout<<"r_u"<<elementResidual_u[i]<<std::endl; */
1662 /* std::cout<<"r_v"<<elementResidual_v[i]<<std::endl; */
1663 /* } */
1664
1665 /* }ebNE */
1666 }
1667
1668 void calculateResidual_supg(arguments_dict& args)
1669 {
1670 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1671 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1672 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1673 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1674 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1675 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1676 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1677 xt::pyarray<double>& h_trial_ref = args.array<double>("h_trial_ref");
1678 xt::pyarray<double>& h_grad_trial_ref = args.array<double>("h_grad_trial_ref");
1679 xt::pyarray<double>& h_test_ref = args.array<double>("h_test_ref");
1680 xt::pyarray<double>& h_grad_test_ref = args.array<double>("h_grad_test_ref");
1681 xt::pyarray<double>& vel_trial_ref = args.array<double>("vel_trial_ref");
1682 xt::pyarray<double>& vel_grad_trial_ref = args.array<double>("vel_grad_trial_ref");
1683 xt::pyarray<double>& vel_test_ref = args.array<double>("vel_test_ref");
1684 xt::pyarray<double>& vel_grad_test_ref = args.array<double>("vel_grad_test_ref");
1685 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1686 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1687 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1688 xt::pyarray<double>& h_trial_trace_ref = args.array<double>("h_trial_trace_ref");
1689 xt::pyarray<double>& h_grad_trial_trace_ref = args.array<double>("h_grad_trial_trace_ref");
1690 xt::pyarray<double>& h_test_trace_ref = args.array<double>("h_test_trace_ref");
1691 xt::pyarray<double>& h_grad_test_trace_ref = args.array<double>("h_grad_test_trace_ref");
1692 xt::pyarray<double>& vel_trial_trace_ref = args.array<double>("vel_trial_trace_ref");
1693 xt::pyarray<double>& vel_grad_trial_trace_ref = args.array<double>("vel_grad_trial_trace_ref");
1694 xt::pyarray<double>& vel_test_trace_ref = args.array<double>("vel_test_trace_ref");
1695 xt::pyarray<double>& vel_grad_test_trace_ref = args.array<double>("vel_grad_test_trace_ref");
1696 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1697 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1698 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1699 int nElements_global = args.scalar<int>("nElements_global");
1700 double useRBLES = args.scalar<double>("useRBLES");
1701 double useMetrics = args.scalar<double>("useMetrics");
1702 double alphaBDF = args.scalar<double>("alphaBDF");
1703 double nu = args.scalar<double>("nu");
1704 double g = args.scalar<double>("g");
1705 double shockCapturingCoefficient = args.scalar<double>("shockCapturingCoefficient");
1706 xt::pyarray<int>& h_l2g = args.array<int>("h_l2g");
1707 xt::pyarray<int>& vel_l2g = args.array<int>("vel_l2g");
1708 xt::pyarray<double>& b_dof = args.array<double>("b_dof");
1709 xt::pyarray<double>& h_dof = args.array<double>("h_dof");
1710 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1711 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
1712 xt::pyarray<double>& h_dof_sge = args.array<double>("h_dof_sge");
1713 xt::pyarray<double>& u_dof_sge = args.array<double>("u_dof_sge");
1714 xt::pyarray<double>& v_dof_sge = args.array<double>("v_dof_sge");
1715 xt::pyarray<double>& q_mass_acc = args.array<double>("q_mass_acc");
1716 xt::pyarray<double>& q_mom_u_acc = args.array<double>("q_mom_u_acc");
1717 xt::pyarray<double>& q_mom_v_acc = args.array<double>("q_mom_v_acc");
1718 xt::pyarray<double>& q_mass_adv = args.array<double>("q_mass_adv");
1719 xt::pyarray<double>& q_mass_acc_beta_bdf = args.array<double>("q_mass_acc_beta_bdf");
1720 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.array<double>("q_mom_u_acc_beta_bdf");
1721 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.array<double>("q_mom_v_acc_beta_bdf");
1722 xt::pyarray<double>& q_velocity_sge = args.array<double>("q_velocity_sge");
1723 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
1724 xt::pyarray<double>& q_numDiff_h = args.array<double>("q_numDiff_h");
1725 xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1726 xt::pyarray<double>& q_numDiff_v = args.array<double>("q_numDiff_v");
1727 xt::pyarray<double>& q_numDiff_h_last = args.array<double>("q_numDiff_h_last");
1728 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1729 xt::pyarray<double>& q_numDiff_v_last = args.array<double>("q_numDiff_v_last");
1730 xt::pyarray<int>& sdInfo_u_u_rowptr = args.array<int>("sdInfo_u_u_rowptr");
1731 xt::pyarray<int>& sdInfo_u_u_colind = args.array<int>("sdInfo_u_u_colind");
1732 xt::pyarray<int>& sdInfo_u_v_rowptr = args.array<int>("sdInfo_u_v_rowptr");
1733 xt::pyarray<int>& sdInfo_u_v_colind = args.array<int>("sdInfo_u_v_colind");
1734 xt::pyarray<int>& sdInfo_v_v_rowptr = args.array<int>("sdInfo_v_v_rowptr");
1735 xt::pyarray<int>& sdInfo_v_v_colind = args.array<int>("sdInfo_v_v_colind");
1736 xt::pyarray<int>& sdInfo_v_u_rowptr = args.array<int>("sdInfo_v_u_rowptr");
1737 xt::pyarray<int>& sdInfo_v_u_colind = args.array<int>("sdInfo_v_u_colind");
1738 int offset_h = args.scalar<int>("offset_h");
1739 int offset_u = args.scalar<int>("offset_u");
1740 int offset_v = args.scalar<int>("offset_v");
1741 int stride_h = args.scalar<int>("stride_h");
1742 int stride_u = args.scalar<int>("stride_u");
1743 int stride_v = args.scalar<int>("stride_v");
1744 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1745 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1746 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1747 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1748 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1749 xt::pyarray<int>& isDOFBoundary_h = args.array<int>("isDOFBoundary_h");
1750 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1751 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
1752 xt::pyarray<int>& isAdvectiveFluxBoundary_h = args.array<int>("isAdvectiveFluxBoundary_h");
1753 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
1754 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.array<int>("isAdvectiveFluxBoundary_v");
1755 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
1756 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.array<int>("isDiffusiveFluxBoundary_v");
1757 xt::pyarray<double>& ebqe_bc_h_ext = args.array<double>("ebqe_bc_h_ext");
1758 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.array<double>("ebqe_bc_flux_mass_ext");
1759 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.array<double>("ebqe_bc_flux_mom_u_adv_ext");
1760 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.array<double>("ebqe_bc_flux_mom_v_adv_ext");
1761 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1762 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.array<double>("ebqe_bc_flux_u_diff_ext");
1763 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
1764 xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
1765 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.array<double>("ebqe_bc_flux_v_diff_ext");
1766 xt::pyarray<double>& q_velocity = args.array<double>("q_velocity");
1767 xt::pyarray<double>& ebqe_velocity = args.array<double>("ebqe_velocity");
1768 xt::pyarray<double>& flux = args.array<double>("flux");
1769 xt::pyarray<double>& elementResidual_h_save = args.array<double>("elementResidual_h");
1770 //
1771 //loop over elements to compute volume integrals and load them into element and global residual
1772 //
1773 double globalConservationError=0.0;
1774 for(int eN=0;eN<nElements_global;eN++)
1775 {
1776 //declare local storage for element residual and initialize
1777 double elementResidual_h[nDOF_test_element],
1778 elementResidual_u[nDOF_test_element],
1779 elementResidual_v[nDOF_test_element];
1780 for (int i=0;i<nDOF_test_element;i++)
1781 {
1782 int eN_i = eN*nDOF_test_element+i;
1783 elementResidual_h_save.data()[eN_i]=0.0;
1784 elementResidual_h[i]=0.0;
1785 elementResidual_u[i]=0.0;
1786 elementResidual_v[i]=0.0;
1787 }//i
1788 //
1789 //loop over quadrature points and compute integrands
1790 //
1791 for(int k=0;k<nQuadraturePoints_element;k++)
1792 {
1793 //compute indices and declare local storage
1794 int eN_k = eN*nQuadraturePoints_element+k,
1795 eN_k_nSpace = eN_k*nSpace,
1796 eN_nDOF_trial_element = eN*nDOF_trial_element;
1797 double b=0.0,h=0.0,u=0.0,v=0.0,h_sge=0.0,u_sge=0.0,v_sge=0.0,
1798 grad_b[nSpace],grad_h[nSpace],grad_u[nSpace],grad_v[nSpace],
1799 mass_acc=0.0,
1800 dmass_acc_h=0.0,
1801 mom_u_acc=0.0,
1802 dmom_u_acc_h=0.0,
1803 dmom_u_acc_u=0.0,
1804 mom_v_acc=0.0,
1805 dmom_v_acc_h=0.0,
1806 dmom_v_acc_v=0.0,
1807 mass_adv[nSpace],
1808 dmass_adv_h[nSpace],
1809 dmass_adv_u[nSpace],
1810 dmass_adv_v[nSpace],
1811 dmass_adv_h_sge[nSpace],
1812 dmass_adv_u_sge[nSpace],
1813 dmass_adv_v_sge[nSpace],
1814 mom_u_adv[nSpace],
1815 dmom_u_adv_h[nSpace],
1816 dmom_u_adv_u[nSpace],
1817 dmom_u_adv_v[nSpace],
1818 dmom_u_adv_h_sge[nSpace],
1819 dmom_u_adv_u_sge[nSpace],
1820 dmom_u_adv_v_sge[nSpace],
1821 mom_v_adv[nSpace],
1822 dmom_v_adv_h[nSpace],
1823 dmom_v_adv_u[nSpace],
1824 dmom_v_adv_v[nSpace],
1825 dmom_v_adv_h_sge[nSpace],
1826 dmom_v_adv_u_sge[nSpace],
1827 dmom_v_adv_v_sge[nSpace],
1828 mom_u_diff_ten[nSpace],
1829 mom_v_diff_ten[nSpace],
1830 mom_uv_diff_ten[1],
1831 mom_vu_diff_ten[1],
1832 mom_u_source=0.0,
1833 dmom_u_source_h=0.0,
1834 mom_v_source=0.0,
1835 dmom_v_source_h=0.0,
1836 mass_acc_t=0.0,
1837 dmass_acc_h_t=0.0,
1838 mom_u_acc_t=0.0,
1839 dmom_u_acc_h_t=0.0,
1840 dmom_u_acc_u_t=0.0,
1841 mom_v_acc_t=0.0,
1842 dmom_v_acc_h_t=0.0,
1843 dmom_v_acc_v_t=0.0,
1844 pdeResidual_h=0.0,
1845 pdeResidual_u=0.0,
1846 pdeResidual_v=0.0,
1847 //mwf switched away from usual VMS
1848 tau_x[9],tau_y[9],
1849 Lhat_x[nDOF_test_element],
1850 Lhat_y[nDOF_test_element],
1851 subgridError_hx=0.0,
1852 subgridError_ux=0.0,
1853 subgridError_vx=0.0,
1854 subgridError_hy=0.0,
1855 subgridError_uy=0.0,
1856 subgridError_vy=0.0,
1857 //
1858 jac[nSpace*nSpace],
1859 jacDet,
1860 jacInv[nSpace*nSpace],
1861 h_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
1862 h_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
1863 h_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
1864 dV,x,y,xt,yt,
1865 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
1866 //mwf added
1867 //mwf todo add as input or calculate based on 'entropy'
1868 double alpha_tau=0.25;
1869 double include_acc_in_strong_mom = 1.0;//omit accumulation terms from momentum strong residual?
1870 //get jacobian, etc for mapping reference element
1871 ck.calculateMapping_element(eN,
1872 k,
1873 mesh_dof.data(),
1874 mesh_l2g.data(),
1875 mesh_trial_ref.data(),
1876 mesh_grad_trial_ref.data(),
1877 jac,
1878 jacDet,
1879 jacInv,
1880 x,y);
1881 ck.calculateMappingVelocity_element(eN,
1882 k,
1883 mesh_velocity_dof.data(),
1884 mesh_l2g.data(),
1885 mesh_trial_ref.data(),
1886 xt,yt);
1887 //xt=0.0;yt=0.0;
1888 //std::cout<<"xt "<<xt<<'\t'<<yt<<'\t'<<std::endl;
1889 //get the physical integration weight
1890 dV = fabs(jacDet)*dV_ref.data()[k];
1891 //get the trial function gradients
1892 ck.gradTrialFromRef(&h_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,h_grad_trial);
1893 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
1894 //get the solution
1895 ck.valFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],b);
1896 ck.valFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h);
1897 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u);
1898 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v);
1899 ck.valFromDOF(h_dof_sge.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h_sge);
1900 ck.valFromDOF(u_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u_sge);
1901 ck.valFromDOF(v_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v_sge);
1902 //get the solution gradients
1903 ck.gradFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_b);
1904 ck.gradFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_h);
1905 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
1906 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
1907 //precalculate test function products with integration weights
1908 for (int j=0;j<nDOF_trial_element;j++)
1909 {
1910 h_test_dV[j] = h_test_ref.data()[k*nDOF_trial_element+j]*dV;
1911 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
1912 for (int I=0;I<nSpace;I++)
1913 {
1914 h_grad_test_dV[j*nSpace+I] = h_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1915 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1916 }
1917 }
1918 //save velocity at quadrature points for other models to use
1919 q_velocity.data()[eN_k_nSpace+0]=u;
1920 q_velocity.data()[eN_k_nSpace+1]=v;
1921 //
1922 //calculate pde coefficients at quadrature points
1923 //
1924 evaluateCoefficients(nu,
1925 g,
1926 grad_b,
1927 h,
1928 u,
1929 v,
1930 mass_acc,
1931 dmass_acc_h,
1932 mom_u_acc,
1933 dmom_u_acc_h,
1934 dmom_u_acc_u,
1935 mom_v_acc,
1936 dmom_v_acc_h,
1937 dmom_v_acc_v,
1938 mass_adv,
1939 dmass_adv_h,
1940 dmass_adv_u,
1941 dmass_adv_v,
1942 mom_u_adv,
1943 dmom_u_adv_h,
1944 dmom_u_adv_u,
1945 dmom_u_adv_v,
1946 mom_v_adv,
1947 dmom_v_adv_h,
1948 dmom_v_adv_u,
1949 dmom_v_adv_v,
1950 mom_u_diff_ten,
1951 mom_v_diff_ten,
1952 mom_uv_diff_ten,
1953 mom_vu_diff_ten,
1954 mom_u_source,
1955 dmom_u_source_h,
1956 mom_v_source,
1957 dmom_v_source_h);
1958 //
1959 //save momentum for time history and velocity for subgrid error
1960 //
1961 q_mass_acc.data()[eN_k] = mass_acc;
1962 q_mom_u_acc.data()[eN_k] = mom_u_acc;
1963 q_mom_v_acc.data()[eN_k] = mom_v_acc;
1964 //subgrid error uses grid scale discharge
1965 q_mass_adv.data()[eN_k_nSpace+0] = h*u;
1966 q_mass_adv.data()[eN_k_nSpace+1] = h*v;
1967
1968 //
1969 //calculate time derivative at quadrature points
1970 //
1971 ck.bdf(alphaBDF,
1972 q_mass_acc_beta_bdf.data()[eN_k],
1973 mass_acc,
1974 dmass_acc_h,
1975 mass_acc_t,
1976 dmass_acc_h_t);
1977 ck.bdfC2(alphaBDF,
1978 q_mom_u_acc_beta_bdf.data()[eN_k],
1979 mom_u_acc,
1980 dmom_u_acc_h,
1981 dmom_u_acc_u,
1982 mom_u_acc_t,
1983 dmom_u_acc_h_t,
1984 dmom_u_acc_u_t);
1985 ck.bdfC2(alphaBDF,
1986 q_mom_v_acc_beta_bdf.data()[eN_k],
1987 mom_v_acc,
1988 dmom_v_acc_h,
1989 dmom_v_acc_v,
1990 mom_v_acc_t,
1991 dmom_v_acc_h_t,
1992 dmom_v_acc_v_t);
1993 //
1994 //calculate subgrid error (strong residual and adjoint)
1995 //
1996 //calculate strong residual
1997 dmass_adv_h_sge[0] = dmass_adv_h[0];
1998 dmass_adv_h_sge[1] = dmass_adv_h[1];
1999 dmass_adv_u_sge[0] = dmass_adv_u[0];
2000 dmass_adv_u_sge[1] = dmass_adv_u[1];
2001 dmass_adv_v_sge[0] = dmass_adv_v[0];
2002 dmass_adv_v_sge[1] = dmass_adv_v[1];
2003 dmom_u_adv_h_sge[0] = dmom_u_adv_h[0];
2004 dmom_u_adv_h_sge[1] = dmom_u_adv_h[1];
2005 dmom_u_adv_u_sge[0] = dmom_u_adv_u[0];
2006 dmom_u_adv_u_sge[1] = dmom_u_adv_u[1];
2007 dmom_u_adv_v_sge[0] = dmom_u_adv_v[0];
2008 dmom_u_adv_v_sge[1] = dmom_u_adv_v[1];
2009 dmom_v_adv_h_sge[0] = dmom_v_adv_h[0];
2010 dmom_v_adv_h_sge[1] = dmom_v_adv_h[1];
2011 dmom_v_adv_u_sge[0] = dmom_v_adv_u[0];
2012 dmom_v_adv_u_sge[1] = dmom_v_adv_u[1];
2013 dmom_v_adv_v_sge[0] = dmom_v_adv_v[0];
2014 dmom_v_adv_v_sge[1] = dmom_v_adv_v[1];
2015
2016 //full linearization, lagged
2017
2018 //mass advective flux
2019 dmass_adv_h_sge[0]=u_sge;
2020 dmass_adv_h_sge[1]=v_sge;
2021
2022 dmass_adv_u_sge[0]=h_sge;
2023 dmass_adv_u_sge[1]=0.0;
2024
2025 dmass_adv_v_sge[0]=0.0;
2026 dmass_adv_v_sge[1]=h_sge;
2027
2028 //u momentum advective flux
2029 dmom_u_adv_h_sge[0]=u_sge*u_sge + g*h_sge;
2030 dmom_u_adv_h_sge[1]=u_sge*v_sge;
2031
2032 dmom_u_adv_u_sge[0]=h_sge*2.0*u_sge;
2033 dmom_u_adv_u_sge[1]=h_sge*v_sge;
2034
2035 dmom_u_adv_v_sge[0]=0.0;
2036 dmom_u_adv_v_sge[1]=h_sge*u_sge;
2037
2038 //v momentum advective_flux
2039 dmom_v_adv_h_sge[0]=v_sge*u_sge;
2040 dmom_v_adv_h_sge[1]=v_sge*v_sge + g*h_sge;
2041
2042 dmom_v_adv_u_sge[0]=h_sge*v_sge;
2043 dmom_v_adv_u_sge[1]=0.0;
2044
2045 dmom_v_adv_v_sge[0]=h_sge*u_sge;
2046 dmom_v_adv_v_sge[1]=h_sge*2.0*v_sge;
2047
2048 pdeResidual_h = ck.Mass_strong(mass_acc_t) +
2049 ck.Advection_strong(dmass_adv_h_sge,grad_h) +
2050 ck.Advection_strong(dmass_adv_u_sge,grad_u) +
2051 ck.Advection_strong(dmass_adv_v_sge,grad_v);
2052
2053 pdeResidual_u = include_acc_in_strong_mom*ck.Mass_strong(mom_u_acc_t) +
2054 ck.Advection_strong(dmom_u_adv_h_sge,grad_h) +
2055 ck.Advection_strong(dmom_u_adv_u_sge,grad_u) +
2056 ck.Advection_strong(dmom_u_adv_v_sge,grad_v) +
2057 ck.Reaction_strong(mom_u_source);
2058
2059 pdeResidual_v = include_acc_in_strong_mom*ck.Mass_strong(mom_v_acc_t) +
2060 ck.Advection_strong(dmom_v_adv_h_sge,grad_h) +
2061 ck.Advection_strong(dmom_v_adv_u_sge,grad_u) +
2062 ck.Advection_strong(dmom_v_adv_v_sge,grad_v) +
2063 ck.Reaction_strong(mom_v_source);
2064 //mwf switch out tau calculations to match adh supg formulation
2065 calculateSubgridError_tau_supg(elementDiameter.data()[eN],
2066 nu,
2067 g,
2068 h_sge,
2069 u_sge,
2070 v_sge,
2071 alpha_tau,
2072 dV,
2073 tau_x,
2074 tau_y,
2075 q_cfl.data()[eN_k]);
2076
2077 subgridError_hx = - tau_x[0*3+0]*pdeResidual_h - tau_x[0*3+1]*pdeResidual_u - tau_x[0*3+2]*pdeResidual_v;
2078 subgridError_ux = - tau_x[1*3+0]*pdeResidual_h - tau_x[1*3+1]*pdeResidual_u - tau_x[1*3+2]*pdeResidual_v;
2079 subgridError_vx = - tau_x[2*3+0]*pdeResidual_h - tau_x[2*3+1]*pdeResidual_u - tau_x[2*3+2]*pdeResidual_v;
2080
2081 subgridError_hy = - tau_y[0*3+0]*pdeResidual_h - tau_y[0*3+1]*pdeResidual_u - tau_y[0*3+2]*pdeResidual_v;
2082 subgridError_uy = - tau_y[1*3+0]*pdeResidual_h - tau_y[1*3+1]*pdeResidual_u - tau_y[1*3+2]*pdeResidual_v;
2083 subgridError_vy = - tau_y[2*3+0]*pdeResidual_h - tau_y[2*3+1]*pdeResidual_u - tau_y[2*3+2]*pdeResidual_v;
2084
2085 //Not really SUPG, but is test function contribution
2086 for (int i=0;i<nDOF_test_element;i++)
2087 {
2088 int i_nSpace = i*nSpace;
2089 Lhat_x[i]= -h_grad_test_dV[i_nSpace];
2090 Lhat_y[i]= -h_grad_test_dV[i_nSpace+1];
2091 }
2092 //mwf end supg tau and test function operator
2093 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v);
2094 /* double */
2095 double grad_norm[2] = {1.0,0.0};
2096 ck.calculateNumericalDiffusion(shockCapturingCoefficient,
2097 elementDiameter.data()[eN],
2098 norm_Rv,
2099 grad_norm,
2100 q_numDiff_u.data()[eN_k]);
2101
2102 /*q_numDiff_u.data()[eN_k] = 0.5*elementDiameter.data()[eN]*norm_Rv/(norm_grad+1.0e-8);*/
2103 q_numDiff_v.data()[eN_k] = q_numDiff_u.data()[eN_k];
2104
2105 /* ck.calculateNumericalDiffusion(1.0, */
2106 /* elementDiameter.data()[eN], */
2107 /* pdeResidual_h, */
2108 /* grad_h, */
2109 /* q_numDiff_h.data()[eN_k]); */
2110
2111 //update element residual
2112
2113 for(int i=0;i<nDOF_test_element;i++)
2114 {
2115 int i_nSpace=i*nSpace;
2116
2117 elementResidual_h[i] += ck.Mass_weak(mass_acc_t,h_test_dV[i]) +
2118 ck.Advection_weak(mass_adv,&h_grad_test_dV[i_nSpace]) +
2119 //mwf changed stabilization terms
2120 ck.SubgridError(subgridError_hx,Lhat_x[i]) +
2121 ck.SubgridError(subgridError_hy,Lhat_y[i]) +
2122 ck.NumericalDiffusion(q_numDiff_h_last.data()[eN_k],grad_h,&h_grad_test_dV[i_nSpace]);
2123
2124 elementResidual_u[i] += ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
2125 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
2126 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_u_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
2127 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
2128 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
2129 //mwf changed stabilization terms
2130 ck.SubgridError(subgridError_ux,Lhat_x[i]) +
2131 ck.SubgridError(subgridError_uy,Lhat_y[i]) +
2132 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]);
2133
2134 elementResidual_v[i] += ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
2135 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
2136 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_v_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
2137 ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
2138 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
2139 //mwf changed stabilization terms
2140 ck.SubgridError(subgridError_vx,Lhat_x[i]) +
2141 ck.SubgridError(subgridError_vy,Lhat_y[i]) +
2142 ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]);
2143 }
2144 }
2145
2146 //load element into global residual and save element residual
2147
2148 for(int i=0;i<nDOF_test_element;i++)
2149 {
2150 int eN_i=eN*nDOF_test_element+i;
2151
2152 elementResidual_h_save.data()[eN_i] += elementResidual_h[i];
2153
2154 globalResidual.data()[offset_h+stride_h*h_l2g.data()[eN_i]]+=elementResidual_h[i];
2155 globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
2156 globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
2157 }
2158 }
2159 }
2160 void calculateJacobian(arguments_dict& args)
2161 {
2162 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
2163 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2164 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
2165 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
2166 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
2167 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
2168 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
2169 xt::pyarray<double>& h_trial_ref = args.array<double>("h_trial_ref");
2170 xt::pyarray<double>& h_grad_trial_ref = args.array<double>("h_grad_trial_ref");
2171 xt::pyarray<double>& h_test_ref = args.array<double>("h_test_ref");
2172 xt::pyarray<double>& h_grad_test_ref = args.array<double>("h_grad_test_ref");
2173 xt::pyarray<double>& vel_trial_ref = args.array<double>("vel_trial_ref");
2174 xt::pyarray<double>& vel_grad_trial_ref = args.array<double>("vel_grad_trial_ref");
2175 xt::pyarray<double>& vel_test_ref = args.array<double>("vel_test_ref");
2176 xt::pyarray<double>& vel_grad_test_ref = args.array<double>("vel_grad_test_ref");
2177 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2178 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2179 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
2180 xt::pyarray<double>& h_trial_trace_ref = args.array<double>("h_trial_trace_ref");
2181 xt::pyarray<double>& h_grad_trial_trace_ref = args.array<double>("h_grad_trial_trace_ref");
2182 xt::pyarray<double>& h_test_trace_ref = args.array<double>("h_test_trace_ref");
2183 xt::pyarray<double>& h_grad_test_trace_ref = args.array<double>("h_grad_test_trace_ref");
2184 xt::pyarray<double>& vel_trial_trace_ref = args.array<double>("vel_trial_trace_ref");
2185 xt::pyarray<double>& vel_grad_trial_trace_ref = args.array<double>("vel_grad_trial_trace_ref");
2186 xt::pyarray<double>& vel_test_trace_ref = args.array<double>("vel_test_trace_ref");
2187 xt::pyarray<double>& vel_grad_test_trace_ref = args.array<double>("vel_grad_test_trace_ref");
2188 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
2189 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
2190 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
2191 int nElements_global = args.scalar<int>("nElements_global");
2192 double useRBLES = args.scalar<double>("useRBLES");
2193 double useMetrics = args.scalar<double>("useMetrics");
2194 double alphaBDF = args.scalar<double>("alphaBDF");
2195 double nu = args.scalar<double>("nu");
2196 double g = args.scalar<double>("g");
2197 xt::pyarray<int>& h_l2g = args.array<int>("h_l2g");
2198 xt::pyarray<int>& vel_l2g = args.array<int>("vel_l2g");
2199 xt::pyarray<double>& b_dof = args.array<double>("b_dof");
2200 xt::pyarray<double>& h_dof = args.array<double>("h_dof");
2201 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
2202 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
2203 xt::pyarray<double>& h_dof_sge = args.array<double>("h_dof_sge");
2204 xt::pyarray<double>& u_dof_sge = args.array<double>("u_dof_sge");
2205 xt::pyarray<double>& v_dof_sge = args.array<double>("v_dof_sge");
2206 xt::pyarray<double>& q_mass_acc_beta_bdf = args.array<double>("q_mass_acc_beta_bdf");
2207 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.array<double>("q_mom_u_acc_beta_bdf");
2208 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.array<double>("q_mom_v_acc_beta_bdf");
2209 xt::pyarray<double>& q_velocity_sge = args.array<double>("q_velocity_sge");
2210 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
2211 xt::pyarray<double>& q_numDiff_h_last = args.array<double>("q_numDiff_h_last");
2212 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
2213 xt::pyarray<double>& q_numDiff_v_last = args.array<double>("q_numDiff_v_last");
2214 xt::pyarray<int>& sdInfo_u_u_rowptr = args.array<int>("sdInfo_u_u_rowptr");
2215 xt::pyarray<int>& sdInfo_u_u_colind = args.array<int>("sdInfo_u_u_colind");
2216 xt::pyarray<int>& sdInfo_u_v_rowptr = args.array<int>("sdInfo_u_v_rowptr");
2217 xt::pyarray<int>& sdInfo_u_v_colind = args.array<int>("sdInfo_u_v_colind");
2218 xt::pyarray<int>& sdInfo_v_v_rowptr = args.array<int>("sdInfo_v_v_rowptr");
2219 xt::pyarray<int>& sdInfo_v_v_colind = args.array<int>("sdInfo_v_v_colind");
2220 xt::pyarray<int>& sdInfo_v_u_rowptr = args.array<int>("sdInfo_v_u_rowptr");
2221 xt::pyarray<int>& sdInfo_v_u_colind = args.array<int>("sdInfo_v_u_colind");
2222 xt::pyarray<int>& csrRowIndeces_h_h = args.array<int>("csrRowIndeces_h_h");
2223 xt::pyarray<int>& csrColumnOffsets_h_h = args.array<int>("csrColumnOffsets_h_h");
2224 xt::pyarray<int>& csrRowIndeces_h_u = args.array<int>("csrRowIndeces_h_u");
2225 xt::pyarray<int>& csrColumnOffsets_h_u = args.array<int>("csrColumnOffsets_h_u");
2226 xt::pyarray<int>& csrRowIndeces_h_v = args.array<int>("csrRowIndeces_h_v");
2227 xt::pyarray<int>& csrColumnOffsets_h_v = args.array<int>("csrColumnOffsets_h_v");
2228 xt::pyarray<int>& csrRowIndeces_u_h = args.array<int>("csrRowIndeces_u_h");
2229 xt::pyarray<int>& csrColumnOffsets_u_h = args.array<int>("csrColumnOffsets_u_h");
2230 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
2231 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
2232 xt::pyarray<int>& csrRowIndeces_u_v = args.array<int>("csrRowIndeces_u_v");
2233 xt::pyarray<int>& csrColumnOffsets_u_v = args.array<int>("csrColumnOffsets_u_v");
2234 xt::pyarray<int>& csrRowIndeces_v_h = args.array<int>("csrRowIndeces_v_h");
2235 xt::pyarray<int>& csrColumnOffsets_v_h = args.array<int>("csrColumnOffsets_v_h");
2236 xt::pyarray<int>& csrRowIndeces_v_u = args.array<int>("csrRowIndeces_v_u");
2237 xt::pyarray<int>& csrColumnOffsets_v_u = args.array<int>("csrColumnOffsets_v_u");
2238 xt::pyarray<int>& csrRowIndeces_v_v = args.array<int>("csrRowIndeces_v_v");
2239 xt::pyarray<int>& csrColumnOffsets_v_v = args.array<int>("csrColumnOffsets_v_v");
2240 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
2241 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2242 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2243 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2244 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2245 xt::pyarray<int>& isDOFBoundary_h = args.array<int>("isDOFBoundary_h");
2246 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
2247 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
2248 xt::pyarray<int>& isAdvectiveFluxBoundary_h = args.array<int>("isAdvectiveFluxBoundary_h");
2249 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
2250 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.array<int>("isAdvectiveFluxBoundary_v");
2251 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
2252 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.array<int>("isDiffusiveFluxBoundary_v");
2253 xt::pyarray<double>& ebqe_bc_h_ext = args.array<double>("ebqe_bc_h_ext");
2254 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.array<double>("ebqe_bc_flux_mass_ext");
2255 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.array<double>("ebqe_bc_flux_mom_u_adv_ext");
2256 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.array<double>("ebqe_bc_flux_mom_v_adv_ext");
2257 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
2258 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.array<double>("ebqe_bc_flux_u_diff_ext");
2259 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
2260 xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
2261 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.array<double>("ebqe_bc_flux_v_diff_ext");
2262 xt::pyarray<int>& csrColumnOffsets_eb_h_h = args.array<int>("csrColumnOffsets_eb_h_h");
2263 xt::pyarray<int>& csrColumnOffsets_eb_h_u = args.array<int>("csrColumnOffsets_eb_h_u");
2264 xt::pyarray<int>& csrColumnOffsets_eb_h_v = args.array<int>("csrColumnOffsets_eb_h_v");
2265 xt::pyarray<int>& csrColumnOffsets_eb_u_h = args.array<int>("csrColumnOffsets_eb_u_h");
2266 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
2267 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.array<int>("csrColumnOffsets_eb_u_v");
2268 xt::pyarray<int>& csrColumnOffsets_eb_v_h = args.array<int>("csrColumnOffsets_eb_v_h");
2269 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.array<int>("csrColumnOffsets_eb_v_u");
2270 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.array<int>("csrColumnOffsets_eb_v_v");
2271 double tauSum=0.0;
2272 //
2273 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2274 //
2275 for(int eN=0;eN<nElements_global;eN++)
2276 {
2277 double elementJacobian_h_h[nDOF_test_element][nDOF_trial_element],
2278 elementJacobian_h_u[nDOF_test_element][nDOF_trial_element],
2279 elementJacobian_h_v[nDOF_test_element][nDOF_trial_element],
2280 elementJacobian_u_h[nDOF_test_element][nDOF_trial_element],
2281 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
2282 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
2283 elementJacobian_v_h[nDOF_test_element][nDOF_trial_element],
2284 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
2285 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element];
2286 for (int i=0;i<nDOF_test_element;i++)
2287 for (int j=0;j<nDOF_trial_element;j++)
2288 {
2289 elementJacobian_h_h[i][j]=0.0;
2290 elementJacobian_h_u[i][j]=0.0;
2291 elementJacobian_h_v[i][j]=0.0;
2292 elementJacobian_u_h[i][j]=0.0;
2293 elementJacobian_u_u[i][j]=0.0;
2294 elementJacobian_u_v[i][j]=0.0;
2295 elementJacobian_v_h[i][j]=0.0;
2296 elementJacobian_v_u[i][j]=0.0;
2297 elementJacobian_v_v[i][j]=0.0;
2298 }
2299 for (int k=0;k<nQuadraturePoints_element;k++)
2300 {
2301 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
2302 eN_k_nSpace = eN_k*nSpace,
2303 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
2304
2305 //declare local storage
2306 double b=0.0,
2307 h=0.0,
2308 u=0.0,
2309 v=0.0,
2310 h_sge=0.0,
2311 u_sge=0.0,
2312 v_sge=0.0,
2313 grad_b[nSpace],
2314 grad_h[nSpace],
2315 grad_u[nSpace],
2316 grad_v[nSpace],
2317 mass_acc=0.0,
2318 dmass_acc_h=0.0,
2319 mom_u_acc=0.0,
2320 dmom_u_acc_h=0.0,
2321 dmom_u_acc_u=0.0,
2322 mom_v_acc=0.0,
2323 dmom_v_acc_h=0.0,
2324 dmom_v_acc_v=0.0,
2325 mass_adv[nSpace],
2326 dmass_adv_h[nSpace],
2327 dmass_adv_u[nSpace],
2328 dmass_adv_v[nSpace],
2329 dmass_adv_h_sge[nSpace],
2330 dmass_adv_u_sge[nSpace],
2331 dmass_adv_v_sge[nSpace],
2332 mom_u_adv[nSpace],
2333 dmom_u_adv_h[nSpace],
2334 dmom_u_adv_u[nSpace],
2335 dmom_u_adv_v[nSpace],
2336 dmom_u_adv_h_sge[nSpace],
2337 dmom_u_adv_u_sge[nSpace],
2338 dmom_u_adv_v_sge[nSpace],
2339 mom_v_adv[nSpace],
2340 dmom_v_adv_h[nSpace],
2341 dmom_v_adv_u[nSpace],
2342 dmom_v_adv_v[nSpace],
2343 dmom_v_adv_h_sge[nSpace],
2344 dmom_v_adv_u_sge[nSpace],
2345 dmom_v_adv_v_sge[nSpace],
2346 mom_u_diff_ten[nSpace],
2347 mom_v_diff_ten[nSpace],
2348 mom_uv_diff_ten[1],
2349 mom_vu_diff_ten[1],
2350 mom_u_source=0.0,
2351 dmom_u_source_h=0.0,
2352 mom_v_source=0.0,
2353 dmom_v_source_h=0.0,
2354 mass_acc_t=0.0,
2355 dmass_acc_h_t=0.0,
2356 mom_u_acc_t=0.0,
2357 dmom_u_acc_h_t=0.0,
2358 dmom_u_acc_u_t=0.0,
2359 mom_v_acc_t=0.0,
2360 dmom_v_acc_h_t=0.0,
2361 dmom_v_acc_v_t=0.0,
2362 tau[9],
2363 pdeResidual_h=0.0,
2364 pdeResidual_u=0.0,
2365 pdeResidual_v=0.0,
2366 dpdeResidual_h_h[nDOF_trial_element],
2367 dpdeResidual_h_u[nDOF_trial_element],
2368 dpdeResidual_h_v[nDOF_trial_element],
2369 dpdeResidual_u_h[nDOF_trial_element],
2370 dpdeResidual_u_u[nDOF_trial_element],
2371 dpdeResidual_u_v[nDOF_trial_element],
2372 dpdeResidual_v_h[nDOF_trial_element],
2373 dpdeResidual_v_u[nDOF_trial_element],
2374 dpdeResidual_v_v[nDOF_trial_element],
2375 Lstar_h_h[nDOF_test_element],
2376 Lstar_u_h[nDOF_test_element],
2377 Lstar_v_h[nDOF_test_element],
2378 Lstar_h_u[nDOF_test_element],
2379 Lstar_u_u[nDOF_test_element],
2380 Lstar_v_u[nDOF_test_element],
2381 Lstar_h_v[nDOF_test_element],
2382 Lstar_u_v[nDOF_test_element],
2383 Lstar_v_v[nDOF_test_element],
2384 subgridError_h=0.0,
2385 subgridError_u=0.0,
2386 subgridError_v=0.0,
2387 dsubgridError_h_h[nDOF_trial_element],
2388 dsubgridError_h_u[nDOF_trial_element],
2389 dsubgridError_h_v[nDOF_trial_element],
2390 dsubgridError_u_h[nDOF_trial_element],
2391 dsubgridError_u_u[nDOF_trial_element],
2392 dsubgridError_u_v[nDOF_trial_element],
2393 dsubgridError_v_h[nDOF_trial_element],
2394 dsubgridError_v_u[nDOF_trial_element],
2395 dsubgridError_v_v[nDOF_trial_element],
2396 tau_h=0.0,tau_h0=0.0,tau_h1=0.0,
2397 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
2398 jac[nSpace*nSpace],
2399 jacDet,
2400 jacInv[nSpace*nSpace],
2401 h_grad_trial[nDOF_trial_element*nSpace],
2402 vel_grad_trial[nDOF_trial_element*nSpace],
2403 dV,
2404 h_test_dV[nDOF_test_element],
2405 vel_test_dV[nDOF_test_element],
2406 h_grad_test_dV[nDOF_test_element*nSpace],
2407 vel_grad_test_dV[nDOF_test_element*nSpace],
2408 x,y,xt,yt,
2409 G[nSpace*nSpace],G_dd_G,tr_G, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
2410 //get jacobian, etc for mapping reference element
2411 ck.calculateMapping_element(eN,
2412 k,
2413 mesh_dof.data(),
2414 mesh_l2g.data(),
2415 mesh_trial_ref.data(),
2416 mesh_grad_trial_ref.data(),
2417 jac,
2418 jacDet,
2419 jacInv,
2420 x,y);
2421 ck.calculateMappingVelocity_element(eN,
2422 k,
2423 mesh_velocity_dof.data(),
2424 mesh_l2g.data(),
2425 mesh_trial_ref.data(),
2426 xt,yt);
2427 //get the physical integration weight
2428 dV = fabs(jacDet)*dV_ref.data()[k];
2429
2430 //get the trial function gradients
2431 ck.gradTrialFromRef(&h_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,h_grad_trial);
2432 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
2433 //get the solution
2434 ck.valFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],b);
2435 ck.valFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h);
2436 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u);
2437 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v);
2438 ck.valFromDOF(h_dof_sge.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h_sge);
2439 ck.valFromDOF(u_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u_sge);
2440 ck.valFromDOF(v_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v_sge);
2441 //get the solution gradients
2442 ck.gradFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_b);
2443 ck.gradFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_h);
2444 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
2445 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
2446 //precalculate test function products with integration weights
2447 for (int j=0;j<nDOF_trial_element;j++)
2448 {
2449 h_test_dV[j] = h_test_ref.data()[k*nDOF_trial_element+j]*dV;
2450 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
2451 for (int I=0;I<nSpace;I++)
2452 {
2453 h_grad_test_dV[j*nSpace+I] = h_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
2454 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin}
2455 }
2456 }
2457 evaluateCoefficients(nu,
2458 g,
2459 grad_b,
2460 h,
2461 u,
2462 v,
2463 mass_acc,
2464 dmass_acc_h,
2465 mom_u_acc,
2466 dmom_u_acc_h,
2467 dmom_u_acc_u,
2468 mom_v_acc,
2469 dmom_v_acc_h,
2470 dmom_v_acc_v,
2471 mass_adv,
2472 dmass_adv_h,
2473 dmass_adv_u,
2474 dmass_adv_v,
2475 mom_u_adv,
2476 dmom_u_adv_h,
2477 dmom_u_adv_u,
2478 dmom_u_adv_v,
2479 mom_v_adv,
2480 dmom_v_adv_h,
2481 dmom_v_adv_u,
2482 dmom_v_adv_v,
2483 mom_u_diff_ten,
2484 mom_v_diff_ten,
2485 mom_uv_diff_ten,
2486 mom_vu_diff_ten,
2487 mom_u_source,
2488 dmom_u_source_h,
2489 mom_v_source,
2490 dmom_v_source_h);
2491 //
2492 //moving mesh
2493 //
2494 /* mass_adv[0] -= MOVING_DOMAIN*mass_acc*xt; */
2495 /* mass_adv[1] -= MOVING_DOMAIN*mass_acc*yt; */
2496
2497 /* dmass_adv_h[0] -= MOVING_DOMAIN*dmass_acc_h*xt; */
2498 /* dmass_adv_h[1] -= MOVING_DOMAIN*dmass_acc_h*yt; */
2499
2500 /* mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*xt; */
2501 /* mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt; */
2502
2503 /* dmom_u_adv_h[0] -= MOVING_DOMAIN*dmom_u_acc_h*xt; */
2504 /* dmom_u_adv_h[1] -= MOVING_DOMAIN*dmom_u_acc_h*yt; */
2505
2506 /* dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*xt; */
2507 /* dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt; */
2508
2509 /* mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*xt; */
2510 /* mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt; */
2511
2512 /* dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_h*xt; */
2513 /* dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_h*yt; */
2514
2515 /* dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*xt; */
2516 /* dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt; */
2517 //
2518 //calculate time derivatives
2519 //
2520 ck.bdf(alphaBDF,
2521 q_mass_acc_beta_bdf.data()[eN_k],
2522 mass_acc,
2523 dmass_acc_h,
2524 mass_acc_t,
2525 dmass_acc_h_t);
2526 ck.bdfC2(alphaBDF,
2527 q_mom_u_acc_beta_bdf.data()[eN_k],
2528 mom_u_acc,
2529 dmom_u_acc_h,
2530 dmom_u_acc_u,
2531 mom_u_acc_t,
2532 dmom_u_acc_h_t,
2533 dmom_u_acc_u_t);
2534 ck.bdfC2(alphaBDF,
2535 q_mom_v_acc_beta_bdf.data()[eN_k],
2536 mom_v_acc,
2537 dmom_v_acc_h,
2538 dmom_v_acc_v,
2539 mom_v_acc_t,
2540 dmom_v_acc_h_t,
2541 dmom_v_acc_v_t);
2542 //
2543 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
2544 //
2545 //
2546 //calculate strong residual
2547 //
2548 dmass_adv_h_sge[0] = dmass_adv_h[0];
2549 dmass_adv_h_sge[1] = dmass_adv_h[1];
2550 dmass_adv_u_sge[0] = dmass_adv_u[0];
2551 dmass_adv_u_sge[1] = dmass_adv_u[1];
2552 dmass_adv_v_sge[0] = dmass_adv_v[0];
2553 dmass_adv_v_sge[1] = dmass_adv_v[1];
2554 dmom_u_adv_h_sge[0] = dmom_u_adv_h[0];
2555 dmom_u_adv_h_sge[1] = dmom_u_adv_h[1];
2556 dmom_u_adv_u_sge[0] = dmom_u_adv_u[0];
2557 dmom_u_adv_u_sge[1] = dmom_u_adv_u[1];
2558 dmom_u_adv_v_sge[0] = dmom_u_adv_v[0];
2559 dmom_u_adv_v_sge[1] = dmom_u_adv_v[1];
2560 dmom_v_adv_h_sge[0] = dmom_v_adv_h[0];
2561 dmom_v_adv_h_sge[1] = dmom_v_adv_h[1];
2562 dmom_v_adv_u_sge[0] = dmom_v_adv_u[0];
2563 dmom_v_adv_u_sge[1] = dmom_v_adv_u[1];
2564 dmom_v_adv_v_sge[0] = dmom_v_adv_v[0];
2565 dmom_v_adv_v_sge[1] = dmom_v_adv_v[1];
2566
2567 /* //mass advective flux */
2568 /* dmass_adv_h_sge[0]=u_sge; */
2569 /* dmass_adv_h_sge[1]=v_sge; */
2570
2571 /* dmass_adv_u_sge[0]=0.0;//h_sge; */
2572 /* dmass_adv_u_sge[1]=0.0; */
2573
2574 /* dmass_adv_v_sge[0]=0.0; */
2575 /* dmass_adv_v_sge[1]=0.0;//h_sge; */
2576
2577 /* //u momentum advective flux */
2578 /* dmom_u_adv_h_sge[0]=0.0;//u_sge*u_sge + g*h_sge; */
2579 /* dmom_u_adv_h_sge[1]=0.0;//u_sge*v_sge; */
2580
2581 /* dmom_u_adv_u_sge[0]=h_sge*u_sge;//h_sge*2.0*u_sge; */
2582 /* dmom_u_adv_u_sge[1]=h_sge*v_sge; */
2583
2584 /* dmom_u_adv_v_sge[0]=0.0; */
2585 /* dmom_u_adv_v_sge[1]=0.0;//h_sge*u_sge; */
2586
2587 /* //v momentum advective_flux */
2588 /* dmom_v_adv_h_sge[0]=0.0;//v_sge*u_sge; */
2589 /* dmom_v_adv_h_sge[1]=0.0;//v_sge*v_sge + g*h_sge; */
2590
2591 /* dmom_v_adv_u_sge[0]=0.0;//h_sge*v_sge; */
2592 /* dmom_v_adv_u_sge[1]=0.0; */
2593
2594 /* dmom_v_adv_v_sge[0]=h_sge*u_sge; */
2595 /* dmom_v_adv_v_sge[1]=h_sge*v_sge;//h_sge*2.0*v_sge; */
2596
2597 //full linearization, lagged
2598
2599 //mass advective flux
2600 dmass_adv_h_sge[0]=u_sge;
2601 dmass_adv_h_sge[1]=v_sge;
2602
2603 dmass_adv_u_sge[0]=h_sge;
2604 dmass_adv_u_sge[1]=0.0;
2605
2606 dmass_adv_v_sge[0]=0.0;
2607 dmass_adv_v_sge[1]=h_sge;
2608
2609 //u momentum advective flux
2610 dmom_u_adv_h_sge[0]=u_sge*u_sge + g*h_sge;
2611 dmom_u_adv_h_sge[1]=u_sge*v_sge;
2612
2613 dmom_u_adv_u_sge[0]=h_sge*2.0*u_sge;
2614 dmom_u_adv_u_sge[1]=h_sge*v_sge;
2615
2616 dmom_u_adv_v_sge[0]=0.0;
2617 dmom_u_adv_v_sge[1]=h_sge*u_sge;
2618
2619 //v momentum advective_flux
2620 dmom_v_adv_h_sge[0]=v_sge*u_sge;
2621 dmom_v_adv_h_sge[1]=v_sge*v_sge + g*h_sge;
2622
2623 dmom_v_adv_u_sge[0]=h_sge*v_sge;
2624 dmom_v_adv_u_sge[1]=0.0;
2625
2626 dmom_v_adv_v_sge[0]=h_sge*u_sge;
2627 dmom_v_adv_v_sge[1]=h_sge*2.0*v_sge;
2628
2629 //calculate the Jacobian of strong residual
2630 for (int j=0;j<nDOF_trial_element;j++)
2631 {
2632 int j_nSpace = j*nSpace;
2633
2634 dpdeResidual_h_h[j]= ck.MassJacobian_strong(dmass_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j])+
2635 ck.AdvectionJacobian_strong(dmass_adv_h_sge,&h_grad_trial[j_nSpace]);
2636
2637 dpdeResidual_h_u[j]=ck.AdvectionJacobian_strong(dmass_adv_u_sge,&vel_grad_trial[j_nSpace]);
2638
2639 dpdeResidual_h_v[j]=ck.AdvectionJacobian_strong(dmass_adv_v_sge,&vel_grad_trial[j_nSpace]);
2640
2641 dpdeResidual_u_h[j]= ck.MassJacobian_strong(dmom_u_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j])+
2642 ck.AdvectionJacobian_strong(dmom_u_adv_h_sge,&h_grad_trial[j_nSpace]) +
2643 ck.ReactionJacobian_strong(dmom_u_source_h,h_trial_ref.data()[k*nDOF_trial_element+j]);
2644
2645 dpdeResidual_u_u[j]= ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j])+
2646 ck.AdvectionJacobian_strong(dmom_u_adv_u_sge,&vel_grad_trial[j_nSpace]);
2647
2648 dpdeResidual_u_v[j]=ck.AdvectionJacobian_strong(dmom_u_adv_v_sge,&vel_grad_trial[j_nSpace]);
2649
2650 dpdeResidual_v_h[j]= ck.MassJacobian_strong(dmom_v_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j])+
2651 ck.AdvectionJacobian_strong(dmom_v_adv_h_sge,&h_grad_trial[j_nSpace])+
2652 ck.ReactionJacobian_strong(dmom_v_source_h,h_trial_ref.data()[k*nDOF_trial_element+j]);
2653
2654 dpdeResidual_v_u[j]=ck.AdvectionJacobian_strong(dmom_v_adv_u_sge,&vel_grad_trial[j_nSpace]);
2655
2656 dpdeResidual_v_v[j]= ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j])+
2657 ck.AdvectionJacobian_strong(dmom_v_adv_v_sge,&vel_grad_trial[j_nSpace]);
2658
2659 }
2660 /* //calculate tau and tau*Res */
2661 calculateSubgridError_tau(elementDiameter.data()[eN],
2662 nu,
2663 g,
2664 h_sge,
2665 u_sge,
2666 v_sge,
2667 tau,
2668 q_cfl.data()[eN_k]);
2669
2670 for (int i=0;i<9;i++)
2671 tauSum += tau[i];
2672
2673 for (int j=0;j<nDOF_trial_element;j++)
2674 {
2675 dsubgridError_h_h[j] = - tau[0*3+0]*dpdeResidual_h_h[j] - tau[0*3+1]*dpdeResidual_u_h[j] - tau[0*3+2]*dpdeResidual_v_h[j];
2676 dsubgridError_h_u[j] = - tau[0*3+0]*dpdeResidual_h_u[j] - tau[0*3+1]*dpdeResidual_u_u[j] - tau[0*3+2]*dpdeResidual_v_u[j];
2677 dsubgridError_h_v[j] = - tau[0*3+0]*dpdeResidual_h_v[j] - tau[0*3+1]*dpdeResidual_u_v[j] - tau[0*3+2]*dpdeResidual_v_v[j];
2678
2679 dsubgridError_u_h[j] = - tau[1*3+0]*dpdeResidual_h_h[j] - tau[1*3+1]*dpdeResidual_u_h[j] - tau[1*3+2]*dpdeResidual_v_h[j];
2680 dsubgridError_u_u[j] = - tau[1*3+0]*dpdeResidual_h_u[j] - tau[1*3+1]*dpdeResidual_u_u[j] - tau[1*3+2]*dpdeResidual_v_u[j];
2681 dsubgridError_u_v[j] = - tau[1*3+0]*dpdeResidual_h_v[j] - tau[1*3+1]*dpdeResidual_u_v[j] - tau[1*3+2]*dpdeResidual_v_v[j];
2682
2683 dsubgridError_v_h[j] = - tau[2*3+0]*dpdeResidual_h_h[j] - tau[2*3+1]*dpdeResidual_u_h[j] - tau[2*3+2]*dpdeResidual_v_h[j];
2684 dsubgridError_v_u[j] = - tau[2*3+0]*dpdeResidual_h_u[j] - tau[2*3+1]*dpdeResidual_u_u[j] - tau[2*3+2]*dpdeResidual_v_u[j];
2685 dsubgridError_v_v[j] = - tau[2*3+0]*dpdeResidual_h_v[j] - tau[2*3+1]*dpdeResidual_u_v[j] - tau[2*3+2]*dpdeResidual_v_v[j];
2686 }
2687 //adjoint times the test functions
2688 for (int i=0;i<nDOF_test_element;i++)
2689 {
2690 int i_nSpace = i*nSpace;
2691 Lstar_h_h[i]=ck.Advection_adjoint(dmass_adv_h_sge,&h_grad_test_dV[i_nSpace]);
2692 Lstar_u_h[i]=ck.Advection_adjoint(dmass_adv_u_sge,&h_grad_test_dV[i_nSpace]);
2693 Lstar_v_h[i]=ck.Advection_adjoint(dmass_adv_v_sge,&h_grad_test_dV[i_nSpace]);
2694
2695 Lstar_h_u[i]=ck.Advection_adjoint(dmom_u_adv_h_sge,&vel_grad_test_dV[i_nSpace]) +
2696 ck.Reaction_adjoint(dmom_u_source_h,vel_test_dV[i]);
2697 Lstar_u_u[i]=ck.Advection_adjoint(dmom_u_adv_u_sge,&vel_grad_test_dV[i_nSpace]);
2698 Lstar_v_u[i]=ck.Advection_adjoint(dmom_u_adv_v_sge,&vel_grad_test_dV[i_nSpace]);
2699
2700 Lstar_h_v[i]=ck.Advection_adjoint(dmom_v_adv_h_sge,&vel_grad_test_dV[i_nSpace])+
2701 ck.Reaction_adjoint(dmom_v_source_h,vel_test_dV[i]);
2702 Lstar_u_v[i]=ck.Advection_adjoint(dmom_v_adv_u_sge,&vel_grad_test_dV[i_nSpace]);
2703 Lstar_v_v[i]=ck.Advection_adjoint(dmom_v_adv_v_sge,&vel_grad_test_dV[i_nSpace]);
2704 }
2705
2706 for(int i=0;i<nDOF_test_element;i++)
2707 {
2708 int i_nSpace = i*nSpace;
2709 for(int j=0;j<nDOF_trial_element;j++)
2710 {
2711 int j_nSpace = j*nSpace;
2712 //h
2713 elementJacobian_h_h[i][j] += ck.MassJacobian_weak(dmass_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j],h_test_dV[i]) +
2714 ck.AdvectionJacobian_weak(dmass_adv_h,h_trial_ref.data()[k*nDOF_trial_element+j],&h_grad_test_dV[i_nSpace]) +
2715 ck.SubgridErrorJacobian(dsubgridError_h_h[j],Lstar_h_h[i]) +
2716 ck.SubgridErrorJacobian(dsubgridError_u_h[j],Lstar_u_h[i]) +
2717 ck.SubgridErrorJacobian(dsubgridError_v_h[j],Lstar_v_h[i]) +
2718 ck.NumericalDiffusionJacobian(q_numDiff_h_last.data()[eN_k],&h_grad_trial[j_nSpace],&h_grad_test_dV[i_nSpace]);
2719
2720 elementJacobian_h_u[i][j] += ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&h_grad_test_dV[i_nSpace]) +
2721 ck.SubgridErrorJacobian(dsubgridError_h_u[j],Lstar_h_h[i]) +
2722 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_h[i]) +
2723 ck.SubgridErrorJacobian(dsubgridError_v_u[j],Lstar_v_h[i]);
2724
2725 elementJacobian_h_v[i][j] += ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&h_grad_test_dV[i_nSpace]) +
2726 ck.SubgridErrorJacobian(dsubgridError_h_v[j],Lstar_h_h[i]) +
2727 ck.SubgridErrorJacobian(dsubgridError_u_v[j],Lstar_u_h[i]) +
2728 ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_h[i]);
2729
2730 //u
2731 elementJacobian_u_h[i][j] += ck.MassJacobian_weak(dmom_u_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
2732 ck.AdvectionJacobian_weak(dmom_u_adv_h,h_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
2733 ck.ReactionJacobian_weak(dmom_u_source_h,h_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
2734 ck.SubgridErrorJacobian(dsubgridError_h_h[j],Lstar_h_u[i]) +
2735 ck.SubgridErrorJacobian(dsubgridError_u_h[j],Lstar_u_u[i]) +
2736 ck.SubgridErrorJacobian(dsubgridError_v_h[j],Lstar_v_u[i]);
2737
2738 elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
2739 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
2740 ck.SimpleDiffusionJacobian_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_u_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
2741 ck.SubgridErrorJacobian(dsubgridError_h_u[j],Lstar_h_u[i]) +
2742 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
2743 ck.SubgridErrorJacobian(dsubgridError_v_u[j],Lstar_v_u[i]) +
2744 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
2745
2746 elementJacobian_u_v[i][j] += ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
2747 ck.SimpleDiffusionJacobian_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
2748 ck.SubgridErrorJacobian(dsubgridError_h_v[j],Lstar_h_u[i]) +
2749 ck.SubgridErrorJacobian(dsubgridError_u_v[j],Lstar_u_u[i]) +
2750 ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_u[i]);
2751
2752 //v
2753 elementJacobian_v_h[i][j] += ck.MassJacobian_weak(dmom_v_acc_h_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
2754 ck.AdvectionJacobian_weak(dmom_v_adv_h,h_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
2755 ck.ReactionJacobian_weak(dmom_v_source_h,h_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
2756 ck.SubgridErrorJacobian(dsubgridError_h_h[j],Lstar_h_v[i]) +
2757 ck.SubgridErrorJacobian(dsubgridError_u_h[j],Lstar_u_v[i]) +
2758 ck.SubgridErrorJacobian(dsubgridError_v_h[j],Lstar_v_v[i]);
2759
2760 elementJacobian_v_u[i][j] += ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
2761 ck.SimpleDiffusionJacobian_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
2762 ck.SubgridErrorJacobian(dsubgridError_h_u[j],Lstar_h_v[i]) +
2763 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_v[i]) +
2764 ck.SubgridErrorJacobian(dsubgridError_v_u[j],Lstar_v_v[i]);
2765
2766 elementJacobian_v_v[i][j] += ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
2767 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
2768 ck.SimpleDiffusionJacobian_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_v_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
2769 ck.SubgridErrorJacobian(dsubgridError_h_v[j],Lstar_h_v[i]) +
2770 ck.SubgridErrorJacobian(dsubgridError_u_v[j],Lstar_u_v[i]) +
2771 ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
2772 ck.NumericalDiffusionJacobian(q_numDiff_v_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
2773 }//j
2774 }//i
2775 }//k
2776 //
2777 //load into element Jacobian into global Jacobian
2778 //
2779 for (int i=0;i<nDOF_test_element;i++)
2780 {
2781 int eN_i = eN*nDOF_test_element+i;
2782 for (int j=0;j<nDOF_trial_element;j++)
2783 {
2784 int eN_i_j = eN_i*nDOF_trial_element+j;
2785 globalJacobian.data()[csrRowIndeces_h_h.data()[eN_i] + csrColumnOffsets_h_h.data()[eN_i_j]] += elementJacobian_h_h[i][j];
2786 globalJacobian.data()[csrRowIndeces_h_u.data()[eN_i] + csrColumnOffsets_h_u.data()[eN_i_j]] += elementJacobian_h_u[i][j];
2787 globalJacobian.data()[csrRowIndeces_h_v.data()[eN_i] + csrColumnOffsets_h_v.data()[eN_i_j]] += elementJacobian_h_v[i][j];
2788
2789 globalJacobian.data()[csrRowIndeces_u_h.data()[eN_i] + csrColumnOffsets_u_h.data()[eN_i_j]] += elementJacobian_u_h[i][j];
2790 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
2791 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_u_v.data()[eN_i_j]] += elementJacobian_u_v[i][j];
2792
2793 globalJacobian.data()[csrRowIndeces_v_h.data()[eN_i] + csrColumnOffsets_v_h.data()[eN_i_j]] += elementJacobian_v_h[i][j];
2794 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_v_u.data()[eN_i_j]] += elementJacobian_v_u[i][j];
2795 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_v_v.data()[eN_i_j]] += elementJacobian_v_v[i][j];
2796 }//j
2797 }//i
2798 }//elements
2799 std::cout<<"tauSum = "<<tauSum<<std::endl;
2800 //
2801 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
2802 //
2803 /* for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) */
2804 /* { */
2805 /* int ebN = exteriorElementBoundariesArray.data()[ebNE], */
2806 /* eN = elementBoundaryElementsArray.data()[ebN*2+0], */
2807 /* eN_nDOF_trial_element = eN*nDOF_trial_element, */
2808 /* ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0]; */
2809 /* double eps_rho,eps_mu; */
2810 /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
2811 /* { */
2812 /* int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb, */
2813 /* ebNE_kb_nSpace = ebNE_kb*nSpace, */
2814 /* ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb, */
2815 /* ebN_local_kb_nSpace = ebN_local_kb*nSpace; */
2816
2817 /* double h_ext=0.0, */
2818 /* u_ext=0.0, */
2819 /* v_ext=0.0, */
2820 /* grad_h_ext[nSpace], */
2821 /* grad_u_ext[nSpace], */
2822 /* grad_v_ext[nSpace], */
2823 /* mom_u_acc_ext=0.0, */
2824 /* dmom_u_acc_u_ext=0.0, */
2825 /* mom_v_acc_ext=0.0, */
2826 /* dmom_v_acc_v_ext=0.0, */
2827 /* mass_adv_ext[nSpace], */
2828 /* dmass_adv_u_ext[nSpace], */
2829 /* dmass_adv_v_ext[nSpace], */
2830 /* mom_u_adv_ext[nSpace], */
2831 /* dmom_u_adv_u_ext[nSpace], */
2832 /* dmom_u_adv_v_ext[nSpace], */
2833 /* mom_v_adv_ext[nSpace], */
2834 /* dmom_v_adv_u_ext[nSpace], */
2835 /* dmom_v_adv_v_ext[nSpace], */
2836 /* mom_u_diff_ten_ext[nSpace], */
2837 /* mom_v_diff_ten_ext[nSpace], */
2838 /* mom_uv_diff_ten_ext[1], */
2839 /* mom_vu_diff_ten_ext[1], */
2840 /* mom_u_source_ext=0.0, */
2841 /* mom_v_source_ext=0.0, */
2842 /* mom_u_ham_ext=0.0, */
2843 /* dmom_u_ham_grad_h_ext[nSpace], */
2844 /* mom_v_ham_ext=0.0, */
2845 /* dmom_v_ham_grad_h_ext[nSpace], */
2846 /* dmom_u_adv_h_ext[nSpace], */
2847 /* dmom_v_adv_h_ext[nSpace], */
2848 /* dflux_mass_u_ext=0.0, */
2849 /* dflux_mass_v_ext=0.0, */
2850 /* dflux_mom_u_adv_h_ext=0.0, */
2851 /* dflux_mom_u_adv_u_ext=0.0, */
2852 /* dflux_mom_u_adv_v_ext=0.0, */
2853 /* dflux_mom_v_adv_h_ext=0.0, */
2854 /* dflux_mom_v_adv_u_ext=0.0, */
2855 /* dflux_mom_v_adv_v_ext=0.0, */
2856 /* bc_h_ext=0.0, */
2857 /* bc_u_ext=0.0, */
2858 /* bc_v_ext=0.0, */
2859 /* bc_mom_u_acc_ext=0.0, */
2860 /* bc_dmom_u_acc_u_ext=0.0, */
2861 /* bc_mom_v_acc_ext=0.0, */
2862 /* bc_dmom_v_acc_v_ext=0.0, */
2863 /* bc_mass_adv_ext[nSpace], */
2864 /* bc_dmass_adv_u_ext[nSpace], */
2865 /* bc_dmass_adv_v_ext[nSpace], */
2866 /* bc_mom_u_adv_ext[nSpace], */
2867 /* bc_dmom_u_adv_u_ext[nSpace], */
2868 /* bc_dmom_u_adv_v_ext[nSpace], */
2869 /* bc_mom_v_adv_ext[nSpace], */
2870 /* bc_dmom_v_adv_u_ext[nSpace], */
2871 /* bc_dmom_v_adv_v_ext[nSpace], */
2872 /* bc_mom_u_diff_ten_ext[nSpace], */
2873 /* bc_mom_v_diff_ten_ext[nSpace], */
2874 /* bc_mom_uv_diff_ten_ext[1], */
2875 /* bc_mom_vu_diff_ten_ext[1], */
2876 /* bc_mom_u_source_ext=0.0, */
2877 /* bc_mom_v_source_ext=0.0, */
2878 /* bc_mom_u_ham_ext=0.0, */
2879 /* bc_dmom_u_ham_grad_h_ext[nSpace], */
2880 /* bc_mom_v_ham_ext=0.0, */
2881 /* bc_dmom_v_ham_grad_h_ext[nSpace], */
2882 /* fluxJacobian_h_h[nDOF_trial_element], */
2883 /* fluxJacobian_h_u[nDOF_trial_element], */
2884 /* fluxJacobian_h_v[nDOF_trial_element], */
2885 /* fluxJacobian_u_h[nDOF_trial_element], */
2886 /* fluxJacobian_u_u[nDOF_trial_element], */
2887 /* fluxJacobian_u_v[nDOF_trial_element], */
2888 /* fluxJacobian_v_h[nDOF_trial_element], */
2889 /* fluxJacobian_v_u[nDOF_trial_element], */
2890 /* fluxJacobian_v_v[nDOF_trial_element], */
2891 /* jac_ext[nSpace*nSpace], */
2892 /* jacDet_ext, */
2893 /* jacInv_ext[nSpace*nSpace], */
2894 /* boundaryJac[nSpace*(nSpace-1)], */
2895 /* metricTensor[(nSpace-1)*(nSpace-1)], */
2896 /* metricTensorDetSqrt, */
2897 /* h_grad_trial_trace[nDOF_trial_element*nSpace], */
2898 /* vel_grad_trial_trace[nDOF_trial_element*nSpace], */
2899 /* dS, */
2900 /* h_test_dS[nDOF_test_element], */
2901 /* vel_test_dS[nDOF_test_element], */
2902 /* normal[3], */
2903 /* x_ext,y_ext,xt_ext,yt_ext,integralScaling, */
2904 /* G[nSpace*nSpace],G_dd_G,tr_G,h_hhi,h_henalty; */
2905 /* ck.calculateMapping_elementBoundary(eN, */
2906 /* ebN_local, */
2907 /* kb, */
2908 /* ebN_local_kb, */
2909 /* mesh_dof.data(), */
2910 /* mesh_l2g.data(), */
2911 /* mesh_trial_trace_ref.data(), */
2912 /* mesh_grad_trial_trace_ref.data(), */
2913 /* boundaryJac_ref.data(), */
2914 /* jac_ext, */
2915 /* jacDet_ext, */
2916 /* jacInv_ext, */
2917 /* boundaryJac, */
2918 /* metricTensor, */
2919 /* metricTensorDetSqrt, */
2920 /* normal_ref.data(), */
2921 /* normal, */
2922 /* x_ext,y_ext); */
2923 /* ck.calculateMappingVelocity_elementBoundary(eN, */
2924 /* ebN_local, */
2925 /* kb, */
2926 /* ebN_local_kb, */
2927 /* mesh_velocity_dof.data(), */
2928 /* mesh_l2g.data(), */
2929 /* mesh_trial_trace_ref.data(), */
2930 /* xt_ext,yt_ext, */
2931 /* normal, */
2932 /* boundaryJac, */
2933 /* metricTensor, */
2934 /* integralScaling); */
2935 /* //xt_ext=0.0;yt_ext=0.0; */
2936 /* //std::cout<<"xt_ext "<<xt_ext<<'\t'<<yt_ext<<'\t'<<std::endl; */
2937 /* dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb]; */
2938 /* ck.calculateG(jacInv_ext,G,G_dd_G,tr_G); */
2939 /* ck.calculateGScale(G,&ebqe_normal_hhi_ext[ebNE_kb_nSpace],h_hhi); */
2940
2941 /* eps_rho = epsFact_rho*(useMetrics*h_hhi+(1.0-useMetrics)*elementDiameter.data()[eN]); */
2942 /* eps_mu = epsFact_mu *(useMetrics*h_hhi+(1.0-useMetrics)*elementDiameter.data()[eN]); */
2943
2944 /* //compute shape and solution information */
2945 /* //shape */
2946 /* ck.gradTrialFromRef(&h_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,h_grad_trial_trace); */
2947 /* ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace); */
2948 /* //solution and gradients */
2949 /* ck.valFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],h_ext); */
2950 /* ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext); */
2951 /* ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext); */
2952 /* ck.gradFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial_trace,grad_h_ext); */
2953 /* ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext); */
2954 /* ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext); */
2955 /* //precalculate test function products with integration weights */
2956 /* for (int j=0;j<nDOF_trial_element;j++) */
2957 /* { */
2958 /* h_test_dS[j] = h_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
2959 /* vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
2960 /* } */
2961 /* // // */
2962 /* // //debugging section for finite element calculations on exterior */
2963 /* // // */
2964 /* // std::cout<<"ebNE = "<<ebNE<<" kb = "<<kb<<std::endl; */
2965 /* // for (int j=0;j<nDOF_trial_element;j++) */
2966 /* // { */
2967 /* // std::cout<<"h_trial_trace["<<j<<"] "<<h_trial_trace_ref.data()[ebN_local_kb*nDOF_trial_element+j]<<std::endl; */
2968 /* // std::cout<<"vel_trial_trace["<<j<<"] "<<vel_trial_trace_ref.data()[ebN_local_kb*nDOF_trial_element+j]<<std::endl; */
2969 /* // std::cout<<"h_test_dS["<<j<<"] "<<h_test_dS[j]<<std::endl; */
2970 /* // std::cout<<"vel_test_dS["<<j<<"] "<<vel_test_dS[j]<<std::endl; */
2971 /* // for (int I=0;I<nSpace;I++) */
2972 /* // { */
2973 /* // std::cout<<"h_grad_trial_trace["<<j<<","<<I<<"] "<<h_grad_trial_trace[j*nSpace+I]<<std::endl; */
2974 /* // std::cout<<"vel_grad_trial_trace["<<j<<","<<I<<"] "<<vel_grad_trial_trace[j*nSpace+I]<<std::endl; */
2975 /* // } */
2976 /* // } */
2977 /* // std::cout<<"h_ext "<<h_ext<<std::endl; */
2978 /* // std::cout<<"u_ext "<<u_ext<<std::endl; */
2979 /* // std::cout<<"v_ext "<<v_ext<<std::endl; */
2980 /* // for(int I=0;I<nSpace;I++) */
2981 /* // { */
2982 /* // std::cout<<"grad_h_ext["<<I<<"] "<<grad_h_ext[I]<<std::endl; */
2983 /* // std::cout<<"grad_u_ext["<<I<<"] "<<grad_u_ext[I]<<std::endl; */
2984 /* // std::cout<<"grad_v_ext["<<I<<"] "<<grad_v_ext[I]<<std::endl; */
2985 /* // } */
2986 /* // */
2987 /* //load the boundary values */
2988 /* // */
2989 /* bc_h_ext = isDOFBoundary_h.data()[ebNE_kb]*ebqe_bc_h_ext.data()[ebNE_kb]+(1-isDOFBoundary_h.data()[ebNE_kb])*h_ext; */
2990 /* bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext; */
2991 /* bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext; */
2992 /* // */
2993 /* //calculate the internal and external trace of the pde coefficients */
2994 /* // */
2995 /* //cek debug */
2996 /* //eps_rho=0.1; */
2997 /* //eps_mu=0.1; */
2998 /* evaluateCoefficients(eps_rho, */
2999 /* eps_mu, */
3000 /* sigma, */
3001 /* rho_0, */
3002 /* nu_0, */
3003 /* rho_1, */
3004 /* nu_1, */
3005 /* g, */
3006 /* ebqe_hhi_ext[ebNE_kb], */
3007 /* &ebqe_normal_hhi_ext[ebNE_kb_nSpace], */
3008 /* ebqe_kappa_hhi_ext[ebNE_kb], */
3009 /* h_ext, */
3010 /* grad_h_ext, */
3011 /* u_ext, */
3012 /* v_ext, */
3013 /* mom_u_acc_ext, */
3014 /* dmom_u_acc_u_ext, */
3015 /* mom_v_acc_ext, */
3016 /* dmom_v_acc_v_ext, */
3017 /* mass_adv_ext, */
3018 /* dmass_adv_u_ext, */
3019 /* dmass_adv_v_ext, */
3020 /* mom_u_adv_ext, */
3021 /* dmom_u_adv_u_ext, */
3022 /* dmom_u_adv_v_ext, */
3023 /* mom_v_adv_ext, */
3024 /* dmom_v_adv_u_ext, */
3025 /* dmom_v_adv_v_ext, */
3026 /* mom_u_diff_ten_ext, */
3027 /* mom_v_diff_ten_ext, */
3028 /* mom_uv_diff_ten_ext, */
3029 /* mom_vu_diff_ten_ext, */
3030 /* mom_u_source_ext, */
3031 /* mom_v_source_ext, */
3032 /* mom_u_ham_ext, */
3033 /* dmom_u_ham_grad_h_ext, */
3034 /* mom_v_ham_ext, */
3035 /* dmom_v_ham_grad_h_ext); */
3036 /* evaluateCoefficients(eps_rho, */
3037 /* eps_mu, */
3038 /* sigma, */
3039 /* rho_0, */
3040 /* nu_0, */
3041 /* rho_1, */
3042 /* nu_1, */
3043 /* g, */
3044 /* ebqe_hhi_ext[ebNE_kb], */
3045 /* &ebqe_normal_hhi_ext[ebNE_kb_nSpace], */
3046 /* ebqe_kappa_hhi_ext[ebNE_kb], */
3047 /* bc_h_ext, */
3048 /* grad_h_ext, //cek shouldn't be used */
3049 /* bc_u_ext, */
3050 /* bc_v_ext, */
3051 /* bc_mom_u_acc_ext, */
3052 /* bc_dmom_u_acc_u_ext, */
3053 /* bc_mom_v_acc_ext, */
3054 /* bc_dmom_v_acc_v_ext, */
3055 /* bc_mass_adv_ext, */
3056 /* bc_dmass_adv_u_ext, */
3057 /* bc_dmass_adv_v_ext, */
3058 /* bc_mom_u_adv_ext, */
3059 /* bc_dmom_u_adv_u_ext, */
3060 /* bc_dmom_u_adv_v_ext, */
3061 /* bc_mom_v_adv_ext, */
3062 /* bc_dmom_v_adv_u_ext, */
3063 /* bc_dmom_v_adv_v_ext, */
3064 /* bc_mom_u_diff_ten_ext, */
3065 /* bc_mom_v_diff_ten_ext, */
3066 /* bc_mom_uv_diff_ten_ext, */
3067 /* bc_mom_vu_diff_ten_ext, */
3068 /* bc_mom_u_source_ext, */
3069 /* bc_mom_v_source_ext, */
3070 /* bc_mom_u_ham_ext, */
3071 /* bc_dmom_u_ham_grad_h_ext, */
3072 /* bc_mom_v_ham_ext, */
3073 /* bc_dmom_v_ham_grad_h_ext); */
3074 /* // */
3075 /* //moving domain */
3076 /* // */
3077 /* mass_adv_ext[0] -= MOVING_DOMAIN*xt_ext; */
3078 /* mass_adv_ext[1] -= MOVING_DOMAIN*yt_ext; */
3079
3080 /* mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext; */
3081 /* mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext; */
3082 /* dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext; */
3083 /* dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext; */
3084
3085 /* mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext; */
3086 /* mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext; */
3087 /* dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext; */
3088 /* dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext; */
3089
3090 /* //moving domain bc's */
3091 /* bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext; */
3092 /* bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext; */
3093
3094 /* bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext; */
3095 /* bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext; */
3096
3097 /* // */
3098 /* //calculate the numerical fluxes */
3099 /* // */
3100 /* exteriorNumericalAdvectiveFluxDerivatives(isDOFBoundary_h.data()[ebNE_kb], */
3101 /* isDOFBoundary_u.data()[ebNE_kb], */
3102 /* isDOFBoundary_v.data()[ebNE_kb], */
3103 /* isAdvectiveFluxBoundary_h.data()[ebNE_kb], */
3104 /* isAdvectiveFluxBoundary_u.data()[ebNE_kb], */
3105 /* isAdvectiveFluxBoundary_v.data()[ebNE_kb], */
3106 /* dmom_u_ham_grad_h_ext[0],//=1/rho */
3107 /* normal, */
3108 /* bc_h_ext, */
3109 /* bc_mass_adv_ext, */
3110 /* bc_mom_u_adv_ext, */
3111 /* bc_mom_v_adv_ext, */
3112 /* ebqe_bc_flux_mass_ext.data()[ebNE_kb], */
3113 /* ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb], */
3114 /* ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb], */
3115 /* h_ext, */
3116 /* mass_adv_ext, */
3117 /* mom_u_adv_ext, */
3118 /* mom_v_adv_ext, */
3119 /* dmass_adv_u_ext, */
3120 /* dmass_adv_v_ext, */
3121 /* dmom_u_adv_h_ext, */
3122 /* dmom_u_adv_u_ext, */
3123 /* dmom_u_adv_v_ext, */
3124 /* dmom_v_adv_h_ext, */
3125 /* dmom_v_adv_u_ext, */
3126 /* dmom_v_adv_v_ext, */
3127 /* dflux_mass_u_ext, */
3128 /* dflux_mass_v_ext, */
3129 /* dflux_mom_u_adv_h_ext, */
3130 /* dflux_mom_u_adv_u_ext, */
3131 /* dflux_mom_u_adv_v_ext, */
3132 /* dflux_mom_v_adv_h_ext, */
3133 /* dflux_mom_v_adv_u_ext, */
3134 /* dflux_mom_v_adv_v_ext); */
3135 /* // */
3136 /* //calculate the flux jacobian */
3137 /* // */
3138 /* ck.calculateGScale(G,normal,h_henalty); */
3139 /* h_henalty = 10.0/h_henalty; */
3140 /* //cek debug, do it the old way */
3141 /* h_henalty = 100.0/elementDiameter.data()[eN]; */
3142 /* for (int j=0;j<nDOF_trial_element;j++) */
3143 /* { */
3144 /* int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j; */
3145 /* //cek debug */
3146 /* //ebqe_henalty_ext[ebNE_kb] = 10.0; */
3147 /* // */
3148 /* //cek todo add full stress on boundaries */
3149
3150 /* fluxJacobian_h_h[j]=0.0; */
3151 /* fluxJacobian_h_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3152 /* fluxJacobian_h_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3153
3154 /* fluxJacobian_u_h[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_h_ext,h_trial_trace_ref.data()[ebN_local_kb_j]); */
3155 /* fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
3156 /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
3157 /* ebqe_hhi_ext[ebNE_kb], */
3158 /* sdInfo_u_u_rowptr.data(), */
3159 /* sdInfo_u_u_colind.data(), */
3160 /* isDOFBoundary_u.data()[ebNE_kb], */
3161 /* normal, */
3162 /* mom_u_diff_ten_ext, */
3163 /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
3164 /* &vel_grad_trial_trace[j_nSpace], */
3165 /* h_henalty);//ebqe_henalty_ext[ebNE_kb]); */
3166 /* fluxJacobian_u_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3167
3168 /* fluxJacobian_v_h[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_h_ext,h_trial_trace_ref.data()[ebN_local_kb_j]); */
3169 /* fluxJacobian_v_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3170 /* fluxJacobian_v_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
3171 /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
3172 /* ebqe_hhi_ext[ebNE_kb], */
3173 /* sdInfo_v_v_rowptr.data(), */
3174 /* sdInfo_v_v_colind.data(), */
3175 /* isDOFBoundary_v.data()[ebNE_kb], */
3176 /* normal, */
3177 /* mom_v_diff_ten_ext, */
3178 /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
3179 /* &vel_grad_trial_trace[j_nSpace], */
3180 /* h_henalty);//ebqe_henalty_ext[ebNE_kb]); */
3181
3182 /* fluxJacobian_h_h[j]=0.0; */
3183 /* fluxJacobian_h_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3184 /* fluxJacobian_h_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3185
3186 /* fluxJacobian_u_h[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_h_ext,h_trial_trace_ref.data()[ebN_local_kb_j]); */
3187 /* fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
3188 /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
3189 /* ebqe_hhi_ext[ebNE_kb], */
3190 /* sdInfo_u_u_rowptr.data(), */
3191 /* sdInfo_u_u_colind.data(), */
3192 /* isDOFBoundary_u.data()[ebNE_kb], */
3193 /* normal, */
3194 /* mom_u_diff_ten_ext, */
3195 /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
3196 /* &vel_grad_trial_trace[j_nSpace], */
3197 /* h_henalty);//ebqe_henalty_ext[ebNE_kb]); */
3198 /* fluxJacobian_u_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3199
3200 /* fluxJacobian_v_h[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_h_ext,h_trial_trace_ref.data()[ebN_local_kb_j]); */
3201 /* fluxJacobian_v_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
3202 /* fluxJacobian_v_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
3203 /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
3204 /* ebqe_hhi_ext[ebNE_kb], */
3205 /* sdInfo_v_v_rowptr.data(), */
3206 /* sdInfo_v_v_colind.data(), */
3207 /* isDOFBoundary_v.data()[ebNE_kb], */
3208 /* normal, */
3209 /* mom_v_diff_ten_ext, */
3210 /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
3211 /* &vel_grad_trial_trace[j_nSpace], */
3212 /* h_henalty);//ebqe_henalty_ext[ebNE_kb]); */
3213 /* // //cek debug */
3214 /* // fluxJacobian_h_h[j]=0.0; */
3215 /* // fluxJacobian_h_u[j]=0.0; */
3216 /* // fluxJacobian_h_v[j]=0.0; */
3217
3218 /* // fluxJacobian_u_h[j]=0.0; */
3219 /* // fluxJacobian_u_u[j]=0.0; */
3220 /* // fluxJacobian_u_v[j]=0.0; */
3221
3222 /* // fluxJacobian_v_h[j]=0.0; */
3223 /* // fluxJacobian_v_u[j]=0.0; */
3224 /* // fluxJacobian_v_v[j]=0.0; */
3225
3226 /* // //cek debug */
3227 /* }//j */
3228 /* // */
3229 /* //update the global Jacobian from the flux Jacobian */
3230 /* // */
3231 /* for (int i=0;i<nDOF_test_element;i++) */
3232 /* { */
3233 /* int eN_i = eN*nDOF_test_element+i; */
3234 /* for (int j=0;j<nDOF_trial_element;j++) */
3235 /* { */
3236 /* int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j; */
3237
3238 /* globalJacobian.data()[csrRowIndeces_h_h.data()[eN_i] + csrColumnOffsets_eb_h_h.data()[ebN_i_j]] += fluxJacobian_h_h[j]*h_test_dS[i]; */
3239 /* globalJacobian.data()[csrRowIndeces_h_u.data()[eN_i] + csrColumnOffsets_eb_h_u.data()[ebN_i_j]] += fluxJacobian_h_u[j]*h_test_dS[i]; */
3240 /* globalJacobian.data()[csrRowIndeces_h_v.data()[eN_i] + csrColumnOffsets_eb_h_v.data()[ebN_i_j]] += fluxJacobian_h_v[j]*h_test_dS[i]; */
3241
3242 /* globalJacobian.data()[csrRowIndeces_u_h.data()[eN_i] + csrColumnOffsets_eb_u_h.data()[ebN_i_j]] += fluxJacobian_u_h[j]*vel_test_dS[i]; */
3243 /* globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*vel_test_dS[i]; */
3244 /* globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_eb_u_v.data()[ebN_i_j]] += fluxJacobian_u_v[j]*vel_test_dS[i]; */
3245
3246 /* globalJacobian.data()[csrRowIndeces_v_h.data()[eN_i] + csrColumnOffsets_eb_v_h.data()[ebN_i_j]] += fluxJacobian_v_h[j]*vel_test_dS[i]; */
3247 /* globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_eb_v_u.data()[ebN_i_j]] += fluxJacobian_v_u[j]*vel_test_dS[i]; */
3248 /* globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_eb_v_v.data()[ebN_i_j]] += fluxJacobian_v_v[j]*vel_test_dS[i]; */
3249
3250 /* }//j */
3251 /* }//i */
3252 /* // //debug */
3253 /* // std::cout<<"flux jacobian ebNE "<<ebNE<<" kb "<<kb<<std::endl; */
3254 /* // for (int i=0;i<nDOF_test_element;i++) */
3255 /* // { */
3256 /* // for (int j=0;j<nDOF_trial_element;j++) */
3257 /* // { */
3258 /* // std::cout<< fluxJacobian_h_h[j]*h_test_dS[i]<<std::endl; */
3259 /* // std::cout<< fluxJacobian_h_u[j]*h_test_dS[i]<<std::endl; */
3260 /* // std::cout<< fluxJacobian_h_v[j]*h_test_dS[i]<<std::endl; */
3261
3262 /* // std::cout<< fluxJacobian_u_h[j]*vel_test_dS[i]<<std::endl; */
3263 /* // std::cout<< fluxJacobian_u_u[j]*vel_test_dS[i]<<std::endl; */
3264 /* // std::cout<< fluxJacobian_u_v[j]*vel_test_dS[i]<<std::endl; */
3265
3266 /* // std::cout<< fluxJacobian_v_h[j]*vel_test_dS[i]<<std::endl; */
3267 /* // std::cout<< fluxJacobian_v_u[j]*vel_test_dS[i]<<std::endl; */
3268 /* // std::cout<< fluxJacobian_v_v[j]*vel_test_dS[i]<<std::endl; */
3269 /* // }//j */
3270 /* // }//i */
3271 /* }//kb */
3272 /* }//ebNE */
3273 }//computeJacobian
3274 void calculateJacobian_supg(arguments_dict& args)
3275 {
3276 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
3277 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
3278 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
3279 xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
3280 double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
3281 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
3282 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
3283 xt::pyarray<double>& h_trial_ref = args.array<double>("h_trial_ref");
3284 xt::pyarray<double>& h_grad_trial_ref = args.array<double>("h_grad_trial_ref");
3285 xt::pyarray<double>& h_test_ref = args.array<double>("h_test_ref");
3286 xt::pyarray<double>& h_grad_test_ref = args.array<double>("h_grad_test_ref");
3287 xt::pyarray<double>& vel_trial_ref = args.array<double>("vel_trial_ref");
3288 xt::pyarray<double>& vel_grad_trial_ref = args.array<double>("vel_grad_trial_ref");
3289 xt::pyarray<double>& vel_test_ref = args.array<double>("vel_test_ref");
3290 xt::pyarray<double>& vel_grad_test_ref = args.array<double>("vel_grad_test_ref");
3291 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
3292 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
3293 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
3294 xt::pyarray<double>& h_trial_trace_ref = args.array<double>("h_trial_trace_ref");
3295 xt::pyarray<double>& h_grad_trial_trace_ref = args.array<double>("h_grad_trial_trace_ref");
3296 xt::pyarray<double>& h_test_trace_ref = args.array<double>("h_test_trace_ref");
3297 xt::pyarray<double>& h_grad_test_trace_ref = args.array<double>("h_grad_test_trace_ref");
3298 xt::pyarray<double>& vel_trial_trace_ref = args.array<double>("vel_trial_trace_ref");
3299 xt::pyarray<double>& vel_grad_trial_trace_ref = args.array<double>("vel_grad_trial_trace_ref");
3300 xt::pyarray<double>& vel_test_trace_ref = args.array<double>("vel_test_trace_ref");
3301 xt::pyarray<double>& vel_grad_test_trace_ref = args.array<double>("vel_grad_test_trace_ref");
3302 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
3303 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
3304 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
3305 int nElements_global = args.scalar<int>("nElements_global");
3306 double useRBLES = args.scalar<double>("useRBLES");
3307 double useMetrics = args.scalar<double>("useMetrics");
3308 double alphaBDF = args.scalar<double>("alphaBDF");
3309 double nu = args.scalar<double>("nu");
3310 double g = args.scalar<double>("g");
3311 xt::pyarray<int>& h_l2g = args.array<int>("h_l2g");
3312 xt::pyarray<int>& vel_l2g = args.array<int>("vel_l2g");
3313 xt::pyarray<double>& b_dof = args.array<double>("b_dof");
3314 xt::pyarray<double>& h_dof = args.array<double>("h_dof");
3315 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
3316 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
3317 xt::pyarray<double>& h_dof_sge = args.array<double>("h_dof_sge");
3318 xt::pyarray<double>& u_dof_sge = args.array<double>("u_dof_sge");
3319 xt::pyarray<double>& v_dof_sge = args.array<double>("v_dof_sge");
3320 xt::pyarray<double>& q_mass_acc_beta_bdf = args.array<double>("q_mass_acc_beta_bdf");
3321 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.array<double>("q_mom_u_acc_beta_bdf");
3322 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.array<double>("q_mom_v_acc_beta_bdf");
3323 xt::pyarray<double>& q_velocity_sge = args.array<double>("q_velocity_sge");
3324 xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
3325 xt::pyarray<double>& q_numDiff_h_last = args.array<double>("q_numDiff_h_last");
3326 xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
3327 xt::pyarray<double>& q_numDiff_v_last = args.array<double>("q_numDiff_v_last");
3328 xt::pyarray<int>& sdInfo_u_u_rowptr = args.array<int>("sdInfo_u_u_rowptr");
3329 xt::pyarray<int>& sdInfo_u_u_colind = args.array<int>("sdInfo_u_u_colind");
3330 xt::pyarray<int>& sdInfo_u_v_rowptr = args.array<int>("sdInfo_u_v_rowptr");
3331 xt::pyarray<int>& sdInfo_u_v_colind = args.array<int>("sdInfo_u_v_colind");
3332 xt::pyarray<int>& sdInfo_v_v_rowptr = args.array<int>("sdInfo_v_v_rowptr");
3333 xt::pyarray<int>& sdInfo_v_v_colind = args.array<int>("sdInfo_v_v_colind");
3334 xt::pyarray<int>& sdInfo_v_u_rowptr = args.array<int>("sdInfo_v_u_rowptr");
3335 xt::pyarray<int>& sdInfo_v_u_colind = args.array<int>("sdInfo_v_u_colind");
3336 xt::pyarray<int>& csrRowIndeces_h_h = args.array<int>("csrRowIndeces_h_h");
3337 xt::pyarray<int>& csrColumnOffsets_h_h = args.array<int>("csrColumnOffsets_h_h");
3338 xt::pyarray<int>& csrRowIndeces_h_u = args.array<int>("csrRowIndeces_h_u");
3339 xt::pyarray<int>& csrColumnOffsets_h_u = args.array<int>("csrColumnOffsets_h_u");
3340 xt::pyarray<int>& csrRowIndeces_h_v = args.array<int>("csrRowIndeces_h_v");
3341 xt::pyarray<int>& csrColumnOffsets_h_v = args.array<int>("csrColumnOffsets_h_v");
3342 xt::pyarray<int>& csrRowIndeces_u_h = args.array<int>("csrRowIndeces_u_h");
3343 xt::pyarray<int>& csrColumnOffsets_u_h = args.array<int>("csrColumnOffsets_u_h");
3344 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
3345 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
3346 xt::pyarray<int>& csrRowIndeces_u_v = args.array<int>("csrRowIndeces_u_v");
3347 xt::pyarray<int>& csrColumnOffsets_u_v = args.array<int>("csrColumnOffsets_u_v");
3348 xt::pyarray<int>& csrRowIndeces_v_h = args.array<int>("csrRowIndeces_v_h");
3349 xt::pyarray<int>& csrColumnOffsets_v_h = args.array<int>("csrColumnOffsets_v_h");
3350 xt::pyarray<int>& csrRowIndeces_v_u = args.array<int>("csrRowIndeces_v_u");
3351 xt::pyarray<int>& csrColumnOffsets_v_u = args.array<int>("csrColumnOffsets_v_u");
3352 xt::pyarray<int>& csrRowIndeces_v_v = args.array<int>("csrRowIndeces_v_v");
3353 xt::pyarray<int>& csrColumnOffsets_v_v = args.array<int>("csrColumnOffsets_v_v");
3354 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
3355 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
3356 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
3357 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
3358 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
3359 xt::pyarray<int>& isDOFBoundary_h = args.array<int>("isDOFBoundary_h");
3360 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
3361 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
3362 xt::pyarray<int>& isAdvectiveFluxBoundary_h = args.array<int>("isAdvectiveFluxBoundary_h");
3363 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
3364 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.array<int>("isAdvectiveFluxBoundary_v");
3365 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
3366 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.array<int>("isDiffusiveFluxBoundary_v");
3367 xt::pyarray<double>& ebqe_bc_h_ext = args.array<double>("ebqe_bc_h_ext");
3368 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.array<double>("ebqe_bc_flux_mass_ext");
3369 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.array<double>("ebqe_bc_flux_mom_u_adv_ext");
3370 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.array<double>("ebqe_bc_flux_mom_v_adv_ext");
3371 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
3372 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.array<double>("ebqe_bc_flux_u_diff_ext");
3373 xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
3374 xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
3375 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.array<double>("ebqe_bc_flux_v_diff_ext");
3376 xt::pyarray<int>& csrColumnOffsets_eb_h_h = args.array<int>("csrColumnOffsets_eb_h_h");
3377 xt::pyarray<int>& csrColumnOffsets_eb_h_u = args.array<int>("csrColumnOffsets_eb_h_u");
3378 xt::pyarray<int>& csrColumnOffsets_eb_h_v = args.array<int>("csrColumnOffsets_eb_h_v");
3379 xt::pyarray<int>& csrColumnOffsets_eb_u_h = args.array<int>("csrColumnOffsets_eb_u_h");
3380 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
3381 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.array<int>("csrColumnOffsets_eb_u_v");
3382 xt::pyarray<int>& csrColumnOffsets_eb_v_h = args.array<int>("csrColumnOffsets_eb_v_h");
3383 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.array<int>("csrColumnOffsets_eb_v_u");
3384 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.array<int>("csrColumnOffsets_eb_v_v");
3385 //
3386 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
3387 //
3388 for(int eN=0;eN<nElements_global;eN++)
3389 {
3390 double elementJacobian_h_h[nDOF_test_element][nDOF_trial_element],
3391 elementJacobian_h_u[nDOF_test_element][nDOF_trial_element],
3392 elementJacobian_h_v[nDOF_test_element][nDOF_trial_element],
3393 elementJacobian_u_h[nDOF_test_element][nDOF_trial_element],
3394 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
3395 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
3396 elementJacobian_v_h[nDOF_test_element][nDOF_trial_element],
3397 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
3398 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element];
3399 for (int i=0;i<nDOF_test_element;i++)
3400 for (int j=0;j<nDOF_trial_element;j++)
3401 {
3402 elementJacobian_h_h[i][j]=0.0;
3403 elementJacobian_h_u[i][j]=0.0;
3404 elementJacobian_h_v[i][j]=0.0;
3405 elementJacobian_u_h[i][j]=0.0;
3406 elementJacobian_u_u[i][j]=0.0;
3407 elementJacobian_u_v[i][j]=0.0;
3408 elementJacobian_v_h[i][j]=0.0;
3409 elementJacobian_v_u[i][j]=0.0;
3410 elementJacobian_v_v[i][j]=0.0;
3411 }
3412 for (int k=0;k<nQuadraturePoints_element;k++)
3413 {
3414 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
3415 eN_k_nSpace = eN_k*nSpace,
3416 eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
3417
3418 //declare local storage
3419 double b=0.0,
3420 h=0.0,
3421 u=0.0,
3422 v=0.0,
3423 h_sge=0.0,
3424 u_sge=0.0,
3425 v_sge=0.0,
3426 grad_b[nSpace],
3427 grad_h[nSpace],
3428 grad_u[nSpace],
3429 grad_v[nSpace],
3430 mass_acc=0.0,
3431 dmass_acc_h=0.0,
3432 mom_u_acc=0.0,
3433 dmom_u_acc_h=0.0,
3434 dmom_u_acc_u=0.0,
3435 mom_v_acc=0.0,
3436 dmom_v_acc_h=0.0,
3437 dmom_v_acc_v=0.0,
3438 mass_adv[nSpace],
3439 dmass_adv_h[nSpace],
3440 dmass_adv_u[nSpace],
3441 dmass_adv_v[nSpace],
3442 dmass_adv_h_sge[nSpace],
3443 dmass_adv_u_sge[nSpace],
3444 dmass_adv_v_sge[nSpace],
3445 mom_u_adv[nSpace],
3446 dmom_u_adv_h[nSpace],
3447 dmom_u_adv_u[nSpace],
3448 dmom_u_adv_v[nSpace],
3449 dmom_u_adv_h_sge[nSpace],
3450 dmom_u_adv_u_sge[nSpace],
3451 dmom_u_adv_v_sge[nSpace],
3452 mom_v_adv[nSpace],
3453 dmom_v_adv_h[nSpace],
3454 dmom_v_adv_u[nSpace],
3455 dmom_v_adv_v[nSpace],
3456 dmom_v_adv_h_sge[nSpace],
3457 dmom_v_adv_u_sge[nSpace],
3458 dmom_v_adv_v_sge[nSpace],
3459 mom_u_diff_ten[nSpace],
3460 mom_v_diff_ten[nSpace],
3461 mom_uv_diff_ten[1],
3462 mom_vu_diff_ten[1],
3463 mom_u_source=0.0,
3464 dmom_u_source_h=0.0,
3465 mom_v_source=0.0,
3466 dmom_v_source_h=0.0,
3467 mass_acc_t=0.0,
3468 dmass_acc_h_t=0.0,
3469 mom_u_acc_t=0.0,
3470 dmom_u_acc_h_t=0.0,
3471 dmom_u_acc_u_t=0.0,
3472 mom_v_acc_t=0.0,
3473 dmom_v_acc_h_t=0.0,
3474 dmom_v_acc_v_t=0.0,
3475 pdeResidual_h=0.0,
3476 pdeResidual_u=0.0,
3477 pdeResidual_v=0.0,
3478 dpdeResidual_h_h[nDOF_trial_element],
3479 dpdeResidual_h_u[nDOF_trial_element],
3480 dpdeResidual_h_v[nDOF_trial_element],
3481 dpdeResidual_u_h[nDOF_trial_element],
3482 dpdeResidual_u_u[nDOF_trial_element],
3483 dpdeResidual_u_v[nDOF_trial_element],
3484 dpdeResidual_v_h[nDOF_trial_element],
3485 dpdeResidual_v_u[nDOF_trial_element],
3486 dpdeResidual_v_v[nDOF_trial_element],
3487 //mwf switched away from usual VMS
3488 tau_x[9],tau_y[9],
3489 Lhat_x[nDOF_test_element],
3490 Lhat_y[nDOF_test_element],
3491 subgridError_hx=0.0,
3492 subgridError_ux=0.0,
3493 subgridError_vx=0.0,
3494 subgridError_hy=0.0,
3495 subgridError_uy=0.0,
3496 subgridError_vy=0.0,
3497 dsubgridError_hx_h[nDOF_trial_element],
3498 dsubgridError_hx_u[nDOF_trial_element],
3499 dsubgridError_hx_v[nDOF_trial_element],
3500 dsubgridError_ux_h[nDOF_trial_element],
3501 dsubgridError_ux_u[nDOF_trial_element],
3502 dsubgridError_ux_v[nDOF_trial_element],
3503 dsubgridError_vx_h[nDOF_trial_element],
3504 dsubgridError_vx_u[nDOF_trial_element],
3505 dsubgridError_vx_v[nDOF_trial_element],
3506 dsubgridError_hy_h[nDOF_trial_element],
3507 dsubgridError_hy_u[nDOF_trial_element],
3508 dsubgridError_hy_v[nDOF_trial_element],
3509 dsubgridError_uy_h[nDOF_trial_element],
3510 dsubgridError_uy_u[nDOF_trial_element],
3511 dsubgridError_uy_v[nDOF_trial_element],
3512 dsubgridError_vy_h[nDOF_trial_element],
3513 dsubgridError_vy_u[nDOF_trial_element],
3514 dsubgridError_vy_v[nDOF_trial_element],
3515 //mwf end changes
3516 jac[nSpace*nSpace],
3517 jacDet,
3518 jacInv[nSpace*nSpace],
3519 h_grad_trial[nDOF_trial_element*nSpace],
3520 vel_grad_trial[nDOF_trial_element*nSpace],
3521 dV,
3522 h_test_dV[nDOF_test_element],
3523 vel_test_dV[nDOF_test_element],
3524 h_grad_test_dV[nDOF_test_element*nSpace],
3525 vel_grad_test_dV[nDOF_test_element*nSpace],
3526 x,y,xt,yt,
3527 G[nSpace*nSpace],G_dd_G,tr_G, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
3528 //mwf added
3529 //mwf todo add as input or calculate based on 'entropy'
3530 double alpha_tau=0.25;
3531 double include_acc_in_strong_mom = 1.0;//omit accumulation terms from momentum strong residual?
3532 //get jacobian, etc for mapping reference element
3533 ck.calculateMapping_element(eN,
3534 k,
3535 mesh_dof.data(),
3536 mesh_l2g.data(),
3537 mesh_trial_ref.data(),
3538 mesh_grad_trial_ref.data(),
3539 jac,
3540 jacDet,
3541 jacInv,
3542 x,y);
3543 ck.calculateMappingVelocity_element(eN,
3544 k,
3545 mesh_velocity_dof.data(),
3546 mesh_l2g.data(),
3547 mesh_trial_ref.data(),
3548 xt,yt);
3549 //get the physical integration weight
3550 dV = fabs(jacDet)*dV_ref.data()[k];
3551
3552 //get the trial function gradients
3553 ck.gradTrialFromRef(&h_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,h_grad_trial);
3554 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
3555 //get the solution
3556 ck.valFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],b);
3557 ck.valFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h);
3558 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u);
3559 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v);
3560 ck.valFromDOF(h_dof_sge.data(),&h_l2g.data()[eN_nDOF_trial_element],&h_trial_ref.data()[k*nDOF_trial_element],h_sge);
3561 ck.valFromDOF(u_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u_sge);
3562 ck.valFromDOF(v_dof_sge.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v_sge);
3563 //get the solution gradients
3564 ck.gradFromDOF(b_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_b);
3565 ck.gradFromDOF(h_dof.data(),&h_l2g.data()[eN_nDOF_trial_element],h_grad_trial,grad_h);
3566 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
3567 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
3568 //precalculate test function products with integration weights
3569 for (int j=0;j<nDOF_trial_element;j++)
3570 {
3571 h_test_dV[j] = h_test_ref.data()[k*nDOF_trial_element+j]*dV;
3572 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
3573 for (int I=0;I<nSpace;I++)
3574 {
3575 h_grad_test_dV[j*nSpace+I] = h_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
3576 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin}
3577 }
3578 }
3579 evaluateCoefficients(nu,
3580 g,
3581 grad_b,
3582 h,
3583 u,
3584 v,
3585 mass_acc,
3586 dmass_acc_h,
3587 mom_u_acc,
3588 dmom_u_acc_h,
3589 dmom_u_acc_u,
3590 mom_v_acc,
3591 dmom_v_acc_h,
3592 dmom_v_acc_v,
3593 mass_adv,
3594 dmass_adv_h,
3595 dmass_adv_u,
3596 dmass_adv_v,
3597 mom_u_adv,
3598 dmom_u_adv_h,
3599 dmom_u_adv_u,
3600 dmom_u_adv_v,
3601 mom_v_adv,
3602 dmom_v_adv_h,
3603 dmom_v_adv_u,
3604 dmom_v_adv_v,
3605 mom_u_diff_ten,
3606 mom_v_diff_ten,
3607 mom_uv_diff_ten,
3608 mom_vu_diff_ten,
3609 mom_u_source,
3610 dmom_u_source_h,
3611 mom_v_source,
3612 dmom_v_source_h);
3613 //
3614 //calculate time derivatives
3615 //
3616 ck.bdf(alphaBDF,
3617 q_mass_acc_beta_bdf.data()[eN_k],
3618 mass_acc,
3619 dmass_acc_h,
3620 mass_acc_t,
3621 dmass_acc_h_t);
3622 ck.bdfC2(alphaBDF,
3623 q_mom_u_acc_beta_bdf.data()[eN_k],
3624 mom_u_acc,
3625 dmom_u_acc_h,
3626 dmom_u_acc_u,
3627 mom_u_acc_t,
3628 dmom_u_acc_h_t,
3629 dmom_u_acc_u_t);
3630 ck.bdfC2(alphaBDF,
3631 q_mom_v_acc_beta_bdf.data()[eN_k],
3632 mom_v_acc,
3633 dmom_v_acc_h,
3634 dmom_v_acc_v,
3635 mom_v_acc_t,
3636 dmom_v_acc_h_t,
3637 dmom_v_acc_v_t);
3638 //
3639 //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
3640 //
3641 //
3642 //calculate strong residual
3643 //
3644 dmass_adv_h_sge[0] = dmass_adv_h[0];
3645 dmass_adv_h_sge[1] = dmass_adv_h[1];
3646 dmass_adv_u_sge[0] = dmass_adv_u[0];
3647 dmass_adv_u_sge[1] = dmass_adv_u[1];
3648 dmass_adv_v_sge[0] = dmass_adv_v[0];
3649 dmass_adv_v_sge[1] = dmass_adv_v[1];
3650 dmom_u_adv_h_sge[0] = dmom_u_adv_h[0];
3651 dmom_u_adv_h_sge[1] = dmom_u_adv_h[1];
3652 dmom_u_adv_u_sge[0] = dmom_u_adv_u[0];
3653 dmom_u_adv_u_sge[1] = dmom_u_adv_u[1];
3654 dmom_u_adv_v_sge[0] = dmom_u_adv_v[0];
3655 dmom_u_adv_v_sge[1] = dmom_u_adv_v[1];
3656 dmom_v_adv_h_sge[0] = dmom_v_adv_h[0];
3657 dmom_v_adv_h_sge[1] = dmom_v_adv_h[1];
3658 dmom_v_adv_u_sge[0] = dmom_v_adv_u[0];
3659 dmom_v_adv_u_sge[1] = dmom_v_adv_u[1];
3660 dmom_v_adv_v_sge[0] = dmom_v_adv_v[0];
3661 dmom_v_adv_v_sge[1] = dmom_v_adv_v[1];
3662
3663 //full linearization, lagged
3664
3665 //mass advective flux
3666 dmass_adv_h_sge[0]=u_sge;
3667 dmass_adv_h_sge[1]=v_sge;
3668
3669 dmass_adv_u_sge[0]=h_sge;
3670 dmass_adv_u_sge[1]=0.0;
3671
3672 dmass_adv_v_sge[0]=0.0;
3673 dmass_adv_v_sge[1]=h_sge;
3674
3675 //u momentum advective flux
3676 dmom_u_adv_h_sge[0]=u_sge*u_sge + g*h_sge;
3677 dmom_u_adv_h_sge[1]=u_sge*v_sge;
3678
3679 dmom_u_adv_u_sge[0]=h_sge*2.0*u_sge;
3680 dmom_u_adv_u_sge[1]=h_sge*v_sge;
3681
3682 dmom_u_adv_v_sge[0]=0.0;
3683 dmom_u_adv_v_sge[1]=h_sge*u_sge;
3684
3685 //v momentum advective_flux
3686 dmom_v_adv_h_sge[0]=v_sge*u_sge;
3687 dmom_v_adv_h_sge[1]=v_sge*v_sge + g*h_sge;
3688
3689 dmom_v_adv_u_sge[0]=h_sge*v_sge;
3690 dmom_v_adv_u_sge[1]=0.0;
3691
3692 dmom_v_adv_v_sge[0]=h_sge*u_sge;
3693 dmom_v_adv_v_sge[1]=h_sge*2.0*v_sge;
3694
3695 //calculate the Jacobian of strong residual
3696 for (int j=0;j<nDOF_trial_element;j++)
3697 {
3698 int j_nSpace = j*nSpace;
3699
3700 dpdeResidual_h_h[j]= ck.MassJacobian_strong(dmass_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j])+
3701 ck.AdvectionJacobian_strong(dmass_adv_h_sge,&h_grad_trial[j_nSpace]);
3702
3703 dpdeResidual_h_u[j]=ck.AdvectionJacobian_strong(dmass_adv_u_sge,&vel_grad_trial[j_nSpace]);
3704
3705 dpdeResidual_h_v[j]=ck.AdvectionJacobian_strong(dmass_adv_v_sge,&vel_grad_trial[j_nSpace]);
3706
3707 dpdeResidual_u_h[j]= include_acc_in_strong_mom*ck.MassJacobian_strong(dmom_u_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j])+
3708 ck.AdvectionJacobian_strong(dmom_u_adv_h_sge,&h_grad_trial[j_nSpace]) +
3709 ck.ReactionJacobian_strong(dmom_u_source_h,h_trial_ref.data()[k*nDOF_trial_element+j]);
3710
3711 dpdeResidual_u_u[j]= include_acc_in_strong_mom*ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j])+
3712 ck.AdvectionJacobian_strong(dmom_u_adv_u_sge,&vel_grad_trial[j_nSpace]);
3713
3714 dpdeResidual_u_v[j]=ck.AdvectionJacobian_strong(dmom_u_adv_v_sge,&vel_grad_trial[j_nSpace]);
3715
3716 dpdeResidual_v_h[j]= include_acc_in_strong_mom*ck.MassJacobian_strong(dmom_v_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j])+
3717 ck.AdvectionJacobian_strong(dmom_v_adv_h_sge,&h_grad_trial[j_nSpace])+
3718 ck.ReactionJacobian_strong(dmom_v_source_h,h_trial_ref.data()[k*nDOF_trial_element+j]);
3719
3720 dpdeResidual_v_u[j]=ck.AdvectionJacobian_strong(dmom_v_adv_u_sge,&vel_grad_trial[j_nSpace]);
3721
3722 dpdeResidual_v_v[j]= include_acc_in_strong_mom*ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j])+
3723 ck.AdvectionJacobian_strong(dmom_v_adv_v_sge,&vel_grad_trial[j_nSpace]);
3724
3725 }
3726 /* //calculate tau and tau*Res */
3727 calculateSubgridError_tau_supg(elementDiameter.data()[eN],
3728 nu,
3729 g,
3730 h_sge,
3731 u_sge,
3732 v_sge,
3733 alpha_tau,
3734 dV,
3735 tau_x,
3736 tau_y,
3737 q_cfl.data()[eN_k]);
3738
3739 for (int j=0;j<nDOF_trial_element;j++)
3740 {
3741 dsubgridError_hx_h[j] = - tau_x[0*3+0]*dpdeResidual_h_h[j] - tau_x[0*3+1]*dpdeResidual_u_h[j] - tau_x[0*3+2]*dpdeResidual_v_h[j];
3742 dsubgridError_hx_u[j] = - tau_x[0*3+0]*dpdeResidual_h_u[j] - tau_x[0*3+1]*dpdeResidual_u_u[j] - tau_x[0*3+2]*dpdeResidual_v_u[j];
3743 dsubgridError_hx_v[j] = - tau_x[0*3+0]*dpdeResidual_h_v[j] - tau_x[0*3+1]*dpdeResidual_u_v[j] - tau_x[0*3+2]*dpdeResidual_v_v[j];
3744
3745 dsubgridError_ux_h[j] = - tau_x[1*3+0]*dpdeResidual_h_h[j] - tau_x[1*3+1]*dpdeResidual_u_h[j] - tau_x[1*3+2]*dpdeResidual_v_h[j];
3746 dsubgridError_ux_u[j] = - tau_x[1*3+0]*dpdeResidual_h_u[j] - tau_x[1*3+1]*dpdeResidual_u_u[j] - tau_x[1*3+2]*dpdeResidual_v_u[j];
3747 dsubgridError_ux_v[j] = - tau_x[1*3+0]*dpdeResidual_h_v[j] - tau_x[1*3+1]*dpdeResidual_u_v[j] - tau_x[1*3+2]*dpdeResidual_v_v[j];
3748
3749 dsubgridError_vx_h[j] = - tau_x[2*3+0]*dpdeResidual_h_h[j] - tau_x[2*3+1]*dpdeResidual_u_h[j] - tau_x[2*3+2]*dpdeResidual_v_h[j];
3750 dsubgridError_vx_u[j] = - tau_x[2*3+0]*dpdeResidual_h_u[j] - tau_x[2*3+1]*dpdeResidual_u_u[j] - tau_x[2*3+2]*dpdeResidual_v_u[j];
3751 dsubgridError_vx_v[j] = - tau_x[2*3+0]*dpdeResidual_h_v[j] - tau_x[2*3+1]*dpdeResidual_u_v[j] - tau_x[2*3+2]*dpdeResidual_v_v[j];
3752 //
3753 dsubgridError_hy_h[j] = - tau_y[0*3+0]*dpdeResidual_h_h[j] - tau_y[0*3+1]*dpdeResidual_u_h[j] - tau_y[0*3+2]*dpdeResidual_v_h[j];
3754 dsubgridError_hy_u[j] = - tau_y[0*3+0]*dpdeResidual_h_u[j] - tau_y[0*3+1]*dpdeResidual_u_u[j] - tau_y[0*3+2]*dpdeResidual_v_u[j];
3755 dsubgridError_hy_v[j] = - tau_y[0*3+0]*dpdeResidual_h_v[j] - tau_y[0*3+1]*dpdeResidual_u_v[j] - tau_y[0*3+2]*dpdeResidual_v_v[j];
3756
3757 dsubgridError_uy_h[j] = - tau_y[1*3+0]*dpdeResidual_h_h[j] - tau_y[1*3+1]*dpdeResidual_u_h[j] - tau_y[1*3+2]*dpdeResidual_v_h[j];
3758 dsubgridError_uy_u[j] = - tau_y[1*3+0]*dpdeResidual_h_u[j] - tau_y[1*3+1]*dpdeResidual_u_u[j] - tau_y[1*3+2]*dpdeResidual_v_u[j];
3759 dsubgridError_uy_v[j] = - tau_y[1*3+0]*dpdeResidual_h_v[j] - tau_y[1*3+1]*dpdeResidual_u_v[j] - tau_y[1*3+2]*dpdeResidual_v_v[j];
3760
3761 dsubgridError_vy_h[j] = - tau_y[2*3+0]*dpdeResidual_h_h[j] - tau_y[2*3+1]*dpdeResidual_u_h[j] - tau_y[2*3+2]*dpdeResidual_v_h[j];
3762 dsubgridError_vy_u[j] = - tau_y[2*3+0]*dpdeResidual_h_u[j] - tau_y[2*3+1]*dpdeResidual_u_u[j] - tau_y[2*3+2]*dpdeResidual_v_u[j];
3763 dsubgridError_vy_v[j] = - tau_y[2*3+0]*dpdeResidual_h_v[j] - tau_y[2*3+1]*dpdeResidual_u_v[j] - tau_y[2*3+2]*dpdeResidual_v_v[j];
3764
3765 }
3766 //Not really SUPG, but is test function contribution
3767 for (int i=0;i<nDOF_test_element;i++)
3768 {
3769 int i_nSpace = i*nSpace;
3770 Lhat_x[i]= -h_grad_test_dV[i_nSpace];
3771 Lhat_y[i]= -h_grad_test_dV[i_nSpace+1];
3772 }
3773
3774 for(int i=0;i<nDOF_test_element;i++)
3775 {
3776 int i_nSpace = i*nSpace;
3777 for(int j=0;j<nDOF_trial_element;j++)
3778 {
3779 int j_nSpace = j*nSpace;
3780 //h
3781 elementJacobian_h_h[i][j] += ck.MassJacobian_weak(dmass_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j],h_test_dV[i]) +
3782 ck.AdvectionJacobian_weak(dmass_adv_h,h_trial_ref.data()[k*nDOF_trial_element+j],&h_grad_test_dV[i_nSpace]) +
3783 //mwf changed stabilization terms
3784 ck.SubgridErrorJacobian(dsubgridError_hx_h[j],Lhat_x[i]) +
3785 ck.SubgridErrorJacobian(dsubgridError_hy_h[j],Lhat_y[i]) +
3786 ck.NumericalDiffusionJacobian(q_numDiff_h_last.data()[eN_k],&h_grad_trial[j_nSpace],&h_grad_test_dV[i_nSpace]);
3787
3788 elementJacobian_h_u[i][j] += ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&h_grad_test_dV[i_nSpace]) +
3789 //mwf changed stabilization terms
3790 ck.SubgridErrorJacobian(dsubgridError_hx_u[j],Lhat_x[i]) +
3791 ck.SubgridErrorJacobian(dsubgridError_hy_u[j],Lhat_y[i]);
3792
3793 elementJacobian_h_v[i][j] += ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&h_grad_test_dV[i_nSpace]) +
3794 //mwf changed stabilization terms
3795 ck.SubgridErrorJacobian(dsubgridError_hx_v[j],Lhat_x[i]) +
3796 ck.SubgridErrorJacobian(dsubgridError_hy_v[j],Lhat_y[i]);
3797
3798 //u
3799 elementJacobian_u_h[i][j] += ck.MassJacobian_weak(dmom_u_acc_h_t,h_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3800 ck.AdvectionJacobian_weak(dmom_u_adv_h,h_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3801 ck.ReactionJacobian_weak(dmom_u_source_h,h_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3802 //mwf changed stabilization terms
3803 ck.SubgridErrorJacobian(dsubgridError_ux_h[j],Lhat_x[i]) +
3804 ck.SubgridErrorJacobian(dsubgridError_uy_h[j],Lhat_y[i]);
3805
3806
3807 elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3808 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3809 ck.SimpleDiffusionJacobian_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_u_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
3810 ck.SubgridErrorJacobian(dsubgridError_ux_u[j],Lhat_x[i]) +
3811 ck.SubgridErrorJacobian(dsubgridError_uy_u[j],Lhat_y[i]) +
3812 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3813
3814 elementJacobian_u_v[i][j] += ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3815 ck.SimpleDiffusionJacobian_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
3816 ck.SubgridErrorJacobian(dsubgridError_ux_v[j],Lhat_x[i]) +
3817 ck.SubgridErrorJacobian(dsubgridError_uy_v[j],Lhat_y[i]);
3818
3819 //v
3820 elementJacobian_v_h[i][j] += ck.MassJacobian_weak(dmom_v_acc_h_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3821 ck.AdvectionJacobian_weak(dmom_v_adv_h,h_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3822 ck.ReactionJacobian_weak(dmom_v_source_h,h_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3823 ck.SubgridErrorJacobian(dsubgridError_vx_h[j],Lhat_x[i]) +
3824 ck.SubgridErrorJacobian(dsubgridError_vy_h[j],Lhat_y[i]);
3825
3826
3827 elementJacobian_v_u[i][j] += ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3828 ck.SimpleDiffusionJacobian_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
3829 ck.SubgridErrorJacobian(dsubgridError_vx_u[j],Lhat_x[i]) +
3830 ck.SubgridErrorJacobian(dsubgridError_vy_u[j],Lhat_y[i]);
3831
3832
3833 elementJacobian_v_v[i][j] += ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3834 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3835 ck.SimpleDiffusionJacobian_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_v_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
3836 ck.SubgridErrorJacobian(dsubgridError_vx_v[j],Lhat_x[i]) +
3837 ck.SubgridErrorJacobian(dsubgridError_vy_v[j],Lhat_y[i]) +
3838 ck.NumericalDiffusionJacobian(q_numDiff_v_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3839 }//j
3840 }//i
3841 }//k
3842 //
3843 //load into element Jacobian into global Jacobian
3844 //
3845 for (int i=0;i<nDOF_test_element;i++)
3846 {
3847 int eN_i = eN*nDOF_test_element+i;
3848 for (int j=0;j<nDOF_trial_element;j++)
3849 {
3850 int eN_i_j = eN_i*nDOF_trial_element+j;
3851 globalJacobian.data()[csrRowIndeces_h_h.data()[eN_i] + csrColumnOffsets_h_h.data()[eN_i_j]] += elementJacobian_h_h[i][j];
3852 globalJacobian.data()[csrRowIndeces_h_u.data()[eN_i] + csrColumnOffsets_h_u.data()[eN_i_j]] += elementJacobian_h_u[i][j];
3853 globalJacobian.data()[csrRowIndeces_h_v.data()[eN_i] + csrColumnOffsets_h_v.data()[eN_i_j]] += elementJacobian_h_v[i][j];
3854
3855 globalJacobian.data()[csrRowIndeces_u_h.data()[eN_i] + csrColumnOffsets_u_h.data()[eN_i_j]] += elementJacobian_u_h[i][j];
3856 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
3857 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_u_v.data()[eN_i_j]] += elementJacobian_u_v[i][j];
3858
3859 globalJacobian.data()[csrRowIndeces_v_h.data()[eN_i] + csrColumnOffsets_v_h.data()[eN_i_j]] += elementJacobian_v_h[i][j];
3860 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_v_u.data()[eN_i_j]] += elementJacobian_v_u[i][j];
3861 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_v_v.data()[eN_i_j]] += elementJacobian_v_v[i][j];
3862 }//j
3863 }//i
3864 }//elements
3865 //
3866 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
3867 //
3868 }//computeJacobian_supg
3869
3870 /* void calculateVelocityAverage(int nExteriorElementBoundaries_global, */
3871 /* int* exteriorElementBoundariesArray, */
3872 /* int nInteriorElementBoundaries_global, */
3873 /* int* interiorElementBoundariesArray, */
3874 /* int* elementBoundaryElementsArray, */
3875 /* int* elementBoundaryLocalElementBoundariesArray, */
3876 /* double* mesh_dof, */
3877 /* int* mesh_l2g, */
3878 /* double* mesh_trial_trace_ref, */
3879 /* double* mesh_grad_trial_trace_ref, */
3880 /* double* normal_ref, */
3881 /* double* boundaryJac_ref, */
3882 /* int* vel_l2g, */
3883 /* double* u_dof, */
3884 /* double* v_dof, */
3885 /* double* vel_trial_trace_ref, */
3886 /* double* ebqe_velocity.data(), */
3887 /* double* velocityAverage) */
3888 /* { */
3889 /* int permutations[nQuadraturePoints_elementBoundary]; */
3890 /* double xArray_left[nQuadraturePoints_elementBoundary*3], */
3891 /* xArray_right[nQuadraturePoints_elementBoundary*3]; */
3892 /* for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) */
3893 /* { */
3894 /* int ebN = exteriorElementBoundariesArray[ebNE]; */
3895 /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
3896 /* { */
3897 /* int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace, */
3898 /* ebNE_kb_nSpace = ebNE*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace; */
3899 /* velocityAverage[ebN_kb_nSpace+0]=ebqe_velocity.data()[ebNE_kb_nSpace+0]; */
3900 /* velocityAverage[ebN_kb_nSpace+1]=ebqe_velocity.data()[ebNE_kb_nSpace+1]; */
3901 /* velocityAverage[ebN_kb_nSpace+2]=ebqe_velocity.data()[ebNE_kb_nSpace+2]; */
3902 /* }//ebNE */
3903 /* } */
3904 /* for (int ebNI = 0; ebNI < nInteriorElementBoundaries_global; ebNI++) */
3905 /* { */
3906 /* int ebN = interiorElementBoundariesArray[ebNI], */
3907 /* left_eN_global = elementBoundaryElementsArray[ebN*2+0], */
3908 /* left_ebN_element = elementBoundaryLocalElementBoundariesArray[ebN*2+0], */
3909 /* right_eN_global = elementBoundaryElementsArray[ebN*2+1], */
3910 /* right_ebN_element = elementBoundaryLocalElementBoundariesArray[ebN*2+1], */
3911 /* left_eN_nDOF_trial_element = left_eN_global*nDOF_trial_element, */
3912 /* right_eN_nDOF_trial_element = right_eN_global*nDOF_trial_element; */
3913 /* double jac[nSpace*nSpace], */
3914 /* jacDet, */
3915 /* jacInv[nSpace*nSpace], */
3916 /* boundaryJac[nSpace*(nSpace-1)], */
3917 /* metricTensor[(nSpace-1)*(nSpace-1)], */
3918 /* metricTensorDetSqrt, */
3919 /* normal[3], */
3920 /* x,y; */
3921 /* //double G[nSpace*nSpace],G_dd_G,tr_G,h_hhi,h_henalty; */
3922
3923 /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
3924 /* { */
3925 /* ck.calculateMapping_elementBoundary(left_eN_global, */
3926 /* left_ebN_element, */
3927 /* kb, */
3928 /* left_ebN_element*kb, */
3929 /* mesh_dof, */
3930 /* mesh_l2g, */
3931 /* mesh_trial_trace_ref, */
3932 /* mesh_grad_trial_trace_ref, */
3933 /* boundaryJac_ref, */
3934 /* jac, */
3935 /* jacDet, */
3936 /* jacInv, */
3937 /* boundaryJac, */
3938 /* metricTensor, */
3939 /* metricTensorDetSqrt, */
3940 /* normal_ref, */
3941 /* normal, */
3942 /* x,y); */
3943 /* xArray_left[kb*3+0] = x; */
3944 /* xArray_left[kb*3+1] = y; */
3945 /* ck.calculateMapping_elementBoundary(right_eN_global, */
3946 /* right_ebN_element, */
3947 /* kb, */
3948 /* right_ebN_element*kb, */
3949 /* mesh_dof, */
3950 /* mesh_l2g, */
3951 /* mesh_trial_trace_ref, */
3952 /* mesh_grad_trial_trace_ref, */
3953 /* boundaryJac_ref, */
3954 /* jac, */
3955 /* jacDet, */
3956 /* jacInv, */
3957 /* boundaryJac, */
3958 /* metricTensor, */
3959 /* metricTensorDetSqrt, */
3960 /* normal_ref, */
3961 /* normal, */
3962 /* x,y); */
3963 /* xArray_right[kb*3+0] = x; */
3964 /* xArray_right[kb*3+1] = y; */
3965 /* } */
3966 /* for (int kb_left=0;kb_left<nQuadraturePoints_elementBoundary;kb_left++) */
3967 /* { */
3968 /* double errorNormMin = 1.0; */
3969 /* for (int kb_right=0;kb_right<nQuadraturePoints_elementBoundary;kb_right++) */
3970 /* { */
3971 /* double errorNorm=0.0; */
3972 /* for (int I=0;I<nSpace;I++) */
3973 /* { */
3974 /* errorNorm += fabs(xArray_left[kb_left*3+I] */
3975 /* - */
3976 /* xArray_right[kb_right*3+I]); */
3977 /* } */
3978 /* if (errorNorm < errorNormMin) */
3979 /* { */
3980 /* permutations[kb_right] = kb_left; */
3981 /* errorNormMin = errorNorm; */
3982 /* } */
3983 /* } */
3984 /* } */
3985 /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
3986 /* { */
3987 /* int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace; */
3988 /* double u_left=0.0, */
3989 /* v_left=0.0, */
3990 /* u_right=0.0, */
3991 /* v_right=0.0; */
3992 /* int left_kb = kb, */
3993 /* right_kb = permutations[kb], */
3994 /* left_ebN_element_kb_nDOF_test_element=left_ebN_element*left_kb*nDOF_test_element, */
3995 /* right_ebN_element_kb_nDOF_test_element=right_ebN_element*right_kb*nDOF_test_element; */
3996 /* // */
3997 /* //calculate the velocity solution at quadrature points on left and right */
3998 /* // */
3999 /* ck.valFromDOF(u_dof,&vel_l2g[left_eN_nDOF_trial_element],&vel_trial_trace_ref[left_ebN_element_kb_nDOF_test_element],u_left); */
4000 /* ck.valFromDOF(v_dof,&vel_l2g[left_eN_nDOF_trial_element],&vel_trial_trace_ref[left_ebN_element_kb_nDOF_test_element],v_left); */
4001 /* // */
4002 /* ck.valFromDOF(u_dof,&vel_l2g[right_eN_nDOF_trial_element],&vel_trial_trace_ref[right_ebN_element_kb_nDOF_test_element],u_right); */
4003 /* ck.valFromDOF(v_dof,&vel_l2g[right_eN_nDOF_trial_element],&vel_trial_trace_ref[right_ebN_element_kb_nDOF_test_element],v_right); */
4004 /* // */
4005 /* velocityAverage[ebN_kb_nSpace+0]=0.5*(u_left + u_right); */
4006 /* velocityAverage[ebN_kb_nSpace+1]=0.5*(v_left + v_right); */
4007 /* }//ebNI */
4008 /* } */
4009 /* } */
4010
4011 };//SW2D
4012
4013 inline SW2D_base* newSW2D(int nSpaceIn,
4014 int nQuadraturePoints_elementIn,
4015 int nDOF_mesh_trial_elementIn,
4016 int nDOF_trial_elementIn,
4017 int nDOF_test_elementIn,
4018 int nQuadraturePoints_elementBoundaryIn,
4019 int CompKernelFlag)
4020 {
4021 return proteus::chooseAndAllocateDiscretization2D<SW2D_base,SW2D,CompKernel>(nSpaceIn,
4022 nQuadraturePoints_elementIn,
4023 nDOF_mesh_trial_elementIn,
4024 nDOF_trial_elementIn,
4025 nDOF_test_elementIn,
4026 nQuadraturePoints_elementBoundaryIn,
4027 CompKernelFlag);
4028 }
4029}//proteus
4030#endif
n
Int n
Definition Headers.h:28
L
Double L
Definition Headers.h:72
u
Double u
Definition Headers.h:89
v
Double v
Definition Headers.h:95
proteus::SW2D_base
Definition SW2D.h:20
proteus::SW2D_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::SW2D_base::calculateJacobian_supg
virtual void calculateJacobian_supg(arguments_dict &args)=0
proteus::SW2D_base::calculateResidual_supg
virtual void calculateResidual_supg(arguments_dict &args)=0
proteus::SW2D_base::~SW2D_base
virtual ~SW2D_base()
Definition SW2D.h:22
proteus::SW2D_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::SW2D::ck
CompKernelType ck
Definition SW2D.h:40
proteus::SW2D::calculateResidual
void calculateResidual(arguments_dict &args)
Definition SW2D.h:731
proteus::SW2D::calculateResidual_supg
void calculateResidual_supg(arguments_dict &args)
Definition SW2D.h:1668
proteus::SW2D::calculateSubgridError_tau_supg
void calculateSubgridError_tau_supg(const double &elementDiameter, const double &nu, const double &g, const double &h, const double &u, const double &v, const double &alpha, const double &area, double tau_x[9], double tau_y[9], double &cfl)
Definition SW2D.h:333
proteus::SW2D::calculateJacobian_supg
void calculateJacobian_supg(arguments_dict &args)
Definition SW2D.h:3274
proteus::SW2D::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isDOFBoundary_h, const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isFluxBoundary_h, const int &isFluxBoundary_u, const int &isFluxBoundary_v, const double n[nSpace], const double &bc_h, const double bc_f_mass[nSpace], const double bc_f_umom[nSpace], const double bc_f_vmom[nSpace], const double &bc_flux_mass, const double &bc_flux_umom, const double &bc_flux_vmom, const double &h, const double f_mass[nSpace], const double f_umom[nSpace], const double f_vmom[nSpace], const double df_mass_dh[nSpace], const double df_mass_du[nSpace], const double df_mass_dv[nSpace], const double df_umom_dh[nSpace], const double df_umom_du[nSpace], const double df_umom_dv[nSpace], const double df_vmom_dh[nSpace], const double df_vmom_du[nSpace], const double df_vmom_dv[nSpace], double &flux_mass, double &flux_umom, double &flux_vmom, double *velocity)
Definition SW2D.h:385
proteus::SW2D::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition SW2D.h:39
proteus::SW2D::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &nu, const double &g, const double &h, const double &u, const double &v, double tau[9], double &cfl)
Definition SW2D.h:165
proteus::SW2D::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition SW2D.h:2160
proteus::SW2D::exteriorNumericalAdvectiveFluxDerivatives
void exteriorNumericalAdvectiveFluxDerivatives(const int &isDOFBoundary_h, const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isFluxBoundary_h, const int &isFluxBoundary_u, const int &isFluxBoundary_v, const double n[nSpace], const double &bc_h, const double bc_f_mass[nSpace], const double bc_f_umom[nSpace], const double bc_f_vmom[nSpace], const double &bc_flux_mass, const double &bc_flux_umom, const double &bc_flux_vmom, const double &h, const double f_mass[nSpace], const double f_umom[nSpace], const double f_vmom[nSpace], const double df_mass_du[nSpace], const double df_mass_dv[nSpace], const double df_umom_dh[nSpace], const double df_umom_du[nSpace], const double df_umom_dv[nSpace], const double df_vmom_dh[nSpace], const double df_vmom_du[nSpace], const double df_vmom_dv[nSpace], double &dflux_mass_dh, double &dflux_mass_du, double &dflux_mass_dv, double &dflux_umom_dh, double &dflux_umom_du, double &dflux_umom_dv, double &dflux_vmom_dh, double &dflux_vmom_du, double &dflux_vmom_dv)
Definition SW2D.h:515
proteus::SW2D::evaluateCoefficients
void evaluateCoefficients(const double nu, const double g, const double grad_b[nSpace], const double &h, const double &u, const double &v, double &mass_acc, double &dmass_acc_h, double &mom_u_acc, double &dmom_u_acc_h, double &dmom_u_acc_u, double &mom_v_acc, double &dmom_v_acc_h, double &dmom_v_acc_v, double mass_adv[nSpace], double dmass_adv_h[nSpace], double dmass_adv_u[nSpace], double dmass_adv_v[nSpace], double mom_u_adv[nSpace], double dmom_u_adv_h[nSpace], double dmom_u_adv_u[nSpace], double dmom_u_adv_v[nSpace], double mom_v_adv[nSpace], double dmom_v_adv_h[nSpace], double dmom_v_adv_u[nSpace], double dmom_v_adv_v[nSpace], double mom_u_diff_ten[nSpace], double mom_v_diff_ten[nSpace], double mom_uv_diff_ten[1], double mom_vu_diff_ten[1], double &mom_u_source, double &dmom_u_source_h, double &mom_v_source, double &dmom_v_source_h)
Definition SW2D.h:56
c
#define c(i)
Definition jf.h:21
proteus
Definition ADR.h:19
proteus::newSW2D
SW2D_base * newSW2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition SW2D.h:4013
proteus::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:309
proteus::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0