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