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