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