proteus  1.8.1
C/C++/Fortran libraries
RANS3PSed2D.h
Go to the documentation of this file.
1 #ifndef RANS3PSed2D_H
2 #define RANS3PSed2D_H
3 #include <cmath>
4 #include <iostream>
5 #include "CompKernel.h"
6 #include "ModelFactory.h"
7 #include "SedClosure.h"
8 #include "ArgumentsDict.h"
9 #include "xtensor-python/pyarray.hpp"
10 
11 namespace py = pybind11;
12 
13 const double DM=0.0;//1-mesh conservation and divergence, 0 - weak div(v) only
14 const double DM2=0.0;//1-point-wise mesh volume strong-residual, 0 - div(v) only
15 const double DM3=1.0;//1-point-wise divergence, 0-point-wise rate of volume change
16 #define DRAG_FAC 1.0
17 #define TURB_FORCE_FAC 1.0
18 
19 namespace proteus
20 {
21  class cppRANS3PSed2D_base
22  {
23  public:
24  virtual ~cppRANS3PSed2D_base(){}
25  virtual void setSedClosure(double aDarcy,
26  double betaForch,
27  double grain,
28  double packFraction,
29  double packMargin,
30  double maxFraction,
31  double frFraction,
32  double sigmaC,
33  double C3e,
34  double C4e,
35  double eR,
36  double fContact,
37  double mContact,
38  double nContact,
39  double angFriction,
40  double vos_limiter,
41  double mu_fr_limiter){}
42 
43  virtual void calculateResidual(arguments_dict& args)=0;
44  virtual void calculateJacobian(arguments_dict& args)=0;
45  virtual void calculateVelocityAverage(arguments_dict& args) = 0;
46  };
47 
48  template<class CompKernelType,
49  int nSpace,
50  int nQuadraturePoints_element,
51  int nDOF_mesh_trial_element,
52  int nDOF_trial_element,
53  int nDOF_test_element,
54  int nQuadraturePoints_elementBoundary>
55  class cppRANS3PSed2D : public cppRANS3PSed2D_base
56  {
57  public:
58  cppHsuSedStress<2> closure;
59  const int nDOF_test_X_trial_element;
60  CompKernelType ck;
61  cppRANS3PSed2D():
62  closure(150.0,
63  0.0,
64  0.0102,
65  0.2,
66  0.01,
67  0.635,
68  0.57,
69  1.1,
70  1.2,
71  1.0,
72  0.8,
73  0.02,
74  2.0,
75  5.0,
76  M_PI/6.,
77  0.05,
78  1.00),
79  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
80  ck()
81  {
82  }
83 
84  void setSedClosure(double aDarcy,
85  double betaForch,
86  double grain,
87  double packFraction,
88  double packMargin,
89  double maxFraction,
90  double frFraction,
91  double sigmaC,
92  double C3e,
93  double C4e,
94  double eR,
95  double fContact,
96  double mContact,
97  double nContact,
98  double angFriction,
99  double vos_limiter,
100  double mu_fr_limiter)
101  {
102  closure = cppHsuSedStress<2>(aDarcy,
103  betaForch,
104  grain,
105  packFraction,
106  packMargin,
107  maxFraction,
108  frFraction,
109  sigmaC,
110  C3e,
111  C4e,
112  eR,
113  fContact,
114  mContact,
115  nContact,
116  angFriction,
117  vos_limiter,
118  mu_fr_limiter);
119  }
120 
121  inline double smoothedHeaviside(double eps, double phi)
122  {
123  double H;
124  if (phi > eps)
125  H=1.0;
126  else if (phi < -eps)
127  H=0.0;
128  else if (phi==0.0)
129  H=0.5;
130  else
131  H = 0.5*(1.0 + phi/eps + sin(M_PI*phi/eps)/M_PI);
132  return H;
133  }
134 
135  inline double smoothedHeaviside_integral(double eps, double phi)
136  {
137  double HI;
138  if (phi > eps)
139  {
140  HI= phi - eps \
141  + 0.5*(eps + 0.5*eps*eps/eps - eps*cos(M_PI*eps/eps)/(M_PI*M_PI)) \
142  - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
143  }
144  else if (phi < -eps)
145  {
146  HI=0.0;
147  }
148  else
149  {
150  HI = 0.5*(phi + 0.5*phi*phi/eps - eps*cos(M_PI*phi/eps)/(M_PI*M_PI)) \
151  - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
152  }
153  return HI;
154  }
155 
156  inline double smoothedDirac(double eps, double phi)
157  {
158  double d;
159  if (phi > eps)
160  d=0.0;
161  else if (phi < -eps)
162  d=0.0;
163  else
164  d = 0.5*(1.0 + cos(M_PI*phi/eps))/eps;
165  return d;
166  }
167 
168  inline
169  void evaluateCoefficients(const double eps_rho,
170  const double eps_mu,
171  const double sigma,
172  const double rho_0,
173  double nu_0,
174  const double rho_1,
175  double nu_1,
176  const double rho_s,
177  const double h_e,
178  const double smagorinskyConstant,
179  const int turbulenceClosureModel,
180  const double g[nSpace],
181  const double useVF,
182  const double& vf,
183  const double& phi,
184  const double n[nSpace],
185  const double& kappa,
186  const double vos,//VRANS specific
187  const double& p,
188  const double grad_p[nSpace],
189  const double grad_u[nSpace],
190  const double grad_v[nSpace],
191  const double grad_w[nSpace],
192  const double& u,
193  const double& v,
194  const double& w,
195  const double& uStar,
196  const double& vStar,
197  const double& wStar,
198  double& eddy_viscosity,
199  double& mom_u_acc,
200  double& dmom_u_acc_u,
201  double& mom_v_acc,
202  double& dmom_v_acc_v,
203  double& mom_w_acc,
204  double& dmom_w_acc_w,
205  double mass_adv[nSpace],
206  double dmass_adv_u[nSpace],
207  double dmass_adv_v[nSpace],
208  double dmass_adv_w[nSpace],
209  double mom_u_adv[nSpace],
210  double dmom_u_adv_u[nSpace],
211  double dmom_u_adv_v[nSpace],
212  double dmom_u_adv_w[nSpace],
213  double mom_v_adv[nSpace],
214  double dmom_v_adv_u[nSpace],
215  double dmom_v_adv_v[nSpace],
216  double dmom_v_adv_w[nSpace],
217  double mom_w_adv[nSpace],
218  double dmom_w_adv_u[nSpace],
219  double dmom_w_adv_v[nSpace],
220  double dmom_w_adv_w[nSpace],
221  double mom_uu_diff_ten[nSpace],
222  double mom_vv_diff_ten[nSpace],
223  double mom_ww_diff_ten[nSpace],
224  double mom_uv_diff_ten[1],
225  double mom_uw_diff_ten[1],
226  double mom_vu_diff_ten[1],
227  double mom_vw_diff_ten[1],
228  double mom_wu_diff_ten[1],
229  double mom_wv_diff_ten[1],
230  double& mom_u_source,
231  double& mom_v_source,
232  double& mom_w_source,
233  double& mom_u_ham,
234  double dmom_u_ham_grad_p[nSpace],
235  double dmom_u_ham_grad_u[nSpace],
236  double& mom_v_ham,
237  double dmom_v_ham_grad_p[nSpace],
238  double dmom_v_ham_grad_v[nSpace],
239  double& mom_w_ham,
240  double dmom_w_ham_grad_p[nSpace],
241  double dmom_w_ham_grad_w[nSpace])
242  {
243  double rho,rho_f,H_rho;
244  H_rho = (1.0-useVF)*smoothedHeaviside(eps_rho,phi) + useVF*fmin(1.0,fmax(0.0,vf));
245  // H_mu = (1.0-useVF)*smoothedHeaviside(eps_mu,phi) + useVF*fmin(1.0,fmax(0.0,vf));
246  // d_mu = (1.0-useVF)*smoothedDirac(eps_mu,phi);
247 
248  //calculate eddy viscosity - NO NEED FOR SOLID PHASE Coefficients
249  /* switch (turbulenceClosureModel)
250  {
251  double norm_S;
252  case 1:
253  {
254  norm_S = sqrt(2.0*(grad_u[0]*grad_u[0] + grad_v[1]*grad_v[1] +
255  0.5*(grad_u[1]+grad_v[0])*(grad_u[1]+grad_v[0])));
256  nu_t0 = smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
257  nu_t1 = smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
258  }
259  case 2:
260  {
261  double re_0,cs_0=0.0,re_1,cs_1=0.0;
262  norm_S = sqrt(2.0*(grad_u[0]*grad_u[0] + grad_v[1]*grad_v[1] +
263  0.5*(grad_u[1]+grad_v[0])*(grad_u[1]+grad_v[0])));
264  re_0 = h_e*h_e*norm_S/nu_0;
265  if (re_0 > 1.0)
266  cs_0=0.027*pow(10.0,-3.23*pow(re_0,-0.92));
267  nu_t0 = cs_0*h_e*h_e*norm_S;
268  re_1 = h_e*h_e*norm_S/nu_1;
269  if (re_1 > 1.0)
270  cs_1=0.027*pow(10.0,-3.23*pow(re_1,-0.92));
271  nu_t1 = cs_1*h_e*h_e*norm_S;
272  }
273  }*/
274 
275  rho = rho_s;
276  rho_f = rho_0*(1.0-H_rho)+rho_1*H_rho;
277  // nu_t= nu_t0*(1.0-H_mu)+nu_t1*H_mu;
278  // nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
279  //nu = nu_0*(1.0-H_mu)+nu_1*H_mu;
280  // nu += nu_t;
281  //mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
282  // mu = rho_0*nu_0*(1.0-H_mu)+rho_1*nu_1*H_mu;
283 
284  // eddy_viscosity = rho_0*nu_t0*(1.0-H_mu)+rho_1*nu_t1*H_mu;
285 
286  // mass (volume accumulation)
287  //..hardwired
288 
289  //u momentum accumulation
290  mom_u_acc=u;//trick for non-conservative form
291  dmom_u_acc_u=rho;
292 
293  //v momentum accumulation
294  mom_v_acc=v;
295  dmom_v_acc_v=rho;
296 
297  /* //w momentum accumulation */
298  /* mom_w_acc=w; */
299  /* dmom_w_acc_w=rho*vos; */
300 
301  //mass advective flux
302  mass_adv[0]=vos*u;
303  mass_adv[1]=vos*v;
304  /* mass_adv[2]=qw; */
305 
306  dmass_adv_u[0]=vos;
307  dmass_adv_u[1]=0.0;
308  /* dmass_adv_u[2]=0.0; */
309 
310  dmass_adv_v[0]=0.0;
311  dmass_adv_v[1]=vos;
312  /* dmass_adv_v[2]=0.0; */
313 
314  /* dmass_adv_w[0]=0.0; */
315  /* dmass_adv_w[1]=0.0; */
316  /* dmass_adv_w[2]=1.0; */
317 
318  //advection switched to non-conservative form but could be used for mesh motion...
319  //u momentum advective flux
320  mom_u_adv[0]=0.0;
321  mom_u_adv[1]=0.0;
322  /* mom_u_adv[2]=0.0; */
323 
324  dmom_u_adv_u[0]=0.0;
325  dmom_u_adv_u[1]=0.0;
326  /* dmom_u_adv_u[2]=0.0; */
327 
328  dmom_u_adv_v[0]=0.0;
329  dmom_u_adv_v[1]=0.0;
330  /* dmom_u_adv_v[2]=0.0; */
331 
332  /* dmom_u_adv_w[0]=0.0; */
333  /* dmom_u_adv_w[1]=0.0; */
334  /* dmom_u_adv_w[2]=0.0; */
335 
336  //v momentum advective_flux
337  mom_v_adv[0]=0.0;
338  mom_v_adv[1]=0.0;
339  /* mom_v_adv[2]=0.0; */
340 
341  dmom_v_adv_u[0]=0.0;
342  dmom_v_adv_u[1]=0.0;
343  /* dmom_v_adv_u[2]=0.0; */
344 
345  /* dmom_v_adv_w[0]=0.0; */
346  /* dmom_v_adv_w[1]=0.0; */
347  /* dmom_v_adv_w[2]=0.0; */
348 
349  dmom_v_adv_v[0]=0.0;
350  dmom_v_adv_v[1]=0.0;
351  /* dmom_v_adv_v[2]=0.0; */
352 
353  /* //w momentum advective_flux */
354  /* mom_w_adv[0]=0.0; */
355  /* mom_w_adv[1]=0.0; */
356  /* mom_w_adv[2]=0.0; */
357 
358  /* dmom_w_adv_u[0]=0.0; */
359  /* dmom_w_adv_u[1]=0.0; */
360  /* dmom_w_adv_u[2]=0.0; */
361 
362  /* dmom_w_adv_v[0]=0.0; */
363  /* dmom_w_adv_v[1]=0.0; */
364  /* dmom_w_adv_v[2]=0.0; */
365 
366  /* dmom_w_adv_w[0]=0.0; */
367  /* dmom_w_adv_w[1]=0.0; */
368  /* dmom_w_adv_w[2]=0.0; */
369 
370  //u momentum diffusion tensor
371  mom_uu_diff_ten[0] =0.0;
372  mom_uu_diff_ten[1] =0.0;
373  /* mom_uu_diff_ten[2] = vos*nu; */
374 
375  mom_uv_diff_ten[0]=0.0;
376 
377  /* mom_uw_diff_ten[0]=vos*nu; */
378 
379  //v momentum diffusion tensor
380  mom_vv_diff_ten[0] =0.0;
381  mom_vv_diff_ten[1] =0.0;
382  /* mom_vv_diff_ten[2] = vos*mu; */
383 
384  mom_vu_diff_ten[0]=0.0;
385 
386  /* mom_vw_diff_ten[0]=vos*mu; */
387 
388  /* //w momentum diffusion tensor */
389  /* mom_ww_diff_ten[0] = vos*mu; */
390  /* mom_ww_diff_ten[1] = vos*mu; */
391  /* mom_ww_diff_ten[2] = vos*2.0*mu; */
392 
393  /* mom_wu_diff_ten[0]=vos*mu; */
394 
395  /* mom_wv_diff_ten[0]=vos*mu; */
396 
397  //momentum sources
398  mom_u_source = -rho*g[0];// - vos*d_mu*sigma*kappa*n[0]/(rho*(norm_n+1.0e-8));
399  mom_v_source = -rho*g[1];// - vos*d_mu*sigma*kappa*n[1]/(rho*(norm_n+1.0e-8));
400  /* mom_w_source = -vos*g[2];// - vos*d_mu*sigma*kappa*n[2]/(rho*(norm_n+1.0e-8)); */
401  if(vos > closure.frFraction_)
402  {
403  mom_u_source += (rho-rho_f)*g[0];
404  mom_v_source += (rho-rho_f)*g[1];
405 
406  }
407 
408  //u momentum Hamiltonian (pressure)
409  mom_u_ham = grad_p[0];
410  dmom_u_ham_grad_p[0]=1.0;
411  dmom_u_ham_grad_p[1]=0.0;
412  /* dmom_u_ham_grad_p[2]=0.0; */
413 
414  //v momentum Hamiltonian (pressure)
415  mom_v_ham = grad_p[1];
416  dmom_v_ham_grad_p[0]=0.0;
417  dmom_v_ham_grad_p[1]=1.0;
418  /* dmom_v_ham_grad_p[2]=0.0; */
419 
420  /* //w momentum Hamiltonian (pressure) */
421  /* mom_w_ham = vos*grad_p[2]; */
422  /* dmom_w_ham_grad_p[0]=0.0; */
423  /* dmom_w_ham_grad_p[1]=0.0; */
424  /* dmom_w_ham_grad_p[2]=vos/rho; */
425 
426  //u momentum Hamiltonian (advection)
427  mom_u_ham += rho*(uStar*grad_u[0]+vStar*grad_u[1]);
428  dmom_u_ham_grad_u[0]= rho*uStar;
429  dmom_u_ham_grad_u[1]= rho*vStar;
430  /* dmom_u_ham_grad_u[2]=vos*wStar; */
431 
432  //v momentum Hamiltonian (advection)
433  mom_v_ham += rho*(uStar*grad_v[0]+vStar*grad_v[1]);
434  dmom_v_ham_grad_v[0]= rho*uStar;
435  dmom_v_ham_grad_v[1]= rho*vStar;
436  /* dmom_v_ham_grad_v[2]=vos*wStar; */
437 
438  /* //w momentum Hamiltonian (advection) */
439  /* mom_w_ham = vos*(uStar*grad_w[0]+vStar*grad_w[1]+wStar*grad_w[2]); */
440  /* dmom_w_ham_grad_w[0]=vos*uStar; */
441  /* dmom_w_ham_grad_w[1]=vos*vStar; */
442  /* dmom_w_ham_grad_w[2]=vos*wStar; */
443  }
444  //VRANS specific
445  inline
446  void updateDarcyForchheimerTerms_Ergun(/* const double linearDragFactor, */
447  /* const double nonlinearDragFactor, */
448  /* const double vos, */
449  /* const double meanGrainSize, */
450  const double alpha,
451  const double beta,
452  const double eps_rho,
453  const double eps_mu,
454  const double rho_0,
455  const double nu_0,
456  const double rho_1,
457  const double nu_1,
458  double nu_t,
459  const double useVF,
460  const double vf,
461  const double phi,
462  const double u,
463  const double v,
464  const double w,
465  const double uStar,
466  const double vStar,
467  const double wStar,
468  const double eps_s,
469  const double vos,
470  const double u_f,
471  const double v_f,
472  const double w_f,
473  const double uStar_f,
474  const double vStar_f,
475  const double wStar_f,
476  double& mom_u_source,
477  double& mom_v_source,
478  double& mom_w_source,
479  double dmom_u_source[nSpace],
480  double dmom_v_source[nSpace],
481  double dmom_w_source[nSpace],
482  double gradC_x,
483  double gradC_y,
484  double gradC_z)
485  {
486  double rhoFluid, nuFluid,H_mu,viscosity;
487  H_mu = (1.0-useVF)*smoothedHeaviside(eps_mu,phi)+useVF*fmin(1.0,fmax(0.0,vf));
488  nuFluid = nu_0*(1.0-H_mu)+nu_1*H_mu;
489  rhoFluid = rho_0*(1.0-H_mu)+rho_1*H_mu;
490  //gco kinematic viscosity used, in sedclosure betaterm is multiplied by fluidDensity
491  viscosity = nuFluid;//mu; gco check
492  //vos is sediment fraction in this case - gco check
493 
494  double solid_velocity[2]={uStar,vStar}, fluid_velocity[2]={uStar_f, vStar_f};
495  double new_beta = closure.betaCoeff(vos,
496  rhoFluid,
497  fluid_velocity,
498  solid_velocity,
499  viscosity)*DRAG_FAC;
500 
501  //new_beta = 254800.0;//hack fall velocity of 0.1 with no pressure gradient
502  double beta2 = 156976.4;//hack, fall velocity of 0.1 with hydrostatic water
503  double one_by_vos = 2.0*vos/(vos*vos + fmax(1.0e-8,vos*vos));
504 
505  //new_beta/=rhoFluid;
506  //std::cout<<"total "<<(1.0-phi_s)*new_beta<<std::endl;
507  mom_u_source += new_beta*((u - u_f) + TURB_FORCE_FAC*DRAG_FAC*one_by_vos*nu_t*gradC_x/closure.sigmaC_) + (1.0-DRAG_FAC)*beta2*(u-u_f);
508 
509  mom_v_source += new_beta*((v - v_f) + TURB_FORCE_FAC*DRAG_FAC*one_by_vos*nu_t*gradC_y/closure.sigmaC_) + (1.0-DRAG_FAC)*beta2*(v-v_f);
510 
511  /* mom_w_source += new_beta*(w-w_s); */
512 
513  dmom_u_source[0] = new_beta + (1.0-DRAG_FAC)*beta2;
514  dmom_u_source[1] = 0.0;
515  /* dmom_u_source[2] = 0.0; */
516 
517  dmom_v_source[0] = 0.0;
518  dmom_v_source[1] = new_beta + (1.0-DRAG_FAC)*beta2;
519  /*dmom_v_source[2] = 0.0;*/
520 
521  /* dmom_w_source[0] = 0.0;
522  dmom_w_source[1] = 0.0;
523  dmom_w_source[2] = new_beta;*/
524  }
525 
526 
527  inline void updatePenaltyForPacking(const double vos,
528  const double u,
529  const double v,
530  const double w,
531  double& mom_u_source,
532  double& mom_v_source,
533  double& mom_w_source,
534  double dmom_u_source[nSpace],
535  double dmom_v_source[nSpace],
536  double dmom_w_source[nSpace])
537 
538  {
539  double meanPack = (closure.maxFraction_ + closure.vos_limiter_)/2.;
540  double epsPack = (closure.maxFraction_ - closure.vos_limiter_)/2.;
541  double dVos = vos - meanPack;
542  double sigma = smoothedHeaviside( epsPack, dVos);
543  double packPenalty = 1e6;
544  mom_u_source += sigma * packPenalty*u;
545  mom_v_source += sigma * packPenalty*v;
546  dmom_u_source[0] += sigma * packPenalty;
547  dmom_v_source[1] += sigma * packPenalty;
548 
549  }
550 
551 
552 
553  inline
554  void updateFrictionalPressure(const double vos,
555  const double grad_vos[nSpace],
556  double& mom_u_source,
557  double& mom_v_source,
558  double& mom_w_source)
559  {
560  double coeff = closure.gradp_friction(vos);
561  double one_by_vos = 2.0*vos/(vos*vos + fmax(1.0e-8,vos*vos));
562 
563  mom_u_source += coeff * grad_vos[0]*one_by_vos;
564  mom_v_source += coeff * grad_vos[1]*one_by_vos;
565  //mom_w_source += coeff * grad_vos[2];
566 
567  }
568 
569  inline
570  void updateFrictionalStress(const double LAG_MU_FR,
571  const double mu_fr_last,
572  double& mu_fr_new,
573  const double vos,
574  const double eps_rho,
575  const double eps_mu,
576  const double rho_0,
577  const double nu_0,
578  const double rho_1,
579  const double nu_1,
580  const double rho_s,
581  const double useVF,
582  const double vf,
583  const double phi,
584  const double grad_u[nSpace],
585  const double grad_v[nSpace],
586  const double grad_w[nSpace],
587  double mom_uu_diff_ten[nSpace],
588  double mom_uv_diff_ten[1],
589  double mom_uw_diff_ten[1],
590  double mom_vv_diff_ten[nSpace],
591  double mom_vu_diff_ten[1],
592  double mom_vw_diff_ten[1],
593  double mom_ww_diff_ten[nSpace],
594  double mom_wu_diff_ten[1],
595  double mom_wv_diff_ten[1])
596  {
597  double one_by_vos = 2.0*vos/(vos*vos + fmax(1.0e-8,vos*vos));
598 
599  mu_fr_new = closure.mu_fr(vos,
600  grad_u[0], grad_u[1], 0.,
601  grad_v[0], grad_v[1], 0.,
602  0., 0. , 0.)*one_by_vos;
603  double mu_fr = LAG_MU_FR*mu_fr_last + (1.0 - LAG_MU_FR)*mu_fr_new;
604 
605  mom_uu_diff_ten[0] += 2. * mu_fr * (2./3.);
606  mom_uu_diff_ten[1] += mu_fr;
607  /*mom_uu_diff_ten[2] += mu_fr; */
608 
609  mom_uv_diff_ten[0] += mu_fr;
610  /*mom_uw_diff_ten[0] += mu_fr; */
611 
612  mom_vv_diff_ten[0] += mu_fr;
613  mom_vv_diff_ten[1] += 2. * mu_fr * (2./3.) ;
614  /*mom_vv_diff_ten[2] += mu_fr; */
615 
616  mom_vu_diff_ten[0] += mu_fr;
617  /*mom_vw_diff_ten[0] += mu_fr; */
618 
619  /*mom_ww_diff_ten[0] += mu_fr; */
620  /*mom_ww_diff_ten[1] += mu_fr; */
621  /*mom_ww_diff_ten[2] += 2. * mu_fr * (2./3.) ; */
622 
623  /*mom_wu_diff_ten[0] += mu_fr; */
624  /*mom_wv_diff_ten[0] += mu_fr; */
625 
626  }
627 
628 
629 
630 
631  inline void calculateSubgridError_tau(const double &hFactor,
632  const double &elementDiameter,
633  const double &dmt,
634  const double &dm,
635  const double df[nSpace],
636  const double &a,
637  const double &pfac,
638  double &tau_v,
639  double &tau_p,
640  double &cfl)
641  {
642  double h, oneByAbsdt, density, viscosity, nrm_df;
643  h = hFactor * elementDiameter;
644  density = dm;
645  viscosity = a;
646  nrm_df = 0.0;
647  for (int I = 0; I < nSpace; I++)
648  {
649  nrm_df += df[I] * df[I];
650  }
651  nrm_df = sqrt(nrm_df);
652  if(density > 1e-8)
653  {
654  cfl = nrm_df / (fabs(h * density));
655  }
656  else
657  {
658  cfl = nrm_df / fabs(h );
659  }
660  oneByAbsdt = fabs(dmt);
661  tau_v = 1.0 / (4.0 * viscosity / (h * h) + 2.0 * nrm_df / h + oneByAbsdt);
662  tau_p = (4.0 * viscosity + 2.0 * nrm_df * h + oneByAbsdt * h * h) / pfac;
663  }
664 
665  inline void calculateSubgridError_tau(const double &Ct_sge,
666  const double &Cd_sge,
667  const double G[nSpace * nSpace],
668  const double &G_dd_G,
669  const double &tr_G,
670  const double &A0,
671  const double Ai[nSpace],
672  const double &Kij,
673  const double &pfac,
674  double &tau_v,
675  double &tau_p,
676  double &q_cfl)
677  {
678  double v_d_Gv = 0.0;
679  for (int I = 0; I < nSpace; I++)
680  for (int J = 0; J < nSpace; J++)
681  v_d_Gv += Ai[I] * G[I * nSpace + J] * Ai[J];
682  tau_v = 1.0 / sqrt(Ct_sge * A0 * A0 + v_d_Gv + Cd_sge * Kij * Kij * G_dd_G + 1.0e-12);
683  tau_p = 1.0 / (pfac * tr_G * tau_v);
684  }
685 
686  inline void calculateSubgridError_tauRes(const double &tau_p,
687  const double &tau_v,
688  const double &pdeResidualP,
689  const double &pdeResidualU,
690  const double &pdeResidualV,
691  const double &pdeResidualW,
692  double &subgridErrorP,
693  double &subgridErrorU,
694  double &subgridErrorV,
695  double &subgridErrorW)
696  {
697  /* GLS pressure */
698  subgridErrorP = -tau_p * pdeResidualP;
699  /* GLS momentum */
700  subgridErrorU = -tau_v * pdeResidualU;
701  subgridErrorV = -tau_v * pdeResidualV;
702  /* subgridErrorW = -tau_v*pdeResidualW; */
703  }
704 
705  inline void calculateSubgridErrorDerivatives_tauRes(const double &tau_p,
706  const double &tau_v,
707  const double dpdeResidualP_du[nDOF_trial_element],
708  const double dpdeResidualP_dv[nDOF_trial_element],
709  const double dpdeResidualP_dw[nDOF_trial_element],
710  const double dpdeResidualU_dp[nDOF_trial_element],
711  const double dpdeResidualU_du[nDOF_trial_element],
712  const double dpdeResidualV_dp[nDOF_trial_element],
713  const double dpdeResidualV_dv[nDOF_trial_element],
714  const double dpdeResidualW_dp[nDOF_trial_element],
715  const double dpdeResidualW_dw[nDOF_trial_element],
716  double dsubgridErrorP_du[nDOF_trial_element],
717  double dsubgridErrorP_dv[nDOF_trial_element],
718  double dsubgridErrorP_dw[nDOF_trial_element],
719  double dsubgridErrorU_dp[nDOF_trial_element],
720  double dsubgridErrorU_du[nDOF_trial_element],
721  double dsubgridErrorV_dp[nDOF_trial_element],
722  double dsubgridErrorV_dv[nDOF_trial_element],
723  double dsubgridErrorW_dp[nDOF_trial_element],
724  double dsubgridErrorW_dw[nDOF_trial_element])
725  {
726  for (int j = 0; j < nDOF_trial_element; j++)
727  {
728  /* GLS pressure */
729  dsubgridErrorP_du[j] = -tau_p * dpdeResidualP_du[j];
730  dsubgridErrorP_dv[j] = -tau_p * dpdeResidualP_dv[j];
731  /* dsubgridErrorP_dw[j] = -tau_p*dpdeResidualP_dw[j]; */
732  /* GLS momentum*/
733  /* u */
734  dsubgridErrorU_dp[j] = -tau_v * dpdeResidualU_dp[j];
735  dsubgridErrorU_du[j] = -tau_v * dpdeResidualU_du[j];
736  /* v */
737  dsubgridErrorV_dp[j] = -tau_v * dpdeResidualV_dp[j];
738  dsubgridErrorV_dv[j] = -tau_v * dpdeResidualV_dv[j];
739  /* /\* w *\/ */
740  /* dsubgridErrorW_dp[j] = -tau_v*dpdeResidualW_dp[j]; */
741  /* dsubgridErrorW_dw[j] = -tau_v*dpdeResidualW_dw[j]; */
742  }
743  }
744 
745  inline
746  void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_p,
747  const int& isDOFBoundary_u,
748  const int& isDOFBoundary_v,
749  const int& isDOFBoundary_w,
750  const int& isFluxBoundary_p,
751  const int& isFluxBoundary_u,
752  const int& isFluxBoundary_v,
753  const int& isFluxBoundary_w,
754  const double& oneByRho,
755  const double& bc_oneByRho,
756  const double n[nSpace],
757  const double& vos,
758  const double& bc_p,
759  const double& bc_u,
760  const double& bc_v,
761  const double& bc_w,
762  const double bc_f_mass[nSpace],
763  const double bc_f_umom[nSpace],
764  const double bc_f_vmom[nSpace],
765  const double bc_f_wmom[nSpace],
766  const double& bc_flux_mass,
767  const double& bc_flux_umom,
768  const double& bc_flux_vmom,
769  const double& bc_flux_wmom,
770  const double& p,
771  const double& u,
772  const double& v,
773  const double& w,
774  const double f_mass[nSpace],
775  const double f_umom[nSpace],
776  const double f_vmom[nSpace],
777  const double f_wmom[nSpace],
778  const double df_mass_du[nSpace],
779  const double df_mass_dv[nSpace],
780  const double df_mass_dw[nSpace],
781  const double df_umom_dp[nSpace],
782  const double df_umom_du[nSpace],
783  const double df_umom_dv[nSpace],
784  const double df_umom_dw[nSpace],
785  const double df_vmom_dp[nSpace],
786  const double df_vmom_du[nSpace],
787  const double df_vmom_dv[nSpace],
788  const double df_vmom_dw[nSpace],
789  const double df_wmom_dp[nSpace],
790  const double df_wmom_du[nSpace],
791  const double df_wmom_dv[nSpace],
792  const double df_wmom_dw[nSpace],
793  double& flux_mass,
794  double& flux_umom,
795  double& flux_vmom,
796  double& flux_wmom,
797  double* velocity_star,
798  double* velocity)
799  {
800  double flowSpeedNormal;
801  flux_mass = 0.0;
802  flux_umom = 0.0;
803  flux_vmom = 0.0;
804  /* flux_wmom = 0.0; */
805  flowSpeedNormal=(n[0]*velocity_star[0] +
806  n[1]*velocity_star[1]);
807  velocity[0] = u;
808  velocity[1] = v;
809  /* velocity[2] = w; */
810  //cek todo, switch back to conservative form weak BC's
811  if (isDOFBoundary_u != 1)
812  {
813  flux_mass += n[0]*f_mass[0];
814  }
815  else
816  {
817  flux_mass += n[0]*f_mass[0];
818  if (flowSpeedNormal < 0.0)
819  {
820  flux_umom+=flowSpeedNormal*(bc_u - u);
821  velocity[0] = bc_u;
822  }
823  }
824  if (isDOFBoundary_v != 1)
825  {
826  flux_mass+=n[1]*f_mass[1];
827  }
828  else
829  {
830  flux_mass+=n[1]*f_mass[1];
831  if (flowSpeedNormal < 0.0)
832  {
833  flux_vmom+=flowSpeedNormal*(bc_v - v);
834  velocity[1] = bc_v;
835  }
836  }
837  /* if (isDOFBoundary_w != 1) */
838  /* { */
839  /* flux_mass+=n[2]*f_mass[2]; */
840  /* } */
841  /* else */
842  /* { */
843  /* flux_mass +=n[2]*f_mass[2]; */
844  /* if (flowSpeedNormal < 0.0) */
845  /* { */
846  /* flux_wmom+=flowSpeedNormal*(bc_w - w); */
847  /* } */
848  /* } */
849  /* if (isDOFBoundary_w != 1) */
850  /* { */
851  /* flux_mass+=n[2]*f_mass[2]; */
852  /* } */
853  /* else */
854  /* { */
855  /* flux_mass +=n[2]*f_mass[2]; */
856  /* if (flowSpeedNormal < 0.0) */
857  /* flux_wmom+=bc_speed*(bc_w - w); */
858  /* } */
859 
860  if (isFluxBoundary_u == 1)
861  {
862  flux_umom = bc_flux_umom;
863  velocity[0] = bc_flux_umom;
864  }
865  if (isFluxBoundary_v == 1)
866  {
867  flux_vmom = bc_flux_vmom;
868  velocity[1] = bc_flux_umom;
869  }
870  /* if (isFluxBoundary_w == 1) */
871  /* { */
872  /* flux_wmom = bc_flux_wmom; */
873  /* } */
874  }
875 
876  inline
877  void exteriorNumericalAdvectiveFluxDerivatives(const int& isDOFBoundary_p,
878  const int& isDOFBoundary_u,
879  const int& isDOFBoundary_v,
880  const int& isDOFBoundary_w,
881  const int& isFluxBoundary_p,
882  const int& isFluxBoundary_u,
883  const int& isFluxBoundary_v,
884  const int& isFluxBoundary_w,
885  const double& oneByRho,
886  const double n[nSpace],
887  const double& vos,
888  const double& bc_p,
889  const double& bc_u,
890  const double& bc_v,
891  const double& bc_w,
892  const double bc_f_mass[nSpace],
893  const double bc_f_umom[nSpace],
894  const double bc_f_vmom[nSpace],
895  const double bc_f_wmom[nSpace],
896  const double& bc_flux_mass,
897  const double& bc_flux_umom,
898  const double& bc_flux_vmom,
899  const double& bc_flux_wmom,
900  const double& p,
901  const double& u,
902  const double& v,
903  const double& w,
904  const double f_mass[nSpace],
905  const double f_umom[nSpace],
906  const double f_vmom[nSpace],
907  const double f_wmom[nSpace],
908  const double df_mass_du[nSpace],
909  const double df_mass_dv[nSpace],
910  const double df_mass_dw[nSpace],
911  const double df_umom_dp[nSpace],
912  const double df_umom_du[nSpace],
913  const double df_umom_dv[nSpace],
914  const double df_umom_dw[nSpace],
915  const double df_vmom_dp[nSpace],
916  const double df_vmom_du[nSpace],
917  const double df_vmom_dv[nSpace],
918  const double df_vmom_dw[nSpace],
919  const double df_wmom_dp[nSpace],
920  const double df_wmom_du[nSpace],
921  const double df_wmom_dv[nSpace],
922  const double df_wmom_dw[nSpace],
923  double& dflux_mass_du,
924  double& dflux_mass_dv,
925  double& dflux_mass_dw,
926  double& dflux_umom_dp,
927  double& dflux_umom_du,
928  double& dflux_umom_dv,
929  double& dflux_umom_dw,
930  double& dflux_vmom_dp,
931  double& dflux_vmom_du,
932  double& dflux_vmom_dv,
933  double& dflux_vmom_dw,
934  double& dflux_wmom_dp,
935  double& dflux_wmom_du,
936  double& dflux_wmom_dv,
937  double& dflux_wmom_dw,
938  double* velocity_star)
939  {
940  double flowSpeedNormal;
941  dflux_mass_du = 0.0;
942  dflux_mass_dv = 0.0;
943  /* dflux_mass_dw = 0.0; */
944 
945  dflux_umom_dp = 0.0;
946  dflux_umom_du = 0.0;
947  dflux_umom_dv = 0.0;
948  /* dflux_umom_dw = 0.0; */
949  dflux_vmom_dp = 0.0;
950  dflux_vmom_du = 0.0;
951  dflux_vmom_dv = 0.0;
952  /* dflux_vmom_dw = 0.0; */
953 
954 
955  dflux_wmom_dp = 0.0;
956  dflux_wmom_du = 0.0;
957  dflux_wmom_dv = 0.0;
958  /* dflux_wmom_dw = 0.0; */
959  flowSpeedNormal=(n[0]*velocity_star[0] +
960  n[1]*velocity_star[1]);
961  if (isDOFBoundary_u != 1)
962  {
963  dflux_mass_du += n[0]*df_mass_du[0];
964  }
965  else
966  {
967  dflux_mass_du += n[0]*df_mass_du[0];
968  if (flowSpeedNormal < 0.0)
969  dflux_umom_du -= flowSpeedNormal;
970  }
971  if (isDOFBoundary_v != 1)
972  {
973  dflux_mass_dv += n[1]*df_mass_dv[1];
974  }
975  else
976  {
977  dflux_mass_dv += n[1]*df_mass_dv[1];
978  if (flowSpeedNormal < 0.0)
979  dflux_vmom_dv -= flowSpeedNormal;
980  }
981  /* if (isDOFBoundary_w != 1) */
982  /* { */
983  /* dflux_mass_dw+=n[2]*df_mass_dw[2]; */
984  /* } */
985  /* else */
986  /* { */
987  /* dflux_mass_dw += n[2]*df_mass_dw[2]; */
988  /* if (flowSpeedNormal < 0.0) */
989  /* dflux_wmom_dw -= flowSpeedNormal; */
990  /* } */
991  /* if (isDOFBoundary_w != 1) */
992  /* { */
993  /* dflux_mass_dw+=n[2]*df_mass_dw[2]; */
994  /* } */
995  /* else */
996  /* { */
997  /* dflux_mass_dw += n[2]*df_mass_dw[2]; */
998  /* if (flowSpeedNormal < 0.0) */
999  /* dflux_wmom_dw += bc_speed; */
1000  /* } */
1001  /* if (isDOFBoundary_p == 1) */
1002  /* { */
1003  /* dflux_umom_dp= -n[0]*oneByRho; */
1004  /* dflux_vmom_dp= -n[1]*oneByRho; */
1005  /* /\* dflux_wmom_dp= -n[2]*oneByRho; *\/ */
1006  /* } */
1007  /* if (isFluxBoundary_p == 1) */
1008  /* { */
1009  /* dflux_mass_du = 0.0; */
1010  /* dflux_mass_dv = 0.0; */
1011  /* /\* dflux_mass_dw = 0.0; *\/ */
1012  /* } */
1013  if (isFluxBoundary_u == 1)
1014  {
1015  dflux_umom_dp = 0.0;
1016  dflux_umom_du = 0.0;
1017  dflux_umom_dv = 0.0;
1018  /* dflux_umom_dw = 0.0; */
1019  }
1020  if (isFluxBoundary_v == 1)
1021  {
1022  dflux_vmom_dp = 0.0;
1023  dflux_vmom_du = 0.0;
1024  dflux_vmom_dv = 0.0;
1025  /* dflux_vmom_dw = 0.0; */
1026  }
1027  /* if (isFluxBoundary_w == 1) */
1028  /* { */
1029  /* dflux_wmom_dp = 0.0; */
1030  /* dflux_wmom_du = 0.0; */
1031  /* dflux_wmom_dv = 0.0; */
1032  /* dflux_wmom_dw = 0.0; */
1033  /* } */
1034  }
1035 
1036  inline
1037  void exteriorNumericalDiffusiveFlux(const double& eps,
1038  const double& phi,
1039  int* rowptr,
1040  int* colind,
1041  const int& isDOFBoundary,
1042  const int& isFluxBoundary,
1043  const double n[nSpace],
1044  double* bc_a,
1045  const double& bc_u,
1046  const double& bc_flux,
1047  double* a,
1048  const double grad_potential[nSpace],
1049  const double& u,
1050  const double& penalty,
1051  double& flux)
1052  {
1053  double diffusiveVelocityComponent_I,penaltyFlux,max_a;
1054  if(isFluxBoundary == 1)
1055  {
1056  flux = bc_flux;
1057  }
1058  else if(isDOFBoundary == 1)
1059  {
1060  flux = 0.0;
1061  max_a=0.0;
1062  for(int I=0;I<nSpace;I++)
1063  {
1064  diffusiveVelocityComponent_I=0.0;
1065  for(int m=rowptr[I];m<rowptr[I+1];m++)
1066  {
1067  diffusiveVelocityComponent_I -= a[m]*grad_potential[colind[m]];
1068  max_a = fmax(max_a,a[m]);
1069  }
1070  flux+= diffusiveVelocityComponent_I*n[I];
1071  }
1072  penaltyFlux = max_a*penalty*(u-bc_u);
1073  flux += penaltyFlux;
1074  //contact line slip
1075  //flux*=(smoothedDirac(eps,0) - smoothedDirac(eps,phi))/smoothedDirac(eps,0);
1076  }
1077  else
1078  {
1079  //std::cerr<<"warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0"<<std::endl;
1080  flux = 0.0;
1081  }
1082  }
1083 
1084  inline
1085  double ExteriorNumericalDiffusiveFluxJacobian(const double& eps,
1086  const double& phi,
1087  int* rowptr,
1088  int* colind,
1089  const int& isDOFBoundary,
1090  const int& isFluxBoundary,
1091  const double n[nSpace],
1092  double* a,
1093  const double& v,
1094  const double grad_v[nSpace],
1095  const double& penalty)
1096  {
1097  double dvel_I,tmp=0.0,max_a=0.0;
1098  if(isFluxBoundary==0 && isDOFBoundary==1)
1099  {
1100  for(int I=0;I<nSpace;I++)
1101  {
1102  dvel_I=0.0;
1103  for(int m=rowptr[I];m<rowptr[I+1];m++)
1104  {
1105  dvel_I -= a[m]*grad_v[colind[m]];
1106  max_a = fmax(max_a,a[m]);
1107  }
1108  tmp += dvel_I*n[I];
1109  }
1110  tmp +=max_a*penalty*v;
1111  //contact line slip
1112  //tmp*=(smoothedDirac(eps,0) - smoothedDirac(eps,phi))/smoothedDirac(eps,0);
1113  }
1114  return tmp;
1115  }
1116 
1117  void calculateResidual(arguments_dict& args)
1118  {
1119  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1120  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1121  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1122  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
1123  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
1124  double PSTAB = args.scalar<double>("PSTAB");
1125  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1126  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1127  xt::pyarray<double>& p_trial_ref = args.array<double>("p_trial_ref");
1128  xt::pyarray<double>& p_grad_trial_ref = args.array<double>("p_grad_trial_ref");
1129  xt::pyarray<double>& p_test_ref = args.array<double>("p_test_ref");
1130  xt::pyarray<double>& p_grad_test_ref = args.array<double>("p_grad_test_ref");
1131  xt::pyarray<double>& q_p = args.array<double>("q_p");
1132  xt::pyarray<double>& q_grad_p = args.array<double>("q_grad_p");
1133  xt::pyarray<double>& ebqe_p = args.array<double>("ebqe_p");
1134  xt::pyarray<double>& ebqe_grad_p = args.array<double>("ebqe_grad_p");
1135  xt::pyarray<double>& vel_trial_ref = args.array<double>("vel_trial_ref");
1136  xt::pyarray<double>& vel_grad_trial_ref = args.array<double>("vel_grad_trial_ref");
1137  xt::pyarray<double>& vel_test_ref = args.array<double>("vel_test_ref");
1138  xt::pyarray<double>& vel_grad_test_ref = args.array<double>("vel_grad_test_ref");
1139  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1140  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1141  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1142  xt::pyarray<double>& p_trial_trace_ref = args.array<double>("p_trial_trace_ref");
1143  xt::pyarray<double>& p_grad_trial_trace_ref = args.array<double>("p_grad_trial_trace_ref");
1144  xt::pyarray<double>& p_test_trace_ref = args.array<double>("p_test_trace_ref");
1145  xt::pyarray<double>& p_grad_test_trace_ref = args.array<double>("p_grad_test_trace_ref");
1146  xt::pyarray<double>& vel_trial_trace_ref = args.array<double>("vel_trial_trace_ref");
1147  xt::pyarray<double>& vel_grad_trial_trace_ref = args.array<double>("vel_grad_trial_trace_ref");
1148  xt::pyarray<double>& vel_test_trace_ref = args.array<double>("vel_test_trace_ref");
1149  xt::pyarray<double>& vel_grad_test_trace_ref = args.array<double>("vel_grad_test_trace_ref");
1150  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1151  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1152  double eb_adjoint_sigma = args.scalar<double>("eb_adjoint_sigma");
1153  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1154  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1155  double hFactor = args.scalar<double>("hFactor");
1156  int nElements_global = args.scalar<int>("nElements_global");
1157  int nElementBoundaries_owned = args.scalar<int>("nElementBoundaries_owned");
1158  double useRBLES = args.scalar<double>("useRBLES");
1159  double useMetrics = args.scalar<double>("useMetrics");
1160  double alphaBDF = args.scalar<double>("alphaBDF");
1161  double epsFact_rho = args.scalar<double>("epsFact_rho");
1162  double epsFact_mu = args.scalar<double>("epsFact_mu");
1163  double sigma = args.scalar<double>("sigma");
1164  double rho_0 = args.scalar<double>("rho_0");
1165  double nu_0 = args.scalar<double>("nu_0");
1166  double rho_1 = args.scalar<double>("rho_1");
1167  double nu_1 = args.scalar<double>("nu_1");
1168  double rho_s = args.scalar<double>("rho_s");
1169  double smagorinskyConstant = args.scalar<double>("smagorinskyConstant");
1170  int turbulenceClosureModel = args.scalar<int>("turbulenceClosureModel");
1171  double Ct_sge = args.scalar<double>("Ct_sge");
1172  double Cd_sge = args.scalar<double>("Cd_sge");
1173  double C_dc = args.scalar<double>("C_dc");
1174  double C_b = args.scalar<double>("C_b");
1175  const xt::pyarray<double>& eps_solid = args.array<double>("eps_solid");
1176  const xt::pyarray<double>& q_velocity_fluid = args.array<double>("q_velocity_fluid");
1177  const xt::pyarray<double>& q_velocityStar_fluid = args.array<double>("q_velocityStar_fluid");
1178  const xt::pyarray<double>& q_vos = args.array<double>("q_vos");
1179  const xt::pyarray<double>& q_dvos_dt = args.array<double>("q_dvos_dt");
1180  const xt::pyarray<double>& q_grad_vos = args.array<double>("q_grad_vos");
1181  const xt::pyarray<double>& q_dragAlpha = args.array<double>("q_dragAlpha");
1182  const xt::pyarray<double>& q_dragBeta = args.array<double>("q_dragBeta");
1183  const xt::pyarray<double>& q_mass_source = args.array<double>("q_mass_source");
1184  const xt::pyarray<double>& q_turb_var_0 = args.array<double>("q_turb_var_0");
1185  const xt::pyarray<double>& q_turb_var_1 = args.array<double>("q_turb_var_1");
1186  const xt::pyarray<double>& q_turb_var_grad_0 = args.array<double>("q_turb_var_grad_0");
1187  xt::pyarray<double>& q_eddy_viscosity = args.array<double>("q_eddy_viscosity");
1188  xt::pyarray<int>& p_l2g = args.array<int>("p_l2g");
1189  xt::pyarray<int>& vel_l2g = args.array<int>("vel_l2g");
1190  xt::pyarray<double>& p_dof = args.array<double>("p_dof");
1191  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1192  xt::pyarray<double>& v_dof = args.array<double>("v_dof");
1193  xt::pyarray<double>& w_dof = args.array<double>("w_dof");
1194  xt::pyarray<double>& g = args.array<double>("g");
1195  const double useVF = args.scalar<double>("useVF");
1196  xt::pyarray<double>& vf = args.array<double>("vf");
1197  xt::pyarray<double>& phi = args.array<double>("phi");
1198  xt::pyarray<double>& normal_phi = args.array<double>("normal_phi");
1199  xt::pyarray<double>& kappa_phi = args.array<double>("kappa_phi");
1200  xt::pyarray<double>& q_mom_u_acc = args.array<double>("q_mom_u_acc");
1201  xt::pyarray<double>& q_mom_v_acc = args.array<double>("q_mom_v_acc");
1202  xt::pyarray<double>& q_mom_w_acc = args.array<double>("q_mom_w_acc");
1203  xt::pyarray<double>& q_mass_adv = args.array<double>("q_mass_adv");
1204  xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.array<double>("q_mom_u_acc_beta_bdf");
1205  xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.array<double>("q_mom_v_acc_beta_bdf");
1206  xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.array<double>("q_mom_w_acc_beta_bdf");
1207  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
1208  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
1209  xt::pyarray<double>& q_velocity_sge = args.array<double>("q_velocity_sge");
1210  xt::pyarray<double>& ebqe_velocity_star = args.array<double>("ebqe_velocity_star");
1211  xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
1212  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1213  xt::pyarray<double>& q_numDiff_v = args.array<double>("q_numDiff_v");
1214  xt::pyarray<double>& q_numDiff_w = args.array<double>("q_numDiff_w");
1215  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1216  xt::pyarray<double>& q_numDiff_v_last = args.array<double>("q_numDiff_v_last");
1217  xt::pyarray<double>& q_numDiff_w_last = args.array<double>("q_numDiff_w_last");
1218  xt::pyarray<int>& sdInfo_u_u_rowptr = args.array<int>("sdInfo_u_u_rowptr");
1219  xt::pyarray<int>& sdInfo_u_u_colind = args.array<int>("sdInfo_u_u_colind");
1220  xt::pyarray<int>& sdInfo_u_v_rowptr = args.array<int>("sdInfo_u_v_rowptr");
1221  xt::pyarray<int>& sdInfo_u_v_colind = args.array<int>("sdInfo_u_v_colind");
1222  xt::pyarray<int>& sdInfo_u_w_rowptr = args.array<int>("sdInfo_u_w_rowptr");
1223  xt::pyarray<int>& sdInfo_u_w_colind = args.array<int>("sdInfo_u_w_colind");
1224  xt::pyarray<int>& sdInfo_v_v_rowptr = args.array<int>("sdInfo_v_v_rowptr");
1225  xt::pyarray<int>& sdInfo_v_v_colind = args.array<int>("sdInfo_v_v_colind");
1226  xt::pyarray<int>& sdInfo_v_u_rowptr = args.array<int>("sdInfo_v_u_rowptr");
1227  xt::pyarray<int>& sdInfo_v_u_colind = args.array<int>("sdInfo_v_u_colind");
1228  xt::pyarray<int>& sdInfo_v_w_rowptr = args.array<int>("sdInfo_v_w_rowptr");
1229  xt::pyarray<int>& sdInfo_v_w_colind = args.array<int>("sdInfo_v_w_colind");
1230  xt::pyarray<int>& sdInfo_w_w_rowptr = args.array<int>("sdInfo_w_w_rowptr");
1231  xt::pyarray<int>& sdInfo_w_w_colind = args.array<int>("sdInfo_w_w_colind");
1232  xt::pyarray<int>& sdInfo_w_u_rowptr = args.array<int>("sdInfo_w_u_rowptr");
1233  xt::pyarray<int>& sdInfo_w_u_colind = args.array<int>("sdInfo_w_u_colind");
1234  xt::pyarray<int>& sdInfo_w_v_rowptr = args.array<int>("sdInfo_w_v_rowptr");
1235  xt::pyarray<int>& sdInfo_w_v_colind = args.array<int>("sdInfo_w_v_colind");
1236  int offset_p = args.scalar<int>("offset_p");
1237  int offset_u = args.scalar<int>("offset_u");
1238  int offset_v = args.scalar<int>("offset_v");
1239  int offset_w = args.scalar<int>("offset_w");
1240  int stride_p = args.scalar<int>("stride_p");
1241  int stride_u = args.scalar<int>("stride_u");
1242  int stride_v = args.scalar<int>("stride_v");
1243  int stride_w = args.scalar<int>("stride_w");
1244  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1245  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1246  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1247  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1248  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1249  xt::pyarray<double>& ebqe_vf_ext = args.array<double>("ebqe_vf_ext");
1250  xt::pyarray<double>& bc_ebqe_vf_ext = args.array<double>("bc_ebqe_vf_ext");
1251  xt::pyarray<double>& ebqe_phi_ext = args.array<double>("ebqe_phi_ext");
1252  xt::pyarray<double>& bc_ebqe_phi_ext = args.array<double>("bc_ebqe_phi_ext");
1253  xt::pyarray<double>& ebqe_normal_phi_ext = args.array<double>("ebqe_normal_phi_ext");
1254  xt::pyarray<double>& ebqe_kappa_phi_ext = args.array<double>("ebqe_kappa_phi_ext");
1255  const xt::pyarray<double>& ebqe_vos_ext = args.array<double>("ebqe_vos_ext");
1256  const xt::pyarray<double>& ebqe_turb_var_0 = args.array<double>("ebqe_turb_var_0");
1257  const xt::pyarray<double>& ebqe_turb_var_1 = args.array<double>("ebqe_turb_var_1");
1258  xt::pyarray<int>& isDOFBoundary_p = args.array<int>("isDOFBoundary_p");
1259  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1260  xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
1261  xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
1262  xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.array<int>("isAdvectiveFluxBoundary_p");
1263  xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
1264  xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.array<int>("isAdvectiveFluxBoundary_v");
1265  xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.array<int>("isAdvectiveFluxBoundary_w");
1266  xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
1267  xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.array<int>("isDiffusiveFluxBoundary_v");
1268  xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.array<int>("isDiffusiveFluxBoundary_w");
1269  xt::pyarray<double>& ebqe_bc_p_ext = args.array<double>("ebqe_bc_p_ext");
1270  xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.array<double>("ebqe_bc_flux_mass_ext");
1271  xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.array<double>("ebqe_bc_flux_mom_u_adv_ext");
1272  xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.array<double>("ebqe_bc_flux_mom_v_adv_ext");
1273  xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.array<double>("ebqe_bc_flux_mom_w_adv_ext");
1274  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1275  xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.array<double>("ebqe_bc_flux_u_diff_ext");
1276  xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
1277  xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
1278  xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.array<double>("ebqe_bc_flux_v_diff_ext");
1279  xt::pyarray<double>& ebqe_bc_w_ext = args.array<double>("ebqe_bc_w_ext");
1280  xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.array<double>("ebqe_bc_flux_w_diff_ext");
1281  xt::pyarray<double>& q_x = args.array<double>("q_x");
1282  xt::pyarray<double>& q_velocity = args.array<double>("q_velocity");
1283  xt::pyarray<double>& ebqe_velocity = args.array<double>("ebqe_velocity");
1284  xt::pyarray<double>& flux = args.array<double>("flux");
1285  xt::pyarray<double>& elementResidual_p_save = args.array<double>("elementResidual_p_save");
1286  xt::pyarray<int>& elementFlags = args.array<int>("elementFlags");
1287  xt::pyarray<int>& boundaryFlags = args.array<int>("boundaryFlags");
1288  xt::pyarray<double>& barycenters = args.array<double>("barycenters");
1289  xt::pyarray<double>& wettedAreas = args.array<double>("wettedAreas");
1290  xt::pyarray<double>& netForces_p = args.array<double>("netForces_p");
1291  xt::pyarray<double>& netForces_v = args.array<double>("netForces_v");
1292  xt::pyarray<double>& netMoments = args.array<double>("netMoments");
1293  xt::pyarray<double>& ncDrag = args.array<double>("ncDrag");
1294  double LAG_MU_FR = args.scalar<double>("LAG_MU_FR");
1295  xt::pyarray<double>& q_mu_fr_last = args.array<double>("q_mu_fr_last");
1296  xt::pyarray<double>& q_mu_fr = args.array<double>("q_mu_fr");
1297  //
1298  //loop over elements to compute volume integrals and load them into element and global residual
1299  //
1300  for(int eN=0;eN<nElements_global;eN++)
1301  {
1302  //declare local storage for element residual and initialize
1303  register double elementResidual_u[nDOF_test_element],
1304  elementResidual_v[nDOF_test_element],
1305  //elementResidual_w[nDOF_test_element],
1306  mom_u_source_i[nDOF_test_element],
1307  mom_v_source_i[nDOF_test_element],
1308  eps_rho,eps_mu;
1309  double mesh_volume_conservation_element=0.0;
1310  for (int i=0;i<nDOF_test_element;i++)
1311  {
1312  int eN_i = eN*nDOF_test_element+i;
1313  elementResidual_p_save[eN_i]=0.0;
1314  elementResidual_u[i]=0.0;
1315  elementResidual_v[i]=0.0;
1316  /* elementResidual_w[i]=0.0; */
1317  mom_u_source_i[i]=0.0;
1318  mom_v_source_i[i]=0.0;
1319  }//i
1320  //
1321  //loop over quadrature points and compute integrands
1322  //
1323  for(int k=0;k<nQuadraturePoints_element;k++)
1324  {
1325  //compute indices and declare local storage
1326  register int eN_k = eN*nQuadraturePoints_element+k,
1327  eN_k_nSpace = eN_k*nSpace,
1328  eN_k_3d = eN_k*3,
1329  eN_nDOF_trial_element = eN*nDOF_trial_element;
1330  register double p=0.0,u=0.0,v=0.0,w=0.0,
1331  grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],grad_vos[nSpace],
1332  mom_u_acc=0.0,
1333  dmom_u_acc_u=0.0,
1334  mom_v_acc=0.0,
1335  dmom_v_acc_v=0.0,
1336  mom_w_acc=0.0,
1337  dmom_w_acc_w=0.0,
1338  mass_adv[nSpace],
1339  dmass_adv_u[nSpace],
1340  dmass_adv_v[nSpace],
1341  dmass_adv_w[nSpace],
1342  mom_u_adv[nSpace],
1343  dmom_u_adv_u[nSpace],
1344  dmom_u_adv_v[nSpace],
1345  dmom_u_adv_w[nSpace],
1346  mom_v_adv[nSpace],
1347  dmom_v_adv_u[nSpace],
1348  dmom_v_adv_v[nSpace],
1349  dmom_v_adv_w[nSpace],
1350  mom_w_adv[nSpace],
1351  dmom_w_adv_u[nSpace],
1352  dmom_w_adv_v[nSpace],
1353  dmom_w_adv_w[nSpace],
1354  mom_uu_diff_ten[nSpace],
1355  mom_vv_diff_ten[nSpace],
1356  mom_ww_diff_ten[nSpace],
1357  mom_uv_diff_ten[1],
1358  mom_uw_diff_ten[1],
1359  mom_vu_diff_ten[1],
1360  mom_vw_diff_ten[1],
1361  mom_wu_diff_ten[1],
1362  mom_wv_diff_ten[1],
1363  mom_u_source=0.0,
1364  mom_v_source=0.0,
1365  mom_w_source=0.0,
1366  mom_u_ham=0.0,
1367  dmom_u_ham_grad_p[nSpace],
1368  dmom_u_ham_grad_u[nSpace],
1369  mom_v_ham=0.0,
1370  dmom_v_ham_grad_p[nSpace],
1371  dmom_v_ham_grad_v[nSpace],
1372  mom_w_ham=0.0,
1373  dmom_w_ham_grad_p[nSpace],
1374  dmom_w_ham_grad_w[nSpace],
1375  mom_u_acc_t=0.0,
1376  dmom_u_acc_u_t=0.0,
1377  mom_v_acc_t=0.0,
1378  dmom_v_acc_v_t=0.0,
1379  pdeResidual_p=0.0,
1380  pdeResidual_u=0.0,
1381  pdeResidual_v=0.0,
1382  pdeResidual_w=0.0,
1383  Lstar_u_u[nDOF_test_element],
1384  Lstar_v_v[nDOF_test_element],
1385  subgridError_p=0.0,
1386  subgridError_u=0.0,
1387  subgridError_v=0.0,
1388  subgridError_w=0.0,
1389  tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
1390  tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
1391  jac[nSpace*nSpace],
1392  jacDet,
1393  jacInv[nSpace*nSpace],
1394  vel_grad_trial[nDOF_trial_element*nSpace],
1395  vel_test_dV[nDOF_trial_element],
1396  vel_grad_test_dV[nDOF_test_element*nSpace],
1397  dV,x,y,z,xt,yt,zt,
1398  //
1399  vos,
1400  //meanGrainSize,
1401  mass_source,
1402  dmom_u_source[nSpace],
1403  dmom_v_source[nSpace],
1404  dmom_w_source[nSpace],
1405  //
1406  G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
1407  //get jacobian, etc for mapping reference element
1408  ck.calculateMapping_element(eN,
1409  k,
1410  mesh_dof.data(),
1411  mesh_l2g.data(),
1412  mesh_trial_ref.data(),
1413  mesh_grad_trial_ref.data(),
1414  jac,
1415  jacDet,
1416  jacInv,
1417  x,y,z);
1418  ck.calculateH_element(eN,
1419  k,
1420  nodeDiametersArray.data(),
1421  mesh_l2g.data(),
1422  mesh_trial_ref.data(),
1423  h_phi);
1424  ck.calculateMappingVelocity_element(eN,
1425  k,
1426  mesh_velocity_dof.data(),
1427  mesh_l2g.data(),
1428  mesh_trial_ref.data(),
1429  xt,yt,zt);
1430  //xt=0.0;yt=0.0;zt=0.0;
1431  //std::cout<<"xt "<<xt<<'\t'<<yt<<'\t'<<zt<<std::endl;
1432  //get the physical integration weight
1433  dV = fabs(jacDet)*dV_ref.data()[k];
1434  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1435  //ck.calculateGScale(G,&normal_phi.data()[eN_k_nSpace],h_phi);
1436 
1437  eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1438  eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1439 
1440  //get the trial function gradients
1441  /* ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial); */
1442  ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
1443  //get the solution
1444  /* ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_ref.data()[k*nDOF_trial_element],p); */
1445  p = q_p.data()[eN_k];
1446  ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u);
1447  ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v);
1448  /* ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],w); */
1449  //get the solution gradients
1450  /* ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial,grad_p); */
1451  for (int I=0;I<nSpace;I++)
1452  grad_p[I] = q_grad_p.data()[eN_k_nSpace + I];
1453  for (int I=0;I<nSpace;I++)
1454  grad_vos[I] = q_grad_vos.data()[eN_k_nSpace + I];
1455  ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
1456  ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
1457  /* ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_w); */
1458  //VRANS
1459  vos = q_vos.data()[eN_k];//sed fraction - gco check
1460  //precalculate test function products with integration weights
1461  for (int j=0;j<nDOF_trial_element;j++)
1462  {
1463  /* p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV; */
1464  vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
1465  for (int I=0;I<nSpace;I++)
1466  {
1467  /* p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin */
1468  vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1469  }
1470  }
1471  //cek hack
1472  double div_mesh_velocity=0.0;
1473  int NDOF_MESH_TRIAL_ELEMENT=3;
1474  for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1475  {
1476  int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
1477  div_mesh_velocity +=
1478  mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+0]*vel_grad_trial[j*nSpace+0] +
1479  mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+1]*vel_grad_trial[j*nSpace+1];
1480  }
1481  mesh_volume_conservation_element += (alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity)*dV;
1482  div_mesh_velocity = DM3*div_mesh_velocity + (1.0-DM3)*alphaBDF*(dV-q_dV_last.data()[eN_k])/dV;
1483 
1484  //meanGrainSize = q_meanGrain[eN_k];
1485  //
1486  q_x.data()[eN_k_3d+0]=x;
1487  q_x.data()[eN_k_3d+1]=y;
1488  /* q_x.data()[eN_k_3d+2]=z; */
1489  /* // */
1490  //calculate pde coefficients at quadrature points
1491  //
1492  evaluateCoefficients(eps_rho,
1493  eps_mu,
1494  sigma,
1495  rho_0,
1496  nu_0,
1497  rho_1,
1498  nu_1,
1499  rho_s,
1500  elementDiameter.data()[eN],
1501  smagorinskyConstant,
1502  turbulenceClosureModel,
1503  g.data(),
1504  useVF,
1505  vf.data()[eN_k],
1506  phi.data()[eN_k],
1507  &normal_phi.data()[eN_k_nSpace],
1508  kappa_phi.data()[eN_k],
1509  //VRANS
1510  vos,
1511  //
1512  p,
1513  grad_p,
1514  grad_u,
1515  grad_v,
1516  grad_w,
1517  u,
1518  v,
1519  w,
1520  q_velocity_sge.data()[eN_k_nSpace+0],
1521  q_velocity_sge.data()[eN_k_nSpace+1],
1522  q_velocity_sge.data()[eN_k_nSpace+1],//hack, shouldn't be used
1523  q_eddy_viscosity.data()[eN_k],
1524  mom_u_acc,
1525  dmom_u_acc_u,
1526  mom_v_acc,
1527  dmom_v_acc_v,
1528  mom_w_acc,
1529  dmom_w_acc_w,
1530  mass_adv,
1531  dmass_adv_u,
1532  dmass_adv_v,
1533  dmass_adv_w,
1534  mom_u_adv,
1535  dmom_u_adv_u,
1536  dmom_u_adv_v,
1537  dmom_u_adv_w,
1538  mom_v_adv,
1539  dmom_v_adv_u,
1540  dmom_v_adv_v,
1541  dmom_v_adv_w,
1542  mom_w_adv,
1543  dmom_w_adv_u,
1544  dmom_w_adv_v,
1545  dmom_w_adv_w,
1546  mom_uu_diff_ten,
1547  mom_vv_diff_ten,
1548  mom_ww_diff_ten,
1549  mom_uv_diff_ten,
1550  mom_uw_diff_ten,
1551  mom_vu_diff_ten,
1552  mom_vw_diff_ten,
1553  mom_wu_diff_ten,
1554  mom_wv_diff_ten,
1555  mom_u_source,
1556  mom_v_source,
1557  mom_w_source,
1558  mom_u_ham,
1559  dmom_u_ham_grad_p,
1560  dmom_u_ham_grad_u,
1561  mom_v_ham,
1562  dmom_v_ham_grad_p,
1563  dmom_v_ham_grad_v,
1564  mom_w_ham,
1565  dmom_w_ham_grad_p,
1566  dmom_w_ham_grad_w);
1567  //VRANS
1568  mass_source = q_mass_source.data()[eN_k];
1569  for (int I=0;I<nSpace;I++)
1570  {
1571  dmom_u_source[I] = 0.0;
1572  dmom_v_source[I] = 0.0;
1573  dmom_w_source[I] = 0.0;
1574  }
1575  updateDarcyForchheimerTerms_Ergun(// linearDragFactor,
1576  // nonlinearDragFactor,
1577  // vos,
1578  // meanGrainSize,
1579  q_dragAlpha.data()[eN_k],
1580  q_dragBeta.data()[eN_k],
1581  eps_rho,
1582  eps_mu,
1583  rho_0,
1584  nu_0,
1585  rho_1,
1586  nu_1,
1587  q_eddy_viscosity.data()[eN_k],
1588  useVF,
1589  vf.data()[eN_k],
1590  phi.data()[eN_k],
1591  u,
1592  v,
1593  w,
1594  q_velocity_sge.data()[eN_k_nSpace+0],
1595  q_velocity_sge.data()[eN_k_nSpace+1],
1596  q_velocity_sge.data()[eN_k_nSpace+1],//hack, shouldn't be used
1597  eps_solid.data()[elementFlags.data()[eN]],
1598  vos,
1599  q_velocity_fluid.data()[eN_k_nSpace+0],
1600  q_velocity_fluid.data()[eN_k_nSpace+1],
1601  q_velocity_fluid.data()[eN_k_nSpace+1],
1602  q_velocityStar_fluid.data()[eN_k_nSpace+0],
1603  q_velocityStar_fluid.data()[eN_k_nSpace+1],
1604  q_velocityStar_fluid.data()[eN_k_nSpace+1],
1605  mom_u_source,
1606  mom_v_source,
1607  mom_w_source,
1608  dmom_u_source,
1609  dmom_v_source,
1610  dmom_w_source,
1611  q_grad_vos.data()[eN_k_nSpace+0],
1612  q_grad_vos.data()[eN_k_nSpace+1],
1613  q_grad_vos.data()[eN_k_nSpace+1]);
1614 
1615  updateFrictionalPressure(vos,
1616  grad_vos,
1617  mom_u_source,
1618  mom_v_source,
1619  mom_w_source);
1620  updateFrictionalStress(LAG_MU_FR,
1621  q_mu_fr_last.data()[eN_k],
1622  q_mu_fr.data()[eN_k],
1623  vos,
1624  eps_rho,
1625  eps_mu,
1626  rho_0,
1627  nu_0,
1628  rho_1,
1629  nu_1,
1630  rho_s,
1631  useVF,
1632  vf.data()[eN_k],
1633  phi.data()[eN_k],
1634  grad_u,
1635  grad_v,
1636  grad_w,
1637  mom_uu_diff_ten,
1638  mom_uv_diff_ten,
1639  mom_uw_diff_ten,
1640  mom_vv_diff_ten,
1641  mom_vu_diff_ten,
1642  mom_vw_diff_ten,
1643  mom_ww_diff_ten,
1644  mom_wu_diff_ten,
1645  mom_wv_diff_ten);
1646  //Turbulence closure model
1647  //
1648  //save momentum for time history and velocity for subgrid error
1649  //
1650  q_mom_u_acc.data()[eN_k] = mom_u_acc;
1651  q_mom_v_acc.data()[eN_k] = mom_v_acc;
1652  /* q_mom_w_acc.data()[eN_k] = mom_w_acc; */
1653  //subgrid error uses grid scale velocity
1654  q_mass_adv.data()[eN_k_nSpace+0] = u;
1655  q_mass_adv.data()[eN_k_nSpace+1] = v;
1656  /* q_mass_adv.data()[eN_k_nSpace+2] = w; */
1657  //
1658  //moving mesh
1659  //
1660  mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*xt;
1661  mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
1662  /* mom_u_adv[2] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*zt; */
1663  dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*xt;
1664  dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
1665  /* dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt; */
1666 
1667  mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*xt;
1668  mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
1669  /* mom_v_adv[2] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*zt; */
1670  dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*xt;
1671  dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
1672  /* dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt; */
1673 
1674  /* mom_w_adv[0] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*xt; */
1675  /* mom_w_adv[1] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*yt; */
1676  /* mom_w_adv[2] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*zt; */
1677  /* dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*xt; */
1678  /* dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt; */
1679  /* dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt; */
1680  //
1681  //calculate time derivative at quadrature points
1682  //
1683  if (q_dV_last.data()[eN_k] <= -100)
1684  q_dV_last.data()[eN_k] = dV;
1685  q_dV.data()[eN_k] = dV;
1686  ck.bdf(alphaBDF,
1687  q_mom_u_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
1688  mom_u_acc,
1689  dmom_u_acc_u,
1690  mom_u_acc_t,
1691  dmom_u_acc_u_t);
1692  ck.bdf(alphaBDF,
1693  q_mom_v_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
1694  mom_v_acc,
1695  dmom_v_acc_v,
1696  mom_v_acc_t,
1697  dmom_v_acc_v_t);
1698 
1699  /* ck.bdf(alphaBDF, */
1700  /* q_mom_w_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV, */
1701  /* mom_w_acc, */
1702  /* dmom_w_acc_w, */
1703  /* mom_w_acc_t, */
1704  /* dmom_w_acc_w_t); */
1705  /* // */
1706 
1707  mom_u_acc_t *= dmom_u_acc_u;
1708  mom_v_acc_t *= dmom_v_acc_v;
1709 
1710  //calculate subgrid error (strong residual and adjoint)
1711  //
1712  //calculate strong residual
1713  pdeResidual_p =
1714  ck.Mass_strong(q_dvos_dt.data()[eN_k]) +
1715  ck.Advection_strong(dmass_adv_u,grad_u) +
1716  ck.Advection_strong(dmass_adv_v,grad_v) +
1717  /* ck.Advection_strong(dmass_adv_w,grad_w) + */
1718  DM2*MOVING_DOMAIN*ck.Reaction_strong(alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity) +
1719  //VRANS
1720  ck.Reaction_strong(mass_source);
1721  //
1722 
1723  dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt);
1724  dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
1725  /* dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt); */
1726 
1727  pdeResidual_u =
1728  ck.Mass_strong(mom_u_acc_t) +
1729  ck.Advection_strong(dmom_adv_sge,grad_u) + //note here and below: same in cons. and non-cons.
1730  ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
1731  ck.Reaction_strong(mom_u_source) -
1732  ck.Reaction_strong(u*div_mesh_velocity);
1733 
1734  pdeResidual_v =
1735  ck.Mass_strong(mom_v_acc_t) +
1736  ck.Advection_strong(dmom_adv_sge,grad_v) +
1737  ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
1738  ck.Reaction_strong(mom_v_source) -
1739  ck.Reaction_strong(v*div_mesh_velocity);
1740 
1741  /* pdeResidual_w = ck.Mass_strong(mom_w_acc_t) + */
1742  /* ck.Advection_strong(dmom_adv_sge,grad_w) + */
1743  /* ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) + */
1744  /* ck.Reaction_strong(mom_w_source) - */
1745  /* ck.Reaction_strong(w*div_mesh_velocity); */
1746 
1747  //calculate tau and tau*Res
1748  //cek debug
1749  double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
1750  calculateSubgridError_tau(hFactor,
1751  elementDiameter.data()[eN],
1752  tmpR,//dmom_u_acc_u_t,
1753  dmom_u_acc_u,
1754  dmom_adv_sge,
1755  mom_uu_diff_ten[1],
1756  dmom_u_ham_grad_p[0],
1757  tau_v0,
1758  tau_p0,
1759  q_cfl.data()[eN_k]);
1760 
1761  calculateSubgridError_tau(Ct_sge,Cd_sge,
1762  G,G_dd_G,tr_G,
1763  tmpR,//dmom_u_acc_u_t,
1764  dmom_adv_sge,
1765  mom_uu_diff_ten[1],
1766  dmom_u_ham_grad_p[0],
1767  tau_v1,
1768  tau_p1,
1769  q_cfl.data()[eN_k]);
1770 
1771  tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
1772  tau_p = PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
1773 
1774  calculateSubgridError_tauRes(tau_p,
1775  tau_v,
1776  pdeResidual_p,
1777  pdeResidual_u,
1778  pdeResidual_v,
1779  pdeResidual_w,
1780  subgridError_p,
1781  subgridError_u,
1782  subgridError_v,
1783  subgridError_w);
1784  // velocity used in adjoint (VMS or RBLES, with or without lagging the grid scale velocity)
1785  dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt + useRBLES*subgridError_u);
1786  dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
1787  /* dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w); */
1788 
1789  mom_u_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+0]);
1790  mom_u_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+0]);
1791  /* mom_u_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge.data()[eN_k_nSpace+0]); */
1792 
1793  // adjoint times the test functions
1794  for (int i=0;i<nDOF_test_element;i++)
1795  {
1796  register int i_nSpace = i*nSpace;
1797  /* Lstar_u_p[i]=ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]); */
1798  /* Lstar_v_p[i]=ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]); */
1799  /* Lstar_w_p[i]=ck.Advection_adjoint(dmass_adv_w,&p_grad_test_dV[i_nSpace]); */
1800  //use the same advection adjoint for all three since we're approximating the linearized adjoint
1801  Lstar_u_u[i]=ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
1802  Lstar_v_v[i]=ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
1803  /* Lstar_w_w[i]=ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]); */
1804  /* Lstar_p_w[i]=ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]); */
1805 
1806  //VRANS account for drag terms, diagonal only here ... decide if need off diagonal terms too
1807  Lstar_u_u[i]+=ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
1808  Lstar_v_v[i]+=ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
1809  /* Lstar_w_w[i]+=ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]); */
1810  //
1811  }
1812 
1813  norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v);// + pdeResidual_w*pdeResidual_w);
1814  q_numDiff_u.data()[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
1815  (1.0-useMetrics)*hFactor*hFactor*elementDiameter.data()[eN]*elementDiameter.data()[eN]);
1816  q_numDiff_v.data()[eN_k] = q_numDiff_u.data()[eN_k];
1817  q_numDiff_w.data()[eN_k] = q_numDiff_u.data()[eN_k];
1818  // Save velocity and its gradient (to be used in other models and to compute errors)
1819  q_velocity.data()[eN_k_nSpace+0]=u;
1820  q_velocity.data()[eN_k_nSpace+1]=v;
1821  /* q_velocity.data()[eN_k_nSpace+2]=w;*/
1822  for(int i=0;i<nDOF_test_element;i++)
1823  {
1824  register int i_nSpace=i*nSpace;
1825  /* std::cout<<"elemRes_mesh "<<mesh_vel[0]<<'\t'<<mesh_vel[2]<<'\t'<<p_test_dV[i]<<'\t'<<(q_dV_last.data()[eN_k]/dV)<<'\t'<<dV<<std::endl; */
1826  /* elementResidual_mesh[i] += ck.Reaction_weak(1.0,p_test_dV[i]) - */
1827  /* ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) - */
1828  /* ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace]); */
1829 
1830  /* elementResidual_p.data()[i] += ck.Mass_weak(-q_dvos_dt.data()[eN_k],p_test_dV[i]) + */
1831  /* ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) + */
1832  /* DM*MOVING_DOMAIN*(ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) - */
1833  /* ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) - */
1834  /* ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])) + */
1835  /* //VRANS */
1836  /* ck.Reaction_weak(mass_source,p_test_dV[i]) + //VRANS source term for wave maker */
1837  /* // */
1838  /* ck.SubgridError(subgridError_u,Lstar_u_p[i]) + */
1839  /* ck.SubgridError(subgridError_v,Lstar_v_p[i]);// + */
1840  /* /\* ck.SubgridError(subgridError_w,Lstar_w_p[i]); *\/ */
1841 
1842  elementResidual_u[i] +=
1843  ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
1844  ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
1845  ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
1846  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]) +
1847  /* ck.Diffusion_weak(sdInfo_u_w_rowptr.data(),sdInfo_u_w_colind.data(),mom_uw_diff_ten,grad_w,&vel_grad_test_dV[i_nSpace]) + */
1848  ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
1849  ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
1850  /* ck.SubgridError(subgridError_p,Lstar_p_u[i]) + */
1851  ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
1852  ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]);
1853  mom_u_source_i[i] += ck.Reaction_weak(mom_u_source,vel_test_dV[i]);
1854 
1855  elementResidual_v[i] +=
1856  ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
1857  ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
1858  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]) +
1859  ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
1860  /* ck.Diffusion_weak(sdInfo_v_w_rowptr.data(),sdInfo_v_w_colind.data(),mom_vw_diff_ten,grad_w,&vel_grad_test_dV[i_nSpace]) + */
1861  ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
1862  ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
1863  /* ck.SubgridError(subgridError_p,Lstar_p_v[i]) + */
1864  ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
1865  ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]);
1866  mom_v_source_i[i] += ck.Reaction_weak(mom_v_source,vel_test_dV[i]);
1867 
1868  /* elementResidual_w[i] +=
1869  ck.Mass_weak(mom_w_acc_t,vel_test_dV[i]) + */
1870  /* ck.Advection_weak(mom_w_adv,&vel_grad_test_dV[i_nSpace]) + */
1871  /* ck.Diffusion_weak(sdInfo_w_u_rowptr.data(),sdInfo_w_u_colind.data(),mom_wu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) + */
1872  /* ck.Diffusion_weak(sdInfo_w_v_rowptr.data(),sdInfo_w_v_colind.data(),mom_wv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) + */
1873  /* ck.Diffusion_weak(sdInfo_w_w_rowptr.data(),sdInfo_w_w_colind.data(),mom_ww_diff_ten,grad_w,&vel_grad_test_dV[i_nSpace]) + */
1874  /* ck.Reaction_weak(mom_w_source,vel_test_dV[i]) + */
1875  /* ck.Hamiltonian_weak(mom_w_ham,vel_test_dV[i]) + */
1876  /* ck.SubgridError(subgridError_p,Lstar_p_w[i]) + */
1877  /* ck.SubgridError(subgridError_w,Lstar_w_w[i]) + */
1878  /* ck.NumericalDiffusion(q_numDiff_w_last.data()[eN_k],grad_w,&vel_grad_test_dV[i_nSpace]); */
1879  }//i
1880  }
1881  //
1882  //load element into global residual and save element residual
1883  //
1884  for(int i=0;i<nDOF_test_element;i++)
1885  {
1886  register int eN_i=eN*nDOF_test_element+i;
1887 
1888  /* elementResidual_p_save[eN_i] += elementResidual_p.data()[i]; */
1889  /* mesh_volume_conservation_element_weak += elementResidual_mesh[i]; */
1890  /* globalResidual.data()[offset_p+stride_p*p_l2g.data()[eN_i]]+=elementResidual_p.data()[i]; */
1891  globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
1892  globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
1893  ncDrag.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=mom_u_source_i[i];
1894  ncDrag.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=mom_v_source_i[i];
1895  /* globalResidual.data()[offset_w+stride_w*vel_l2g.data()[eN_i]]+=elementResidual_w[i]; */
1896  }//i
1897  /* mesh_volume_conservation += mesh_volume_conservation_element; */
1898  /* mesh_volume_conservation_weak += mesh_volume_conservation_element_weak; */
1899  /* mesh_volume_conservation_err_max=fmax(mesh_volume_conservation_err_max,fabs(mesh_volume_conservation_element)); */
1900  /* mesh_volume_conservation_err_max_weak=fmax(mesh_volume_conservation_err_max_weak,fabs(mesh_volume_conservation_element_weak)); */
1901  }//elements
1902  //
1903  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
1904  //
1905  //ebNE is the Exterior element boundary INdex
1906  //ebN is the element boundary INdex
1907  //eN is the element index
1908  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1909  {
1910  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
1911  eN = elementBoundaryElementsArray.data()[ebN*2+0],
1912  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1913  eN_nDOF_trial_element = eN*nDOF_trial_element;
1914  register double elementResidual_u[nDOF_test_element],
1915  elementResidual_v[nDOF_test_element],
1916  //elementResidual_w[nDOF_test_element],
1917  eps_rho,eps_mu;
1918  for (int i=0;i<nDOF_test_element;i++)
1919  {
1920  elementResidual_u[i]=0.0;
1921  elementResidual_v[i]=0.0;
1922  /* elementResidual_w[i]=0.0; */
1923  }
1924  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1925  {
1926  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1927  ebNE_kb_nSpace = ebNE_kb*nSpace,
1928  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1929  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1930  register double p_ext=0.0,
1931  u_ext=0.0,
1932  v_ext=0.0,
1933  w_ext=0.0,
1934  grad_p_ext[nSpace],
1935  grad_u_ext[nSpace],
1936  grad_v_ext[nSpace],
1937  grad_w_ext[nSpace],
1938  mom_u_acc_ext=0.0,
1939  dmom_u_acc_u_ext=0.0,
1940  mom_v_acc_ext=0.0,
1941  dmom_v_acc_v_ext=0.0,
1942  mom_w_acc_ext=0.0,
1943  dmom_w_acc_w_ext=0.0,
1944  mass_adv_ext[nSpace],
1945  dmass_adv_u_ext[nSpace],
1946  dmass_adv_v_ext[nSpace],
1947  dmass_adv_w_ext[nSpace],
1948  mom_u_adv_ext[nSpace],
1949  dmom_u_adv_u_ext[nSpace],
1950  dmom_u_adv_v_ext[nSpace],
1951  dmom_u_adv_w_ext[nSpace],
1952  mom_v_adv_ext[nSpace],
1953  dmom_v_adv_u_ext[nSpace],
1954  dmom_v_adv_v_ext[nSpace],
1955  dmom_v_adv_w_ext[nSpace],
1956  mom_w_adv_ext[nSpace],
1957  dmom_w_adv_u_ext[nSpace],
1958  dmom_w_adv_v_ext[nSpace],
1959  dmom_w_adv_w_ext[nSpace],
1960  mom_uu_diff_ten_ext[nSpace],
1961  mom_vv_diff_ten_ext[nSpace],
1962  mom_ww_diff_ten_ext[nSpace],
1963  mom_uv_diff_ten_ext[1],
1964  mom_uw_diff_ten_ext[1],
1965  mom_vu_diff_ten_ext[1],
1966  mom_vw_diff_ten_ext[1],
1967  mom_wu_diff_ten_ext[1],
1968  mom_wv_diff_ten_ext[1],
1969  mom_u_source_ext=0.0,
1970  mom_v_source_ext=0.0,
1971  mom_w_source_ext=0.0,
1972  mom_u_ham_ext=0.0,
1973  dmom_u_ham_grad_p_ext[nSpace],
1974  dmom_u_ham_grad_u_ext[nSpace],
1975  mom_v_ham_ext=0.0,
1976  dmom_v_ham_grad_p_ext[nSpace],
1977  dmom_v_ham_grad_v_ext[nSpace],
1978  mom_w_ham_ext=0.0,
1979  dmom_w_ham_grad_p_ext[nSpace],
1980  dmom_w_ham_grad_w_ext[nSpace],
1981  dmom_u_adv_p_ext[nSpace],
1982  dmom_v_adv_p_ext[nSpace],
1983  dmom_w_adv_p_ext[nSpace],
1984  flux_mass_ext=0.0,
1985  flux_mom_u_adv_ext=0.0,
1986  flux_mom_v_adv_ext=0.0,
1987  flux_mom_w_adv_ext=0.0,
1988  flux_mom_uu_diff_ext=0.0,
1989  flux_mom_uv_diff_ext=0.0,
1990  flux_mom_uw_diff_ext=0.0,
1991  flux_mom_vu_diff_ext=0.0,
1992  flux_mom_vv_diff_ext=0.0,
1993  flux_mom_vw_diff_ext=0.0,
1994  bc_p_ext=0.0,
1995  bc_u_ext=0.0,
1996  bc_v_ext=0.0,
1997  bc_w_ext=0.0,
1998  bc_mom_u_acc_ext=0.0,
1999  bc_dmom_u_acc_u_ext=0.0,
2000  bc_mom_v_acc_ext=0.0,
2001  bc_dmom_v_acc_v_ext=0.0,
2002  bc_mom_w_acc_ext=0.0,
2003  bc_dmom_w_acc_w_ext=0.0,
2004  bc_mass_adv_ext[nSpace],
2005  bc_dmass_adv_u_ext[nSpace],
2006  bc_dmass_adv_v_ext[nSpace],
2007  bc_dmass_adv_w_ext[nSpace],
2008  bc_mom_u_adv_ext[nSpace],
2009  bc_dmom_u_adv_u_ext[nSpace],
2010  bc_dmom_u_adv_v_ext[nSpace],
2011  bc_dmom_u_adv_w_ext[nSpace],
2012  bc_mom_v_adv_ext[nSpace],
2013  bc_dmom_v_adv_u_ext[nSpace],
2014  bc_dmom_v_adv_v_ext[nSpace],
2015  bc_dmom_v_adv_w_ext[nSpace],
2016  bc_mom_w_adv_ext[nSpace],
2017  bc_dmom_w_adv_u_ext[nSpace],
2018  bc_dmom_w_adv_v_ext[nSpace],
2019  bc_dmom_w_adv_w_ext[nSpace],
2020  bc_mom_uu_diff_ten_ext[nSpace],
2021  bc_mom_vv_diff_ten_ext[nSpace],
2022  bc_mom_ww_diff_ten_ext[nSpace],
2023  bc_mom_uv_diff_ten_ext[1],
2024  bc_mom_uw_diff_ten_ext[1],
2025  bc_mom_vu_diff_ten_ext[1],
2026  bc_mom_vw_diff_ten_ext[1],
2027  bc_mom_wu_diff_ten_ext[1],
2028  bc_mom_wv_diff_ten_ext[1],
2029  bc_mom_u_source_ext=0.0,
2030  bc_mom_v_source_ext=0.0,
2031  bc_mom_w_source_ext=0.0,
2032  bc_mom_u_ham_ext=0.0,
2033  bc_dmom_u_ham_grad_p_ext[nSpace],
2034  bc_dmom_u_ham_grad_u_ext[nSpace],
2035  bc_mom_v_ham_ext=0.0,
2036  bc_dmom_v_ham_grad_p_ext[nSpace],
2037  bc_dmom_v_ham_grad_v_ext[nSpace],
2038  bc_mom_w_ham_ext=0.0,
2039  bc_dmom_w_ham_grad_p_ext[nSpace],
2040  bc_dmom_w_ham_grad_w_ext[nSpace],
2041  jac_ext[nSpace*nSpace],
2042  jacDet_ext,
2043  jacInv_ext[nSpace*nSpace],
2044  boundaryJac[nSpace*(nSpace-1)],
2045  metricTensor[(nSpace-1)*(nSpace-1)],
2046  metricTensorDetSqrt,
2047  dS,vel_test_dS[nDOF_test_element],
2048  vel_grad_trial_trace[nDOF_trial_element*nSpace],
2049  vel_grad_test_dS[nDOF_trial_element*nSpace],
2050  normal[2],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
2051  //VRANS
2052  vos_ext,
2053  //
2054  G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
2055  force_x,force_y,force_p_x,force_p_y,force_v_x,force_v_y,r_x,r_y;
2056  //compute information about mapping from reference element to physical element
2057  ck.calculateMapping_elementBoundary(eN,
2058  ebN_local,
2059  kb,
2060  ebN_local_kb,
2061  mesh_dof.data(),
2062  mesh_l2g.data(),
2063  mesh_trial_trace_ref.data(),
2064  mesh_grad_trial_trace_ref.data(),
2065  boundaryJac_ref.data(),
2066  jac_ext,
2067  jacDet_ext,
2068  jacInv_ext,
2069  boundaryJac,
2070  metricTensor,
2071  metricTensorDetSqrt,
2072  normal_ref.data(),
2073  normal,
2074  x_ext,y_ext,z_ext);
2075  ck.calculateMappingVelocity_elementBoundary(eN,
2076  ebN_local,
2077  kb,
2078  ebN_local_kb,
2079  mesh_velocity_dof.data(),
2080  mesh_l2g.data(),
2081  mesh_trial_trace_ref.data(),
2082  xt_ext,yt_ext,zt_ext,
2083  normal,
2084  boundaryJac,
2085  metricTensor,
2086  integralScaling);
2087  //xt_ext=0.0;yt_ext=0.0;zt_ext=0.0;
2088  //std::cout<<"xt_ext "<<xt_ext<<'\t'<<yt_ext<<'\t'<<zt_ext<<std::endl;
2089  //std::cout<<"x_ext "<<x_ext<<'\t'<<y_ext<<'\t'<<z_ext<<std::endl;
2090  //std::cout<<"integralScaling - metricTensorDetSrt ==============================="<<integralScaling-metricTensorDetSqrt<<std::endl;
2091  /* std::cout<<"metricTensorDetSqrt "<<metricTensorDetSqrt */
2092  /* <<"dS_ref.data()[kb]"<<dS_ref.data()[kb]<<std::endl; */
2093  //dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];//cek need to test effect on accuracy
2094  dS = metricTensorDetSqrt*dS_ref.data()[kb];
2095  //get the metric tensor
2096  //cek todo use symmetry
2097  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
2098  ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
2099 
2100  eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2101  eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2102 
2103  //compute shape and solution information
2104  //shape
2105  /* ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace); */
2106  ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
2107  //cek hack use trial ck.gradTrialFromRef(&vel_grad_test_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_test_trace);
2108  //solution and gradients
2109  /* ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_ext); */
2110  p_ext = ebqe_p.data()[ebNE_kb];
2111  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);
2112  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);
2113  /* ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext); */
2114  /* ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext); */
2115  for (int I=0;I<nSpace;I++)
2116  grad_p_ext[I] = ebqe_grad_p.data()[ebNE_kb_nSpace + I];
2117  ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
2118  ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
2119  /* ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext); */
2120  //precalculate test function products with integration weights
2121  for (int j=0;j<nDOF_trial_element;j++)
2122  {
2123  /* p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
2124  vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2125  for (int I=0;I<nSpace;I++)
2126  vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
2127  }
2128  bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
2129  //note, our convention is that bc values at moving boundaries are relative to boundary velocity so we add it here
2130  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*(ebqe_bc_u_ext.data()[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
2131  bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*(ebqe_bc_v_ext.data()[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
2132  /* bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*(ebqe_bc_w_ext.data()[ebNE_kb] + MOVING_DOMAIN*zt_ext) + (1-isDOFBoundary_w.data()[ebNE_kb])*w_ext; */
2133  //VRANS
2134  vos_ext = ebqe_vos_ext.data()[ebNE_kb];//sed fraction - gco check
2135 
2136  //
2137  //calculate the pde coefficients using the solution and the boundary values for the solution
2138  //
2139  double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.); //not interested in saving boundary eddy viscosity for now
2140  evaluateCoefficients(eps_rho,
2141  eps_mu,
2142  sigma,
2143  rho_0,
2144  nu_0,
2145  rho_1,
2146  nu_1,
2147  rho_s,
2148  elementDiameter.data()[eN],
2149  smagorinskyConstant,
2150  turbulenceClosureModel,
2151  g.data(),
2152  useVF,
2153  ebqe_vf_ext.data()[ebNE_kb],
2154  ebqe_phi_ext.data()[ebNE_kb],
2155  &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
2156  ebqe_kappa_phi_ext.data()[ebNE_kb],
2157  //VRANS
2158  vos_ext,
2159  //
2160  p_ext,
2161  grad_p_ext,
2162  grad_u_ext,
2163  grad_v_ext,
2164  grad_w_ext,
2165  u_ext,
2166  v_ext,
2167  w_ext,
2168  ebqe_velocity_star.data()[ebNE_kb_nSpace+0],
2169  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],
2170  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],//hack,not used
2171  eddy_viscosity_ext,
2172  mom_u_acc_ext,
2173  dmom_u_acc_u_ext,
2174  mom_v_acc_ext,
2175  dmom_v_acc_v_ext,
2176  mom_w_acc_ext,
2177  dmom_w_acc_w_ext,
2178  mass_adv_ext,
2179  dmass_adv_u_ext,
2180  dmass_adv_v_ext,
2181  dmass_adv_w_ext,
2182  mom_u_adv_ext,
2183  dmom_u_adv_u_ext,
2184  dmom_u_adv_v_ext,
2185  dmom_u_adv_w_ext,
2186  mom_v_adv_ext,
2187  dmom_v_adv_u_ext,
2188  dmom_v_adv_v_ext,
2189  dmom_v_adv_w_ext,
2190  mom_w_adv_ext,
2191  dmom_w_adv_u_ext,
2192  dmom_w_adv_v_ext,
2193  dmom_w_adv_w_ext,
2194  mom_uu_diff_ten_ext,
2195  mom_vv_diff_ten_ext,
2196  mom_ww_diff_ten_ext,
2197  mom_uv_diff_ten_ext,
2198  mom_uw_diff_ten_ext,
2199  mom_vu_diff_ten_ext,
2200  mom_vw_diff_ten_ext,
2201  mom_wu_diff_ten_ext,
2202  mom_wv_diff_ten_ext,
2203  mom_u_source_ext,
2204  mom_v_source_ext,
2205  mom_w_source_ext,
2206  mom_u_ham_ext,
2207  dmom_u_ham_grad_p_ext,
2208  dmom_u_ham_grad_u_ext,
2209  mom_v_ham_ext,
2210  dmom_v_ham_grad_p_ext,
2211  dmom_v_ham_grad_v_ext,
2212  mom_w_ham_ext,
2213  dmom_w_ham_grad_p_ext,
2214  dmom_w_ham_grad_w_ext);
2215  evaluateCoefficients(eps_rho,
2216  eps_mu,
2217  sigma,
2218  rho_0,
2219  nu_0,
2220  rho_1,
2221  nu_1,
2222  rho_s,
2223  elementDiameter.data()[eN],
2224  smagorinskyConstant,
2225  turbulenceClosureModel,
2226  g.data(),
2227  useVF,
2228  bc_ebqe_vf_ext.data()[ebNE_kb],
2229  bc_ebqe_phi_ext.data()[ebNE_kb],
2230  &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
2231  ebqe_kappa_phi_ext.data()[ebNE_kb],
2232  //VRANS
2233  vos_ext,
2234  //
2235  p_ext,
2236  grad_p_ext,
2237  grad_u_ext,
2238  grad_v_ext,
2239  grad_w_ext,
2240  bc_u_ext,
2241  bc_v_ext,
2242  bc_w_ext,
2243  ebqe_velocity_star.data()[ebNE_kb_nSpace+0],
2244  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],
2245  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],//hack,not used
2246  bc_eddy_viscosity_ext,
2247  bc_mom_u_acc_ext,
2248  bc_dmom_u_acc_u_ext,
2249  bc_mom_v_acc_ext,
2250  bc_dmom_v_acc_v_ext,
2251  bc_mom_w_acc_ext,
2252  bc_dmom_w_acc_w_ext,
2253  bc_mass_adv_ext,
2254  bc_dmass_adv_u_ext,
2255  bc_dmass_adv_v_ext,
2256  bc_dmass_adv_w_ext,
2257  bc_mom_u_adv_ext,
2258  bc_dmom_u_adv_u_ext,
2259  bc_dmom_u_adv_v_ext,
2260  bc_dmom_u_adv_w_ext,
2261  bc_mom_v_adv_ext,
2262  bc_dmom_v_adv_u_ext,
2263  bc_dmom_v_adv_v_ext,
2264  bc_dmom_v_adv_w_ext,
2265  bc_mom_w_adv_ext,
2266  bc_dmom_w_adv_u_ext,
2267  bc_dmom_w_adv_v_ext,
2268  bc_dmom_w_adv_w_ext,
2269  bc_mom_uu_diff_ten_ext,
2270  bc_mom_vv_diff_ten_ext,
2271  bc_mom_ww_diff_ten_ext,
2272  bc_mom_uv_diff_ten_ext,
2273  bc_mom_uw_diff_ten_ext,
2274  bc_mom_vu_diff_ten_ext,
2275  bc_mom_vw_diff_ten_ext,
2276  bc_mom_wu_diff_ten_ext,
2277  bc_mom_wv_diff_ten_ext,
2278  bc_mom_u_source_ext,
2279  bc_mom_v_source_ext,
2280  bc_mom_w_source_ext,
2281  bc_mom_u_ham_ext,
2282  bc_dmom_u_ham_grad_p_ext,
2283  bc_dmom_u_ham_grad_u_ext,
2284  bc_mom_v_ham_ext,
2285  bc_dmom_v_ham_grad_p_ext,
2286  bc_dmom_v_ham_grad_v_ext,
2287  bc_mom_w_ham_ext,
2288  bc_dmom_w_ham_grad_p_ext,
2289  bc_dmom_w_ham_grad_w_ext);
2290 
2291  //Turbulence closure model
2292 
2293  //
2294  //moving domain
2295  //
2296  mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
2297  mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
2298  /* mom_u_adv_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*zt_ext; */
2299  dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
2300  dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
2301  /* dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext; */
2302 
2303  mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
2304  mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
2305  /* mom_v_adv_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*zt_ext; */
2306  dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
2307  dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
2308  /* dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext; */
2309 
2310  /* mom_w_adv_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*xt_ext; */
2311  /* mom_w_adv_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*yt_ext; */
2312  /* mom_w_adv_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*zt_ext; */
2313  /* dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext; */
2314  /* dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext; */
2315  /* dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext; */
2316 
2317  //bc's
2318  // mql. CHECK.
2319  bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
2320  bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
2321  /* bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*zt_ext; */
2322 
2323  bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
2324  bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
2325  /* bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*zt_ext; */
2326 
2327  /* bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*bc_mom_w_acc_ext*xt_ext; */
2328  /* bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*bc_mom_w_acc_ext*yt_ext; */
2329  /* bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*bc_mom_w_acc_ext*zt_ext; */
2330  //
2331  //calculate the numerical fluxes
2332  //
2333  ck.calculateGScale(G,normal,h_penalty);
2334  penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
2335  exteriorNumericalAdvectiveFlux(isDOFBoundary_p.data()[ebNE_kb],
2336  isDOFBoundary_u.data()[ebNE_kb],
2337  isDOFBoundary_v.data()[ebNE_kb],
2338  isDOFBoundary_w.data()[ebNE_kb],
2339  isAdvectiveFluxBoundary_p.data()[ebNE_kb],
2340  isAdvectiveFluxBoundary_u.data()[ebNE_kb],
2341  isAdvectiveFluxBoundary_v.data()[ebNE_kb],
2342  isAdvectiveFluxBoundary_w.data()[ebNE_kb],
2343  dmom_u_ham_grad_p_ext[0],//=1/rho,
2344  bc_dmom_u_ham_grad_p_ext[0],//=1/bc_rho,
2345  normal,
2346  dmom_u_acc_u_ext,
2347  bc_p_ext,
2348  bc_u_ext,
2349  bc_v_ext,
2350  bc_w_ext,
2351  bc_mass_adv_ext,
2352  bc_mom_u_adv_ext,
2353  bc_mom_v_adv_ext,
2354  bc_mom_w_adv_ext,
2355  ebqe_bc_flux_mass_ext.data()[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),//BC is relative mass flux
2356  ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
2357  ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
2358  ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
2359  p_ext,
2360  u_ext,
2361  v_ext,
2362  w_ext,
2363  mass_adv_ext,
2364  mom_u_adv_ext,
2365  mom_v_adv_ext,
2366  mom_w_adv_ext,
2367  dmass_adv_u_ext,
2368  dmass_adv_v_ext,
2369  dmass_adv_w_ext,
2370  dmom_u_adv_p_ext,
2371  dmom_u_adv_u_ext,
2372  dmom_u_adv_v_ext,
2373  dmom_u_adv_w_ext,
2374  dmom_v_adv_p_ext,
2375  dmom_v_adv_u_ext,
2376  dmom_v_adv_v_ext,
2377  dmom_v_adv_w_ext,
2378  dmom_w_adv_p_ext,
2379  dmom_w_adv_u_ext,
2380  dmom_w_adv_v_ext,
2381  dmom_w_adv_w_ext,
2382  flux_mass_ext,
2383  flux_mom_u_adv_ext,
2384  flux_mom_v_adv_ext,
2385  flux_mom_w_adv_ext,
2386  &ebqe_velocity_star.data()[ebNE_kb_nSpace],
2387  &ebqe_velocity.data()[ebNE_kb_nSpace]);
2388  exteriorNumericalDiffusiveFlux(eps_rho,
2389  ebqe_phi_ext.data()[ebNE_kb],
2390  sdInfo_u_u_rowptr.data(),
2391  sdInfo_u_u_colind.data(),
2392  isDOFBoundary_u.data()[ebNE_kb],
2393  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2394  normal,
2395  bc_mom_uu_diff_ten_ext,
2396  bc_u_ext,
2397  ebqe_bc_flux_u_diff_ext.data()[ebNE_kb],
2398  mom_uu_diff_ten_ext,
2399  grad_u_ext,
2400  u_ext,
2401  penalty,//ebqe_penalty_ext.data()[ebNE_kb],
2402  flux_mom_uu_diff_ext);
2403  exteriorNumericalDiffusiveFlux(eps_rho,
2404  ebqe_phi_ext.data()[ebNE_kb],
2405  sdInfo_u_v_rowptr.data(),
2406  sdInfo_u_v_colind.data(),
2407  isDOFBoundary_v.data()[ebNE_kb],
2408  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2409  normal,
2410  bc_mom_uv_diff_ten_ext,
2411  bc_v_ext,
2412  0.0,//assume all of the flux gets applied in diagonal component
2413  mom_uv_diff_ten_ext,
2414  grad_v_ext,
2415  v_ext,
2416  penalty,//ebqe_penalty_ext.data()[ebNE_kb],
2417  flux_mom_uv_diff_ext);
2418  /* exteriorNumericalDiffusiveFlux(eps_rho, */
2419  /* ebqe_phi_ext.data()[ebNE_kb], */
2420  /* sdInfo_u_w_rowptr.data(), */
2421  /* sdInfo_u_w_colind.data(), */
2422  /* isDOFBoundary_w.data()[ebNE_kb], */
2423  /* isDiffusiveFluxBoundary_u.data()[ebNE_kb], */
2424  /* normal, */
2425  /* bc_mom_uw_diff_ten_ext, */
2426  /* bc_w_ext, */
2427  /* 0.0,//see above */
2428  /* mom_uw_diff_ten_ext, */
2429  /* grad_w_ext, */
2430  /* w_ext, */
2431  /* penalty,//ebqe_penalty_ext.data()[ebNE_kb], */
2432  /* flux_mom_uw_diff_ext); */
2433  exteriorNumericalDiffusiveFlux(eps_rho,
2434  ebqe_phi_ext.data()[ebNE_kb],
2435  sdInfo_v_u_rowptr.data(),
2436  sdInfo_v_u_colind.data(),
2437  isDOFBoundary_u.data()[ebNE_kb],
2438  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2439  normal,
2440  bc_mom_vu_diff_ten_ext,
2441  bc_u_ext,
2442  0.0,//see above
2443  mom_vu_diff_ten_ext,
2444  grad_u_ext,
2445  u_ext,
2446  penalty,//ebqe_penalty_ext.data()[ebNE_kb],
2447  flux_mom_vu_diff_ext);
2448  exteriorNumericalDiffusiveFlux(eps_rho,
2449  ebqe_phi_ext.data()[ebNE_kb],
2450  sdInfo_v_v_rowptr.data(),
2451  sdInfo_v_v_colind.data(),
2452  isDOFBoundary_v.data()[ebNE_kb],
2453  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2454  normal,
2455  bc_mom_vv_diff_ten_ext,
2456  bc_v_ext,
2457  ebqe_bc_flux_v_diff_ext.data()[ebNE_kb],
2458  mom_vv_diff_ten_ext,
2459  grad_v_ext,
2460  v_ext,
2461  penalty,//ebqe_penalty_ext.data()[ebNE_kb],
2462  flux_mom_vv_diff_ext);
2463  /* exteriorNumericalDiffusiveFlux(eps_rho, */
2464  /* ebqe_phi_ext.data()[ebNE_kb], */
2465  /* sdInfo_v_w_rowptr.data(), */
2466  /* sdInfo_v_w_colind.data(), */
2467  /* isDOFBoundary_w.data()[ebNE_kb], */
2468  /* isDiffusiveFluxBoundary_v.data()[ebNE_kb], */
2469  /* normal, */
2470  /* bc_mom_vw_diff_ten_ext, */
2471  /* bc_w_ext, */
2472  /* 0.0,//see above */
2473  /* mom_vw_diff_ten_ext, */
2474  /* grad_w_ext, */
2475  /* w_ext, */
2476  /* penalty,//ebqe_penalty_ext.data()[ebNE_kb], */
2477  /* flux_mom_vw_diff_ext); */
2478  /* exteriorNumericalDiffusiveFlux(eps_rho, */
2479  /* ebqe_phi_ext.data()[ebNE_kb], */
2480  /* sdInfo_w_u_rowptr.data(), */
2481  /* sdInfo_w_u_colind.data(), */
2482  /* isDOFBoundary_u.data()[ebNE_kb], */
2483  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
2484  /* normal, */
2485  /* bc_mom_wu_diff_ten_ext, */
2486  /* bc_u_ext, */
2487  /* 0.0,//see above */
2488  /* mom_wu_diff_ten_ext, */
2489  /* grad_u_ext, */
2490  /* u_ext, */
2491  /* penalty,//ebqe_penalty_ext.data()[ebNE_kb], */
2492  /* flux_mom_wu_diff_ext); */
2493  /* exteriorNumericalDiffusiveFlux(eps_rho, */
2494  /* ebqe_phi_ext.data()[ebNE_kb], */
2495  /* sdInfo_w_v_rowptr.data(), */
2496  /* sdInfo_w_v_colind.data(), */
2497  /* isDOFBoundary_v.data()[ebNE_kb], */
2498  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
2499  /* normal, */
2500  /* bc_mom_wv_diff_ten_ext, */
2501  /* bc_v_ext, */
2502  /* 0.0,//see above */
2503  /* mom_wv_diff_ten_ext, */
2504  /* grad_v_ext, */
2505  /* v_ext, */
2506  /* penalty,//ebqe_penalty_ext.data()[ebNE_kb], */
2507  /* flux_mom_wv_diff_ext); */
2508  /* exteriorNumericalDiffusiveFlux(eps_rho, */
2509  /* ebqe_phi_ext.data()[ebNE_kb], */
2510  /* sdInfo_w_w_rowptr.data(), */
2511  /* sdInfo_w_w_colind.data(), */
2512  /* isDOFBoundary_w.data()[ebNE_kb], */
2513  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
2514  /* normal, */
2515  /* bc_mom_ww_diff_ten_ext, */
2516  /* bc_w_ext, */
2517  /* ebqe_bc_flux_w_diff_ext.data()[ebNE_kb], */
2518  /* mom_ww_diff_ten_ext, */
2519  /* grad_w_ext, */
2520  /* w_ext, */
2521  /* penalty,//ebqe_penalty_ext.data()[ebNE_kb], */
2522  /* flux_mom_ww_diff_ext); */
2523  flux.data()[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
2524  /* std::cout<<"external u,v,u_n " */
2525  /* <<ebqe_velocity.data()[ebNE_kb_nSpace+0]<<'\t' */
2526  /* <<ebqe_velocity.data()[ebNE_kb_nSpace+1]<<'\t' */
2527  /* <<flux.data()[ebN*nQuadraturePoints_elementBoundary+kb]<<std::endl; */
2528  //
2529  //integrate the net force and moment on flagged boundaries
2530  //
2531  if (ebN < nElementBoundaries_owned)
2532  {
2533  force_v_x = (flux_mom_u_adv_ext + flux_mom_uu_diff_ext + flux_mom_uv_diff_ext + flux_mom_uw_diff_ext)/dmom_u_ham_grad_p_ext[0];//same as *rho
2534  force_v_y = (flux_mom_v_adv_ext + flux_mom_vu_diff_ext + flux_mom_vv_diff_ext + flux_mom_vw_diff_ext)/dmom_u_ham_grad_p_ext[0];
2535  //force_v_z = (flux_mom_wu_diff_ext + flux_mom_wv_diff_ext + flux_mom_ww_diff_ext)/dmom_u_ham_grad_p_ext[0];
2536 
2537  force_p_x = p_ext*normal[0];
2538  force_p_y = p_ext*normal[1];
2539  //force_p_z = p_ext*normal[2];
2540 
2541  force_x = force_p_x + force_v_x;
2542  force_y = force_p_y + force_v_y;
2543  //force_z = force_p_z + force_v_z;
2544 
2545  r_x = x_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+0];
2546  r_y = y_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+1];
2547  //r_z = z_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+2];
2548 
2549  wettedAreas.data()[boundaryFlags.data()[ebN]] += dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2550 
2551  netForces_p.data()[3*boundaryFlags.data()[ebN]+0] += force_p_x*dS;
2552  netForces_p.data()[3*boundaryFlags.data()[ebN]+1] += force_p_y*dS;
2553  //netForces_p.data()[3*boundaryFlags.data()[ebN]+2] += force_p_z*dS;
2554 
2555  netForces_v.data()[3*boundaryFlags.data()[ebN]+0] += force_v_x*dS;
2556  netForces_v.data()[3*boundaryFlags.data()[ebN]+1] += force_v_y*dS;
2557  //netForces_v.data()[3*boundaryFlags.data()[ebN]+2] += force_v_z*dS;
2558 
2559  //netMoments.data()[3*boundaryFlags.data()[ebN]+0] += (r_y*force_z - r_z*force_y)*dS;
2560  //netMoments.data()[3*boundaryFlags.data()[ebN]+1] += (r_z*force_x - r_x*force_z)*dS;
2561  netMoments.data()[3*boundaryFlags.data()[ebN]+2] += (r_x*force_y - r_y*force_x)*dS;
2562  }
2563  //
2564  //update residuals
2565  //
2566  for (int i=0;i<nDOF_test_element;i++)
2567  {
2568  /* elementResidual_mesh[i] -= ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),p_test_dS[i]); */
2569  /* elementResidual_p.data()[i] += ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]); */
2570  /* elementResidual_p.data()[i] -= DM*ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),p_test_dS[i]); */
2571  /* globalConservationError += ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]); */
2572  elementResidual_u[i] +=
2573  ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
2574  ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
2575  ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
2576  ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2577  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2578  eb_adjoint_sigma,
2579  u_ext,
2580  bc_u_ext,
2581  normal,
2582  sdInfo_u_u_rowptr.data(),
2583  sdInfo_u_u_colind.data(),
2584  mom_uu_diff_ten_ext,
2585  &vel_grad_test_dS[i*nSpace])+
2586  ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2587  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2588  eb_adjoint_sigma,
2589  v_ext,
2590  bc_v_ext,
2591  normal,
2592  sdInfo_u_v_rowptr.data(),
2593  sdInfo_u_v_colind.data(),
2594  mom_uv_diff_ten_ext,
2595  &vel_grad_test_dS[i*nSpace]);//+
2596  /* ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb], */
2597  /* isDiffusiveFluxBoundary_u.data()[ebNE_kb], */
2598  /* eb_adjoint_sigma, */
2599  /* w_ext, */
2600  /* bc_w_ext, */
2601  /* normal, */
2602  /* sdInfo_u_w_rowptr.data(), */
2603  /* sdInfo_u_w_colind.data(), */
2604  /* mom_uw_diff_ten_ext, */
2605  /* &vel_grad_test_dS[i*nSpace]); */
2606  elementResidual_v[i] += ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
2607  ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
2608  ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
2609  ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2610  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2611  eb_adjoint_sigma,
2612  u_ext,
2613  bc_u_ext,
2614  normal,
2615  sdInfo_v_u_rowptr.data(),
2616  sdInfo_v_u_colind.data(),
2617  mom_vu_diff_ten_ext,
2618  &vel_grad_test_dS[i*nSpace])+
2619  ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2620  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2621  eb_adjoint_sigma,
2622  v_ext,
2623  bc_v_ext,
2624  normal,
2625  sdInfo_v_v_rowptr.data(),
2626  sdInfo_v_v_colind.data(),
2627  mom_vv_diff_ten_ext,
2628  &vel_grad_test_dS[i*nSpace]);//+
2629  /* ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb], */
2630  /* isDiffusiveFluxBoundary_v.data()[ebNE_kb], */
2631  /* eb_adjoint_sigma, */
2632  /* w_ext, */
2633  /* bc_w_ext, */
2634  /* normal, */
2635  /* sdInfo_v_w_rowptr.data(), */
2636  /* sdInfo_v_w_colind.data(), */
2637  /* mom_vw_diff_ten_ext, */
2638  /* &vel_grad_test_dS[i*nSpace]); */
2639 
2640  /* elementResidual_w[i] += ck.ExteriorElementBoundaryFlux(flux_mom_w_adv_ext,vel_test_dS[i]) + */
2641  /* ck.ExteriorElementBoundaryFlux(flux_mom_wu_diff_ext,vel_test_dS[i])+ */
2642  /* ck.ExteriorElementBoundaryFlux(flux_mom_wv_diff_ext,vel_test_dS[i])+ */
2643  /* ck.ExteriorElementBoundaryFlux(flux_mom_ww_diff_ext,vel_test_dS[i])+ */
2644  /* ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb], */
2645  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
2646  /* eb_adjoint_sigma, */
2647  /* u_ext, */
2648  /* bc_u_ext, */
2649  /* normal, */
2650  /* sdInfo_w_u_rowptr.data(), */
2651  /* sdInfo_w_u_colind.data(), */
2652  /* mom_wu_diff_ten_ext, */
2653  /* &vel_grad_test_dS[i*nSpace])+ */
2654  /* ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb], */
2655  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
2656  /* eb_adjoint_sigma, */
2657  /* v_ext, */
2658  /* bc_v_ext, */
2659  /* normal, */
2660  /* sdInfo_w_v_rowptr.data(), */
2661  /* sdInfo_w_v_colind.data(), */
2662  /* mom_wv_diff_ten_ext, */
2663  /* &vel_grad_test_dS[i*nSpace])+ */
2664  /* ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb], */
2665  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
2666  /* eb_adjoint_sigma, */
2667  /* w_ext, */
2668  /* bc_w_ext, */
2669  /* normal, */
2670  /* sdInfo_w_w_rowptr.data(), */
2671  /* sdInfo_w_w_colind.data(), */
2672  /* mom_ww_diff_ten_ext, */
2673  /* &vel_grad_test_dS[i*nSpace]); */
2674  }//i
2675  }//kb
2676  //
2677  //update the element and global residual storage
2678  //
2679  for (int i=0;i<nDOF_test_element;i++)
2680  {
2681  int eN_i = eN*nDOF_test_element+i;
2682 
2683  /* elementResidual_p_save[eN_i] += elementResidual_p.data()[i]; */
2684  /* mesh_volume_conservation_weak += elementResidual_mesh[i]; */
2685  /* globalResidual.data()[offset_p+stride_p*p_l2g.data()[eN_i]]+=elementResidual_p.data()[i]; */
2686  globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
2687  globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
2688  /* globalResidual.data()[offset_w+stride_w*vel_l2g.data()[eN_i]]+=elementResidual_w[i]; */
2689  }//i
2690  }//ebNE
2691  /* std::cout<<"mesh volume conservation = "<<mesh_volume_conservation<<std::endl; */
2692  /* std::cout<<"mesh volume conservation weak = "<<mesh_volume_conservation_weak<<std::endl; */
2693  /* std::cout<<"mesh volume conservation err max= "<<mesh_volume_conservation_err_max<<std::endl; */
2694  /* std::cout<<"mesh volume conservation err max weak = "<<mesh_volume_conservation_err_max_weak<<std::endl; */
2695  }
2696 
2697  void calculateJacobian(arguments_dict& args)
2698  {
2699  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
2700  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2701  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
2702  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
2703  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
2704  double PSTAB = args.scalar<double>("PSTAB");
2705  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
2706  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
2707  xt::pyarray<double>& p_trial_ref = args.array<double>("p_trial_ref");
2708  xt::pyarray<double>& p_grad_trial_ref = args.array<double>("p_grad_trial_ref");
2709  xt::pyarray<double>& p_test_ref = args.array<double>("p_test_ref");
2710  xt::pyarray<double>& p_grad_test_ref = args.array<double>("p_grad_test_ref");
2711  xt::pyarray<double>& q_p = args.array<double>("q_p");
2712  xt::pyarray<double>& q_grad_p = args.array<double>("q_grad_p");
2713  xt::pyarray<double>& ebqe_p = args.array<double>("ebqe_p");
2714  xt::pyarray<double>& ebqe_grad_p = args.array<double>("ebqe_grad_p");
2715  xt::pyarray<double>& vel_trial_ref = args.array<double>("vel_trial_ref");
2716  xt::pyarray<double>& vel_grad_trial_ref = args.array<double>("vel_grad_trial_ref");
2717  xt::pyarray<double>& vel_test_ref = args.array<double>("vel_test_ref");
2718  xt::pyarray<double>& vel_grad_test_ref = args.array<double>("vel_grad_test_ref");
2719  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2720  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2721  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
2722  xt::pyarray<double>& p_trial_trace_ref = args.array<double>("p_trial_trace_ref");
2723  xt::pyarray<double>& p_grad_trial_trace_ref = args.array<double>("p_grad_trial_trace_ref");
2724  xt::pyarray<double>& p_test_trace_ref = args.array<double>("p_test_trace_ref");
2725  xt::pyarray<double>& p_grad_test_trace_ref = args.array<double>("p_grad_test_trace_ref");
2726  xt::pyarray<double>& vel_trial_trace_ref = args.array<double>("vel_trial_trace_ref");
2727  xt::pyarray<double>& vel_grad_trial_trace_ref = args.array<double>("vel_grad_trial_trace_ref");
2728  xt::pyarray<double>& vel_test_trace_ref = args.array<double>("vel_test_trace_ref");
2729  xt::pyarray<double>& vel_grad_test_trace_ref = args.array<double>("vel_grad_test_trace_ref");
2730  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
2731  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
2732  double eb_adjoint_sigma = args.scalar<double>("eb_adjoint_sigma");
2733  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
2734  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
2735  double hFactor = args.scalar<double>("hFactor");
2736  int nElements_global = args.scalar<int>("nElements_global");
2737  double useRBLES = args.scalar<double>("useRBLES");
2738  double useMetrics = args.scalar<double>("useMetrics");
2739  double alphaBDF = args.scalar<double>("alphaBDF");
2740  double epsFact_rho = args.scalar<double>("epsFact_rho");
2741  double epsFact_mu = args.scalar<double>("epsFact_mu");
2742  double sigma = args.scalar<double>("sigma");
2743  double rho_0 = args.scalar<double>("rho_0");
2744  double nu_0 = args.scalar<double>("nu_0");
2745  double rho_1 = args.scalar<double>("rho_1");
2746  double nu_1 = args.scalar<double>("nu_1");
2747  double rho_s = args.scalar<double>("rho_s");
2748  double smagorinskyConstant = args.scalar<double>("smagorinskyConstant");
2749  int turbulenceClosureModel = args.scalar<int>("turbulenceClosureModel");
2750  double Ct_sge = args.scalar<double>("Ct_sge");
2751  double Cd_sge = args.scalar<double>("Cd_sge");
2752  double C_dg = args.scalar<double>("C_dg");
2753  double C_b = args.scalar<double>("C_b");
2754  const xt::pyarray<double>& eps_solid = args.array<double>("eps_solid");
2755  const xt::pyarray<double>& q_velocity_fluid = args.array<double>("q_velocity_fluid");
2756  const xt::pyarray<double>& q_velocityStar_fluid = args.array<double>("q_velocityStar_fluid");
2757  const xt::pyarray<double>& q_vos = args.array<double>("q_vos");
2758  const xt::pyarray<double>& q_dvos_dt = args.array<double>("q_dvos_dt");
2759  const xt::pyarray<double>& q_grad_vos = args.array<double>("q_grad_vos");
2760  const xt::pyarray<double>& q_dragAlpha = args.array<double>("q_dragAlpha");
2761  const xt::pyarray<double>& q_dragBeta = args.array<double>("q_dragBeta");
2762  const xt::pyarray<double>& q_mass_source = args.array<double>("q_mass_source");
2763  const xt::pyarray<double>& q_turb_var_0 = args.array<double>("q_turb_var_0");
2764  const xt::pyarray<double>& q_turb_var_1 = args.array<double>("q_turb_var_1");
2765  const xt::pyarray<double>& q_turb_var_grad_0 = args.array<double>("q_turb_var_grad_0");
2766  xt::pyarray<int>& p_l2g = args.array<int>("p_l2g");
2767  xt::pyarray<int>& vel_l2g = args.array<int>("vel_l2g");
2768  xt::pyarray<double>& p_dof = args.array<double>("p_dof");
2769  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
2770  xt::pyarray<double>& v_dof = args.array<double>("v_dof");
2771  xt::pyarray<double>& w_dof = args.array<double>("w_dof");
2772  xt::pyarray<double>& g = args.array<double>("g");
2773  const double useVF = args.scalar<double>("useVF");
2774  xt::pyarray<double>& vf = args.array<double>("vf");
2775  xt::pyarray<double>& phi = args.array<double>("phi");
2776  xt::pyarray<double>& normal_phi = args.array<double>("normal_phi");
2777  xt::pyarray<double>& kappa_phi = args.array<double>("kappa_phi");
2778  xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.array<double>("q_mom_u_acc_beta_bdf");
2779  xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.array<double>("q_mom_v_acc_beta_bdf");
2780  xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.array<double>("q_mom_w_acc_beta_bdf");
2781  xt::pyarray<double>& q_dV = args.array<double>("q_dV");
2782  xt::pyarray<double>& q_dV_last = args.array<double>("q_dV_last");
2783  xt::pyarray<double>& q_velocity_sge = args.array<double>("q_velocity_sge");
2784  xt::pyarray<double>& ebqe_velocity_star = args.array<double>("ebqe_velocity_star");
2785  xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
2786  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
2787  xt::pyarray<double>& q_numDiff_v_last = args.array<double>("q_numDiff_v_last");
2788  xt::pyarray<double>& q_numDiff_w_last = args.array<double>("q_numDiff_w_last");
2789  xt::pyarray<int>& sdInfo_u_u_rowptr = args.array<int>("sdInfo_u_u_rowptr");
2790  xt::pyarray<int>& sdInfo_u_u_colind = args.array<int>("sdInfo_u_u_colind");
2791  xt::pyarray<int>& sdInfo_u_v_rowptr = args.array<int>("sdInfo_u_v_rowptr");
2792  xt::pyarray<int>& sdInfo_u_v_colind = args.array<int>("sdInfo_u_v_colind");
2793  xt::pyarray<int>& sdInfo_u_w_rowptr = args.array<int>("sdInfo_u_w_rowptr");
2794  xt::pyarray<int>& sdInfo_u_w_colind = args.array<int>("sdInfo_u_w_colind");
2795  xt::pyarray<int>& sdInfo_v_v_rowptr = args.array<int>("sdInfo_v_v_rowptr");
2796  xt::pyarray<int>& sdInfo_v_v_colind = args.array<int>("sdInfo_v_v_colind");
2797  xt::pyarray<int>& sdInfo_v_u_rowptr = args.array<int>("sdInfo_v_u_rowptr");
2798  xt::pyarray<int>& sdInfo_v_u_colind = args.array<int>("sdInfo_v_u_colind");
2799  xt::pyarray<int>& sdInfo_v_w_rowptr = args.array<int>("sdInfo_v_w_rowptr");
2800  xt::pyarray<int>& sdInfo_v_w_colind = args.array<int>("sdInfo_v_w_colind");
2801  xt::pyarray<int>& sdInfo_w_w_rowptr = args.array<int>("sdInfo_w_w_rowptr");
2802  xt::pyarray<int>& sdInfo_w_w_colind = args.array<int>("sdInfo_w_w_colind");
2803  xt::pyarray<int>& sdInfo_w_u_rowptr = args.array<int>("sdInfo_w_u_rowptr");
2804  xt::pyarray<int>& sdInfo_w_u_colind = args.array<int>("sdInfo_w_u_colind");
2805  xt::pyarray<int>& sdInfo_w_v_rowptr = args.array<int>("sdInfo_w_v_rowptr");
2806  xt::pyarray<int>& sdInfo_w_v_colind = args.array<int>("sdInfo_w_v_colind");
2807  xt::pyarray<int>& csrRowIndeces_p_p = args.array<int>("csrRowIndeces_p_p");
2808  xt::pyarray<int>& csrColumnOffsets_p_p = args.array<int>("csrColumnOffsets_p_p");
2809  xt::pyarray<int>& csrRowIndeces_p_u = args.array<int>("csrRowIndeces_p_u");
2810  xt::pyarray<int>& csrColumnOffsets_p_u = args.array<int>("csrColumnOffsets_p_u");
2811  xt::pyarray<int>& csrRowIndeces_p_v = args.array<int>("csrRowIndeces_p_v");
2812  xt::pyarray<int>& csrColumnOffsets_p_v = args.array<int>("csrColumnOffsets_p_v");
2813  xt::pyarray<int>& csrRowIndeces_p_w = args.array<int>("csrRowIndeces_p_w");
2814  xt::pyarray<int>& csrColumnOffsets_p_w = args.array<int>("csrColumnOffsets_p_w");
2815  xt::pyarray<int>& csrRowIndeces_u_p = args.array<int>("csrRowIndeces_u_p");
2816  xt::pyarray<int>& csrColumnOffsets_u_p = args.array<int>("csrColumnOffsets_u_p");
2817  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
2818  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
2819  xt::pyarray<int>& csrRowIndeces_u_v = args.array<int>("csrRowIndeces_u_v");
2820  xt::pyarray<int>& csrColumnOffsets_u_v = args.array<int>("csrColumnOffsets_u_v");
2821  xt::pyarray<int>& csrRowIndeces_u_w = args.array<int>("csrRowIndeces_u_w");
2822  xt::pyarray<int>& csrColumnOffsets_u_w = args.array<int>("csrColumnOffsets_u_w");
2823  xt::pyarray<int>& csrRowIndeces_v_p = args.array<int>("csrRowIndeces_v_p");
2824  xt::pyarray<int>& csrColumnOffsets_v_p = args.array<int>("csrColumnOffsets_v_p");
2825  xt::pyarray<int>& csrRowIndeces_v_u = args.array<int>("csrRowIndeces_v_u");
2826  xt::pyarray<int>& csrColumnOffsets_v_u = args.array<int>("csrColumnOffsets_v_u");
2827  xt::pyarray<int>& csrRowIndeces_v_v = args.array<int>("csrRowIndeces_v_v");
2828  xt::pyarray<int>& csrColumnOffsets_v_v = args.array<int>("csrColumnOffsets_v_v");
2829  xt::pyarray<int>& csrRowIndeces_v_w = args.array<int>("csrRowIndeces_v_w");
2830  xt::pyarray<int>& csrColumnOffsets_v_w = args.array<int>("csrColumnOffsets_v_w");
2831  xt::pyarray<int>& csrRowIndeces_w_p = args.array<int>("csrRowIndeces_w_p");
2832  xt::pyarray<int>& csrColumnOffsets_w_p = args.array<int>("csrColumnOffsets_w_p");
2833  xt::pyarray<int>& csrRowIndeces_w_u = args.array<int>("csrRowIndeces_w_u");
2834  xt::pyarray<int>& csrColumnOffsets_w_u = args.array<int>("csrColumnOffsets_w_u");
2835  xt::pyarray<int>& csrRowIndeces_w_v = args.array<int>("csrRowIndeces_w_v");
2836  xt::pyarray<int>& csrColumnOffsets_w_v = args.array<int>("csrColumnOffsets_w_v");
2837  xt::pyarray<int>& csrRowIndeces_w_w = args.array<int>("csrRowIndeces_w_w");
2838  xt::pyarray<int>& csrColumnOffsets_w_w = args.array<int>("csrColumnOffsets_w_w");
2839  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
2840  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2841  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2842  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2843  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2844  xt::pyarray<double>& ebqe_vf_ext = args.array<double>("ebqe_vf_ext");
2845  xt::pyarray<double>& bc_ebqe_vf_ext = args.array<double>("bc_ebqe_vf_ext");
2846  xt::pyarray<double>& ebqe_phi_ext = args.array<double>("ebqe_phi_ext");
2847  xt::pyarray<double>& bc_ebqe_phi_ext = args.array<double>("bc_ebqe_phi_ext");
2848  xt::pyarray<double>& ebqe_normal_phi_ext = args.array<double>("ebqe_normal_phi_ext");
2849  xt::pyarray<double>& ebqe_kappa_phi_ext = args.array<double>("ebqe_kappa_phi_ext");
2850  const xt::pyarray<double>& ebqe_vos_ext = args.array<double>("ebqe_vos_ext");
2851  const xt::pyarray<double>& ebqe_turb_var_0 = args.array<double>("ebqe_turb_var_0");
2852  const xt::pyarray<double>& ebqe_turb_var_1 = args.array<double>("ebqe_turb_var_1");
2853  xt::pyarray<int>& isDOFBoundary_p = args.array<int>("isDOFBoundary_p");
2854  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
2855  xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
2856  xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
2857  xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.array<int>("isAdvectiveFluxBoundary_p");
2858  xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
2859  xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.array<int>("isAdvectiveFluxBoundary_v");
2860  xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.array<int>("isAdvectiveFluxBoundary_w");
2861  xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
2862  xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.array<int>("isDiffusiveFluxBoundary_v");
2863  xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.array<int>("isDiffusiveFluxBoundary_w");
2864  xt::pyarray<double>& ebqe_bc_p_ext = args.array<double>("ebqe_bc_p_ext");
2865  xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.array<double>("ebqe_bc_flux_mass_ext");
2866  xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.array<double>("ebqe_bc_flux_mom_u_adv_ext");
2867  xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.array<double>("ebqe_bc_flux_mom_v_adv_ext");
2868  xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.array<double>("ebqe_bc_flux_mom_w_adv_ext");
2869  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
2870  xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.array<double>("ebqe_bc_flux_u_diff_ext");
2871  xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
2872  xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
2873  xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.array<double>("ebqe_bc_flux_v_diff_ext");
2874  xt::pyarray<double>& ebqe_bc_w_ext = args.array<double>("ebqe_bc_w_ext");
2875  xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.array<double>("ebqe_bc_flux_w_diff_ext");
2876  xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.array<int>("csrColumnOffsets_eb_p_p");
2877  xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.array<int>("csrColumnOffsets_eb_p_u");
2878  xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.array<int>("csrColumnOffsets_eb_p_v");
2879  xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.array<int>("csrColumnOffsets_eb_p_w");
2880  xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.array<int>("csrColumnOffsets_eb_u_p");
2881  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
2882  xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.array<int>("csrColumnOffsets_eb_u_v");
2883  xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.array<int>("csrColumnOffsets_eb_u_w");
2884  xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.array<int>("csrColumnOffsets_eb_v_p");
2885  xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.array<int>("csrColumnOffsets_eb_v_u");
2886  xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.array<int>("csrColumnOffsets_eb_v_v");
2887  xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.array<int>("csrColumnOffsets_eb_v_w");
2888  xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.array<int>("csrColumnOffsets_eb_w_p");
2889  xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.array<int>("csrColumnOffsets_eb_w_u");
2890  xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.array<int>("csrColumnOffsets_eb_w_v");
2891  xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.array<int>("csrColumnOffsets_eb_w_w");
2892  xt::pyarray<int>& elementFlags = args.array<int>("elementFlags");
2893  double LAG_MU_FR = args.scalar<double>("LAG_MU_FR");
2894  xt::pyarray<double>& q_mu_fr_last = args.array<double>("q_mu_fr_last");
2895  xt::pyarray<double>& q_mu_fr = args.array<double>("q_mu_fr");
2896  //
2897  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2898  //
2899  for(int eN=0;eN<nElements_global;eN++)
2900  {
2901  register double eps_rho,eps_mu;
2902 
2903  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
2904  elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
2905  elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
2906  elementJacobian_v_v[nDOF_test_element][nDOF_trial_element];
2907  for (int i=0;i<nDOF_test_element;i++)
2908  for (int j=0;j<nDOF_trial_element;j++)
2909 
2910  {
2911  elementJacobian_u_u[i][j]=0.0;
2912  elementJacobian_u_v[i][j]=0.0;
2913  elementJacobian_v_u[i][j]=0.0;
2914  elementJacobian_v_v[i][j]=0.0;
2915  }
2916 
2917  for (int k=0;k<nQuadraturePoints_element;k++)
2918  {
2919  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
2920  eN_k_nSpace = eN_k*nSpace,
2921  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
2922 
2923  //declare local storage
2924  register double p=0.0,u=0.0,v=0.0,w=0.0,
2925  grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],grad_vos[nSpace],
2926  mom_u_acc=0.0,
2927  dmom_u_acc_u=0.0,
2928  mom_v_acc=0.0,
2929  dmom_v_acc_v=0.0,
2930  mom_w_acc=0.0,
2931  dmom_w_acc_w=0.0,
2932  mass_adv[nSpace],
2933  dmass_adv_u[nSpace],
2934  dmass_adv_v[nSpace],
2935  dmass_adv_w[nSpace],
2936  mom_u_adv[nSpace],
2937  dmom_u_adv_u[nSpace],
2938  dmom_u_adv_v[nSpace],
2939  dmom_u_adv_w[nSpace],
2940  mom_v_adv[nSpace],
2941  dmom_v_adv_u[nSpace],
2942  dmom_v_adv_v[nSpace],
2943  dmom_v_adv_w[nSpace],
2944  mom_w_adv[nSpace],
2945  dmom_w_adv_u[nSpace],
2946  dmom_w_adv_v[nSpace],
2947  dmom_w_adv_w[nSpace],
2948  mom_uu_diff_ten[nSpace],
2949  mom_vv_diff_ten[nSpace],
2950  mom_ww_diff_ten[nSpace],
2951  mom_uv_diff_ten[1],
2952  mom_uw_diff_ten[1],
2953  mom_vu_diff_ten[1],
2954  mom_vw_diff_ten[1],
2955  mom_wu_diff_ten[1],
2956  mom_wv_diff_ten[1],
2957  mom_u_source=0.0,
2958  mom_v_source=0.0,
2959  mom_w_source=0.0,
2960  mom_u_ham=0.0,
2961  dmom_u_ham_grad_p[nSpace],
2962  dmom_u_ham_grad_u[nSpace],
2963  mom_v_ham=0.0,
2964  dmom_v_ham_grad_p[nSpace],
2965  dmom_v_ham_grad_v[nSpace],
2966  mom_w_ham=0.0,
2967  dmom_w_ham_grad_p[nSpace],
2968  dmom_w_ham_grad_w[nSpace],
2969  mom_u_acc_t=0.0,
2970  dmom_u_acc_u_t=0.0,
2971  mom_v_acc_t=0.0,
2972  dmom_v_acc_v_t=0.0,
2973  pdeResidual_p=0.0,
2974  pdeResidual_u=0.0,
2975  pdeResidual_v=0.0,
2976  pdeResidual_w=0.0,
2977  dpdeResidual_p_u[nDOF_trial_element],dpdeResidual_p_v[nDOF_trial_element],dpdeResidual_p_w[nDOF_trial_element],
2978  dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_trial_element],
2979  dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_trial_element],
2980  dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_trial_element],
2981  Lstar_u_u[nDOF_test_element],
2982  Lstar_v_v[nDOF_test_element],
2983  subgridError_p=0.0,
2984  subgridError_u=0.0,
2985  subgridError_v=0.0,
2986  subgridError_w=0.0,
2987  dsubgridError_p_u[nDOF_trial_element],
2988  dsubgridError_p_v[nDOF_trial_element],
2989  dsubgridError_p_w[nDOF_trial_element],
2990  dsubgridError_u_p[nDOF_trial_element],
2991  dsubgridError_u_u[nDOF_trial_element],
2992  dsubgridError_v_p[nDOF_trial_element],
2993  dsubgridError_v_v[nDOF_trial_element],
2994  dsubgridError_w_p[nDOF_trial_element],
2995  dsubgridError_w_w[nDOF_trial_element],
2996  tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
2997  tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
2998  jac[nSpace*nSpace],
2999  jacDet,
3000  jacInv[nSpace*nSpace],
3001  p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
3002  dV,
3003  vel_test_dV[nDOF_test_element],
3004  vel_grad_test_dV[nDOF_test_element*nSpace],
3005  x,y,z,xt,yt,zt,
3006  //VRANS
3007  vos,
3008  //meanGrainSize,
3009  dmom_u_source[nSpace],
3010  dmom_v_source[nSpace],
3011  dmom_w_source[nSpace],
3012  mass_source,
3013  //
3014  G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
3015  //get jacobian, etc for mapping reference element
3016  ck.calculateMapping_element(eN,
3017  k,
3018  mesh_dof.data(),
3019  mesh_l2g.data(),
3020  mesh_trial_ref.data(),
3021  mesh_grad_trial_ref.data(),
3022  jac,
3023  jacDet,
3024  jacInv,
3025  x,y,z);
3026  ck.calculateH_element(eN,
3027  k,
3028  nodeDiametersArray.data(),
3029  mesh_l2g.data(),
3030  mesh_trial_ref.data(),
3031  h_phi);
3032  ck.calculateMappingVelocity_element(eN,
3033  k,
3034  mesh_velocity_dof.data(),
3035  mesh_l2g.data(),
3036  mesh_trial_ref.data(),
3037  xt,yt,zt);
3038  //xt=0.0;yt=0.0;zt=0.0;
3039  //std::cout<<"xt "<<xt<<'\t'<<yt<<'\t'<<zt<<std::endl;
3040  //get the physical integration weight
3041  dV = fabs(jacDet)*dV_ref.data()[k];
3042  ck.calculateG(jacInv,G,G_dd_G,tr_G);
3043  //ck.calculateGScale(G,&normal_phi.data()[eN_k_nSpace],h_phi);
3044 
3045  eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3046  eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3047 
3048  //get the trial function gradients
3049  /* ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial); */
3050  ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
3051  //get the solution
3052  /* ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_ref.data()[k*nDOF_trial_element],p); */
3053  p = q_p.data()[eN_k];
3054  ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],u);
3055  ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],v);
3056  /* ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],w); */
3057  //get the solution gradients
3058  /* ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial,grad_p); */
3059  for (int I=0;I<nSpace;I++)
3060  grad_p[I] = q_grad_p.data()[eN_k_nSpace+I];
3061  for (int I=0;I<nSpace;I++)
3062  grad_vos[I] = q_grad_vos.data()[eN_k_nSpace + I];
3063  ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
3064  ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
3065  /* ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_w); */
3066  //VRANS
3067  vos = q_vos.data()[eN_k];//sed fraction - gco check
3068  //precalculate test function products with integration weights
3069  for (int j=0;j<nDOF_trial_element;j++)
3070  {
3071  /* p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV; */
3072  vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
3073  for (int I=0;I<nSpace;I++)
3074  {
3075  /* p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin */
3076  vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin}
3077  }
3078  }
3079  //cek hack
3080  double div_mesh_velocity=0.0;
3081  int NDOF_MESH_TRIAL_ELEMENT=3;
3082  for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
3083  {
3084  int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
3085  div_mesh_velocity +=
3086  mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+0]*vel_grad_trial[j*2+0] +
3087  mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+1]*vel_grad_trial[j*2+1];
3088  }
3089  div_mesh_velocity = DM3*div_mesh_velocity + (1.0-DM3)*alphaBDF*(dV-q_dV_last.data()[eN_k])/dV;
3090  //
3091 
3092  //
3093  //
3094  //calculate pde coefficients and derivatives at quadrature points
3095  //
3096  double eddy_viscosity(0.);//not really interested in saving eddy_viscosity in jacobian
3097  evaluateCoefficients(eps_rho,
3098  eps_mu,
3099  sigma,
3100  rho_0,
3101  nu_0,
3102  rho_1,
3103  nu_1,
3104  rho_s,
3105  elementDiameter.data()[eN],
3106  smagorinskyConstant,
3107  turbulenceClosureModel,
3108  g.data(),
3109  useVF,
3110  vf.data()[eN_k],
3111  phi.data()[eN_k],
3112  &normal_phi.data()[eN_k_nSpace],
3113  kappa_phi.data()[eN_k],
3114  //VRANS
3115  vos,
3116  //
3117  p,
3118  grad_p,
3119  grad_u,
3120  grad_v,
3121  grad_w,
3122  u,
3123  v,
3124  w,
3125  q_velocity_sge.data()[eN_k_nSpace+0],
3126  q_velocity_sge.data()[eN_k_nSpace+1],
3127  q_velocity_sge.data()[eN_k_nSpace+1],//hack, shouldn't be used
3128  eddy_viscosity,
3129  mom_u_acc,
3130  dmom_u_acc_u,
3131  mom_v_acc,
3132  dmom_v_acc_v,
3133  mom_w_acc,
3134  dmom_w_acc_w,
3135  mass_adv,
3136  dmass_adv_u,
3137  dmass_adv_v,
3138  dmass_adv_w,
3139  mom_u_adv,
3140  dmom_u_adv_u,
3141  dmom_u_adv_v,
3142  dmom_u_adv_w,
3143  mom_v_adv,
3144  dmom_v_adv_u,
3145  dmom_v_adv_v,
3146  dmom_v_adv_w,
3147  mom_w_adv,
3148  dmom_w_adv_u,
3149  dmom_w_adv_v,
3150  dmom_w_adv_w,
3151  mom_uu_diff_ten,
3152  mom_vv_diff_ten,
3153  mom_ww_diff_ten,
3154  mom_uv_diff_ten,
3155  mom_uw_diff_ten,
3156  mom_vu_diff_ten,
3157  mom_vw_diff_ten,
3158  mom_wu_diff_ten,
3159  mom_wv_diff_ten,
3160  mom_u_source,
3161  mom_v_source,
3162  mom_w_source,
3163  mom_u_ham,
3164  dmom_u_ham_grad_p,
3165  dmom_u_ham_grad_u,
3166  mom_v_ham,
3167  dmom_v_ham_grad_p,
3168  dmom_v_ham_grad_v,
3169  mom_w_ham,
3170  dmom_w_ham_grad_p,
3171  dmom_w_ham_grad_w);
3172  //VRANS
3173  mass_source = q_mass_source.data()[eN_k];
3174  for (int I=0;I<nSpace;I++)
3175  {
3176  dmom_u_source[I] = 0.0;
3177  dmom_v_source[I] = 0.0;
3178  dmom_w_source[I] = 0.0;
3179  }
3180  updateDarcyForchheimerTerms_Ergun(// linearDragFactor,
3181  // nonlinearDragFactor,
3182  // vos,
3183  // meanGrainSize,
3184  q_dragAlpha.data()[eN_k],
3185  q_dragBeta.data()[eN_k],
3186  eps_rho,
3187  eps_mu,
3188  rho_0,
3189  nu_0,
3190  rho_1,
3191  nu_1,
3192  eddy_viscosity,
3193  useVF,
3194  vf.data()[eN_k],
3195  phi.data()[eN_k],
3196  u,
3197  v,
3198  w,
3199  q_velocity_sge.data()[eN_k_nSpace+0],
3200  q_velocity_sge.data()[eN_k_nSpace+1],
3201  q_velocity_sge.data()[eN_k_nSpace+1],//hack, shouldn't be used
3202  eps_solid.data()[elementFlags.data()[eN]],
3203  vos,
3204  q_velocity_fluid.data()[eN_k_nSpace+0],
3205  q_velocity_fluid.data()[eN_k_nSpace+1],
3206  q_velocity_fluid.data()[eN_k_nSpace+1],
3207  q_velocityStar_fluid.data()[eN_k_nSpace+0],
3208  q_velocityStar_fluid.data()[eN_k_nSpace+1],
3209  q_velocityStar_fluid.data()[eN_k_nSpace+1],
3210  mom_u_source,
3211  mom_v_source,
3212  mom_w_source,
3213  dmom_u_source,
3214  dmom_v_source,
3215  dmom_w_source,
3216  q_grad_vos.data()[eN_k_nSpace+0],
3217  q_grad_vos.data()[eN_k_nSpace+1],
3218  q_grad_vos.data()[eN_k_nSpace+1]);
3219 
3220  double mu_fr_tmp=0.0;
3221  updateFrictionalPressure(vos,
3222  grad_vos,
3223  mom_u_source,
3224  mom_v_source,
3225  mom_w_source);
3226  updateFrictionalStress(LAG_MU_FR,
3227  q_mu_fr_last.data()[eN_k],
3228  mu_fr_tmp,
3229  vos,
3230  eps_rho,
3231  eps_mu,
3232  rho_0,
3233  nu_0,
3234  rho_1,
3235  nu_1,
3236  rho_s,
3237  useVF,
3238  vf.data()[eN_k],
3239  phi.data()[eN_k],
3240  grad_u,
3241  grad_v,
3242  grad_w,
3243  mom_uu_diff_ten,
3244  mom_uv_diff_ten,
3245  mom_uw_diff_ten,
3246  mom_vv_diff_ten,
3247  mom_vu_diff_ten,
3248  mom_vw_diff_ten,
3249  mom_ww_diff_ten,
3250  mom_wu_diff_ten,
3251  mom_wv_diff_ten);
3252  //
3253  //moving mesh
3254  //
3255  mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*xt;
3256  mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
3257  /* mom_u_adv[2] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*zt; */
3258  dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*xt;
3259  dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
3260  /* dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt; */
3261 
3262  mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*xt;
3263  mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
3264  /* mom_v_adv[2] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*zt; */
3265  dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*xt;
3266  dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
3267  /* dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt; */
3268 
3269  /* mom_w_adv[0] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*xt; */
3270  /* mom_w_adv[1] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*yt; */
3271  /* mom_w_adv[2] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*zt; */
3272  /* dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*xt; */
3273  /* dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt; */
3274  /* dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt; */
3275  //
3276  //calculate time derivatives
3277  //
3278  ck.bdf(alphaBDF,
3279  q_mom_u_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
3280  mom_u_acc,
3281  dmom_u_acc_u,
3282  mom_u_acc_t,
3283  dmom_u_acc_u_t);
3284  ck.bdf(alphaBDF,
3285  q_mom_v_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
3286  mom_v_acc,
3287  dmom_v_acc_v,
3288  mom_v_acc_t,
3289  dmom_v_acc_v_t);
3290  /* ck.bdf(alphaBDF, */
3291  /* q_mom_w_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV, */
3292  /* mom_w_acc, */
3293  /* dmom_w_acc_w, */
3294  /* mom_w_acc_t, */
3295  /* dmom_w_acc_w_t); */
3296  //
3297  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
3298  //
3299 
3300  mom_u_acc_t *= dmom_u_acc_u;
3301  mom_v_acc_t *= dmom_v_acc_v;
3302 
3303  //
3304  dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt);
3305  dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
3306  /* dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt); */
3307  //
3308  //calculate strong residual
3309  //
3310  pdeResidual_p =
3311  ck.Mass_strong(q_dvos_dt.data()[eN_k]) +
3312  ck.Advection_strong(dmass_adv_u,grad_u) +
3313  ck.Advection_strong(dmass_adv_v,grad_v) +
3314  /* ck.Advection_strong(dmass_adv_w,grad_w) + */
3315  DM2*MOVING_DOMAIN*ck.Reaction_strong(alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity) +
3316  //VRANS
3317  ck.Reaction_strong(mass_source);
3318  //
3319 
3320  pdeResidual_u =
3321  ck.Mass_strong(mom_u_acc_t) +
3322  ck.Advection_strong(dmom_adv_sge,grad_u) +
3323  ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
3324  ck.Reaction_strong(mom_u_source) -
3325  ck.Reaction_strong(u*div_mesh_velocity);
3326 
3327  pdeResidual_v =
3328  ck.Mass_strong(mom_v_acc_t) +
3329  ck.Advection_strong(dmom_adv_sge,grad_v) +
3330  ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
3331  ck.Reaction_strong(mom_v_source) -
3332  ck.Reaction_strong(v*div_mesh_velocity);
3333 
3334  /* pdeResidual_w =
3335  ck.Mass_strong(mom_w_acc_t) + */
3336  /* ck.Advection_strong(dmom_adv_sge,grad_w) + */
3337  /* ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) + */
3338  /* ck.Reaction_strong(mom_w_source) - */
3339  /* ck.Reaction_strong(w*div_mesh_velocity); */
3340 
3341  //calculate the Jacobian of strong residual
3342  for (int j=0;j<nDOF_trial_element;j++)
3343  {
3344  register int j_nSpace = j*nSpace;
3345  dpdeResidual_p_u[j]=ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial[j_nSpace]);
3346  dpdeResidual_p_v[j]=ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial[j_nSpace]);
3347  /* dpdeResidual_p_w[j]=ck.AdvectionJacobian_strong(dmass_adv_w,&vel_grad_trial[j_nSpace]); */
3348 
3349  dpdeResidual_u_p[j]=ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace]);
3350  dpdeResidual_u_u[j]=ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3351  ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
3352  ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref.data()[k*nDOF_trial_element+j]);
3353 
3354  dpdeResidual_v_p[j]=ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace]);
3355  dpdeResidual_v_v[j]=ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3356  ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
3357  ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref.data()[k*nDOF_trial_element+j]);
3358 
3359  /* dpdeResidual_w_p[j]=ck.HamiltonianJacobian_strong(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace]); */
3360  /* dpdeResidual_w_w[j]=ck.MassJacobian_strong(dmom_w_acc_w_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) + */
3361  /* ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
3362  ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref.data()[k*nDOF_trial_element+j]); */
3363 
3364  //VRANS account for drag terms, diagonal only here ... decide if need off diagonal terms too
3365  dpdeResidual_u_u[j]+= ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3366  dpdeResidual_v_v[j]+= ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3367  /* dpdeResidual_w_w[j]+= ck.ReactionJacobian_strong(dmom_w_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j]); */
3368  //
3369  }
3370  //calculate tau and tau*Res
3371  //cek debug
3372  double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
3373  calculateSubgridError_tau(hFactor,
3374  elementDiameter.data()[eN],
3375  tmpR,//dmom_u_acc_u_t,
3376  dmom_u_acc_u,
3377  dmom_adv_sge,
3378  mom_uu_diff_ten[1],
3379  dmom_u_ham_grad_p[0],
3380  tau_v0,
3381  tau_p0,
3382  q_cfl.data()[eN_k]);
3383 
3384  calculateSubgridError_tau(Ct_sge,Cd_sge,
3385  G,G_dd_G,tr_G,
3386  tmpR,//dmom_u_acc_u_t,
3387  dmom_adv_sge,
3388  mom_uu_diff_ten[1],
3389  dmom_u_ham_grad_p[0],
3390  tau_v1,
3391  tau_p1,
3392  q_cfl.data()[eN_k]);
3393 
3394 
3395  tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
3396  tau_p = PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
3397  calculateSubgridError_tauRes(tau_p,
3398  tau_v,
3399  pdeResidual_p,
3400  pdeResidual_u,
3401  pdeResidual_v,
3402  pdeResidual_w,
3403  subgridError_p,
3404  subgridError_u,
3405  subgridError_v,
3406  subgridError_w);
3407 
3408  calculateSubgridErrorDerivatives_tauRes(tau_p,
3409  tau_v,
3410  dpdeResidual_p_u,
3411  dpdeResidual_p_v,
3412  dpdeResidual_p_w,
3413  dpdeResidual_u_p,
3414  dpdeResidual_u_u,
3415  dpdeResidual_v_p,
3416  dpdeResidual_v_v,
3417  dpdeResidual_w_p,
3418  dpdeResidual_w_w,
3419  dsubgridError_p_u,
3420  dsubgridError_p_v,
3421  dsubgridError_p_w,
3422  dsubgridError_u_p,
3423  dsubgridError_u_u,
3424  dsubgridError_v_p,
3425  dsubgridError_v_v,
3426  dsubgridError_w_p,
3427  dsubgridError_w_w);
3428  // velocity used in adjoint (VMS or RBLES, with or without lagging the grid scale velocity)
3429  dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt + useRBLES*subgridError_u);
3430  dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
3431  /* dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w); */
3432 
3433  //calculate the adjoint times the test functions
3434  for (int i=0;i<nDOF_test_element;i++)
3435  {
3436  register int i_nSpace = i*nSpace;
3437  Lstar_u_u[i]=ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3438  Lstar_v_v[i]=ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3439  //VRANS account for drag terms, diagonal only here ... decide if need off diagonal terms too
3440  Lstar_u_u[i]+=ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
3441  Lstar_v_v[i]+=ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
3442  }
3443 
3444  // Assumes non-lagged subgrid velocity
3445  dmom_u_adv_u[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3446  dmom_u_adv_u[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3447  /* dmom_u_adv_u[2] += dmom_u_acc_u*(useRBLES*subgridError_w); */
3448 
3449  dmom_v_adv_v[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3450  dmom_v_adv_v[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3451  /* dmom_v_adv_v[2] += dmom_u_acc_u*(useRBLES*subgridError_w); */
3452 
3453  /* dmom_w_adv_w[0] += dmom_u_acc_u*(useRBLES*subgridError_u); */
3454  /* dmom_w_adv_w[1] += dmom_u_acc_u*(useRBLES*subgridError_v); */
3455  /* dmom_w_adv_w[2] += dmom_u_acc_u*(useRBLES*subgridError_w); */
3456 
3457 
3458  //cek todo add RBLES terms consistent to residual modifications or ignore the partials w.r.t the additional RBLES terms
3459  for(int i=0;i<nDOF_test_element;i++)
3460  {
3461  register int i_nSpace = i*nSpace;
3462  for(int j=0;j<nDOF_trial_element;j++)
3463  {
3464  register int j_nSpace = j*nSpace;
3465  /* elementJacobian_p_p[i][j] += ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_p[i]) + */
3466  /* ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_p[i]);// + */
3467  /* /\* ck.SubgridErrorJacobian(dsubgridError_w_p[j],Lstar_w_p[i]); *\/ */
3468 
3469  /* elementJacobian_p_u[i][j] += ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + */
3470  /* ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_p[i]); */
3471  /* elementJacobian_p_v[i][j] += ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + */
3472  /* ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_p[i]); */
3473  /* elementJacobian_p_w[i][j] += ck.AdvectionJacobian_weak(dmass_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + */
3474  /* ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_p[i]); */
3475 
3476  /* elementJacobian_u_p[i][j] += ck.HamiltonianJacobian_weak(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) + */
3477  /* ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_u[i]); */
3478  elementJacobian_u_u[i][j] +=
3479  ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3480  ck.HamiltonianJacobian_weak(dmom_u_ham_grad_u,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
3481  ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3482  ck.SimpleDiffusionJacobian_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
3483  //VRANS
3484  ck.ReactionJacobian_weak(dmom_u_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3485  //
3486  //ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_u[i]) +
3487  ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
3488  ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3489  elementJacobian_u_v[i][j] +=
3490  ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3491  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]) +
3492  //VRANS
3493  ck.ReactionJacobian_weak(dmom_u_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]);// +
3494  //
3495  //ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_u[i]);
3496  /* elementJacobian_u_w[i][j] += ck.AdvectionJacobian_weak(dmom_u_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) + */
3497  /* ck.SimpleDiffusionJacobian_weak(sdInfo_u_w_rowptr.data(),sdInfo_u_w_colind.data(),mom_uw_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) + */
3498  /* //VRANS */
3499  /* ck.ReactionJacobian_weak(dmom_u_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) + */
3500  /* // */
3501  /* ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_u[i]); */
3502 
3503  /* elementJacobian_v_p[i][j] += ck.HamiltonianJacobian_weak(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) + */
3504  /* ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_v[i]); */
3505  elementJacobian_v_u[i][j] +=
3506  ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3507  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]) +
3508  //VRANS
3509  ck.ReactionJacobian_weak(dmom_v_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]);// +
3510  //
3511  //ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_v[i]);
3512  elementJacobian_v_v[i][j] +=
3513  ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3514  ck.HamiltonianJacobian_weak(dmom_v_ham_grad_v,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
3515  ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3516  ck.SimpleDiffusionJacobian_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
3517  //VRANS
3518  ck.ReactionJacobian_weak(dmom_v_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3519  //
3520  //ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_v[i]) +
3521  ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
3522  ck.NumericalDiffusionJacobian(q_numDiff_v_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3523  /* elementJacobian_v_w[i][j] += ck.AdvectionJacobian_weak(dmom_v_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) + */
3524  /* ck.SimpleDiffusionJacobian_weak(sdInfo_v_w_rowptr.data(),sdInfo_v_w_colind.data(),mom_vw_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) + */
3525  /* //VRANS */
3526  /* ck.ReactionJacobian_weak(dmom_v_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) + */
3527  /* // */
3528  /* ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_v[i]); */
3529 
3530  /* elementJacobian_w_p[i][j] += ck.HamiltonianJacobian_weak(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) + */
3531  /* ck.SubgridErrorJacobian(dsubgridError_w_p[j],Lstar_w_w[i]); */
3532  /* elementJacobian_w_u[i][j] += ck.AdvectionJacobian_weak(dmom_w_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) + */
3533  /* ck.SimpleDiffusionJacobian_weak(sdInfo_w_u_rowptr.data(),sdInfo_w_u_colind.data(),mom_wu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) + */
3534  /* //VRANS */
3535  /* ck.ReactionJacobian_weak(dmom_w_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) + */
3536  /* // */
3537  /* ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_w[i]); */
3538  /* elementJacobian_w_v[i][j] += ck.AdvectionJacobian_weak(dmom_w_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) + */
3539  /* ck.SimpleDiffusionJacobian_weak(sdInfo_w_v_rowptr.data(),sdInfo_w_v_colind.data(),mom_wv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) + */
3540  /* //VRANS */
3541  /* ck.ReactionJacobian_weak(dmom_w_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) + */
3542  /* // */
3543  /* ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_w[i]); */
3544  /* elementJacobian_w_w[i][j] += ck.MassJacobian_weak(dmom_w_acc_w_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) + */
3545  /* ck.HamiltonianJacobian_weak(dmom_w_ham_grad_w,&vel_grad_trial[j_nSpace],vel_test_dV[i]) + */
3546  /* ck.AdvectionJacobian_weak(dmom_w_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) + */
3547  /* ck.SimpleDiffusionJacobian_weak(sdInfo_w_w_rowptr.data(),sdInfo_w_w_colind.data(),mom_ww_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) + */
3548  /* //VRANS */
3549  /* ck.ReactionJacobian_weak(dmom_w_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) + */
3550  /* // */
3551  /* ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_w[i]) + */
3552  /* ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_w[i]) + */
3553  /* ck.NumericalDiffusionJacobian(q_numDiff_w_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]); */
3554  }//j
3555  }//i
3556  }//k
3557  //
3558  //load into element Jacobian into global Jacobian
3559  //
3560  for (int i=0;i<nDOF_test_element;i++)
3561  {
3562  register int eN_i = eN*nDOF_test_element+i;
3563  for (int j=0;j<nDOF_trial_element;j++)
3564  {
3565  register int eN_i_j = eN_i*nDOF_trial_element+j;
3566  /* globalJacobian.data()[csrRowIndeces_p_p[eN_i] + csrColumnOffsets_p_p[eN_i_j]] += elementJacobian_p_p[i][j]; */
3567  /* globalJacobian.data()[csrRowIndeces_p_u[eN_i] + csrColumnOffsets_p_u[eN_i_j]] += elementJacobian_p_u[i][j]; */
3568  /* globalJacobian.data()[csrRowIndeces_p_v[eN_i] + csrColumnOffsets_p_v[eN_i_j]] += elementJacobian_p_v[i][j]; */
3569  /* globalJacobian.data()[csrRowIndeces_p_w[eN_i] + csrColumnOffsets_p_w[eN_i_j]] += elementJacobian_p_w[i][j]; */
3570 
3571  /* globalJacobian.data()[csrRowIndeces_u_p[eN_i] + csrColumnOffsets_u_p[eN_i_j]] += elementJacobian_u_p[i][j]; */
3572  globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
3573  globalJacobian.data()[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_u_v[eN_i_j]] += elementJacobian_u_v[i][j];
3574  /* globalJacobian.data()[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_u_w[eN_i_j]] += elementJacobian_u_w[i][j]; */
3575 
3576  /* globalJacobian.data()[csrRowIndeces_v_p[eN_i] + csrColumnOffsets_v_p[eN_i_j]] += elementJacobian_v_p[i][j]; */
3577  globalJacobian.data()[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_v_u[eN_i_j]] += elementJacobian_v_u[i][j];
3578  globalJacobian.data()[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_v_v[eN_i_j]] += elementJacobian_v_v[i][j];
3579  /* globalJacobian.data()[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_v_w[eN_i_j]] += elementJacobian_v_w[i][j]; */
3580 
3581  /* globalJacobian.data()[csrRowIndeces_w_p[eN_i] + csrColumnOffsets_w_p[eN_i_j]] += elementJacobian_w_p[i][j]; */
3582  /* globalJacobian.data()[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_w_u[eN_i_j]] += elementJacobian_w_u[i][j]; */
3583  /* globalJacobian.data()[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_w_v[eN_i_j]] += elementJacobian_w_v[i][j]; */
3584  /* globalJacobian.data()[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_w_w[eN_i_j]] += elementJacobian_w_w[i][j]; */
3585  }//j
3586  }//i
3587  }//elements
3588  //
3589  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
3590  //
3591  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3592  {
3593  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
3594  eN = elementBoundaryElementsArray.data()[ebN*2+0],
3595  eN_nDOF_trial_element = eN*nDOF_trial_element,
3596  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
3597  register double eps_rho,eps_mu;
3598  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3599  {
3600  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3601  ebNE_kb_nSpace = ebNE_kb*nSpace,
3602  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3603  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3604 
3605  register double p_ext=0.0,
3606  u_ext=0.0,
3607  v_ext=0.0,
3608  w_ext=0.0,
3609  grad_p_ext[nSpace],
3610  grad_u_ext[nSpace],
3611  grad_v_ext[nSpace],
3612  grad_w_ext[nSpace],
3613  mom_u_acc_ext=0.0,
3614  dmom_u_acc_u_ext=0.0,
3615  mom_v_acc_ext=0.0,
3616  dmom_v_acc_v_ext=0.0,
3617  mom_w_acc_ext=0.0,
3618  dmom_w_acc_w_ext=0.0,
3619  mass_adv_ext[nSpace],
3620  dmass_adv_u_ext[nSpace],
3621  dmass_adv_v_ext[nSpace],
3622  dmass_adv_w_ext[nSpace],
3623  mom_u_adv_ext[nSpace],
3624  dmom_u_adv_u_ext[nSpace],
3625  dmom_u_adv_v_ext[nSpace],
3626  dmom_u_adv_w_ext[nSpace],
3627  mom_v_adv_ext[nSpace],
3628  dmom_v_adv_u_ext[nSpace],
3629  dmom_v_adv_v_ext[nSpace],
3630  dmom_v_adv_w_ext[nSpace],
3631  mom_w_adv_ext[nSpace],
3632  dmom_w_adv_u_ext[nSpace],
3633  dmom_w_adv_v_ext[nSpace],
3634  dmom_w_adv_w_ext[nSpace],
3635  mom_uu_diff_ten_ext[nSpace],
3636  mom_vv_diff_ten_ext[nSpace],
3637  mom_ww_diff_ten_ext[nSpace],
3638  mom_uv_diff_ten_ext[1],
3639  mom_uw_diff_ten_ext[1],
3640  mom_vu_diff_ten_ext[1],
3641  mom_vw_diff_ten_ext[1],
3642  mom_wu_diff_ten_ext[1],
3643  mom_wv_diff_ten_ext[1],
3644  mom_u_source_ext=0.0,
3645  mom_v_source_ext=0.0,
3646  mom_w_source_ext=0.0,
3647  mom_u_ham_ext=0.0,
3648  dmom_u_ham_grad_p_ext[nSpace],
3649  dmom_u_ham_grad_u_ext[nSpace],
3650  mom_v_ham_ext=0.0,
3651  dmom_v_ham_grad_p_ext[nSpace],
3652  dmom_v_ham_grad_v_ext[nSpace],
3653  mom_w_ham_ext=0.0,
3654  dmom_w_ham_grad_p_ext[nSpace],
3655  dmom_w_ham_grad_w_ext[nSpace],
3656  dmom_u_adv_p_ext[nSpace],
3657  dmom_v_adv_p_ext[nSpace],
3658  dmom_w_adv_p_ext[nSpace],
3659  dflux_mass_u_ext=0.0,
3660  dflux_mass_v_ext=0.0,
3661  dflux_mass_w_ext=0.0,
3662  dflux_mom_u_adv_p_ext=0.0,
3663  dflux_mom_u_adv_u_ext=0.0,
3664  dflux_mom_u_adv_v_ext=0.0,
3665  dflux_mom_u_adv_w_ext=0.0,
3666  dflux_mom_v_adv_p_ext=0.0,
3667  dflux_mom_v_adv_u_ext=0.0,
3668  dflux_mom_v_adv_v_ext=0.0,
3669  dflux_mom_v_adv_w_ext=0.0,
3670  dflux_mom_w_adv_p_ext=0.0,
3671  dflux_mom_w_adv_u_ext=0.0,
3672  dflux_mom_w_adv_v_ext=0.0,
3673  dflux_mom_w_adv_w_ext=0.0,
3674  bc_p_ext=0.0,
3675  bc_u_ext=0.0,
3676  bc_v_ext=0.0,
3677  bc_w_ext=0.0,
3678  bc_mom_u_acc_ext=0.0,
3679  bc_dmom_u_acc_u_ext=0.0,
3680  bc_mom_v_acc_ext=0.0,
3681  bc_dmom_v_acc_v_ext=0.0,
3682  bc_mom_w_acc_ext=0.0,
3683  bc_dmom_w_acc_w_ext=0.0,
3684  bc_mass_adv_ext[nSpace],
3685  bc_dmass_adv_u_ext[nSpace],
3686  bc_dmass_adv_v_ext[nSpace],
3687  bc_dmass_adv_w_ext[nSpace],
3688  bc_mom_u_adv_ext[nSpace],
3689  bc_dmom_u_adv_u_ext[nSpace],
3690  bc_dmom_u_adv_v_ext[nSpace],
3691  bc_dmom_u_adv_w_ext[nSpace],
3692  bc_mom_v_adv_ext[nSpace],
3693  bc_dmom_v_adv_u_ext[nSpace],
3694  bc_dmom_v_adv_v_ext[nSpace],
3695  bc_dmom_v_adv_w_ext[nSpace],
3696  bc_mom_w_adv_ext[nSpace],
3697  bc_dmom_w_adv_u_ext[nSpace],
3698  bc_dmom_w_adv_v_ext[nSpace],
3699  bc_dmom_w_adv_w_ext[nSpace],
3700  bc_mom_uu_diff_ten_ext[nSpace],
3701  bc_mom_vv_diff_ten_ext[nSpace],
3702  bc_mom_ww_diff_ten_ext[nSpace],
3703  bc_mom_uv_diff_ten_ext[1],
3704  bc_mom_uw_diff_ten_ext[1],
3705  bc_mom_vu_diff_ten_ext[1],
3706  bc_mom_vw_diff_ten_ext[1],
3707  bc_mom_wu_diff_ten_ext[1],
3708  bc_mom_wv_diff_ten_ext[1],
3709  bc_mom_u_source_ext=0.0,
3710  bc_mom_v_source_ext=0.0,
3711  bc_mom_w_source_ext=0.0,
3712  bc_mom_u_ham_ext=0.0,
3713  bc_dmom_u_ham_grad_p_ext[nSpace],
3714  bc_dmom_u_ham_grad_u_ext[nSpace],
3715  bc_mom_v_ham_ext=0.0,
3716  bc_dmom_v_ham_grad_p_ext[nSpace],
3717  bc_dmom_v_ham_grad_v_ext[nSpace],
3718  bc_mom_w_ham_ext=0.0,
3719  bc_dmom_w_ham_grad_p_ext[nSpace],
3720  bc_dmom_w_ham_grad_w_ext[nSpace],
3721  fluxJacobian_u_u[nDOF_trial_element],
3722  fluxJacobian_u_v[nDOF_trial_element],
3723  fluxJacobian_v_u[nDOF_trial_element],
3724  fluxJacobian_v_v[nDOF_trial_element],
3725  jac_ext[nSpace*nSpace],
3726  jacDet_ext,
3727  jacInv_ext[nSpace*nSpace],
3728  boundaryJac[nSpace*(nSpace-1)],
3729  metricTensor[(nSpace-1)*(nSpace-1)],
3730  metricTensorDetSqrt,
3731  vel_grad_trial_trace[nDOF_trial_element*nSpace],
3732  dS,
3733  vel_test_dS[nDOF_test_element],
3734  normal[2],
3735  x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3736  vel_grad_test_dS[nDOF_trial_element*nSpace],
3737  //VRANS
3738  vos_ext,
3739  //
3740  G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
3741  ck.calculateMapping_elementBoundary(eN,
3742  ebN_local,
3743  kb,
3744  ebN_local_kb,
3745  mesh_dof.data(),
3746  mesh_l2g.data(),
3747  mesh_trial_trace_ref.data(),
3748  mesh_grad_trial_trace_ref.data(),
3749  boundaryJac_ref.data(),
3750  jac_ext,
3751  jacDet_ext,
3752  jacInv_ext,
3753  boundaryJac,
3754  metricTensor,
3755  metricTensorDetSqrt,
3756  normal_ref.data(),
3757  normal,
3758  x_ext,y_ext,z_ext);
3759  ck.calculateMappingVelocity_elementBoundary(eN,
3760  ebN_local,
3761  kb,
3762  ebN_local_kb,
3763  mesh_velocity_dof.data(),
3764  mesh_l2g.data(),
3765  mesh_trial_trace_ref.data(),
3766  xt_ext,yt_ext,zt_ext,
3767  normal,
3768  boundaryJac,
3769  metricTensor,
3770  integralScaling);
3771  //dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
3772  dS = metricTensorDetSqrt*dS_ref.data()[kb];
3773  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3774  ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
3775 
3776  eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3777  eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3778 
3779  //compute shape and solution information
3780  //shape
3781  /* ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace); */
3782  ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3783  //solution and gradients
3784  /* ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_ext); */
3785  p_ext = ebqe_p.data()[ebNE_kb];
3786  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);
3787  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);
3788  /* ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext); */
3789  /* ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext); */
3790  for (int I=0;I<nSpace;I++)
3791  grad_p_ext[I] = ebqe_grad_p.data()[ebNE_kb_nSpace+I];
3792  ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3793  ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3794  /* ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext); */
3795  //precalculate test function products with integration weights
3796  for (int j=0;j<nDOF_trial_element;j++)
3797  {
3798  /* p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
3799  vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
3800  for (int I=0;I<nSpace;I++)
3801  vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
3802  }
3803  //
3804  //load the boundary values
3805  //
3806  bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
3807  //bc values at moving boundaries are specified relative to boundary motion so we need to add it here
3808  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*(ebqe_bc_u_ext.data()[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
3809  bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*(ebqe_bc_v_ext.data()[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
3810  /* bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*(ebqe_bc_w_ext.data()[ebNE_kb] + MOVING_DOMAIN*zt_ext) + (1-isDOFBoundary_w.data()[ebNE_kb])*w_ext; */
3811  //VRANS
3812  vos_ext = ebqe_vos_ext.data()[ebNE_kb];//sed fraction - gco check
3813 
3814  //
3815  //calculate the internal and external trace of the pde coefficients
3816  //
3817  double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);//not interested in saving boundary eddy viscosity for now
3818  evaluateCoefficients(eps_rho,
3819  eps_mu,
3820  sigma,
3821  rho_0,
3822  nu_0,
3823  rho_1,
3824  nu_1,
3825  rho_s,
3826  elementDiameter.data()[eN],
3827  smagorinskyConstant,
3828  turbulenceClosureModel,
3829  g.data(),
3830  useVF,
3831  ebqe_vf_ext.data()[ebNE_kb],
3832  ebqe_phi_ext.data()[ebNE_kb],
3833  &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
3834  ebqe_kappa_phi_ext.data()[ebNE_kb],
3835  //VRANS
3836  vos_ext,
3837  //
3838  p_ext,
3839  grad_p_ext,
3840  grad_u_ext,
3841  grad_v_ext,
3842  grad_w_ext,
3843  u_ext,
3844  v_ext,
3845  w_ext,
3846  ebqe_velocity_star.data()[ebNE_kb_nSpace+0],
3847  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],
3848  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],//hack,not used
3849  eddy_viscosity_ext,
3850  mom_u_acc_ext,
3851  dmom_u_acc_u_ext,
3852  mom_v_acc_ext,
3853  dmom_v_acc_v_ext,
3854  mom_w_acc_ext,
3855  dmom_w_acc_w_ext,
3856  mass_adv_ext,
3857  dmass_adv_u_ext,
3858  dmass_adv_v_ext,
3859  dmass_adv_w_ext,
3860  mom_u_adv_ext,
3861  dmom_u_adv_u_ext,
3862  dmom_u_adv_v_ext,
3863  dmom_u_adv_w_ext,
3864  mom_v_adv_ext,
3865  dmom_v_adv_u_ext,
3866  dmom_v_adv_v_ext,
3867  dmom_v_adv_w_ext,
3868  mom_w_adv_ext,
3869  dmom_w_adv_u_ext,
3870  dmom_w_adv_v_ext,
3871  dmom_w_adv_w_ext,
3872  mom_uu_diff_ten_ext,
3873  mom_vv_diff_ten_ext,
3874  mom_ww_diff_ten_ext,
3875  mom_uv_diff_ten_ext,
3876  mom_uw_diff_ten_ext,
3877  mom_vu_diff_ten_ext,
3878  mom_vw_diff_ten_ext,
3879  mom_wu_diff_ten_ext,
3880  mom_wv_diff_ten_ext,
3881  mom_u_source_ext,
3882  mom_v_source_ext,
3883  mom_w_source_ext,
3884  mom_u_ham_ext,
3885  dmom_u_ham_grad_p_ext,
3886  dmom_u_ham_grad_u_ext,
3887  mom_v_ham_ext,
3888  dmom_v_ham_grad_p_ext,
3889  dmom_v_ham_grad_v_ext,
3890  mom_w_ham_ext,
3891  dmom_w_ham_grad_p_ext,
3892  dmom_w_ham_grad_w_ext);
3893  evaluateCoefficients(eps_rho,
3894  eps_mu,
3895  sigma,
3896  rho_0,
3897  nu_0,
3898  rho_1,
3899  nu_1,
3900  rho_s,
3901  elementDiameter.data()[eN],
3902  smagorinskyConstant,
3903  turbulenceClosureModel,
3904  g.data(),
3905  useVF,
3906  bc_ebqe_vf_ext.data()[ebNE_kb],
3907  bc_ebqe_phi_ext.data()[ebNE_kb],
3908  &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
3909  ebqe_kappa_phi_ext.data()[ebNE_kb],
3910  //VRANS
3911  vos_ext,
3912  //
3913  bc_p_ext,
3914  grad_p_ext,
3915  grad_u_ext,
3916  grad_v_ext,
3917  grad_w_ext,
3918  bc_u_ext,
3919  bc_v_ext,
3920  bc_w_ext,
3921  ebqe_velocity_star.data()[ebNE_kb_nSpace+0],
3922  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],
3923  ebqe_velocity_star.data()[ebNE_kb_nSpace+1],//hack,not used
3924  bc_eddy_viscosity_ext,
3925  bc_mom_u_acc_ext,
3926  bc_dmom_u_acc_u_ext,
3927  bc_mom_v_acc_ext,
3928  bc_dmom_v_acc_v_ext,
3929  bc_mom_w_acc_ext,
3930  bc_dmom_w_acc_w_ext,
3931  bc_mass_adv_ext,
3932  bc_dmass_adv_u_ext,
3933  bc_dmass_adv_v_ext,
3934  bc_dmass_adv_w_ext,
3935  bc_mom_u_adv_ext,
3936  bc_dmom_u_adv_u_ext,
3937  bc_dmom_u_adv_v_ext,
3938  bc_dmom_u_adv_w_ext,
3939  bc_mom_v_adv_ext,
3940  bc_dmom_v_adv_u_ext,
3941  bc_dmom_v_adv_v_ext,
3942  bc_dmom_v_adv_w_ext,
3943  bc_mom_w_adv_ext,
3944  bc_dmom_w_adv_u_ext,
3945  bc_dmom_w_adv_v_ext,
3946  bc_dmom_w_adv_w_ext,
3947  bc_mom_uu_diff_ten_ext,
3948  bc_mom_vv_diff_ten_ext,
3949  bc_mom_ww_diff_ten_ext,
3950  bc_mom_uv_diff_ten_ext,
3951  bc_mom_uw_diff_ten_ext,
3952  bc_mom_vu_diff_ten_ext,
3953  bc_mom_vw_diff_ten_ext,
3954  bc_mom_wu_diff_ten_ext,
3955  bc_mom_wv_diff_ten_ext,
3956  bc_mom_u_source_ext,
3957  bc_mom_v_source_ext,
3958  bc_mom_w_source_ext,
3959  bc_mom_u_ham_ext,
3960  bc_dmom_u_ham_grad_p_ext,
3961  bc_dmom_u_ham_grad_u_ext,
3962  bc_mom_v_ham_ext,
3963  bc_dmom_v_ham_grad_p_ext,
3964  bc_dmom_v_ham_grad_v_ext,
3965  bc_mom_w_ham_ext,
3966  bc_dmom_w_ham_grad_p_ext,
3967  bc_dmom_w_ham_grad_w_ext);
3968 
3969  //
3970  //moving domain
3971  //
3972  mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
3973  mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
3974  /* mom_u_adv_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*zt_ext; */
3975  dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
3976  dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
3977  /* dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext; */
3978 
3979  mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
3980  mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
3981  /* mom_v_adv_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*zt_ext; */
3982  dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
3983  dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
3984  /* dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext; */
3985 
3986  /* mom_w_adv_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*xt_ext; */
3987  /* mom_w_adv_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*yt_ext; */
3988  /* mom_w_adv_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*zt_ext; */
3989  /* dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext; */
3990  /* dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext; */
3991  /* dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext; */
3992 
3993  //moving domain bc's
3994  // mql. CHECK.
3995  bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
3996  bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
3997  /* bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*zt_ext; */
3998 
3999  bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
4000  bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
4001  /* bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*zt_ext; */
4002 
4003  /* bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*bc_mom_w_acc_ext*xt_ext; */
4004  /* bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*bc_mom_w_acc_ext*yt_ext; */
4005  /* bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*bc_mom_w_acc_ext*zt_ext; */
4006  //
4007  //calculate the numerical fluxes
4008  //
4009  exteriorNumericalAdvectiveFluxDerivatives(isDOFBoundary_p.data()[ebNE_kb],
4010  isDOFBoundary_u.data()[ebNE_kb],
4011  isDOFBoundary_v.data()[ebNE_kb],
4012  isDOFBoundary_w.data()[ebNE_kb],
4013  isAdvectiveFluxBoundary_p.data()[ebNE_kb],
4014  isAdvectiveFluxBoundary_u.data()[ebNE_kb],
4015  isAdvectiveFluxBoundary_v.data()[ebNE_kb],
4016  isAdvectiveFluxBoundary_w.data()[ebNE_kb],
4017  dmom_u_ham_grad_p_ext[0],//=1/rho
4018  normal,
4019  dmom_u_acc_u_ext,
4020  bc_p_ext,
4021  bc_u_ext,
4022  bc_v_ext,
4023  bc_w_ext,
4024  bc_mass_adv_ext,
4025  bc_mom_u_adv_ext,
4026  bc_mom_v_adv_ext,
4027  bc_mom_w_adv_ext,
4028  ebqe_bc_flux_mass_ext.data()[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),//bc is relative mass flux
4029  ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
4030  ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
4031  ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
4032  p_ext,
4033  u_ext,
4034  v_ext,
4035  w_ext,
4036  mass_adv_ext,
4037  mom_u_adv_ext,
4038  mom_v_adv_ext,
4039  mom_w_adv_ext,
4040  dmass_adv_u_ext,
4041  dmass_adv_v_ext,
4042  dmass_adv_w_ext,
4043  dmom_u_adv_p_ext,
4044  dmom_u_adv_u_ext,
4045  dmom_u_adv_v_ext,
4046  dmom_u_adv_w_ext,
4047  dmom_v_adv_p_ext,
4048  dmom_v_adv_u_ext,
4049  dmom_v_adv_v_ext,
4050  dmom_v_adv_w_ext,
4051  dmom_w_adv_p_ext,
4052  dmom_w_adv_u_ext,
4053  dmom_w_adv_v_ext,
4054  dmom_w_adv_w_ext,
4055  dflux_mass_u_ext,
4056  dflux_mass_v_ext,
4057  dflux_mass_w_ext,
4058  dflux_mom_u_adv_p_ext,
4059  dflux_mom_u_adv_u_ext,
4060  dflux_mom_u_adv_v_ext,
4061  dflux_mom_u_adv_w_ext,
4062  dflux_mom_v_adv_p_ext,
4063  dflux_mom_v_adv_u_ext,
4064  dflux_mom_v_adv_v_ext,
4065  dflux_mom_v_adv_w_ext,
4066  dflux_mom_w_adv_p_ext,
4067  dflux_mom_w_adv_u_ext,
4068  dflux_mom_w_adv_v_ext,
4069  dflux_mom_w_adv_w_ext,
4070  &ebqe_velocity_star.data()[ebNE_kb_nSpace]);
4071  //
4072  //calculate the flux jacobian
4073  //
4074  ck.calculateGScale(G,normal,h_penalty);
4075  penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
4076  for (int j=0;j<nDOF_trial_element;j++)
4077  {
4078  register int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
4079  /* fluxJacobian_p_p[j]=0.0; */
4080  /* fluxJacobian_p_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
4081  /* fluxJacobian_p_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
4082  /* fluxJacobian_p_w[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]); */
4083 
4084  /* fluxJacobian_u_p[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]); */
4085  fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4086  ExteriorNumericalDiffusiveFluxJacobian(eps_rho,
4087  ebqe_phi_ext.data()[ebNE_kb],
4088  sdInfo_u_u_rowptr.data(),
4089  sdInfo_u_u_colind.data(),
4090  isDOFBoundary_u.data()[ebNE_kb],
4091  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4092  normal,
4093  mom_uu_diff_ten_ext,
4094  vel_trial_trace_ref.data()[ebN_local_kb_j],
4095  &vel_grad_trial_trace[j_nSpace],
4096  penalty);//ebqe_penalty_ext.data()[ebNE_kb]);
4097  fluxJacobian_u_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4098  ExteriorNumericalDiffusiveFluxJacobian(eps_rho,
4099  ebqe_phi_ext.data()[ebNE_kb],
4100  sdInfo_u_v_rowptr.data(),
4101  sdInfo_u_v_colind.data(),
4102  isDOFBoundary_v.data()[ebNE_kb],
4103  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4104  normal,
4105  mom_uv_diff_ten_ext,
4106  vel_trial_trace_ref.data()[ebN_local_kb_j],
4107  &vel_grad_trial_trace[j_nSpace],
4108  penalty);//ebqe_penalty_ext.data()[ebNE_kb]);
4109  /* fluxJacobian_u_w[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
4110  /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
4111  /* ebqe_phi_ext.data()[ebNE_kb], */
4112  /* sdInfo_u_w_rowptr.data(), */
4113  /* sdInfo_u_w_colind.data(), */
4114  /* isDOFBoundary_w.data()[ebNE_kb], */
4115  /* isDiffusiveFluxBoundary_u.data()[ebNE_kb], */
4116  /* normal, */
4117  /* mom_uw_diff_ten_ext, */
4118  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4119  /* &vel_grad_trial_trace[j_nSpace], */
4120  /* penalty);//ebqe_penalty_ext.data()[ebNE_kb]); */
4121 
4122  /* fluxJacobian_v_p[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]); */
4123  fluxJacobian_v_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4124  ExteriorNumericalDiffusiveFluxJacobian(eps_rho,
4125  ebqe_phi_ext.data()[ebNE_kb],
4126  sdInfo_v_u_rowptr.data(),
4127  sdInfo_v_u_colind.data(),
4128  isDOFBoundary_u.data()[ebNE_kb],
4129  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4130  normal,
4131  mom_vu_diff_ten_ext,
4132  vel_trial_trace_ref.data()[ebN_local_kb_j],
4133  &vel_grad_trial_trace[j_nSpace],
4134  penalty);//ebqe_penalty_ext.data()[ebNE_kb]);
4135  fluxJacobian_v_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4136  ExteriorNumericalDiffusiveFluxJacobian(eps_rho,
4137  ebqe_phi_ext.data()[ebNE_kb],
4138  sdInfo_v_v_rowptr.data(),
4139  sdInfo_v_v_colind.data(),
4140  isDOFBoundary_v.data()[ebNE_kb],
4141  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4142  normal,
4143  mom_vv_diff_ten_ext,
4144  vel_trial_trace_ref.data()[ebN_local_kb_j],
4145  &vel_grad_trial_trace[j_nSpace],
4146  penalty);//ebqe_penalty_ext.data()[ebNE_kb]);
4147  /* fluxJacobian_v_w[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
4148  /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
4149  /* ebqe_phi_ext.data()[ebNE_kb], */
4150  /* sdInfo_v_w_rowptr.data(), */
4151  /* sdInfo_v_w_colind.data(), */
4152  /* isDOFBoundary_w.data()[ebNE_kb], */
4153  /* isDiffusiveFluxBoundary_v.data()[ebNE_kb], */
4154  /* normal, */
4155  /* mom_vw_diff_ten_ext, */
4156  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4157  /* &vel_grad_trial_trace[j_nSpace], */
4158  /* penalty);//ebqe_penalty_ext.data()[ebNE_kb]); */
4159 
4160  /* fluxJacobian_w_p[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]); */
4161  /* fluxJacobian_w_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
4162  /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
4163  /* ebqe_phi_ext.data()[ebNE_kb], */
4164  /* sdInfo_w_u_rowptr.data(), */
4165  /* sdInfo_w_u_colind.data(), */
4166  /* isDOFBoundary_u.data()[ebNE_kb], */
4167  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
4168  /* normal, */
4169  /* mom_wu_diff_ten_ext, */
4170  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4171  /* &vel_grad_trial_trace[j_nSpace], */
4172  /* penalty);//ebqe_penalty_ext.data()[ebNE_kb]); */
4173  /* fluxJacobian_w_v[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
4174  /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
4175  /* ebqe_phi_ext.data()[ebNE_kb], */
4176  /* sdInfo_w_v_rowptr.data(), */
4177  /* sdInfo_w_v_colind.data(), */
4178  /* isDOFBoundary_v.data()[ebNE_kb], */
4179  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
4180  /* normal, */
4181  /* mom_wv_diff_ten_ext, */
4182  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4183  /* &vel_grad_trial_trace[j_nSpace], */
4184  /* penalty);//ebqe_penalty_ext.data()[ebNE_kb]); */
4185  /* fluxJacobian_w_w[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) + */
4186  /* ExteriorNumericalDiffusiveFluxJacobian(eps_rho, */
4187  /* ebqe_phi_ext.data()[ebNE_kb], */
4188  /* sdInfo_w_w_rowptr.data(), */
4189  /* sdInfo_w_w_colind.data(), */
4190  /* isDOFBoundary_w.data()[ebNE_kb], */
4191  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
4192  /* normal, */
4193  /* mom_ww_diff_ten_ext, */
4194  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4195  /* &vel_grad_trial_trace[j_nSpace], */
4196  /* penalty);//ebqe_penalty_ext.data()[ebNE_kb]); */
4197  }//j
4198  //
4199  //update the global Jacobian from the flux Jacobian
4200  //
4201  for (int i=0;i<nDOF_test_element;i++)
4202  {
4203  register int eN_i = eN*nDOF_test_element+i;
4204  for (int j=0;j<nDOF_trial_element;j++)
4205  {
4206  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
4207 
4208  /* globalJacobian.data()[csrRowIndeces_p_p[eN_i] + csrColumnOffsets_eb_p_p[ebN_i_j]] += fluxJacobian_p_p[j]*p_test_dS[i]; */
4209  /* globalJacobian.data()[csrRowIndeces_p_u[eN_i] + csrColumnOffsets_eb_p_u[ebN_i_j]] += fluxJacobian_p_u[j]*p_test_dS[i]; */
4210  /* globalJacobian.data()[csrRowIndeces_p_v[eN_i] + csrColumnOffsets_eb_p_v[ebN_i_j]] += fluxJacobian_p_v[j]*p_test_dS[i]; */
4211  /* globalJacobian.data()[csrRowIndeces_p_w[eN_i] + csrColumnOffsets_eb_p_w[ebN_i_j]] += fluxJacobian_p_w[j]*p_test_dS[i]; */
4212 
4213  /* globalJacobian.data()[csrRowIndeces_u_p[eN_i] + csrColumnOffsets_eb_u_p[ebN_i_j]] += fluxJacobian_u_p[j]*vel_test_dS[i]; */
4214  globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
4215  fluxJacobian_u_u[j]*vel_test_dS[i]+
4216  ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4217  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4218  eb_adjoint_sigma,
4219  vel_trial_trace_ref.data()[ebN_local_kb_j],
4220  normal,
4221  sdInfo_u_u_rowptr.data(),
4222  sdInfo_u_u_colind.data(),
4223  mom_uu_diff_ten_ext,
4224  &vel_grad_test_dS[i*nSpace]);
4225  globalJacobian.data()[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] +=
4226  fluxJacobian_u_v[j]*vel_test_dS[i]+
4227  ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4228  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4229  eb_adjoint_sigma,
4230  vel_trial_trace_ref.data()[ebN_local_kb_j],
4231  normal,
4232  sdInfo_u_v_rowptr.data(),
4233  sdInfo_u_v_colind.data(),
4234  mom_uv_diff_ten_ext,
4235  &vel_grad_test_dS[i*nSpace]);
4236  /* globalJacobian.data()[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_eb_u_w[ebN_i_j]] += fluxJacobian_u_w[j]*vel_test_dS[i]+ */
4237  /* ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb], */
4238  /* isDiffusiveFluxBoundary_u.data()[ebNE_kb], */
4239  /* eb_adjoint_sigma, */
4240  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4241  /* normal, */
4242  /* sdInfo_u_w_rowptr.data(), */
4243  /* sdInfo_u_w_colind.data(), */
4244  /* mom_uw_diff_ten_ext, */
4245  /* &vel_grad_test_dS[i*nSpace]); */
4246 
4247  /* globalJacobian.data()[csrRowIndeces_v_p[eN_i] + csrColumnOffsets_eb_v_p[ebN_i_j]] += fluxJacobian_v_p[j]*vel_test_dS[i]; */
4248  globalJacobian.data()[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] +=
4249  fluxJacobian_v_u[j]*vel_test_dS[i]+
4250  ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4251  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4252  eb_adjoint_sigma,
4253  vel_trial_trace_ref.data()[ebN_local_kb_j],
4254  normal,
4255  sdInfo_v_u_rowptr.data(),
4256  sdInfo_v_u_colind.data(),
4257  mom_vu_diff_ten_ext,
4258  &vel_grad_test_dS[i*nSpace]);
4259  globalJacobian.data()[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
4260  fluxJacobian_v_v[j]*vel_test_dS[i]+
4261  ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4262  isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4263  eb_adjoint_sigma,
4264  vel_trial_trace_ref.data()[ebN_local_kb_j],
4265  normal,
4266  sdInfo_v_v_rowptr.data(),
4267  sdInfo_v_v_colind.data(),
4268  mom_vv_diff_ten_ext,
4269  &vel_grad_test_dS[i*nSpace]);
4270  /* globalJacobian.data()[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_eb_v_w[ebN_i_j]] += fluxJacobian_v_w[j]*vel_test_dS[i]+ */
4271  /* ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb], */
4272  /* isDiffusiveFluxBoundary_v.data()[ebNE_kb], */
4273  /* eb_adjoint_sigma, */
4274  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4275  /* normal, */
4276  /* sdInfo_v_w_rowptr.data(), */
4277  /* sdInfo_v_w_colind.data(), */
4278  /* mom_vw_diff_ten_ext, */
4279  /* &vel_grad_test_dS[i*nSpace]); */
4280  /* globalJacobian.data()[csrRowIndeces_w_p[eN_i] + csrColumnOffsets_eb_w_p[ebN_i_j]] += fluxJacobian_w_p[j]*vel_test_dS[i]; */
4281  /* globalJacobian.data()[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_eb_w_u[ebN_i_j]] += fluxJacobian_w_u[j]*vel_test_dS[i]+ */
4282  /* ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb], */
4283  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
4284  /* eb_adjoint_sigma, */
4285  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4286  /* normal, */
4287  /* sdInfo_w_u_rowptr.data(), */
4288  /* sdInfo_w_u_colind.data(), */
4289  /* mom_wu_diff_ten_ext, */
4290  /* &vel_grad_test_dS[i*nSpace]); */
4291  /* globalJacobian.data()[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_eb_w_v[ebN_i_j]] += fluxJacobian_w_v[j]*vel_test_dS[i]+ */
4292  /* ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb], */
4293  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
4294  /* eb_adjoint_sigma, */
4295  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4296  /* normal, */
4297  /* sdInfo_w_v_rowptr.data(), */
4298  /* sdInfo_w_v_colind.data(), */
4299  /* mom_wv_diff_ten_ext, */
4300  /* &vel_grad_test_dS[i*nSpace]); */
4301  /* globalJacobian.data()[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] += fluxJacobian_w_w[j]*vel_test_dS[i]+ */
4302  /* ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb], */
4303  /* isDiffusiveFluxBoundary_w.data()[ebNE_kb], */
4304  /* eb_adjoint_sigma, */
4305  /* vel_trial_trace_ref.data()[ebN_local_kb_j], */
4306  /* normal, */
4307  /* sdInfo_w_w_rowptr.data(), */
4308  /* sdInfo_w_w_colind.data(), */
4309  /* mom_ww_diff_ten_ext, */
4310  /* &vel_grad_test_dS[i*nSpace]); */
4311  }//j
4312  }//i
4313  }//kb
4314  }//ebNE
4315  }//computeJacobian
4316 
4317  void calculateVelocityAverage(arguments_dict& args)
4318  {
4319  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
4320  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
4321  int nInteriorElementBoundaries_global = args.scalar<int>("nInteriorElementBoundaries_global");
4322  xt::pyarray<int>& interiorElementBoundariesArray = args.array<int>("interiorElementBoundariesArray");
4323  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
4324  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
4325  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
4326  xt::pyarray<double>& mesh_velocity_dof = args.array<double>("mesh_velocity_dof");
4327  double MOVING_DOMAIN = args.scalar<double>("MOVING_DOMAIN");
4328  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
4329  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
4330  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
4331  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
4332  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
4333  xt::pyarray<int>& vel_l2g = args.array<int>("vel_l2g");
4334  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
4335  xt::pyarray<double>& v_dof = args.array<double>("v_dof");
4336  xt::pyarray<double>& w_dof = args.array<double>("w_dof");
4337  xt::pyarray<double>& vos_dof = args.array<double>("vos_dof");
4338  xt::pyarray<double>& vel_trial_trace_ref = args.array<double>("vel_trial_trace_ref");
4339  xt::pyarray<double>& ebqe_velocity = args.array<double>("ebqe_velocity");
4340  xt::pyarray<double>& velocityAverage = args.array<double>("velocityAverage");
4341  int permutations[nQuadraturePoints_elementBoundary];
4342  double xArray_left[nQuadraturePoints_elementBoundary*2],
4343  xArray_right[nQuadraturePoints_elementBoundary*2];
4344  for (int i=0;i<nQuadraturePoints_elementBoundary;i++)
4345  permutations[i]=i;//just to initialize
4346  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
4347  {
4348  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
4349  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4350  {
4351  register int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace,
4352  ebNE_kb_nSpace = ebNE*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace;
4353  velocityAverage.data()[ebN_kb_nSpace+0]=ebqe_velocity.data()[ebNE_kb_nSpace+0];
4354  velocityAverage.data()[ebN_kb_nSpace+1]=ebqe_velocity.data()[ebNE_kb_nSpace+1];
4355  }//ebNE
4356  }
4357  for (int ebNI = 0; ebNI < nInteriorElementBoundaries_global; ebNI++)
4358  {
4359  register int ebN = interiorElementBoundariesArray.data()[ebNI],
4360  left_eN_global = elementBoundaryElementsArray.data()[ebN*2+0],
4361  left_ebN_element = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
4362  right_eN_global = elementBoundaryElementsArray.data()[ebN*2+1],
4363  right_ebN_element = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+1],
4364  left_eN_nDOF_trial_element = left_eN_global*nDOF_trial_element,
4365  right_eN_nDOF_trial_element = right_eN_global*nDOF_trial_element;
4366  double jac[nSpace*nSpace],
4367  jacDet,
4368  jacInv[nSpace*nSpace],
4369  boundaryJac[nSpace*(nSpace-1)],
4370  metricTensor[(nSpace-1)*(nSpace-1)],
4371  metricTensorDetSqrt,
4372  normal[2],
4373  x,y,z,
4374  xt,yt,zt,integralScaling;
4375 
4376  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4377  {
4378  ck.calculateMapping_elementBoundary(left_eN_global,
4379  left_ebN_element,
4380  kb,
4381  left_ebN_element*nQuadraturePoints_elementBoundary+kb,
4382  mesh_dof.data(),
4383  mesh_l2g.data(),
4384  mesh_trial_trace_ref.data(),
4385  mesh_grad_trial_trace_ref.data(),
4386  boundaryJac_ref.data(),
4387  jac,
4388  jacDet,
4389  jacInv,
4390  boundaryJac,
4391  metricTensor,
4392  metricTensorDetSqrt,
4393  normal_ref.data(),
4394  normal,
4395  x,y,z);
4396  xArray_left[kb*2+0] = x;
4397  xArray_left[kb*2+1] = y;
4398  /* xArray_left[kb*3+2] = z; */
4399  ck.calculateMapping_elementBoundary(right_eN_global,
4400  right_ebN_element,
4401  kb,
4402  right_ebN_element*nQuadraturePoints_elementBoundary+kb,
4403  mesh_dof.data(),
4404  mesh_l2g.data(),
4405  mesh_trial_trace_ref.data(),
4406  mesh_grad_trial_trace_ref.data(),
4407  boundaryJac_ref.data(),
4408  jac,
4409  jacDet,
4410  jacInv,
4411  boundaryJac,
4412  metricTensor,
4413  metricTensorDetSqrt,
4414  normal_ref.data(),
4415  normal,
4416  x,y,z);
4417  ck.calculateMappingVelocity_elementBoundary(left_eN_global,
4418  left_ebN_element,
4419  kb,
4420  left_ebN_element*nQuadraturePoints_elementBoundary+kb,
4421  mesh_velocity_dof.data(),
4422  mesh_l2g.data(),
4423  mesh_trial_trace_ref.data(),
4424  xt,yt,zt,
4425  normal,
4426  boundaryJac,
4427  metricTensor,
4428  integralScaling);
4429  xArray_right[kb*2+0] = x;
4430  xArray_right[kb*2+1] = y;
4431  /* xArray_right[kb*3+2] = z; */
4432  }
4433  for (int kb_left=0;kb_left<nQuadraturePoints_elementBoundary;kb_left++)
4434  {
4435  double errorNormMin = 1.0;
4436  for (int kb_right=0;kb_right<nQuadraturePoints_elementBoundary;kb_right++)
4437  {
4438  double errorNorm=0.0;
4439  for (int I=0;I<nSpace;I++)
4440  {
4441  errorNorm += fabs(xArray_left[kb_left*2+I]
4442  -
4443  xArray_right[kb_right*2+I]);
4444  }
4445  if (errorNorm < errorNormMin)
4446  {
4447  permutations[kb_right] = kb_left;
4448  errorNormMin = errorNorm;
4449  }
4450  }
4451  }
4452  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4453  {
4454  register int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace;
4455  register double u_left=0.0,
4456  v_left=0.0,
4457  u_right=0.0,
4458  v_right=0.0,
4459  vos_left=0.0,
4460  vos_right=0.0;
4461  register int left_kb = kb,
4462  right_kb = permutations[kb],
4463  left_ebN_element_kb_nDOF_test_element=(left_ebN_element*nQuadraturePoints_elementBoundary+left_kb)*nDOF_test_element,
4464  right_ebN_element_kb_nDOF_test_element=(right_ebN_element*nQuadraturePoints_elementBoundary+right_kb)*nDOF_test_element;
4465  //
4466  //calculate the velocity solution at quadrature points on left and right
4467  //
4468  ck.valFromDOF(vos_dof.data(),&vel_l2g.data()[left_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[left_ebN_element_kb_nDOF_test_element],vos_left);
4469  ck.valFromDOF(u_dof.data(),&vel_l2g.data()[left_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[left_ebN_element_kb_nDOF_test_element],u_left);
4470  ck.valFromDOF(v_dof.data(),&vel_l2g.data()[left_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[left_ebN_element_kb_nDOF_test_element],v_left);
4471  /* ck.valFromDOF(w_dof.data(),&vel_l2g.data()[left_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[left_ebN_element_kb_nDOF_test_element],w_left); */
4472  //
4473  ck.valFromDOF(vos_dof.data(),&vel_l2g.data()[right_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[right_ebN_element_kb_nDOF_test_element],vos_right);
4474  ck.valFromDOF(u_dof.data(),&vel_l2g.data()[right_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[right_ebN_element_kb_nDOF_test_element],u_right);
4475  ck.valFromDOF(v_dof.data(),&vel_l2g.data()[right_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[right_ebN_element_kb_nDOF_test_element],v_right);
4476  /* ck.valFromDOF(w_dof.data(),&vel_l2g.data()[right_eN_nDOF_trial_element],&vel_trial_trace_ref.data()[right_ebN_element_kb_nDOF_test_element],w_right); */
4477  //
4478  velocityAverage.data()[ebN_kb_nSpace+0]=0.5*(u_left + u_right);
4479  velocityAverage.data()[ebN_kb_nSpace+1]=0.5*(v_left + v_right);
4480  /* velocityAverage.data()[ebN_kb_nSpace+2]=0.5*(w_left + w_right); */
4481  }//ebNI
4482  }
4483  }
4484  };//RANS3PSed2D
4485 
4486 
4487  inline cppRANS3PSed2D_base* newRANS3PSed2D(int nSpaceIn,
4488  int nQuadraturePoints_elementIn,
4489  int nDOF_mesh_trial_elementIn,
4490  int nDOF_trial_elementIn,
4491  int nDOF_test_elementIn,
4492  int nQuadraturePoints_elementBoundaryIn,
4493  int CompKernelFlag,
4494  double aDarcy,
4495  double betaForch,
4496  double grain,
4497  double packFraction,
4498  double packMargin,
4499  double maxFraction,
4500  double frFraction,
4501  double sigmaC,
4502  double C3e,
4503  double C4e,
4504  double eR,
4505  double fContact,
4506  double mContact,
4507  double nContact,
4508  double angFriction,
4509  double vos_limiter,
4510  double mu_fr_limiter)
4511  {
4512  cppRANS3PSed2D_base* rvalue = proteus::chooseAndAllocateDiscretization2D
4513  <cppRANS3PSed2D_base,cppRANS3PSed2D,CompKernel>
4514  (nSpaceIn,
4515  nQuadraturePoints_elementIn,
4516  nDOF_mesh_trial_elementIn,
4517  nDOF_trial_elementIn,
4518  nDOF_test_elementIn,
4519  nQuadraturePoints_elementBoundaryIn,
4520  CompKernelFlag);
4521  rvalue->setSedClosure(aDarcy,
4522  betaForch,
4523  grain,
4524  packFraction,
4525  packMargin,
4526  maxFraction,
4527  frFraction,
4528  sigmaC,
4529  C3e,
4530  C4e,
4531  eR,
4532  fContact,
4533  mContact,
4534  nContact,
4535  angFriction,
4536  vos_limiter,
4537  mu_fr_limiter);
4538  return rvalue;
4539  }
4540 } //proteus
4541 
4542 #endif
proteus::cppRANS3PSed2D_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::cppRANS3PSed2D::updateFrictionalPressure
void updateFrictionalPressure(const double vos, const double grad_vos[nSpace], double &mom_u_source, double &mom_v_source, double &mom_w_source)
Definition: RANS3PSed2D.h:554
proteus::newRANS3PSed2D
cppRANS3PSed2D_base * newRANS3PSed2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag, double aDarcy, double betaForch, double grain, double packFraction, double packMargin, double maxFraction, double frFraction, double sigmaC, double C3e, double C4e, double eR, double fContact, double mContact, double nContact, double angFriction, double vos_limiter, double mu_fr_limiter)
Definition: RANS3PSed2D.h:4487
proteus::cppRANS3PSed2D::calculateSubgridError_tauRes
void calculateSubgridError_tauRes(const double &tau_p, const double &tau_v, const double &pdeResidualP, const double &pdeResidualU, const double &pdeResidualV, const double &pdeResidualW, double &subgridErrorP, double &subgridErrorU, double &subgridErrorV, double &subgridErrorW)
Definition: RANS3PSed2D.h:686
proteus::cppRANS3PSed2D::updateDarcyForchheimerTerms_Ergun
void updateDarcyForchheimerTerms_Ergun(const double alpha, const double beta, const double eps_rho, const double eps_mu, const double rho_0, const double nu_0, const double rho_1, const double nu_1, double nu_t, const double useVF, const double vf, const double phi, const double u, const double v, const double w, const double uStar, const double vStar, const double wStar, const double eps_s, const double vos, const double u_f, const double v_f, const double w_f, const double uStar_f, const double vStar_f, const double wStar_f, double &mom_u_source, double &mom_v_source, double &mom_w_source, double dmom_u_source[nSpace], double dmom_v_source[nSpace], double dmom_w_source[nSpace], double gradC_x, double gradC_y, double gradC_z)
Definition: RANS3PSed2D.h:446
w
#define w(x)
Definition: jf.h:22
DM
const double DM
Definition: RANS3PSed2D.h:13
HI
#define HI
Definition: Headers.h:4
DM3
const double DM3
Definition: RANS3PSed2D.h:15
proteus::cppRANS3PSed2D::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isDOFBoundary_p, const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isDOFBoundary_w, const int &isFluxBoundary_p, const int &isFluxBoundary_u, const int &isFluxBoundary_v, const int &isFluxBoundary_w, const double &oneByRho, const double &bc_oneByRho, const double n[nSpace], const double &vos, const double &bc_p, const double &bc_u, const double &bc_v, const double &bc_w, const double bc_f_mass[nSpace], const double bc_f_umom[nSpace], const double bc_f_vmom[nSpace], const double bc_f_wmom[nSpace], const double &bc_flux_mass, const double &bc_flux_umom, const double &bc_flux_vmom, const double &bc_flux_wmom, const double &p, const double &u, const double &v, const double &w, const double f_mass[nSpace], const double f_umom[nSpace], const double f_vmom[nSpace], const double f_wmom[nSpace], const double df_mass_du[nSpace], const double df_mass_dv[nSpace], const double df_mass_dw[nSpace], const double df_umom_dp[nSpace], const double df_umom_du[nSpace], const double df_umom_dv[nSpace], const double df_umom_dw[nSpace], const double df_vmom_dp[nSpace], const double df_vmom_du[nSpace], const double df_vmom_dv[nSpace], const double df_vmom_dw[nSpace], const double df_wmom_dp[nSpace], const double df_wmom_du[nSpace], const double df_wmom_dv[nSpace], const double df_wmom_dw[nSpace], double &flux_mass, double &flux_umom, double &flux_vmom, double &flux_wmom, double *velocity_star, double *velocity)
Definition: RANS3PSed2D.h:746
proteus::cppRANS3PSed2D::calculateSubgridError_tau
void calculateSubgridError_tau(const double &hFactor, const double &elementDiameter, const double &dmt, const double &dm, const double df[nSpace], const double &a, const double &pfac, double &tau_v, double &tau_p, double &cfl)
Definition: RANS3PSed2D.h:631
coeff
Double * coeff
Definition: Headers.h:42
proteus::cppRANS3PSed2D::smoothedHeaviside_integral
double smoothedHeaviside_integral(double eps, double phi)
Definition: RANS3PSed2D.h:135
proteus::cppRANS3PSed2D::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: RANS3PSed2D.h:1117
proteus::cppRANS3PSed2D_base
Definition: RANS3PSed2D.h:22
proteus::cppRANS3PSed2D::calculateSubgridError_tau
void calculateSubgridError_tau(const double &Ct_sge, const double &Cd_sge, const double G[nSpace *nSpace], const double &G_dd_G, const double &tr_G, const double &A0, const double Ai[nSpace], const double &Kij, const double &pfac, double &tau_v, double &tau_p, double &q_cfl)
Definition: RANS3PSed2D.h:665
proteus::cppRANS3PSed2D::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: RANS3PSed2D.h:59
n
Int n
Definition: Headers.h:28
proteus::cppRANS3PSed2D_base::~cppRANS3PSed2D_base
virtual ~cppRANS3PSed2D_base()
Definition: RANS3PSed2D.h:24
DRAG_FAC
#define DRAG_FAC
Definition: RANS3PSed2D.h:16
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
proteus::cppRANS3PSed2D::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: RANS3PSed2D.h:2697
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::cppHsuSedStress::frFraction_
double frFraction_
Definition: SedClosure.h:930
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
SedClosure.h
proteus::cppRANS3PSed2D::evaluateCoefficients
void evaluateCoefficients(const double eps_rho, const double eps_mu, const double sigma, const double rho_0, double nu_0, const double rho_1, double nu_1, const double rho_s, const double h_e, const double smagorinskyConstant, const int turbulenceClosureModel, const double g[nSpace], const double useVF, const double &vf, const double &phi, const double n[nSpace], const double &kappa, const double vos, const double &p, const double grad_p[nSpace], const double grad_u[nSpace], const double grad_v[nSpace], const double grad_w[nSpace], const double &u, const double &v, const double &w, const double &uStar, const double &vStar, const double &wStar, double &eddy_viscosity, double &mom_u_acc, double &dmom_u_acc_u, double &mom_v_acc, double &dmom_v_acc_v, double &mom_w_acc, double &dmom_w_acc_w, double mass_adv[nSpace], double dmass_adv_u[nSpace], double dmass_adv_v[nSpace], double dmass_adv_w[nSpace], double mom_u_adv[nSpace], double dmom_u_adv_u[nSpace], double dmom_u_adv_v[nSpace], double dmom_u_adv_w[nSpace], double mom_v_adv[nSpace], double dmom_v_adv_u[nSpace], double dmom_v_adv_v[nSpace], double dmom_v_adv_w[nSpace], double mom_w_adv[nSpace], double dmom_w_adv_u[nSpace], double dmom_w_adv_v[nSpace], double dmom_w_adv_w[nSpace], double mom_uu_diff_ten[nSpace], double mom_vv_diff_ten[nSpace], double mom_ww_diff_ten[nSpace], double mom_uv_diff_ten[1], double mom_uw_diff_ten[1], double mom_vu_diff_ten[1], double mom_vw_diff_ten[1], double mom_wu_diff_ten[1], double mom_wv_diff_ten[1], double &mom_u_source, double &mom_v_source, double &mom_w_source, double &mom_u_ham, double dmom_u_ham_grad_p[nSpace], double dmom_u_ham_grad_u[nSpace], double &mom_v_ham, double dmom_v_ham_grad_p[nSpace], double dmom_v_ham_grad_v[nSpace], double &mom_w_ham, double dmom_w_ham_grad_p[nSpace], double dmom_w_ham_grad_w[nSpace])
Definition: RANS3PSed2D.h:169
proteus::cppRANS3PSed2D_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
H
Double H
Definition: Headers.h:65
DM2
const double DM2
Definition: RANS3PSed2D.h:14
nu_0
double nu_0
Definition: ErrorResidualMethod.cpp:22
v
Double v
Definition: Headers.h:95
proteus::cppHsuSedStress::maxFraction_
double maxFraction_
Definition: SedClosure.h:931
nu_1
double nu_1
Definition: ErrorResidualMethod.cpp:22
proteus::cppRANS3PSed2D::updateFrictionalStress
void updateFrictionalStress(const double LAG_MU_FR, const double mu_fr_last, double &mu_fr_new, const double vos, const double eps_rho, const double eps_mu, const double rho_0, const double nu_0, const double rho_1, const double nu_1, const double rho_s, const double useVF, const double vf, const double phi, const double grad_u[nSpace], const double grad_v[nSpace], const double grad_w[nSpace], double mom_uu_diff_ten[nSpace], double mom_uv_diff_ten[1], double mom_uw_diff_ten[1], double mom_vv_diff_ten[nSpace], double mom_vu_diff_ten[1], double mom_vw_diff_ten[1], double mom_ww_diff_ten[nSpace], double mom_wu_diff_ten[1], double mom_wv_diff_ten[1])
Definition: RANS3PSed2D.h:570
proteus::cppHsuSedStress::betaCoeff
double betaCoeff(double sedF, double rhoFluid, const double uFluid[nSpace], const double uSolid[nSpace], double nu)
Definition: SedClosure.h:71
proteus::cppRANS3PSed2D::ExteriorNumericalDiffusiveFluxJacobian
double ExteriorNumericalDiffusiveFluxJacobian(const double &eps, const double &phi, int *rowptr, int *colind, const int &isDOFBoundary, const int &isFluxBoundary, const double n[nSpace], double *a, const double &v, const double grad_v[nSpace], const double &penalty)
Definition: RANS3PSed2D.h:1085
proteus::cppRANS3PSed2D::smoothedDirac
double smoothedDirac(double eps, double phi)
Definition: RANS3PSed2D.h:156
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
TURB_FORCE_FAC
#define TURB_FORCE_FAC
Definition: RANS3PSed2D.h:17
proteus::phi
double phi(const double &g, const double &h, const double &hL, const double &hR, const double &uL, const double &uR)
Definition: SW2DCV.h:62
proteus::cppHsuSedStress::gradp_friction
double gradp_friction(double sedF)
Definition: SedClosure.h:609
proteus::cppRANS3PSed2D::calculateSubgridErrorDerivatives_tauRes
void calculateSubgridErrorDerivatives_tauRes(const double &tau_p, const double &tau_v, const double dpdeResidualP_du[nDOF_trial_element], const double dpdeResidualP_dv[nDOF_trial_element], const double dpdeResidualP_dw[nDOF_trial_element], const double dpdeResidualU_dp[nDOF_trial_element], const double dpdeResidualU_du[nDOF_trial_element], const double dpdeResidualV_dp[nDOF_trial_element], const double dpdeResidualV_dv[nDOF_trial_element], const double dpdeResidualW_dp[nDOF_trial_element], const double dpdeResidualW_dw[nDOF_trial_element], double dsubgridErrorP_du[nDOF_trial_element], double dsubgridErrorP_dv[nDOF_trial_element], double dsubgridErrorP_dw[nDOF_trial_element], double dsubgridErrorU_dp[nDOF_trial_element], double dsubgridErrorU_du[nDOF_trial_element], double dsubgridErrorV_dp[nDOF_trial_element], double dsubgridErrorV_dv[nDOF_trial_element], double dsubgridErrorW_dp[nDOF_trial_element], double dsubgridErrorW_dw[nDOF_trial_element])
Definition: RANS3PSed2D.h:705
xt
Definition: AddedMass.cpp:7
rho_1
double rho_1
Definition: ErrorResidualMethod.cpp:22
proteus::cppRANS3PSed2D::exteriorNumericalDiffusiveFlux
void exteriorNumericalDiffusiveFlux(const double &eps, const double &phi, int *rowptr, int *colind, const int &isDOFBoundary, const int &isFluxBoundary, const double n[nSpace], double *bc_a, const double &bc_u, const double &bc_flux, double *a, const double grad_potential[nSpace], const double &u, const double &penalty, double &flux)
Definition: RANS3PSed2D.h:1037
proteus::cppRANS3PSed2D
Definition: RANS3PSed2D.h:56
proteus::cppHsuSedStress::mu_fr
double mu_fr(double sedF, double du_dx, double du_dy, double du_dz, double dv_dx, double dv_dy, double dv_dz, double dw_dx, double dw_dy, double dw_dz)
Definition: SedClosure.h:626
proteus::cppHsuSedStress::sigmaC_
double sigmaC_
Definition: SedClosure.h:932
CompKernel
Definition: CompKernel.h:1018
proteus::cppRANS3PSed2D::setSedClosure
void setSedClosure(double aDarcy, double betaForch, double grain, double packFraction, double packMargin, double maxFraction, double frFraction, double sigmaC, double C3e, double C4e, double eR, double fContact, double mContact, double nContact, double angFriction, double vos_limiter, double mu_fr_limiter)
Definition: RANS3PSed2D.h:84
proteus::cppRANS3PSed2D::smoothedHeaviside
double smoothedHeaviside(double eps, double phi)
Definition: RANS3PSed2D.h:121
proteus::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: MixedModelFactory.h:309
proteus::cppHsuSedStress::vos_limiter_
double vos_limiter_
Definition: SedClosure.h:943
proteus::cppRANS3PSed2D::exteriorNumericalAdvectiveFluxDerivatives
void exteriorNumericalAdvectiveFluxDerivatives(const int &isDOFBoundary_p, const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isDOFBoundary_w, const int &isFluxBoundary_p, const int &isFluxBoundary_u, const int &isFluxBoundary_v, const int &isFluxBoundary_w, const double &oneByRho, const double n[nSpace], const double &vos, const double &bc_p, const double &bc_u, const double &bc_v, const double &bc_w, const double bc_f_mass[nSpace], const double bc_f_umom[nSpace], const double bc_f_vmom[nSpace], const double bc_f_wmom[nSpace], const double &bc_flux_mass, const double &bc_flux_umom, const double &bc_flux_vmom, const double &bc_flux_wmom, const double &p, const double &u, const double &v, const double &w, const double f_mass[nSpace], const double f_umom[nSpace], const double f_vmom[nSpace], const double f_wmom[nSpace], const double df_mass_du[nSpace], const double df_mass_dv[nSpace], const double df_mass_dw[nSpace], const double df_umom_dp[nSpace], const double df_umom_du[nSpace], const double df_umom_dv[nSpace], const double df_umom_dw[nSpace], const double df_vmom_dp[nSpace], const double df_vmom_du[nSpace], const double df_vmom_dv[nSpace], const double df_vmom_dw[nSpace], const double df_wmom_dp[nSpace], const double df_wmom_du[nSpace], const double df_wmom_dv[nSpace], const double df_wmom_dw[nSpace], double &dflux_mass_du, double &dflux_mass_dv, double &dflux_mass_dw, double &dflux_umom_dp, double &dflux_umom_du, double &dflux_umom_dv, double &dflux_umom_dw, double &dflux_vmom_dp, double &dflux_vmom_du, double &dflux_vmom_dv, double &dflux_vmom_dw, double &dflux_wmom_dp, double &dflux_wmom_du, double &dflux_wmom_dv, double &dflux_wmom_dw, double *velocity_star)
Definition: RANS3PSed2D.h:877
proteus::cppRANS3PSed2D_base::calculateVelocityAverage
virtual void calculateVelocityAverage(arguments_dict &args)=0
proteus
Definition: ADR.h:17
proteus::cppRANS3PSed2D::closure
cppHsuSedStress< 2 > closure
Definition: RANS3PSed2D.h:58
proteus::cppHsuSedStress< 2 >
proteus::cppRANS3PSed2D::cppRANS3PSed2D
cppRANS3PSed2D()
Definition: RANS3PSed2D.h:61
proteus::cppRANS3PSed2D::updatePenaltyForPacking
void updatePenaltyForPacking(const double vos, const double u, const double v, const double w, double &mom_u_source, double &mom_v_source, double &mom_w_source, double dmom_u_source[nSpace], double dmom_v_source[nSpace], double dmom_w_source[nSpace])
Definition: RANS3PSed2D.h:527
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::cppRANS3PSed2D::calculateVelocityAverage
void calculateVelocityAverage(arguments_dict &args)
Definition: RANS3PSed2D.h:4317
rho_0
double rho_0
Definition: ErrorResidualMethod.cpp:22
proteus::cppRANS3PSed2D::ck
CompKernelType ck
Definition: RANS3PSed2D.h:60
proteus::cppRANS3PSed2D_base::setSedClosure
virtual void setSedClosure(double aDarcy, double betaForch, double grain, double packFraction, double packMargin, double maxFraction, double frFraction, double sigmaC, double C3e, double C4e, double eR, double fContact, double mContact, double nContact, double angFriction, double vos_limiter, double mu_fr_limiter)
Definition: RANS3PSed2D.h:25
ArgumentsDict.h
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20