proteus  1.8.1
C/C++/Fortran libraries
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 
10 namespace 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?
17 namespace 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  register 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  register int eN_k = eN*nQuadraturePoints_element+k,
858  eN_k_nSpace = eN_k*nSpace,
859  eN_nDOF_trial_element = eN*nDOF_trial_element;
860  register 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  register 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  register 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  register 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  /* register 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  /* register 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  /* register 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  /* register 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  register 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  register int eN_k = eN*nQuadraturePoints_element+k,
1795  eN_k_nSpace = eN_k*nSpace,
1796  eN_nDOF_trial_element = eN*nDOF_trial_element;
1797  register 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  register 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  register 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  register 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  register 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  register 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  register 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  register 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  register int i_nSpace = i*nSpace;
2709  for(int j=0;j<nDOF_trial_element;j++)
2710  {
2711  register 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  register int eN_i = eN*nDOF_test_element+i;
2782  for (int j=0;j<nDOF_trial_element;j++)
2783  {
2784  register 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  /* register 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  /* register double eps_rho,eps_mu; */
2810  /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
2811  /* { */
2812  /* register 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  /* register 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  /* register 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  /* register int eN_i = eN*nDOF_test_element+i; */
3234  /* for (int j=0;j<nDOF_trial_element;j++) */
3235  /* { */
3236  /* register 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  register 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  register 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  register 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  register 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  register int i_nSpace = i*nSpace;
3777  for(int j=0;j<nDOF_trial_element;j++)
3778  {
3779  register 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  register int eN_i = eN*nDOF_test_element+i;
3848  for (int j=0;j<nDOF_trial_element;j++)
3849  {
3850  register 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  /* register int ebN = exteriorElementBoundariesArray[ebNE]; */
3895  /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
3896  /* { */
3897  /* register 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  /* register 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  /* register int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace; */
3988  /* register double u_left=0.0, */
3989  /* v_left=0.0, */
3990  /* u_right=0.0, */
3991  /* v_right=0.0; */
3992  /* register 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
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_base::calculateResidual_supg
virtual void calculateResidual_supg(arguments_dict &args)=0
proteus::SW2D::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: SW2D.h:731
proteus::SW2D::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: SW2D.h:2160
proteus::SW2D_base::calculateJacobian_supg
virtual void calculateJacobian_supg(arguments_dict &args)=0
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
L
Double L
Definition: Headers.h:72
n
Int n
Definition: Headers.h:28
proteus::SW2D_base::~SW2D_base
virtual ~SW2D_base()
Definition: SW2D.h:22
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
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::ck
CompKernelType ck
Definition: SW2D.h:40
proteus::SW2D_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
v
Double v
Definition: Headers.h:95
proteus::SW2D::calculateResidual_supg
void calculateResidual_supg(arguments_dict &args)
Definition: SW2D.h:1668
u
Double u
Definition: Headers.h:89
c
Double c
Definition: Headers.h:54
proteus::SW2D::calculateJacobian_supg
void calculateJacobian_supg(arguments_dict &args)
Definition: SW2D.h:3274
xt
Definition: AddedMass.cpp:7
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
Definition: ADR.h:17
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::arguments_dict
Definition: ArgumentsDict.h:70
proteus::SW2D::SW2D
SW2D()
Definition: SW2D.h:41
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
proteus::SW2D_base
Definition: SW2D.h:20
proteus::SW2D_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
ArgumentsDict.h
proteus::SW2D
Definition: SW2D.h:37
Generated on Fri Jul 1 2022 11:17:46 for proteus by doxygen 1.8.20