8 #include "xtensor-python/pyarray.hpp"
10 namespace py = pybind11;
26 template<
class CompKernelType,
28 int nQuadraturePoints_element,
29 int nDOF_mesh_trial_element,
30 int nDOF_trial_element,
31 int nDOF_test_element,
32 int nQuadraturePoints_elementBoundary>
65 H = 0.5*(1.0 +
phi/eps + sin(M_PI*
phi/eps)/M_PI);
75 + 0.5*(eps + 0.5*eps*eps/eps - eps*cos(M_PI*eps/eps)/(M_PI*M_PI)) \
76 - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
84 HI = 0.5*(
phi + 0.5*
phi*
phi/eps - eps*cos(M_PI*
phi/eps)/(M_PI*M_PI)) \
85 - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
98 d = 0.5*(1.0 + cos(M_PI*
phi/eps))/eps;
108 delta= 1.0/8.0*(5.0+2.0*
r-sqrt(-7.0-12.0*
r-4*
r*
r));
110 delta= 1.0/8.0*(3.0+2.0*
r+sqrt(1.0-4.0*
r-4.0*
r*
r));
112 delta= 1.0/8.0*(3.0-2.0*
r+sqrt(1.0+4.0*
r-4.0*
r*
r));
114 delta= 1.0/8.0*(5.0-2.0*
r-sqrt(-7.0+12.0*
r-4.0*
r*
r));
122 const int nBeamElements,
123 const int beam_quadOrder,
124 const double beam_Cd,
125 const double* beamRadius,
132 const double *Beam_h,
133 const double *dV_beam,
137 const double eps_rho,
141 double& mom_u_source,
142 double& mom_v_source,
143 double& mom_w_source,
147 double rho, H_rho, delt,
vel;
154 for(
int I=0; I<nBeams; I++)
156 for(
int k=0;k<nBeamElements; k++)
158 for(
int l=0;l<beam_quadOrder; l++)
160 delt=
delta_h((x-xq[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l])/Beam_h[I*nBeamElements+k])*
delta_h((y-yq[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l])/Beam_h[I*nBeamElements+k])*
delta_h((
z-zq[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l])/Beam_h[I*nBeamElements+k])/(Beam_h[I*nBeamElements+k]*Beam_h[I*nBeamElements+k]*Beam_h[I*nBeamElements+k]);
161 mom_u_source += beam_Cd*rho*beamRadius[I]*
u*
vel*delt*dV_beam[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l];
162 mom_v_source += beam_Cd*rho*beamRadius[I]*
v*
vel*delt*dV_beam[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l];
163 mom_w_source += beam_Cd*rho*beamRadius[I]*
w*
vel*delt*dV_beam[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l];
171 netBeamDrag[0]+= dV*sqrt(mom_u_source*mom_u_source+mom_v_source*mom_v_source+mom_w_source*mom_w_source);
175 const int nBeamElements,
176 const int beam_quadOrder,
177 const double beam_Cd,
178 const double* beamRadius,
185 const double *Beam_h,
189 const double eps_rho,
199 double rho, H_rho, delt,
vel,h_save, buoy;
209 for(
int I=0; I<nBeams; I++)
213 for(
int k=0;k<nBeamElements; k++)
215 for(
int l=0;l<beam_quadOrder; l++)
218 h_save=0.25*Beam_h[I*nBeamElements+k];
220 h_save = 0.5*Beam_h[I*nBeamElements+k];
222 h_save = Beam_h[I*nBeamElements+k];
223 delt=
delta_h((x-xq[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l])/h_save)*
delta_h((y-yq[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l])/h_save)*
delta_h((
z-zq[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l])/h_save)/(h_save*h_save*h_save);
227 q1[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l] += beam_Cd*rho*beamRadius[I]*
u*
vel*delt*dV;
228 q2[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l] += beam_Cd*rho*beamRadius[I]*
v*
vel*delt*dV;
229 q3[I*nBeamElements*beam_quadOrder+beam_quadOrder*k+l] += (beam_Cd*rho*beamRadius[I]*
w*
vel+buoy)*delt*dV;
244 const double smagorinskyConstant,
245 const int turbulenceClosureModel,
246 const double g[nSpace],
250 const double n[nSpace],
252 const double porosity,
254 const double grad_p[nSpace],
255 const double grad_u[nSpace],
256 const double grad_v[nSpace],
257 const double grad_w[nSpace],
262 double& dmom_u_acc_u,
264 double& dmom_v_acc_v,
266 double& dmom_w_acc_w,
267 double mass_adv[nSpace],
268 double dmass_adv_u[nSpace],
269 double dmass_adv_v[nSpace],
270 double dmass_adv_w[nSpace],
271 double mom_u_adv[nSpace],
272 double dmom_u_adv_u[nSpace],
273 double dmom_u_adv_v[nSpace],
274 double dmom_u_adv_w[nSpace],
275 double mom_v_adv[nSpace],
276 double dmom_v_adv_u[nSpace],
277 double dmom_v_adv_v[nSpace],
278 double dmom_v_adv_w[nSpace],
279 double mom_w_adv[nSpace],
280 double dmom_w_adv_u[nSpace],
281 double dmom_w_adv_v[nSpace],
282 double dmom_w_adv_w[nSpace],
283 double mom_uu_diff_ten[nSpace],
284 double mom_vv_diff_ten[nSpace],
285 double mom_ww_diff_ten[nSpace],
286 double mom_uv_diff_ten[1],
287 double mom_uw_diff_ten[1],
288 double mom_vu_diff_ten[1],
289 double mom_vw_diff_ten[1],
290 double mom_wu_diff_ten[1],
291 double mom_wv_diff_ten[1],
292 double& mom_u_source,
293 double& mom_v_source,
294 double& mom_w_source,
296 double dmom_u_ham_grad_p[nSpace],
298 double dmom_v_ham_grad_p[nSpace],
300 double dmom_w_ham_grad_p[nSpace],
301 const double& dragBeam1,
302 const double& dragBeam2,
303 const double& dragBeam3)
305 double rho,nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n;
312 switch (turbulenceClosureModel)
317 norm_S = sqrt(2.0*(grad_u[0]*grad_u[0] + grad_v[1]*grad_v[1] + grad_w[2]*grad_w[2] +
318 0.5*(grad_u[1]+grad_v[0])*(grad_u[1]+grad_v[0]) +
319 0.5*(grad_u[2]+grad_w[0])*(grad_u[2]+grad_w[0]) +
320 0.5*(grad_v[2]+grad_w[1])*(grad_v[2]+grad_w[1])));
321 nu_0 += smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
322 nu_1 += smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
326 double re_0,cs_0,re_1,cs_1;
327 norm_S = sqrt(2.0*(grad_u[0]*grad_u[0] + grad_v[1]*grad_v[1] + grad_w[2]*grad_w[2] +
328 0.5*(grad_u[1]+grad_v[0])*(grad_u[1]+grad_v[0]) +
329 0.5*(grad_u[2]+grad_w[0])*(grad_u[2]+grad_w[0]) +
330 0.5*(grad_v[2]+grad_w[1])*(grad_v[2]+grad_w[1])));
331 re_0 = h_e*h_e*norm_S/
nu_0;
332 cs_0=0.027*pow(10.0,-3.23*pow(re_0,-0.92));
333 nu_0 += cs_0*h_e*h_e*norm_S;
334 re_1 = h_e*h_e*norm_S/
nu_1;
335 cs_1=0.027*pow(10.0,-3.23*pow(re_1,-0.92));
336 nu_1 += cs_1*h_e*h_e*norm_S;
344 #ifdef COMPRESSIBLE_FORM
346 mom_u_acc=porosity*rho*
u;
347 dmom_u_acc_u=porosity*rho;
350 mom_v_acc=porosity*rho*
v;
351 dmom_v_acc_v=porosity*rho;
354 mom_w_acc=porosity*rho*
w;
355 dmom_w_acc_w=porosity*rho;
359 mass_adv[0]=porosity*
u;
360 mass_adv[1]=porosity*
v;
361 mass_adv[2]=porosity*
w;
363 dmass_adv_u[0]=porosity;
368 dmass_adv_v[1]=porosity;
373 dmass_adv_w[2]=porosity;
376 mom_u_adv[0]=porosity*rho*
u*
u;
377 mom_u_adv[1]=porosity*rho*
u*
v;
378 mom_u_adv[2]=porosity*rho*
u*
w;
380 dmom_u_adv_u[0]=porosity*rho*2.0*
u;
381 dmom_u_adv_u[1]=porosity*rho*
v;
382 dmom_u_adv_u[2]=porosity*rho*
w;
385 dmom_u_adv_v[1]=porosity*rho*
u;
390 dmom_u_adv_w[2]=porosity*rho*
u;
393 mom_v_adv[0]=porosity*rho*
v*
u;
394 mom_v_adv[1]=porosity*rho*
v*
v;
395 mom_v_adv[2]=porosity*rho*
v*
w;
397 dmom_v_adv_u[0]=porosity*rho*
v;
403 dmom_v_adv_w[2]=porosity*rho*
v;
405 dmom_v_adv_v[0]=porosity*rho*
u;
406 dmom_v_adv_v[1]=porosity*rho*2.0*
v;
407 dmom_v_adv_v[2]=porosity*rho*
w;
410 mom_w_adv[0]=porosity*rho*
w*
u;
411 mom_w_adv[1]=porosity*rho*
w*
v;
412 mom_w_adv[2]=porosity*rho*
w*
w;
414 dmom_w_adv_u[0]=porosity*rho*
w;
419 dmom_w_adv_v[1]=porosity*rho*
w;
422 dmom_w_adv_w[0]=porosity*rho*
u;
423 dmom_w_adv_w[1]=porosity*rho*
v;
424 dmom_w_adv_w[2]=porosity*rho*2.0*
w;
427 mom_uu_diff_ten[0] = 2.0*porosity*mu;
428 mom_uu_diff_ten[1] = porosity*mu;
429 mom_uu_diff_ten[2] = porosity*mu;
431 mom_uv_diff_ten[0]=porosity*mu;
433 mom_uw_diff_ten[0]=porosity*mu;
436 mom_vv_diff_ten[0] = porosity*mu;
437 mom_vv_diff_ten[1] = 2.0*porosity*mu;
438 mom_vv_diff_ten[2] = porosity*mu;
440 mom_vu_diff_ten[0]=porosity*mu;
442 mom_vw_diff_ten[0]=porosity*mu;
445 mom_ww_diff_ten[0] = porosity*mu;
446 mom_ww_diff_ten[1] = porosity*mu;
447 mom_ww_diff_ten[2] = 2.0*porosity*mu;
449 mom_wu_diff_ten[0]=porosity*mu;
451 mom_wv_diff_ten[0]=porosity*mu;
454 norm_n = sqrt(
n[0]*
n[0]+
n[1]*
n[1]+
n[2]*
n[2]);
455 mom_u_source = dragBeam1-porosity*rho*g[0];
456 mom_v_source = dragBeam2-porosity*rho*g[1];
457 mom_w_source = dragBeam3-porosity*rho*g[2];
460 mom_u_ham = porosity*grad_p[0];
461 dmom_u_ham_grad_p[0]=porosity;
462 dmom_u_ham_grad_p[1]=0.0;
463 dmom_u_ham_grad_p[2]=0.0;
466 mom_v_ham = porosity*grad_p[1];
467 dmom_v_ham_grad_p[0]=0.0;
468 dmom_v_ham_grad_p[1]=porosity;
469 dmom_v_ham_grad_p[2]=0.0;
472 mom_w_ham = porosity*grad_p[2];
473 dmom_w_ham_grad_p[0]=0.0;
474 dmom_w_ham_grad_p[1]=0.0;
475 dmom_w_ham_grad_p[2]=porosity;
478 mom_u_acc=porosity*
u;
479 dmom_u_acc_u=porosity;
482 mom_v_acc=porosity*
v;
483 dmom_v_acc_v=porosity;
486 mom_w_acc=porosity*
w;
487 dmom_w_acc_w=porosity;
491 mass_adv[0]=porosity*
u;
492 mass_adv[1]=porosity*
v;
493 mass_adv[2]=porosity*
w;
495 dmass_adv_u[0]=porosity;
500 dmass_adv_v[1]=porosity;
505 dmass_adv_w[2]=porosity;
508 mom_u_adv[0]=porosity*
u*
u;
509 mom_u_adv[1]=porosity*
u*
v;
510 mom_u_adv[2]=porosity*
u*
w;
512 dmom_u_adv_u[0]=2.0*porosity*
u;
513 dmom_u_adv_u[1]=porosity*
v;
514 dmom_u_adv_u[2]=porosity*
w;
517 dmom_u_adv_v[1]=porosity*
u;
522 dmom_u_adv_w[2]=porosity*
u;
525 mom_v_adv[0]=porosity*
v*
u;
526 mom_v_adv[1]=porosity*
v*
v;
527 mom_v_adv[2]=porosity*
v*
w;
529 dmom_v_adv_u[0]=porosity*
v;
535 dmom_v_adv_w[2]=porosity*
v;
537 dmom_v_adv_v[0]=porosity*
u;
538 dmom_v_adv_v[1]=2.0*porosity*
v;
539 dmom_v_adv_v[2]=porosity*
w;
542 mom_w_adv[0]=porosity*
w*
u;
543 mom_w_adv[1]=porosity*
w*
v;
544 mom_w_adv[2]=porosity*
w*
w;
546 dmom_w_adv_u[0]=porosity*
w;
551 dmom_w_adv_v[1]=porosity*
w;
554 dmom_w_adv_w[0]=porosity*
u;
555 dmom_w_adv_w[1]=porosity*
v;
556 dmom_w_adv_w[2]=2.0*porosity*
w;
559 mom_uu_diff_ten[0] = 2.0*porosity*nu;
560 mom_uu_diff_ten[1] = porosity*nu;
561 mom_uu_diff_ten[2] = porosity*nu;
563 mom_uv_diff_ten[0]=porosity*nu;
565 mom_uw_diff_ten[0]=porosity*nu;
568 mom_vv_diff_ten[0] = porosity*nu;
569 mom_vv_diff_ten[1] = 2.0*porosity*nu;
570 mom_vv_diff_ten[2] = porosity*nu;
572 mom_vu_diff_ten[0]=porosity*nu;
574 mom_vw_diff_ten[0]=porosity*nu;
577 mom_ww_diff_ten[0] = porosity*nu;
578 mom_ww_diff_ten[1] = porosity*nu;
579 mom_ww_diff_ten[2] = 2.0*porosity*nu;
581 mom_wu_diff_ten[0]=porosity*nu;
583 mom_wv_diff_ten[0]=porosity*nu;
586 norm_n = sqrt(
n[0]*
n[0]+
n[1]*
n[1]+
n[2]*
n[2]);
587 mom_u_source = dragBeam1/rho-porosity*g[0];
588 mom_v_source = dragBeam2/rho-porosity*g[1];
589 mom_w_source = dragBeam3/rho-porosity*g[2];
592 mom_u_ham = porosity*grad_p[0]/rho;
593 dmom_u_ham_grad_p[0]=porosity/rho;
594 dmom_u_ham_grad_p[1]=0.0;
595 dmom_u_ham_grad_p[2]=0.0;
598 mom_v_ham = porosity*grad_p[1]/rho;
599 dmom_v_ham_grad_p[0]=0.0;
600 dmom_v_ham_grad_p[1]=porosity/rho;
601 dmom_v_ham_grad_p[2]=0.0;
604 mom_w_ham = porosity*grad_p[2]/rho;
605 dmom_w_ham_grad_p[0]=0.0;
606 dmom_w_ham_grad_p[1]=0.0;
607 dmom_w_ham_grad_p[2]=porosity/rho;
618 const double eps_rho,
635 double& mom_u_source,
636 double& mom_v_source,
637 double& mom_w_source,
638 double dmom_u_source[nSpace],
639 double dmom_v_source[nSpace],
640 double dmom_w_source[nSpace])
642 double mu,nu,H_mu,uc,duc_du,duc_dv,duc_dw,viscosity,H_s;
646 #ifdef COMPRESSIBLE_FORM
654 duc_du =
u/(uc+1.0e-12);
655 duc_dv =
v/(uc+1.0e-12);
656 duc_dw =
w/(uc+1.0e-12);
658 mom_u_source += H_s*(viscosity*alpha + beta*uc)*(
u-u_s);
659 mom_v_source += H_s*(viscosity*alpha + beta*uc)*(
v-v_s);
660 mom_w_source += H_s*(viscosity*alpha + beta*uc)*(
w-w_s);
662 dmom_u_source[0] = H_s*(viscosity*alpha + beta*(uc +
u*duc_du));
663 dmom_u_source[1] = H_s*beta*
u*duc_dv;
664 dmom_u_source[2] = H_s*beta*
u*duc_dw;
666 dmom_v_source[0] = H_s*beta*
v*duc_du;
667 dmom_v_source[1] = H_s*(viscosity*alpha + beta*(uc +
v*duc_dv));
668 dmom_v_source[2] = H_s*beta*
w*duc_dw;
670 dmom_w_source[0] = H_s*beta*
w*duc_du;
671 dmom_w_source[1] = H_s*beta*
w*duc_dv;
672 dmom_w_source[2] = H_s*(viscosity*alpha + beta*(uc +
w*duc_dw));
677 const double eps_rho,
686 const double porosity,
687 const double eddy_visc_coef_0,
688 const double turb_var_0,
689 const double turb_var_1,
690 const double turb_grad_0[nSpace],
691 double mom_uu_diff_ten[nSpace],
692 double mom_vv_diff_ten[nSpace],
693 double mom_ww_diff_ten[nSpace],
694 double mom_uv_diff_ten[1],
695 double mom_uw_diff_ten[1],
696 double mom_vu_diff_ten[1],
697 double mom_vw_diff_ten[1],
698 double mom_wu_diff_ten[1],
699 double mom_wv_diff_ten[1],
700 double& mom_u_source,
701 double& mom_v_source,
702 double& mom_w_source)
710 assert (turbulenceClosureModel >=3);
711 double rho,nu,H_mu,eddy_viscosity,nu_t=0.0,nu_t_keps =0.0, nu_t_komega=0.0;
712 double isKEpsilon = 1.0;
713 if (turbulenceClosureModel == 4)
719 const double twoThirds = 2.0/3.0;
const double div_zero = 1.0e-2*fmin(
nu_0,
nu_1);
720 mom_u_source += twoThirds*turb_grad_0[0];
721 mom_v_source += twoThirds*turb_grad_0[1];
722 mom_w_source += twoThirds*turb_grad_0[2];
726 nu_t_keps = eddy_visc_coef_0*turb_var_0*turb_var_0/(fabs(turb_var_1) + div_zero);
728 nu_t_komega = turb_var_0/(fabs(turb_var_1) + div_zero);
730 nu_t = isKEpsilon*nu_t_keps + (1.0-isKEpsilon)*nu_t_komega;
737 nu_t = fmax(nu_t,1.0e-4*nu);
739 nu_t = fmin(nu_t,1.0e6*nu);
741 #ifdef COMPRESSIBLE_FORM
742 eddy_viscosity = nu_t*rho;
744 mom_uu_diff_ten[0] += 2.0*porosity*eddy_viscosity;
745 mom_uu_diff_ten[1] += porosity*eddy_viscosity;
746 mom_uu_diff_ten[2] += porosity*eddy_viscosity;
748 mom_uv_diff_ten[0] +=porosity*eddy_viscosity;
750 mom_uw_diff_ten[0] +=porosity*eddy_viscosity;
753 mom_vv_diff_ten[0] += porosity*eddy_viscosity;
754 mom_vv_diff_ten[1] += 2.0*porosity*eddy_viscosity;
755 mom_vv_diff_ten[2] += porosity*eddy_viscosity;
757 mom_vu_diff_ten[0] += porosity*eddy_viscosity;
759 mom_vw_diff_ten[0] += porosity*eddy_viscosity;
762 mom_ww_diff_ten[0] += porosity*eddy_viscosity;
763 mom_ww_diff_ten[1] += porosity*eddy_viscosity;
764 mom_ww_diff_ten[2] += 2.0*porosity*eddy_viscosity;
766 mom_wu_diff_ten[0] += porosity*eddy_viscosity;
768 mom_wv_diff_ten[0] += porosity*eddy_viscosity;
771 eddy_viscosity = nu_t;
773 mom_uu_diff_ten[0] += 2.0*porosity*eddy_viscosity;
774 mom_uu_diff_ten[1] += porosity*eddy_viscosity;
775 mom_uu_diff_ten[2] += porosity*eddy_viscosity;
777 mom_uv_diff_ten[0]+=porosity*eddy_viscosity;
779 mom_uw_diff_ten[0]+=porosity*eddy_viscosity;
782 mom_vv_diff_ten[0] += porosity*eddy_viscosity;
783 mom_vv_diff_ten[1] += 2.0*porosity*eddy_viscosity;
784 mom_vv_diff_ten[2] += porosity*eddy_viscosity;
786 mom_vu_diff_ten[0]+=porosity*eddy_viscosity;
788 mom_vw_diff_ten[0]+=porosity*eddy_viscosity;
791 mom_ww_diff_ten[0] += porosity*eddy_viscosity;
792 mom_ww_diff_ten[1] += porosity*eddy_viscosity;
793 mom_ww_diff_ten[2] += 2.0*porosity*eddy_viscosity;
795 mom_wu_diff_ten[0]+=porosity*eddy_viscosity;
797 mom_wv_diff_ten[0]+=porosity*eddy_viscosity;
805 const double& elementDiameter,
808 const double df[nSpace],
815 double h,oneByAbsdt,density,viscosity,nrm_df;
816 h = hFactor*elementDiameter;
820 for(
int I=0;I<nSpace;I++)
822 nrm_df = sqrt(nrm_df);
823 cfl = nrm_df/(h*density);
824 oneByAbsdt = fabs(dmt);
825 tau_v = 1.0/(4.0*viscosity/(h*h) + 2.0*nrm_df/h + oneByAbsdt);
826 tau_p = (4.0*viscosity + 2.0*nrm_df*h + oneByAbsdt*h*h)/pfac;
832 const double& Cd_sge,
833 const double G[nSpace*nSpace],
834 const double& G_dd_G,
837 const double Ai[nSpace],
845 for(
int I=0;I<nSpace;I++)
846 for (
int J=0;J<nSpace;J++)
847 v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
849 tau_v = 1.0/sqrt(Ct_sge*A0*A0 + v_d_Gv + Cd_sge*Kij*Kij*G_dd_G + 1.0e-12);
850 tau_p = 1.0/(pfac*tr_G*tau_v);
856 const double& pdeResidualP,
857 const double& pdeResidualU,
858 const double& pdeResidualV,
859 const double& pdeResidualW,
860 double& subgridErrorP,
861 double& subgridErrorU,
862 double& subgridErrorV,
863 double& subgridErrorW)
866 subgridErrorP = -tau_p*pdeResidualP;
868 subgridErrorU = -tau_v*pdeResidualU;
869 subgridErrorV = -tau_v*pdeResidualV;
870 subgridErrorW = -tau_v*pdeResidualW;
876 const double dpdeResidualP_du[nDOF_trial_element],
877 const double dpdeResidualP_dv[nDOF_trial_element],
878 const double dpdeResidualP_dw[nDOF_trial_element],
879 const double dpdeResidualU_dp[nDOF_trial_element],
880 const double dpdeResidualU_du[nDOF_trial_element],
881 const double dpdeResidualV_dp[nDOF_trial_element],
882 const double dpdeResidualV_dv[nDOF_trial_element],
883 const double dpdeResidualW_dp[nDOF_trial_element],
884 const double dpdeResidualW_dw[nDOF_trial_element],
885 double dsubgridErrorP_du[nDOF_trial_element],
886 double dsubgridErrorP_dv[nDOF_trial_element],
887 double dsubgridErrorP_dw[nDOF_trial_element],
888 double dsubgridErrorU_dp[nDOF_trial_element],
889 double dsubgridErrorU_du[nDOF_trial_element],
890 double dsubgridErrorV_dp[nDOF_trial_element],
891 double dsubgridErrorV_dv[nDOF_trial_element],
892 double dsubgridErrorW_dp[nDOF_trial_element],
893 double dsubgridErrorW_dw[nDOF_trial_element])
895 for (
int j=0;j<nDOF_trial_element;j++)
898 dsubgridErrorP_du[j] = -tau_p*dpdeResidualP_du[j];
899 dsubgridErrorP_dv[j] = -tau_p*dpdeResidualP_dv[j];
900 dsubgridErrorP_dw[j] = -tau_p*dpdeResidualP_dw[j];
903 dsubgridErrorU_dp[j] = -tau_v*dpdeResidualU_dp[j];
904 dsubgridErrorU_du[j] = -tau_v*dpdeResidualU_du[j];
906 dsubgridErrorV_dp[j] = -tau_v*dpdeResidualV_dp[j];
907 dsubgridErrorV_dv[j] = -tau_v*dpdeResidualV_dv[j];
909 dsubgridErrorW_dp[j] = -tau_v*dpdeResidualW_dp[j];
910 dsubgridErrorW_dw[j] = -tau_v*dpdeResidualW_dw[j];
916 const int& isDOFBoundary_u,
917 const int& isDOFBoundary_v,
918 const int& isDOFBoundary_w,
919 const int& isFluxBoundary_p,
920 const int& isFluxBoundary_u,
921 const int& isFluxBoundary_v,
922 const int& isFluxBoundary_w,
923 const double& oneByRho,
924 const double& bc_oneByRho,
925 const double n[nSpace],
927 const double bc_f_mass[nSpace],
928 const double bc_f_umom[nSpace],
929 const double bc_f_vmom[nSpace],
930 const double bc_f_wmom[nSpace],
931 const double& bc_flux_mass,
932 const double& bc_flux_umom,
933 const double& bc_flux_vmom,
934 const double& bc_flux_wmom,
936 const double f_mass[nSpace],
937 const double f_umom[nSpace],
938 const double f_vmom[nSpace],
939 const double f_wmom[nSpace],
940 const double df_mass_du[nSpace],
941 const double df_mass_dv[nSpace],
942 const double df_mass_dw[nSpace],
943 const double df_umom_dp[nSpace],
944 const double df_umom_du[nSpace],
945 const double df_umom_dv[nSpace],
946 const double df_umom_dw[nSpace],
947 const double df_vmom_dp[nSpace],
948 const double df_vmom_du[nSpace],
949 const double df_vmom_dv[nSpace],
950 const double df_vmom_dw[nSpace],
951 const double df_wmom_dp[nSpace],
952 const double df_wmom_du[nSpace],
953 const double df_wmom_dv[nSpace],
954 const double df_wmom_dw[nSpace],
961 double flowDirection;
966 flowDirection=
n[0]*f_mass[0]+
n[1]*f_mass[1]+
n[2]*f_mass[2];
967 if (isDOFBoundary_u != 1)
969 flux_mass +=
n[0]*f_mass[0];
970 velocity[0] = f_mass[0];
971 if (flowDirection >= 0.0)
973 flux_umom +=
n[0]*f_umom[0];
974 flux_vmom +=
n[0]*f_vmom[0];
975 flux_wmom +=
n[0]*f_wmom[0];
980 flux_mass +=
n[0]*f_mass[0];
981 velocity[0] = f_mass[0];
982 if (flowDirection >= 0.0)
984 flux_umom +=
n[0]*f_umom[0];
985 flux_vmom +=
n[0]*f_vmom[0];
986 flux_wmom +=
n[0]*f_wmom[0];
990 flux_umom+=
n[0]*bc_f_umom[0];
991 flux_vmom+=
n[0]*bc_f_vmom[0];
992 flux_wmom+=
n[0]*bc_f_wmom[0];
995 if (isDOFBoundary_v != 1)
997 flux_mass+=
n[1]*f_mass[1];
998 velocity[1] = f_mass[1];
999 if (flowDirection >= 0.0)
1001 flux_umom+=
n[1]*f_umom[1];
1002 flux_vmom+=
n[1]*f_vmom[1];
1003 flux_wmom+=
n[1]*f_wmom[1];
1008 flux_mass+=
n[1]*f_mass[1];
1009 velocity[1] = f_mass[1];
1010 if (flowDirection >= 0.0)
1012 flux_umom+=
n[1]*f_umom[1];
1013 flux_vmom+=
n[1]*f_vmom[1];
1014 flux_wmom+=
n[1]*f_wmom[1];
1018 flux_umom+=
n[1]*bc_f_umom[1];
1019 flux_vmom+=
n[1]*bc_f_vmom[1];
1020 flux_wmom+=
n[1]*bc_f_wmom[1];
1023 if (isDOFBoundary_w != 1)
1025 flux_mass+=
n[2]*f_mass[2];
1026 velocity[2] = f_mass[2];
1027 if (flowDirection >= 0.0)
1029 flux_umom+=
n[2]*f_umom[2];
1030 flux_vmom+=
n[2]*f_vmom[2];
1031 flux_wmom+=
n[2]*f_wmom[2];
1036 flux_mass +=
n[2]*f_mass[2];
1037 velocity[2] = f_mass[2];
1038 if (flowDirection >= 0.0)
1040 flux_umom+=
n[2]*f_umom[2];
1041 flux_vmom+=
n[2]*f_vmom[2];
1042 flux_wmom+=
n[2]*f_wmom[2];
1046 flux_umom+=
n[2]*bc_f_umom[2];
1047 flux_vmom+=
n[2]*bc_f_vmom[2];
1048 flux_wmom+=
n[2]*bc_f_wmom[2];
1051 if (isDOFBoundary_p == 1)
1053 flux_umom+=
n[0]*(bc_p*bc_oneByRho-p*oneByRho);
1054 flux_vmom+=
n[1]*(bc_p*bc_oneByRho-p*oneByRho);
1055 flux_wmom+=
n[2]*(bc_p*bc_oneByRho-p*oneByRho);
1057 if (isFluxBoundary_p == 1)
1059 velocity[0] += (bc_flux_mass - flux_mass)*
n[0];
1060 velocity[1] += (bc_flux_mass - flux_mass)*
n[1];
1061 velocity[2] += (bc_flux_mass - flux_mass)*
n[2];
1062 flux_mass = bc_flux_mass;
1064 if (isFluxBoundary_u == 1)
1066 flux_umom = bc_flux_umom;
1068 if (isFluxBoundary_v == 1)
1070 flux_vmom = bc_flux_vmom;
1072 if (isFluxBoundary_w == 1)
1074 flux_wmom = bc_flux_wmom;
1090 const int& isDOFBoundary_u,
1091 const int& isDOFBoundary_v,
1092 const int& isDOFBoundary_w,
1093 const int& isFluxBoundary_p,
1094 const int& isFluxBoundary_u,
1095 const int& isFluxBoundary_v,
1096 const int& isFluxBoundary_w,
1097 const double& oneByRho,
1098 const double n[nSpace],
1100 const double bc_f_mass[nSpace],
1101 const double bc_f_umom[nSpace],
1102 const double bc_f_vmom[nSpace],
1103 const double bc_f_wmom[nSpace],
1104 const double& bc_flux_mass,
1105 const double& bc_flux_umom,
1106 const double& bc_flux_vmom,
1107 const double& bc_flux_wmom,
1109 const double f_mass[nSpace],
1110 const double f_umom[nSpace],
1111 const double f_vmom[nSpace],
1112 const double f_wmom[nSpace],
1113 const double df_mass_du[nSpace],
1114 const double df_mass_dv[nSpace],
1115 const double df_mass_dw[nSpace],
1116 const double df_umom_dp[nSpace],
1117 const double df_umom_du[nSpace],
1118 const double df_umom_dv[nSpace],
1119 const double df_umom_dw[nSpace],
1120 const double df_vmom_dp[nSpace],
1121 const double df_vmom_du[nSpace],
1122 const double df_vmom_dv[nSpace],
1123 const double df_vmom_dw[nSpace],
1124 const double df_wmom_dp[nSpace],
1125 const double df_wmom_du[nSpace],
1126 const double df_wmom_dv[nSpace],
1127 const double df_wmom_dw[nSpace],
1128 double& dflux_mass_du,
1129 double& dflux_mass_dv,
1130 double& dflux_mass_dw,
1131 double& dflux_umom_dp,
1132 double& dflux_umom_du,
1133 double& dflux_umom_dv,
1134 double& dflux_umom_dw,
1135 double& dflux_vmom_dp,
1136 double& dflux_vmom_du,
1137 double& dflux_vmom_dv,
1138 double& dflux_vmom_dw,
1139 double& dflux_wmom_dp,
1140 double& dflux_wmom_du,
1141 double& dflux_wmom_dv,
1142 double& dflux_wmom_dw)
1144 double flowDirection;
1145 dflux_mass_du = 0.0;
1146 dflux_mass_dv = 0.0;
1147 dflux_mass_dw = 0.0;
1149 dflux_umom_dp = 0.0;
1150 dflux_umom_du = 0.0;
1151 dflux_umom_dv = 0.0;
1152 dflux_umom_dw = 0.0;
1154 dflux_vmom_dp = 0.0;
1155 dflux_vmom_du = 0.0;
1156 dflux_vmom_dv = 0.0;
1157 dflux_vmom_dw = 0.0;
1159 dflux_wmom_dp = 0.0;
1160 dflux_wmom_du = 0.0;
1161 dflux_wmom_dv = 0.0;
1162 dflux_wmom_dw = 0.0;
1164 flowDirection=
n[0]*f_mass[0]+
n[1]*f_mass[1]+
n[2]*f_mass[2];
1165 if (isDOFBoundary_u != 1)
1167 dflux_mass_du +=
n[0]*df_mass_du[0];
1168 if (flowDirection >= 0.0)
1170 dflux_umom_du +=
n[0]*df_umom_du[0];
1171 dflux_vmom_du +=
n[0]*df_vmom_du[0];
1172 dflux_vmom_dv +=
n[0]*df_vmom_dv[0];
1173 dflux_wmom_du +=
n[0]*df_wmom_du[0];
1174 dflux_wmom_dw +=
n[0]*df_wmom_dw[0];
1180 dflux_mass_du +=
n[0]*df_mass_du[0];
1181 if (flowDirection >= 0.0)
1183 dflux_umom_du +=
n[0]*df_umom_du[0];
1184 dflux_vmom_du +=
n[0]*df_vmom_du[0];
1185 dflux_vmom_dv +=
n[0]*df_vmom_dv[0];
1186 dflux_wmom_du +=
n[0]*df_wmom_du[0];
1187 dflux_wmom_dw +=
n[0]*df_wmom_dw[0];
1191 if (isDOFBoundary_v != 1)
1192 dflux_vmom_dv +=
n[0]*df_vmom_dv[0];
1193 if (isDOFBoundary_w != 1)
1194 dflux_wmom_dw +=
n[0]*df_wmom_dw[0];
1197 if (isDOFBoundary_v != 1)
1199 dflux_mass_dv +=
n[1]*df_mass_dv[1];
1200 if (flowDirection >= 0.0)
1202 dflux_umom_du +=
n[1]*df_umom_du[1];
1203 dflux_umom_dv +=
n[1]*df_umom_dv[1];
1204 dflux_vmom_dv +=
n[1]*df_vmom_dv[1];
1205 dflux_wmom_dw +=
n[1]*df_wmom_dw[1];
1206 dflux_wmom_dv +=
n[1]*df_wmom_dv[1];
1212 dflux_mass_dv +=
n[1]*df_mass_dv[1];
1213 if (flowDirection >= 0.0)
1215 dflux_umom_du +=
n[1]*df_umom_du[1];
1216 dflux_umom_dv +=
n[1]*df_umom_dv[1];
1217 dflux_vmom_dv +=
n[1]*df_vmom_dv[1];
1218 dflux_wmom_dw +=
n[1]*df_wmom_dw[1];
1219 dflux_wmom_dv +=
n[1]*df_wmom_dv[1];
1223 if (isDOFBoundary_u != 1)
1224 dflux_umom_du +=
n[1]*df_umom_du[1];
1225 if (isDOFBoundary_w != 1)
1226 dflux_wmom_dw +=
n[1]*df_wmom_dw[1];
1229 if (isDOFBoundary_w != 1)
1231 dflux_mass_dw+=
n[2]*df_mass_dw[2];
1232 if (flowDirection >= 0.0)
1234 dflux_umom_du +=
n[2]*df_umom_du[2];
1235 dflux_umom_dw +=
n[2]*df_umom_dw[2];
1236 dflux_vmom_dv +=
n[2]*df_vmom_dv[2];
1237 dflux_vmom_dw +=
n[2]*df_vmom_dw[2];
1238 dflux_wmom_dw +=
n[2]*df_wmom_dw[2];
1244 dflux_mass_dw +=
n[2]*df_mass_dw[2];
1245 if (flowDirection >= 0.0)
1247 dflux_umom_du +=
n[2]*df_umom_du[2];
1248 dflux_umom_dw +=
n[2]*df_umom_dw[2];
1249 dflux_vmom_dv +=
n[2]*df_vmom_dv[2];
1250 dflux_vmom_dw +=
n[2]*df_vmom_dw[2];
1251 dflux_wmom_dw +=
n[2]*df_wmom_dw[2];
1255 if (isDOFBoundary_u != 1)
1256 dflux_umom_du +=
n[2]*df_umom_du[2];
1257 if (isDOFBoundary_v != 1)
1258 dflux_vmom_dv +=
n[2]*df_vmom_dv[2];
1261 if (isDOFBoundary_p == 1)
1263 dflux_umom_dp= -
n[0]*oneByRho;
1264 dflux_vmom_dp= -
n[1]*oneByRho;
1265 dflux_wmom_dp= -
n[2]*oneByRho;
1267 if (isFluxBoundary_p == 1)
1269 dflux_mass_du = 0.0;
1270 dflux_mass_dv = 0.0;
1271 dflux_mass_dw = 0.0;
1273 if (isFluxBoundary_u == 1)
1275 dflux_umom_dp = 0.0;
1276 dflux_umom_du = 0.0;
1277 dflux_umom_dv = 0.0;
1278 dflux_umom_dw = 0.0;
1280 if (isFluxBoundary_v == 1)
1282 dflux_vmom_dp = 0.0;
1283 dflux_vmom_du = 0.0;
1284 dflux_vmom_dv = 0.0;
1285 dflux_vmom_dw = 0.0;
1287 if (isFluxBoundary_w == 1)
1289 dflux_wmom_dp = 0.0;
1290 dflux_wmom_du = 0.0;
1291 dflux_wmom_dv = 0.0;
1292 dflux_wmom_dw = 0.0;
1301 const int& isDOFBoundary,
1302 const int& isFluxBoundary,
1303 const double n[nSpace],
1306 const double& bc_flux,
1308 const double grad_potential[nSpace],
1310 const double& penalty,
1313 double diffusiveVelocityComponent_I,penaltyFlux,max_a;
1314 if(isFluxBoundary == 1)
1318 else if(isDOFBoundary == 1)
1322 for(
int I=0;I<nSpace;I++)
1324 diffusiveVelocityComponent_I=0.0;
1325 for(
int m=rowptr[I];m<rowptr[I+1];m++)
1327 diffusiveVelocityComponent_I -= a[m]*grad_potential[colind[m]];
1328 max_a = fmax(max_a,a[m]);
1330 flux+= diffusiveVelocityComponent_I*
n[I];
1332 penaltyFlux = max_a*penalty*(
u-bc_u);
1333 flux += penaltyFlux;
1350 const int& isDOFBoundary,
1351 const int& isFluxBoundary,
1352 const double n[nSpace],
1355 const double grad_v[nSpace],
1356 const double& penalty)
1358 double dvel_I,tmp=0.0,max_a=0.0;
1359 if(isFluxBoundary==0 && isDOFBoundary==1)
1361 for(
int I=0;I<nSpace;I++)
1364 for(
int m=rowptr[I];m<rowptr[I+1];m++)
1366 dvel_I -= a[m]*grad_v[colind[m]];
1367 max_a = fmax(max_a,a[m]);
1371 tmp +=max_a*penalty*
v;
1380 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1381 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1382 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1383 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1384 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1385 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1386 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1387 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
1388 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
1389 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
1390 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
1391 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
1392 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
1393 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
1394 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
1395 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1396 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1397 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1398 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
1399 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
1400 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
1401 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
1402 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
1403 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
1404 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
1405 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
1406 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1407 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1408 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
1409 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1410 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1411 double hFactor = args.
scalar<
double>(
"hFactor");
1412 int nElements_global = args.
scalar<
int>(
"nElements_global");
1413 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
1414 double useRBLES = args.
scalar<
double>(
"useRBLES");
1415 double useMetrics = args.
scalar<
double>(
"useMetrics");
1416 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1417 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
1418 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
1419 double sigma = args.
scalar<
double>(
"sigma");
1424 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
1425 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
1426 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
1427 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
1428 double C_dc = args.
scalar<
double>(
"C_dc");
1429 double C_b = args.
scalar<
double>(
"C_b");
1430 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
1431 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
1432 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
1433 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
1434 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
1435 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
1436 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
1437 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
1438 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
1439 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
1440 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
1441 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
1442 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
1443 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1444 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
1445 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
1446 xt::pyarray<double>& g = args.
array<
double>(
"g");
1447 const double useVF = args.
scalar<
double>(
"useVF");
1448 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
1449 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
1450 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
1451 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
1452 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
1453 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
1454 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
1455 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
1456 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
1457 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
1458 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
1459 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
1460 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
1461 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1462 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
1463 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
1464 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1465 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
1466 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
1467 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
1468 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
1469 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
1470 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
1471 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
1472 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
1473 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
1474 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
1475 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
1476 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
1477 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
1478 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
1479 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
1480 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
1481 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
1482 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
1483 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
1484 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
1485 int offset_p = args.
scalar<
int>(
"offset_p");
1486 int offset_u = args.
scalar<
int>(
"offset_u");
1487 int offset_v = args.
scalar<
int>(
"offset_v");
1488 int offset_w = args.
scalar<
int>(
"offset_w");
1489 int stride_p = args.
scalar<
int>(
"stride_p");
1490 int stride_u = args.
scalar<
int>(
"stride_u");
1491 int stride_v = args.
scalar<
int>(
"stride_v");
1492 int stride_w = args.
scalar<
int>(
"stride_w");
1493 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1494 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1495 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1496 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1497 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1498 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
1499 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
1500 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
1501 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
1502 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
1503 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
1504 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
1505 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
1506 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
1507 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
1508 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1509 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
1510 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
1511 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
1512 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
1513 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
1514 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
1515 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
1516 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
1517 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
1518 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
1519 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
1520 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
1521 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
1522 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
1523 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1524 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
1525 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1526 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
1527 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
1528 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
1529 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
1530 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
1531 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
1532 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
1533 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
1534 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
1535 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
1536 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
1537 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
1538 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
1539 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
1540 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
1541 xt::pyarray<double>& q_dragBeam1 = args.
array<
double>(
"q_dragBeam1");
1542 xt::pyarray<double>& q_dragBeam2 = args.
array<
double>(
"q_dragBeam2");
1543 xt::pyarray<double>& q_dragBeam3 = args.
array<
double>(
"q_dragBeam3");
1544 xt::pyarray<double>& ebqe_dragBeam1 = args.
array<
double>(
"ebqe_dragBeam1");
1545 xt::pyarray<double>& ebqe_dragBeam2 = args.
array<
double>(
"ebqe_dragBeam2");
1546 xt::pyarray<double>& ebqe_dragBeam3 = args.
array<
double>(
"ebqe_dragBeam3");
1550 double globalConservationError=0.0;
1551 for(
int eN=0;eN<nElements_global;eN++)
1554 register double elementResidual_p[nDOF_test_element],
1555 elementResidual_u[nDOF_test_element],
1556 elementResidual_v[nDOF_test_element],
1557 elementResidual_w[nDOF_test_element],
1559 for (
int i=0;i<nDOF_test_element;i++)
1561 int eN_i = eN*nDOF_test_element+i;
1562 elementResidual_p_save.data()[eN_i]=0.0;
1563 elementResidual_p[i]=0.0;
1564 elementResidual_u[i]=0.0;
1565 elementResidual_v[i]=0.0;
1566 elementResidual_w[i]=0.0;
1571 for(
int k=0;k<nQuadraturePoints_element;k++)
1574 register int eN_k = eN*nQuadraturePoints_element+k,
1575 eN_k_nSpace = eN_k*nSpace,
1576 eN_nDOF_trial_element = eN*nDOF_trial_element;
1577 register double p=0.0,
u=0.0,
v=0.0,
w=0.0,
1578 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
1586 dmass_adv_u[nSpace],
1587 dmass_adv_v[nSpace],
1588 dmass_adv_w[nSpace],
1590 dmom_u_adv_u[nSpace],
1591 dmom_u_adv_v[nSpace],
1592 dmom_u_adv_w[nSpace],
1594 dmom_v_adv_u[nSpace],
1595 dmom_v_adv_v[nSpace],
1596 dmom_v_adv_w[nSpace],
1598 dmom_w_adv_u[nSpace],
1599 dmom_w_adv_v[nSpace],
1600 dmom_w_adv_w[nSpace],
1601 mom_uu_diff_ten[nSpace],
1602 mom_vv_diff_ten[nSpace],
1603 mom_ww_diff_ten[nSpace],
1614 dmom_u_ham_grad_p[nSpace],
1616 dmom_v_ham_grad_p[nSpace],
1618 dmom_w_ham_grad_p[nSpace],
1629 Lstar_u_p[nDOF_test_element],
1630 Lstar_v_p[nDOF_test_element],
1631 Lstar_w_p[nDOF_test_element],
1632 Lstar_u_u[nDOF_test_element],
1633 Lstar_v_v[nDOF_test_element],
1634 Lstar_w_w[nDOF_test_element],
1635 Lstar_p_u[nDOF_test_element],
1636 Lstar_p_v[nDOF_test_element],
1637 Lstar_p_w[nDOF_test_element],
1642 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
1643 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
1646 jacInv[nSpace*nSpace],
1647 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
1648 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
1649 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
1655 dmom_u_source[nSpace],
1656 dmom_v_source[nSpace],
1657 dmom_w_source[nSpace],
1659 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
1661 ck.calculateMapping_element(eN,
1665 mesh_trial_ref.data(),
1666 mesh_grad_trial_ref.data(),
1671 ck.calculateH_element(eN,
1673 nodeDiametersArray.data(),
1675 mesh_trial_ref.data(),
1678 ck.calculateMappingVelocity_element(eN,
1680 mesh_velocity_dof.data(),
1682 mesh_trial_ref.data(),
1687 dV = fabs(jacDet)*dV_ref.data()[k];
1688 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1691 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1692 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1695 ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial);
1696 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
1698 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_ref.data()[k*nDOF_trial_element],p);
1699 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
u);
1700 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
v);
1701 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
w);
1703 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial,grad_p);
1704 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
1705 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
1706 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_w);
1708 for (
int j=0;j<nDOF_trial_element;j++)
1710 p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV;
1711 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
1712 for (
int I=0;I<nSpace;I++)
1714 p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;
1715 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
1719 porosity = q_porosity.data()[eN_k];
1723 q_velocity.data()[eN_k_nSpace+0]=
u;
1724 q_velocity.data()[eN_k_nSpace+1]=
v;
1725 q_velocity.data()[eN_k_nSpace+2]=
w;
1726 q_x.data()[eN_k_nSpace+0]=x;
1727 q_x.data()[eN_k_nSpace+1]=y;
1728 q_x.data()[eN_k_nSpace+2]=
z;
1739 elementDiameter.data()[eN],
1740 smagorinskyConstant,
1741 turbulenceClosureModel,
1746 &normal_phi.data()[eN_k_nSpace],
1747 kappa_phi.data()[eN_k],
1799 q_dragBeam1.data()[eN_k],
1800 q_dragBeam2.data()[eN_k],
1801 q_dragBeam3.data()[eN_k]);
1803 mass_source = q_mass_source.data()[eN_k];
1809 q_dragAlpha.data()[eN_k],
1810 q_dragBeta.data()[eN_k],
1823 eps_solid.data()[0],
1824 phi_solid.data()[eN_k],
1825 q_velocity_solid.data()[eN_k_nSpace+0],
1826 q_velocity_solid.data()[eN_k_nSpace+1],
1827 q_velocity_solid.data()[eN_k_nSpace+2],
1836 if (turbulenceClosureModel >= 3)
1838 const double c_mu = 0.09;
1851 q_turb_var_0.data()[eN_k],
1852 q_turb_var_1.data()[eN_k],
1853 &q_turb_var_grad_0.data()[eN_k_nSpace],
1871 q_mom_u_acc.data()[eN_k] = mom_u_acc;
1872 q_mom_v_acc.data()[eN_k] = mom_v_acc;
1873 q_mom_w_acc.data()[eN_k] = mom_w_acc;
1875 q_mass_adv.data()[eN_k_nSpace+0] =
u;
1876 q_mass_adv.data()[eN_k_nSpace+1] =
v;
1877 q_mass_adv.data()[eN_k_nSpace+2] =
w;
1882 mass_adv[0] -= MOVING_DOMAIN*
xt;
1883 mass_adv[1] -= MOVING_DOMAIN*yt;
1884 mass_adv[2] -= MOVING_DOMAIN*zt;
1886 mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*
xt;
1887 mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt;
1888 mom_u_adv[2] -= MOVING_DOMAIN*mom_u_acc*zt;
1889 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
1890 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
1891 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
1893 mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*
xt;
1894 mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt;
1895 mom_v_adv[2] -= MOVING_DOMAIN*mom_v_acc*zt;
1896 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
1897 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
1898 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
1900 mom_w_adv[0] -= MOVING_DOMAIN*mom_w_acc*
xt;
1901 mom_w_adv[1] -= MOVING_DOMAIN*mom_w_acc*yt;
1902 mom_w_adv[2] -= MOVING_DOMAIN*mom_w_acc*zt;
1903 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
1904 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
1905 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
1910 q_mom_u_acc_beta_bdf.data()[eN_k],
1916 q_mom_v_acc_beta_bdf.data()[eN_k],
1922 q_mom_w_acc_beta_bdf.data()[eN_k],
1931 pdeResidual_p =
ck.Advection_strong(dmass_adv_u,grad_u) +
1933 ck.Reaction_strong(mass_source) +
1935 ck.Advection_strong(dmass_adv_v,grad_v) +
1936 ck.Advection_strong(dmass_adv_w,grad_w);
1938 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
1939 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
1940 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
1942 pdeResidual_u =
ck.Mass_strong(mom_u_acc_t) +
1943 ck.Advection_strong(dmom_adv_sge,grad_u) +
1944 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
1945 ck.Reaction_strong(mom_u_source);
1947 pdeResidual_v =
ck.Mass_strong(mom_v_acc_t) +
1948 ck.Advection_strong(dmom_adv_sge,grad_v) +
1949 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
1950 ck.Reaction_strong(mom_v_source);
1952 pdeResidual_w =
ck.Mass_strong(mom_w_acc_t) +
1953 ck.Advection_strong(dmom_adv_sge,grad_w) +
1954 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
1955 ck.Reaction_strong(mom_w_source);
1959 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
1961 elementDiameter.data()[eN],
1966 dmom_u_ham_grad_p[0],
1969 q_cfl.data()[eN_k]);
1976 dmom_u_ham_grad_p[0],
1979 q_cfl.data()[eN_k]);
1981 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
1982 tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0;
1995 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
1996 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
1997 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
1999 mom_u_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+0]);
2000 mom_u_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+0]);
2001 mom_u_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge.data()[eN_k_nSpace+0]);
2003 mom_v_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+1]);
2004 mom_v_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+1]);
2005 mom_v_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge.data()[eN_k_nSpace+1]);
2007 mom_w_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+2]);
2008 mom_w_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+2]);
2009 mom_w_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge.data()[eN_k_nSpace+2]);
2012 for (
int i=0;i<nDOF_test_element;i++)
2014 register int i_nSpace = i*nSpace;
2015 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
2016 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
2017 Lstar_w_p[i]=
ck.Advection_adjoint(dmass_adv_w,&p_grad_test_dV[i_nSpace]);
2019 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2020 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2021 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2022 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2023 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2024 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2027 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
2028 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
2029 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
2033 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v + pdeResidual_w*pdeResidual_w);
2034 q_numDiff_u.data()[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
2035 (1.0-useMetrics)*hFactor*hFactor*elementDiameter.data()[eN]*elementDiameter.data()[eN]);
2036 q_numDiff_v.data()[eN_k] = q_numDiff_u.data()[eN_k];
2037 q_numDiff_w.data()[eN_k] = q_numDiff_u.data()[eN_k];
2041 for(
int i=0;i<nDOF_test_element;i++)
2043 register int i_nSpace=i*nSpace;
2045 elementResidual_p[i] +=
ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) +
2047 ck.Reaction_weak(mass_source,p_test_dV[i]) +
2049 ck.SubgridError(subgridError_u,Lstar_u_p[i]) +
2050 ck.SubgridError(subgridError_v,Lstar_v_p[i]) +
2051 ck.SubgridError(subgridError_w,Lstar_w_p[i]);
2053 elementResidual_u[i] +=
ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
2054 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
2055 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]) +
2056 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]) +
2057 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]) +
2058 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
2059 ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
2060 ck.SubgridError(subgridError_p,Lstar_p_u[i]) +
2061 ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
2062 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]);
2064 elementResidual_v[i] +=
ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
2065 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
2066 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]) +
2067 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]) +
2068 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]) +
2069 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
2070 ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
2071 ck.SubgridError(subgridError_p,Lstar_p_v[i]) +
2072 ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
2073 ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]);
2075 elementResidual_w[i] +=
ck.Mass_weak(mom_w_acc_t,vel_test_dV[i]) +
2076 ck.Advection_weak(mom_w_adv,&vel_grad_test_dV[i_nSpace]) +
2077 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]) +
2078 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]) +
2079 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]) +
2080 ck.Reaction_weak(mom_w_source,vel_test_dV[i]) +
2081 ck.Hamiltonian_weak(mom_w_ham,vel_test_dV[i]) +
2082 ck.SubgridError(subgridError_p,Lstar_p_w[i]) +
2083 ck.SubgridError(subgridError_w,Lstar_w_w[i]) +
2084 ck.NumericalDiffusion(q_numDiff_w_last.data()[eN_k],grad_w,&vel_grad_test_dV[i_nSpace]);
2090 for(
int i=0;i<nDOF_test_element;i++)
2092 register int eN_i=eN*nDOF_test_element+i;
2094 elementResidual_p_save.data()[eN_i] += elementResidual_p[i];
2096 globalResidual.data()[offset_p+stride_p*p_l2g.data()[eN_i]]+=elementResidual_p[i];
2097 globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
2098 globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
2099 globalResidual.data()[offset_w+stride_w*vel_l2g.data()[eN_i]]+=elementResidual_w[i];
2108 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2110 register int ebN = exteriorElementBoundariesArray.data()[ebNE],
2111 eN = elementBoundaryElementsArray.data()[ebN*2+0],
2112 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
2113 eN_nDOF_trial_element = eN*nDOF_trial_element;
2114 register double elementResidual_p[nDOF_test_element],
2115 elementResidual_u[nDOF_test_element],
2116 elementResidual_v[nDOF_test_element],
2117 elementResidual_w[nDOF_test_element],
2119 for (
int i=0;i<nDOF_test_element;i++)
2121 elementResidual_p[i]=0.0;
2122 elementResidual_u[i]=0.0;
2123 elementResidual_v[i]=0.0;
2124 elementResidual_w[i]=0.0;
2126 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
2128 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
2129 ebNE_kb_nSpace = ebNE_kb*nSpace,
2130 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
2131 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
2132 register double p_ext=0.0,
2141 dmom_u_acc_u_ext=0.0,
2143 dmom_v_acc_v_ext=0.0,
2145 dmom_w_acc_w_ext=0.0,
2146 mass_adv_ext[nSpace],
2147 dmass_adv_u_ext[nSpace],
2148 dmass_adv_v_ext[nSpace],
2149 dmass_adv_w_ext[nSpace],
2150 mom_u_adv_ext[nSpace],
2151 dmom_u_adv_u_ext[nSpace],
2152 dmom_u_adv_v_ext[nSpace],
2153 dmom_u_adv_w_ext[nSpace],
2154 mom_v_adv_ext[nSpace],
2155 dmom_v_adv_u_ext[nSpace],
2156 dmom_v_adv_v_ext[nSpace],
2157 dmom_v_adv_w_ext[nSpace],
2158 mom_w_adv_ext[nSpace],
2159 dmom_w_adv_u_ext[nSpace],
2160 dmom_w_adv_v_ext[nSpace],
2161 dmom_w_adv_w_ext[nSpace],
2162 mom_uu_diff_ten_ext[nSpace],
2163 mom_vv_diff_ten_ext[nSpace],
2164 mom_ww_diff_ten_ext[nSpace],
2165 mom_uv_diff_ten_ext[1],
2166 mom_uw_diff_ten_ext[1],
2167 mom_vu_diff_ten_ext[1],
2168 mom_vw_diff_ten_ext[1],
2169 mom_wu_diff_ten_ext[1],
2170 mom_wv_diff_ten_ext[1],
2171 mom_u_source_ext=0.0,
2172 mom_v_source_ext=0.0,
2173 mom_w_source_ext=0.0,
2175 dmom_u_ham_grad_p_ext[nSpace],
2177 dmom_v_ham_grad_p_ext[nSpace],
2179 dmom_w_ham_grad_p_ext[nSpace],
2180 dmom_u_adv_p_ext[nSpace],
2181 dmom_v_adv_p_ext[nSpace],
2182 dmom_w_adv_p_ext[nSpace],
2184 flux_mom_u_adv_ext=0.0,
2185 flux_mom_v_adv_ext=0.0,
2186 flux_mom_w_adv_ext=0.0,
2187 flux_mom_uu_diff_ext=0.0,
2188 flux_mom_uv_diff_ext=0.0,
2189 flux_mom_uw_diff_ext=0.0,
2190 flux_mom_vu_diff_ext=0.0,
2191 flux_mom_vv_diff_ext=0.0,
2192 flux_mom_vw_diff_ext=0.0,
2193 flux_mom_wu_diff_ext=0.0,
2194 flux_mom_wv_diff_ext=0.0,
2195 flux_mom_ww_diff_ext=0.0,
2200 bc_mom_u_acc_ext=0.0,
2201 bc_dmom_u_acc_u_ext=0.0,
2202 bc_mom_v_acc_ext=0.0,
2203 bc_dmom_v_acc_v_ext=0.0,
2204 bc_mom_w_acc_ext=0.0,
2205 bc_dmom_w_acc_w_ext=0.0,
2206 bc_mass_adv_ext[nSpace],
2207 bc_dmass_adv_u_ext[nSpace],
2208 bc_dmass_adv_v_ext[nSpace],
2209 bc_dmass_adv_w_ext[nSpace],
2210 bc_mom_u_adv_ext[nSpace],
2211 bc_dmom_u_adv_u_ext[nSpace],
2212 bc_dmom_u_adv_v_ext[nSpace],
2213 bc_dmom_u_adv_w_ext[nSpace],
2214 bc_mom_v_adv_ext[nSpace],
2215 bc_dmom_v_adv_u_ext[nSpace],
2216 bc_dmom_v_adv_v_ext[nSpace],
2217 bc_dmom_v_adv_w_ext[nSpace],
2218 bc_mom_w_adv_ext[nSpace],
2219 bc_dmom_w_adv_u_ext[nSpace],
2220 bc_dmom_w_adv_v_ext[nSpace],
2221 bc_dmom_w_adv_w_ext[nSpace],
2222 bc_mom_uu_diff_ten_ext[nSpace],
2223 bc_mom_vv_diff_ten_ext[nSpace],
2224 bc_mom_ww_diff_ten_ext[nSpace],
2225 bc_mom_uv_diff_ten_ext[1],
2226 bc_mom_uw_diff_ten_ext[1],
2227 bc_mom_vu_diff_ten_ext[1],
2228 bc_mom_vw_diff_ten_ext[1],
2229 bc_mom_wu_diff_ten_ext[1],
2230 bc_mom_wv_diff_ten_ext[1],
2231 bc_mom_u_source_ext=0.0,
2232 bc_mom_v_source_ext=0.0,
2233 bc_mom_w_source_ext=0.0,
2234 bc_mom_u_ham_ext=0.0,
2235 bc_dmom_u_ham_grad_p_ext[nSpace],
2236 bc_mom_v_ham_ext=0.0,
2237 bc_dmom_v_ham_grad_p_ext[nSpace],
2238 bc_mom_w_ham_ext=0.0,
2239 bc_dmom_w_ham_grad_p_ext[nSpace],
2240 jac_ext[nSpace*nSpace],
2242 jacInv_ext[nSpace*nSpace],
2243 boundaryJac[nSpace*(nSpace-1)],
2244 metricTensor[(nSpace-1)*(nSpace-1)],
2245 metricTensorDetSqrt,
2246 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
2247 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
2248 vel_grad_test_dS[nDOF_trial_element*nSpace],
2249 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
2253 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
2254 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;
2256 ck.calculateMapping_elementBoundary(eN,
2262 mesh_trial_trace_ref.data(),
2263 mesh_grad_trial_trace_ref.data(),
2264 boundaryJac_ref.data(),
2270 metricTensorDetSqrt,
2274 ck.calculateMappingVelocity_elementBoundary(eN,
2278 mesh_velocity_dof.data(),
2280 mesh_trial_trace_ref.data(),
2281 xt_ext,yt_ext,zt_ext,
2291 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
2294 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
2295 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
2297 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2298 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2302 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
2303 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
2306 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);
2307 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);
2308 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);
2309 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);
2310 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
2311 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
2312 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
2313 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
2315 for (
int j=0;j<nDOF_trial_element;j++)
2317 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2318 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2319 for (
int I=0;I<nSpace;I++)
2320 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
2322 bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
2323 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
2324 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
2325 bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext;
2327 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
2338 elementDiameter.data()[eN],
2339 smagorinskyConstant,
2340 turbulenceClosureModel,
2343 ebqe_vf_ext.data()[ebNE_kb],
2344 ebqe_phi_ext.data()[ebNE_kb],
2345 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
2346 ebqe_kappa_phi_ext.data()[ebNE_kb],
2380 mom_uu_diff_ten_ext,
2381 mom_vv_diff_ten_ext,
2382 mom_ww_diff_ten_ext,
2383 mom_uv_diff_ten_ext,
2384 mom_uw_diff_ten_ext,
2385 mom_vu_diff_ten_ext,
2386 mom_vw_diff_ten_ext,
2387 mom_wu_diff_ten_ext,
2388 mom_wv_diff_ten_ext,
2393 dmom_u_ham_grad_p_ext,
2395 dmom_v_ham_grad_p_ext,
2397 dmom_w_ham_grad_p_ext,
2398 ebqe_dragBeam1.data()[ebNE_kb],
2399 ebqe_dragBeam2.data()[ebNE_kb],
2400 ebqe_dragBeam3.data()[ebNE_kb]);
2408 elementDiameter.data()[eN],
2409 smagorinskyConstant,
2410 turbulenceClosureModel,
2413 bc_ebqe_vf_ext.data()[ebNE_kb],
2414 bc_ebqe_phi_ext.data()[ebNE_kb],
2415 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
2416 ebqe_kappa_phi_ext.data()[ebNE_kb],
2429 bc_dmom_u_acc_u_ext,
2431 bc_dmom_v_acc_v_ext,
2433 bc_dmom_w_acc_w_ext,
2439 bc_dmom_u_adv_u_ext,
2440 bc_dmom_u_adv_v_ext,
2441 bc_dmom_u_adv_w_ext,
2443 bc_dmom_v_adv_u_ext,
2444 bc_dmom_v_adv_v_ext,
2445 bc_dmom_v_adv_w_ext,
2447 bc_dmom_w_adv_u_ext,
2448 bc_dmom_w_adv_v_ext,
2449 bc_dmom_w_adv_w_ext,
2450 bc_mom_uu_diff_ten_ext,
2451 bc_mom_vv_diff_ten_ext,
2452 bc_mom_ww_diff_ten_ext,
2453 bc_mom_uv_diff_ten_ext,
2454 bc_mom_uw_diff_ten_ext,
2455 bc_mom_vu_diff_ten_ext,
2456 bc_mom_vw_diff_ten_ext,
2457 bc_mom_wu_diff_ten_ext,
2458 bc_mom_wv_diff_ten_ext,
2459 bc_mom_u_source_ext,
2460 bc_mom_v_source_ext,
2461 bc_mom_w_source_ext,
2463 bc_dmom_u_ham_grad_p_ext,
2465 bc_dmom_v_ham_grad_p_ext,
2467 bc_dmom_w_ham_grad_p_ext,
2468 ebqe_dragBeam1.data()[ebNE_kb],
2469 ebqe_dragBeam2.data()[ebNE_kb],
2470 ebqe_dragBeam3.data()[ebNE_kb]);
2473 if (turbulenceClosureModel >= 3)
2475 const double turb_var_grad_0_dummy[3] = {0.,0.,0.};
2476 const double c_mu = 0.09;
2485 ebqe_vf_ext.data()[ebNE_kb],
2486 ebqe_phi_ext.data()[ebNE_kb],
2489 ebqe_turb_var_0.data()[ebNE_kb],
2490 ebqe_turb_var_1.data()[ebNE_kb],
2491 turb_var_grad_0_dummy,
2492 mom_uu_diff_ten_ext,
2493 mom_vv_diff_ten_ext,
2494 mom_ww_diff_ten_ext,
2495 mom_uv_diff_ten_ext,
2496 mom_uw_diff_ten_ext,
2497 mom_vu_diff_ten_ext,
2498 mom_vw_diff_ten_ext,
2499 mom_wu_diff_ten_ext,
2500 mom_wv_diff_ten_ext,
2513 bc_ebqe_vf_ext.data()[ebNE_kb],
2514 bc_ebqe_phi_ext.data()[ebNE_kb],
2517 ebqe_turb_var_0.data()[ebNE_kb],
2518 ebqe_turb_var_1.data()[ebNE_kb],
2519 turb_var_grad_0_dummy,
2520 bc_mom_uu_diff_ten_ext,
2521 bc_mom_vv_diff_ten_ext,
2522 bc_mom_ww_diff_ten_ext,
2523 bc_mom_uv_diff_ten_ext,
2524 bc_mom_uw_diff_ten_ext,
2525 bc_mom_vu_diff_ten_ext,
2526 bc_mom_vw_diff_ten_ext,
2527 bc_mom_wu_diff_ten_ext,
2528 bc_mom_wv_diff_ten_ext,
2529 bc_mom_u_source_ext,
2530 bc_mom_v_source_ext,
2531 bc_mom_w_source_ext);
2538 mass_adv_ext[0] -= MOVING_DOMAIN*xt_ext;
2539 mass_adv_ext[1] -= MOVING_DOMAIN*yt_ext;
2540 mass_adv_ext[2] -= MOVING_DOMAIN*zt_ext;
2542 mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext;
2543 mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext;
2544 mom_u_adv_ext[2] -= MOVING_DOMAIN*mom_u_acc_ext*zt_ext;
2545 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
2546 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
2547 dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext;
2549 mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext;
2550 mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext;
2551 mom_v_adv_ext[2] -= MOVING_DOMAIN*mom_v_acc_ext*zt_ext;
2552 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
2553 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
2554 dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext;
2556 mom_w_adv_ext[0] -= MOVING_DOMAIN*mom_w_acc_ext*xt_ext;
2557 mom_w_adv_ext[1] -= MOVING_DOMAIN*mom_w_acc_ext*yt_ext;
2558 mom_w_adv_ext[2] -= MOVING_DOMAIN*mom_w_acc_ext*zt_ext;
2559 dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext;
2560 dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext;
2561 dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext;
2564 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext;
2565 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext;
2566 bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*bc_mom_u_acc_ext*zt_ext;
2568 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext;
2569 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext;
2570 bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*bc_mom_v_acc_ext*zt_ext;
2572 bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*bc_mom_w_acc_ext*xt_ext;
2573 bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*bc_mom_w_acc_ext*yt_ext;
2574 bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*bc_mom_w_acc_ext*zt_ext;
2578 ck.calculateGScale(G,normal,h_penalty);
2579 penalty = useMetrics*C_b*h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
2581 isDOFBoundary_u.data()[ebNE_kb],
2582 isDOFBoundary_v.data()[ebNE_kb],
2583 isDOFBoundary_w.data()[ebNE_kb],
2584 isAdvectiveFluxBoundary_p.data()[ebNE_kb],
2585 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
2586 isAdvectiveFluxBoundary_v.data()[ebNE_kb],
2587 isAdvectiveFluxBoundary_w.data()[ebNE_kb],
2588 dmom_u_ham_grad_p_ext[0],
2589 bc_dmom_u_ham_grad_p_ext[0],
2596 ebqe_bc_flux_mass_ext.data()[ebNE_kb],
2597 ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
2598 ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
2599 ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
2624 &ebqe_velocity.data()[ebNE_kb_nSpace]);
2626 ebqe_phi_ext.data()[ebNE_kb],
2627 sdInfo_u_u_rowptr.data(),
2628 sdInfo_u_u_colind.data(),
2629 isDOFBoundary_u.data()[ebNE_kb],
2630 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2632 bc_mom_uu_diff_ten_ext,
2634 ebqe_bc_flux_u_diff_ext.data()[ebNE_kb],
2635 mom_uu_diff_ten_ext,
2639 flux_mom_uu_diff_ext);
2641 ebqe_phi_ext.data()[ebNE_kb],
2642 sdInfo_u_v_rowptr.data(),
2643 sdInfo_u_v_colind.data(),
2644 isDOFBoundary_v.data()[ebNE_kb],
2645 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2647 bc_mom_uv_diff_ten_ext,
2650 mom_uv_diff_ten_ext,
2654 flux_mom_uv_diff_ext);
2656 ebqe_phi_ext.data()[ebNE_kb],
2657 sdInfo_u_w_rowptr.data(),
2658 sdInfo_u_w_colind.data(),
2659 isDOFBoundary_w.data()[ebNE_kb],
2660 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2662 bc_mom_uw_diff_ten_ext,
2665 mom_uw_diff_ten_ext,
2669 flux_mom_uw_diff_ext);
2671 ebqe_phi_ext.data()[ebNE_kb],
2672 sdInfo_v_u_rowptr.data(),
2673 sdInfo_v_u_colind.data(),
2674 isDOFBoundary_u.data()[ebNE_kb],
2675 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2677 bc_mom_vu_diff_ten_ext,
2680 mom_vu_diff_ten_ext,
2684 flux_mom_vu_diff_ext);
2686 ebqe_phi_ext.data()[ebNE_kb],
2687 sdInfo_v_v_rowptr.data(),
2688 sdInfo_v_v_colind.data(),
2689 isDOFBoundary_v.data()[ebNE_kb],
2690 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2692 bc_mom_vv_diff_ten_ext,
2694 ebqe_bc_flux_v_diff_ext.data()[ebNE_kb],
2695 mom_vv_diff_ten_ext,
2699 flux_mom_vv_diff_ext);
2701 ebqe_phi_ext.data()[ebNE_kb],
2702 sdInfo_v_w_rowptr.data(),
2703 sdInfo_v_w_colind.data(),
2704 isDOFBoundary_w.data()[ebNE_kb],
2705 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2707 bc_mom_vw_diff_ten_ext,
2710 mom_vw_diff_ten_ext,
2714 flux_mom_vw_diff_ext);
2716 ebqe_phi_ext.data()[ebNE_kb],
2717 sdInfo_w_u_rowptr.data(),
2718 sdInfo_w_u_colind.data(),
2719 isDOFBoundary_u.data()[ebNE_kb],
2720 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2722 bc_mom_wu_diff_ten_ext,
2725 mom_wu_diff_ten_ext,
2729 flux_mom_wu_diff_ext);
2731 ebqe_phi_ext.data()[ebNE_kb],
2732 sdInfo_w_v_rowptr.data(),
2733 sdInfo_w_v_colind.data(),
2734 isDOFBoundary_v.data()[ebNE_kb],
2735 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2737 bc_mom_wv_diff_ten_ext,
2740 mom_wv_diff_ten_ext,
2744 flux_mom_wv_diff_ext);
2746 ebqe_phi_ext.data()[ebNE_kb],
2747 sdInfo_w_w_rowptr.data(),
2748 sdInfo_w_w_colind.data(),
2749 isDOFBoundary_w.data()[ebNE_kb],
2750 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2752 bc_mom_ww_diff_ten_ext,
2754 ebqe_bc_flux_w_diff_ext.data()[ebNE_kb],
2755 mom_ww_diff_ten_ext,
2759 flux_mom_ww_diff_ext);
2760 flux.data()[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
2764 if (ebN < nElementBoundaries_owned)
2766 force_v_x = flux_mom_uu_diff_ext + flux_mom_uv_diff_ext + flux_mom_uw_diff_ext;
2767 force_v_y = flux_mom_vu_diff_ext + flux_mom_vv_diff_ext + flux_mom_vw_diff_ext;
2768 force_v_z = flux_mom_wu_diff_ext + flux_mom_wv_diff_ext + flux_mom_ww_diff_ext;
2770 force_p_x = p_ext*normal[0];
2771 force_p_y = p_ext*normal[1];
2772 force_p_z = p_ext*normal[2];
2774 force_x = force_p_x + force_v_x;
2775 force_y = force_p_y + force_v_y;
2776 force_z = force_p_z + force_v_z;
2778 r_x = x_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+0];
2779 r_y = y_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+1];
2780 r_z = z_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+2];
2782 wettedAreas.data()[boundaryFlags.data()[ebN]] += dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2784 netForces_p.data()[3*boundaryFlags.data()[ebN]+0] += force_p_x*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2785 netForces_p.data()[3*boundaryFlags.data()[ebN]+1] += force_p_y*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2786 netForces_p.data()[3*boundaryFlags.data()[ebN]+2] += force_p_z*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2788 netForces_v.data()[3*boundaryFlags.data()[ebN]+0] += force_v_x*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2789 netForces_v.data()[3*boundaryFlags.data()[ebN]+1] += force_v_y*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2790 netForces_v.data()[3*boundaryFlags.data()[ebN]+2] += force_v_z*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2792 netMoments.data()[3*boundaryFlags.data()[ebN]+0] += (r_y*force_z - r_z*force_y)*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2793 netMoments.data()[3*boundaryFlags.data()[ebN]+1] += (r_z*force_x - r_x*force_z)*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2794 netMoments.data()[3*boundaryFlags.data()[ebN]+2] += (r_x*force_y - r_y*force_x)*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2799 for (
int i=0;i<nDOF_test_element;i++)
2801 elementResidual_p[i] +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
2802 globalConservationError +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
2804 elementResidual_u[i] +=
ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
2805 ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
2806 ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
2807 ck.ExteriorElementBoundaryFlux(flux_mom_uw_diff_ext,vel_test_dS[i])+
2808 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2809 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2814 sdInfo_u_u_rowptr.data(),
2815 sdInfo_u_u_colind.data(),
2816 mom_uu_diff_ten_ext,
2817 &vel_grad_test_dS[i*nSpace])+
2818 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2819 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2824 sdInfo_u_v_rowptr.data(),
2825 sdInfo_u_v_colind.data(),
2826 mom_uv_diff_ten_ext,
2827 &vel_grad_test_dS[i*nSpace])+
2828 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb],
2829 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2834 sdInfo_u_w_rowptr.data(),
2835 sdInfo_u_w_colind.data(),
2836 mom_uw_diff_ten_ext,
2837 &vel_grad_test_dS[i*nSpace]);
2838 elementResidual_v[i] +=
ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
2839 ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
2840 ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
2841 ck.ExteriorElementBoundaryFlux(flux_mom_vw_diff_ext,vel_test_dS[i])+
2842 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2843 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2848 sdInfo_v_u_rowptr.data(),
2849 sdInfo_v_u_colind.data(),
2850 mom_vu_diff_ten_ext,
2851 &vel_grad_test_dS[i*nSpace])+
2852 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2853 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2858 sdInfo_v_v_rowptr.data(),
2859 sdInfo_v_v_colind.data(),
2860 mom_vv_diff_ten_ext,
2861 &vel_grad_test_dS[i*nSpace])+
2862 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb],
2863 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2868 sdInfo_v_w_rowptr.data(),
2869 sdInfo_v_w_colind.data(),
2870 mom_vw_diff_ten_ext,
2871 &vel_grad_test_dS[i*nSpace]);
2873 elementResidual_w[i] +=
ck.ExteriorElementBoundaryFlux(flux_mom_w_adv_ext,vel_test_dS[i]) +
2874 ck.ExteriorElementBoundaryFlux(flux_mom_wu_diff_ext,vel_test_dS[i])+
2875 ck.ExteriorElementBoundaryFlux(flux_mom_wv_diff_ext,vel_test_dS[i])+
2876 ck.ExteriorElementBoundaryFlux(flux_mom_ww_diff_ext,vel_test_dS[i])+
2877 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2878 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2883 sdInfo_w_u_rowptr.data(),
2884 sdInfo_w_u_colind.data(),
2885 mom_wu_diff_ten_ext,
2886 &vel_grad_test_dS[i*nSpace])+
2887 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2888 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2893 sdInfo_w_v_rowptr.data(),
2894 sdInfo_w_v_colind.data(),
2895 mom_wv_diff_ten_ext,
2896 &vel_grad_test_dS[i*nSpace])+
2897 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb],
2898 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2903 sdInfo_w_w_rowptr.data(),
2904 sdInfo_w_w_colind.data(),
2905 mom_ww_diff_ten_ext,
2906 &vel_grad_test_dS[i*nSpace]);
2912 for (
int i=0;i<nDOF_test_element;i++)
2914 int eN_i = eN*nDOF_test_element+i;
2916 elementResidual_p_save.data()[eN_i] += elementResidual_p[i];
2918 globalResidual.data()[offset_p+stride_p*p_l2g.data()[eN_i]]+=elementResidual_p[i];
2919 globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
2920 globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
2921 globalResidual.data()[offset_w+stride_w*vel_l2g.data()[eN_i]]+=elementResidual_w[i];
2928 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
2929 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
2930 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
2931 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
2932 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
2933 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
2934 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
2935 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
2936 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
2937 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
2938 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
2939 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
2940 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
2941 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
2942 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
2943 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
2944 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
2945 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
2946 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
2947 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
2948 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
2949 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
2950 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
2951 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
2952 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
2953 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
2954 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
2955 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
2956 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
2957 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
2958 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
2959 double hFactor = args.
scalar<
double>(
"hFactor");
2960 int nElements_global = args.
scalar<
int>(
"nElements_global");
2961 double useRBLES = args.
scalar<
double>(
"useRBLES");
2962 double useMetrics = args.
scalar<
double>(
"useMetrics");
2963 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
2964 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
2965 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
2966 double sigma = args.
scalar<
double>(
"sigma");
2971 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
2972 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
2973 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
2974 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
2975 double C_dg = args.
scalar<
double>(
"C_dg");
2976 double C_b = args.
scalar<
double>(
"C_b");
2977 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
2978 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
2979 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
2980 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
2981 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
2982 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
2983 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
2984 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
2985 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
2986 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
2987 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
2988 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
2989 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
2990 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
2991 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
2992 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
2993 xt::pyarray<double>& g = args.
array<
double>(
"g");
2994 const double useVF = args.
scalar<
double>(
"useVF");
2995 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
2996 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
2997 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
2998 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
2999 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
3000 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
3001 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
3002 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
3003 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
3004 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
3005 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
3006 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
3007 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
3008 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
3009 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
3010 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
3011 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
3012 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
3013 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
3014 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
3015 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
3016 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
3017 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
3018 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
3019 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
3020 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
3021 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
3022 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
3023 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
3024 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
3025 xt::pyarray<int>& csrRowIndeces_p_p = args.
array<
int>(
"csrRowIndeces_p_p");
3026 xt::pyarray<int>& csrColumnOffsets_p_p = args.
array<
int>(
"csrColumnOffsets_p_p");
3027 xt::pyarray<int>& csrRowIndeces_p_u = args.
array<
int>(
"csrRowIndeces_p_u");
3028 xt::pyarray<int>& csrColumnOffsets_p_u = args.
array<
int>(
"csrColumnOffsets_p_u");
3029 xt::pyarray<int>& csrRowIndeces_p_v = args.
array<
int>(
"csrRowIndeces_p_v");
3030 xt::pyarray<int>& csrColumnOffsets_p_v = args.
array<
int>(
"csrColumnOffsets_p_v");
3031 xt::pyarray<int>& csrRowIndeces_p_w = args.
array<
int>(
"csrRowIndeces_p_w");
3032 xt::pyarray<int>& csrColumnOffsets_p_w = args.
array<
int>(
"csrColumnOffsets_p_w");
3033 xt::pyarray<int>& csrRowIndeces_u_p = args.
array<
int>(
"csrRowIndeces_u_p");
3034 xt::pyarray<int>& csrColumnOffsets_u_p = args.
array<
int>(
"csrColumnOffsets_u_p");
3035 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
3036 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
3037 xt::pyarray<int>& csrRowIndeces_u_v = args.
array<
int>(
"csrRowIndeces_u_v");
3038 xt::pyarray<int>& csrColumnOffsets_u_v = args.
array<
int>(
"csrColumnOffsets_u_v");
3039 xt::pyarray<int>& csrRowIndeces_u_w = args.
array<
int>(
"csrRowIndeces_u_w");
3040 xt::pyarray<int>& csrColumnOffsets_u_w = args.
array<
int>(
"csrColumnOffsets_u_w");
3041 xt::pyarray<int>& csrRowIndeces_v_p = args.
array<
int>(
"csrRowIndeces_v_p");
3042 xt::pyarray<int>& csrColumnOffsets_v_p = args.
array<
int>(
"csrColumnOffsets_v_p");
3043 xt::pyarray<int>& csrRowIndeces_v_u = args.
array<
int>(
"csrRowIndeces_v_u");
3044 xt::pyarray<int>& csrColumnOffsets_v_u = args.
array<
int>(
"csrColumnOffsets_v_u");
3045 xt::pyarray<int>& csrRowIndeces_v_v = args.
array<
int>(
"csrRowIndeces_v_v");
3046 xt::pyarray<int>& csrColumnOffsets_v_v = args.
array<
int>(
"csrColumnOffsets_v_v");
3047 xt::pyarray<int>& csrRowIndeces_v_w = args.
array<
int>(
"csrRowIndeces_v_w");
3048 xt::pyarray<int>& csrColumnOffsets_v_w = args.
array<
int>(
"csrColumnOffsets_v_w");
3049 xt::pyarray<int>& csrRowIndeces_w_p = args.
array<
int>(
"csrRowIndeces_w_p");
3050 xt::pyarray<int>& csrColumnOffsets_w_p = args.
array<
int>(
"csrColumnOffsets_w_p");
3051 xt::pyarray<int>& csrRowIndeces_w_u = args.
array<
int>(
"csrRowIndeces_w_u");
3052 xt::pyarray<int>& csrColumnOffsets_w_u = args.
array<
int>(
"csrColumnOffsets_w_u");
3053 xt::pyarray<int>& csrRowIndeces_w_v = args.
array<
int>(
"csrRowIndeces_w_v");
3054 xt::pyarray<int>& csrColumnOffsets_w_v = args.
array<
int>(
"csrColumnOffsets_w_v");
3055 xt::pyarray<int>& csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
3056 xt::pyarray<int>& csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
3057 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
3058 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
3059 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
3060 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
3061 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
3062 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
3063 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
3064 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
3065 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
3066 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
3067 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
3068 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
3069 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
3070 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
3071 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
3072 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
3073 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
3074 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
3075 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
3076 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
3077 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
3078 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
3079 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
3080 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
3081 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
3082 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
3083 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
3084 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
3085 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
3086 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
3087 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
3088 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
3089 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
3090 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
3091 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
3092 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
3093 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
3094 xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.
array<
int>(
"csrColumnOffsets_eb_p_p");
3095 xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.
array<
int>(
"csrColumnOffsets_eb_p_u");
3096 xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.
array<
int>(
"csrColumnOffsets_eb_p_v");
3097 xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.
array<
int>(
"csrColumnOffsets_eb_p_w");
3098 xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.
array<
int>(
"csrColumnOffsets_eb_u_p");
3099 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
3100 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.
array<
int>(
"csrColumnOffsets_eb_u_v");
3101 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.
array<
int>(
"csrColumnOffsets_eb_u_w");
3102 xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.
array<
int>(
"csrColumnOffsets_eb_v_p");
3103 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.
array<
int>(
"csrColumnOffsets_eb_v_u");
3104 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.
array<
int>(
"csrColumnOffsets_eb_v_v");
3105 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.
array<
int>(
"csrColumnOffsets_eb_v_w");
3106 xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.
array<
int>(
"csrColumnOffsets_eb_w_p");
3107 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.
array<
int>(
"csrColumnOffsets_eb_w_u");
3108 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.
array<
int>(
"csrColumnOffsets_eb_w_v");
3109 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.
array<
int>(
"csrColumnOffsets_eb_w_w");
3110 xt::pyarray<double>& q_dragBeam1 = args.
array<
double>(
"q_dragBeam1");
3111 xt::pyarray<double>& q_dragBeam2 = args.
array<
double>(
"q_dragBeam2");
3112 xt::pyarray<double>& q_dragBeam3 = args.
array<
double>(
"q_dragBeam3");
3113 xt::pyarray<double>& ebqe_dragBeam1 = args.
array<
double>(
"ebqe_dragBeam1");
3114 xt::pyarray<double>& ebqe_dragBeam2 = args.
array<
double>(
"ebqe_dragBeam2");
3115 xt::pyarray<double>& ebqe_dragBeam3 = args.
array<
double>(
"ebqe_dragBeam3");
3119 for(
int eN=0;eN<nElements_global;eN++)
3121 register double eps_rho,eps_mu;
3123 register double elementJacobian_p_p[nDOF_test_element][nDOF_trial_element],
3124 elementJacobian_p_u[nDOF_test_element][nDOF_trial_element],
3125 elementJacobian_p_v[nDOF_test_element][nDOF_trial_element],
3126 elementJacobian_p_w[nDOF_test_element][nDOF_trial_element],
3127 elementJacobian_u_p[nDOF_test_element][nDOF_trial_element],
3128 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
3129 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
3130 elementJacobian_u_w[nDOF_test_element][nDOF_trial_element],
3131 elementJacobian_v_p[nDOF_test_element][nDOF_trial_element],
3132 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
3133 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element],
3134 elementJacobian_v_w[nDOF_test_element][nDOF_trial_element],
3135 elementJacobian_w_p[nDOF_test_element][nDOF_trial_element],
3136 elementJacobian_w_u[nDOF_test_element][nDOF_trial_element],
3137 elementJacobian_w_v[nDOF_test_element][nDOF_trial_element],
3138 elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
3139 for (
int i=0;i<nDOF_test_element;i++)
3140 for (
int j=0;j<nDOF_trial_element;j++)
3142 elementJacobian_p_p[i][j]=0.0;
3143 elementJacobian_p_u[i][j]=0.0;
3144 elementJacobian_p_v[i][j]=0.0;
3145 elementJacobian_p_w[i][j]=0.0;
3146 elementJacobian_u_p[i][j]=0.0;
3147 elementJacobian_u_u[i][j]=0.0;
3148 elementJacobian_u_v[i][j]=0.0;
3149 elementJacobian_u_w[i][j]=0.0;
3150 elementJacobian_v_p[i][j]=0.0;
3151 elementJacobian_v_u[i][j]=0.0;
3152 elementJacobian_v_v[i][j]=0.0;
3153 elementJacobian_v_w[i][j]=0.0;
3154 elementJacobian_w_p[i][j]=0.0;
3155 elementJacobian_w_u[i][j]=0.0;
3156 elementJacobian_w_v[i][j]=0.0;
3157 elementJacobian_w_w[i][j]=0.0;
3159 for (
int k=0;k<nQuadraturePoints_element;k++)
3161 int eN_k = eN*nQuadraturePoints_element+k,
3162 eN_k_nSpace = eN_k*nSpace,
3163 eN_nDOF_trial_element = eN*nDOF_trial_element;
3166 register double p=0.0,
u=0.0,
v=0.0,
w=0.0,
3167 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
3175 dmass_adv_u[nSpace],
3176 dmass_adv_v[nSpace],
3177 dmass_adv_w[nSpace],
3179 dmom_u_adv_u[nSpace],
3180 dmom_u_adv_v[nSpace],
3181 dmom_u_adv_w[nSpace],
3183 dmom_v_adv_u[nSpace],
3184 dmom_v_adv_v[nSpace],
3185 dmom_v_adv_w[nSpace],
3187 dmom_w_adv_u[nSpace],
3188 dmom_w_adv_v[nSpace],
3189 dmom_w_adv_w[nSpace],
3190 mom_uu_diff_ten[nSpace],
3191 mom_vv_diff_ten[nSpace],
3192 mom_ww_diff_ten[nSpace],
3203 dmom_u_ham_grad_p[nSpace],
3205 dmom_v_ham_grad_p[nSpace],
3207 dmom_w_ham_grad_p[nSpace],
3218 dpdeResidual_p_u[nDOF_trial_element],dpdeResidual_p_v[nDOF_trial_element],dpdeResidual_p_w[nDOF_trial_element],
3219 dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_trial_element],
3220 dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_trial_element],
3221 dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_trial_element],
3222 Lstar_u_p[nDOF_test_element],
3223 Lstar_v_p[nDOF_test_element],
3224 Lstar_w_p[nDOF_test_element],
3225 Lstar_u_u[nDOF_test_element],
3226 Lstar_v_v[nDOF_test_element],
3227 Lstar_w_w[nDOF_test_element],
3228 Lstar_p_u[nDOF_test_element],
3229 Lstar_p_v[nDOF_test_element],
3230 Lstar_p_w[nDOF_test_element],
3235 dsubgridError_p_u[nDOF_trial_element],
3236 dsubgridError_p_v[nDOF_trial_element],
3237 dsubgridError_p_w[nDOF_trial_element],
3238 dsubgridError_u_p[nDOF_trial_element],
3239 dsubgridError_u_u[nDOF_trial_element],
3240 dsubgridError_v_p[nDOF_trial_element],
3241 dsubgridError_v_v[nDOF_trial_element],
3242 dsubgridError_w_p[nDOF_trial_element],
3243 dsubgridError_w_w[nDOF_trial_element],
3244 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
3245 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
3248 jacInv[nSpace*nSpace],
3249 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
3251 p_test_dV[nDOF_test_element],vel_test_dV[nDOF_test_element],
3252 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
3257 dmom_u_source[nSpace],
3258 dmom_v_source[nSpace],
3259 dmom_w_source[nSpace],
3262 G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
3264 ck.calculateMapping_element(eN,
3268 mesh_trial_ref.data(),
3269 mesh_grad_trial_ref.data(),
3274 ck.calculateH_element(eN,
3276 nodeDiametersArray.data(),
3278 mesh_trial_ref.data(),
3280 ck.calculateMappingVelocity_element(eN,
3282 mesh_velocity_dof.data(),
3284 mesh_trial_ref.data(),
3289 dV = fabs(jacDet)*dV_ref.data()[k];
3290 ck.calculateG(jacInv,G,G_dd_G,tr_G);
3293 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3294 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3297 ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial);
3298 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
3300 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_ref.data()[k*nDOF_trial_element],p);
3301 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
u);
3302 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
v);
3303 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
w);
3305 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial,grad_p);
3306 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
3307 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
3308 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_w);
3310 for (
int j=0;j<nDOF_trial_element;j++)
3312 p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV;
3313 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
3314 for (
int I=0;I<nSpace;I++)
3316 p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;
3317 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
3322 porosity = q_porosity.data()[eN_k];
3334 elementDiameter.data()[eN],
3335 smagorinskyConstant,
3336 turbulenceClosureModel,
3341 &normal_phi.data()[eN_k_nSpace],
3342 kappa_phi.data()[eN_k],
3394 q_dragBeam1.data()[eN_k],
3395 q_dragBeam2.data()[eN_k],
3396 q_dragBeam3.data()[eN_k]);
3398 mass_source = q_mass_source.data()[eN_k];
3404 q_dragAlpha.data()[eN_k],
3405 q_dragBeta.data()[eN_k],
3418 eps_solid.data()[0],
3419 phi_solid.data()[eN_k],
3420 q_velocity_solid.data()[eN_k_nSpace+0],
3421 q_velocity_solid.data()[eN_k_nSpace+1],
3422 q_velocity_solid.data()[eN_k_nSpace+2],
3430 if (turbulenceClosureModel >= 3)
3432 const double c_mu = 0.09;
3445 q_turb_var_0.data()[eN_k],
3446 q_turb_var_1.data()[eN_k],
3447 &q_turb_var_grad_0.data()[eN_k_nSpace],
3466 mass_adv[0] -= MOVING_DOMAIN*
xt;
3467 mass_adv[1] -= MOVING_DOMAIN*yt;
3468 mass_adv[2] -= MOVING_DOMAIN*zt;
3470 mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*
xt;
3471 mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt;
3472 mom_u_adv[2] -= MOVING_DOMAIN*mom_u_acc*zt;
3473 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
3474 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
3475 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
3477 mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*
xt;
3478 mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt;
3479 mom_v_adv[2] -= MOVING_DOMAIN*mom_v_acc*zt;
3480 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
3481 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
3482 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
3484 mom_w_adv[0] -= MOVING_DOMAIN*mom_w_acc*
xt;
3485 mom_w_adv[1] -= MOVING_DOMAIN*mom_w_acc*yt;
3486 mom_w_adv[2] -= MOVING_DOMAIN*mom_w_acc*zt;
3487 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
3488 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
3489 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
3494 q_mom_u_acc_beta_bdf.data()[eN_k],
3500 q_mom_v_acc_beta_bdf.data()[eN_k],
3506 q_mom_w_acc_beta_bdf.data()[eN_k],
3514 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
3515 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
3516 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
3520 pdeResidual_p =
ck.Advection_strong(dmass_adv_u,grad_u) +
3521 ck.Advection_strong(dmass_adv_v,grad_v) +
3522 ck.Advection_strong(dmass_adv_w,grad_w);
3524 pdeResidual_u =
ck.Mass_strong(mom_u_acc_t) +
3525 ck.Advection_strong(dmom_adv_sge,grad_u) +
3526 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
3527 ck.Reaction_strong(mom_u_source);
3529 pdeResidual_v =
ck.Mass_strong(mom_v_acc_t) +
3530 ck.Advection_strong(dmom_adv_sge,grad_v) +
3531 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
3532 ck.Reaction_strong(mom_v_source);
3534 pdeResidual_w =
ck.Mass_strong(mom_w_acc_t) +
3535 ck.Advection_strong(dmom_adv_sge,grad_w) +
3536 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
3537 ck.Reaction_strong(mom_w_source);
3540 for (
int j=0;j<nDOF_trial_element;j++)
3542 register int j_nSpace = j*nSpace;
3543 dpdeResidual_p_u[j]=
ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial[j_nSpace]);
3544 dpdeResidual_p_v[j]=
ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial[j_nSpace]);
3545 dpdeResidual_p_w[j]=
ck.AdvectionJacobian_strong(dmass_adv_w,&vel_grad_trial[j_nSpace]);
3547 dpdeResidual_u_p[j]=
ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace]);
3548 dpdeResidual_u_u[j]=
ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3549 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]);
3551 dpdeResidual_v_p[j]=
ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace]);
3552 dpdeResidual_v_v[j]=
ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3553 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]);
3555 dpdeResidual_w_p[j]=
ck.HamiltonianJacobian_strong(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace]);
3556 dpdeResidual_w_w[j]=
ck.MassJacobian_strong(dmom_w_acc_w_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3557 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]);
3560 dpdeResidual_u_u[j]+=
ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3561 dpdeResidual_v_v[j]+=
ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3562 dpdeResidual_w_w[j]+=
ck.ReactionJacobian_strong(dmom_w_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3567 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
3569 elementDiameter.data()[eN],
3574 dmom_u_ham_grad_p[0],
3577 q_cfl.data()[eN_k]);
3584 dmom_u_ham_grad_p[0],
3587 q_cfl.data()[eN_k]);
3590 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
3591 tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0;
3625 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
3626 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
3627 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
3630 for (
int i=0;i<nDOF_test_element;i++)
3632 register int i_nSpace = i*nSpace;
3633 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
3634 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
3635 Lstar_w_p[i]=
ck.Advection_adjoint(dmass_adv_w,&p_grad_test_dV[i_nSpace]);
3636 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3637 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3638 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3639 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
3640 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
3641 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
3643 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
3644 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
3645 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
3649 dmom_u_adv_u[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3650 dmom_u_adv_u[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3651 dmom_u_adv_u[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
3653 dmom_v_adv_v[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3654 dmom_v_adv_v[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3655 dmom_v_adv_v[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
3657 dmom_w_adv_w[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3658 dmom_w_adv_w[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3659 dmom_w_adv_w[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
3663 for(
int i=0;i<nDOF_test_element;i++)
3665 register int i_nSpace = i*nSpace;
3666 for(
int j=0;j<nDOF_trial_element;j++)
3668 register int j_nSpace = j*nSpace;
3669 elementJacobian_p_p[i][j] +=
ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_p[i]) +
3670 ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_p[i]) +
3671 ck.SubgridErrorJacobian(dsubgridError_w_p[j],Lstar_w_p[i]);
3673 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]) +
3674 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_p[i]);
3675 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]) +
3676 ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_p[i]);
3677 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]) +
3678 ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_p[i]);
3680 elementJacobian_u_p[i][j] +=
ck.HamiltonianJacobian_weak(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) +
3681 ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_u[i]);
3682 elementJacobian_u_u[i][j] +=
ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3683 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3684 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]) +
3686 ck.ReactionJacobian_weak(dmom_u_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3688 ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_u[i]) +
3689 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
3690 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3691 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]) +
3692 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]) +
3694 ck.ReactionJacobian_weak(dmom_u_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3696 ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_u[i]);
3697 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]) +
3698 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]) +
3700 ck.ReactionJacobian_weak(dmom_u_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3702 ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_u[i]);
3704 elementJacobian_v_p[i][j] +=
ck.HamiltonianJacobian_weak(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) +
3705 ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_v[i]);
3706 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]) +
3707 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]) +
3709 ck.ReactionJacobian_weak(dmom_v_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3711 ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_v[i]);
3712 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]) +
3713 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3714 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]) +
3716 ck.ReactionJacobian_weak(dmom_v_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3718 ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_v[i]) +
3719 ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
3720 ck.NumericalDiffusionJacobian(q_numDiff_v_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3721 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]) +
3722 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]) +
3724 ck.ReactionJacobian_weak(dmom_v_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3726 ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_v[i]);
3728 elementJacobian_w_p[i][j] +=
ck.HamiltonianJacobian_weak(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) +
3729 ck.SubgridErrorJacobian(dsubgridError_w_p[j],Lstar_w_w[i]);
3730 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]) +
3731 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]) +
3733 ck.ReactionJacobian_weak(dmom_w_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3735 ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_w[i]);
3736 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]) +
3737 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]) +
3739 ck.ReactionJacobian_weak(dmom_w_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3741 ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_w[i]);
3742 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]) +
3743 ck.AdvectionJacobian_weak(dmom_w_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3744 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]) +
3746 ck.ReactionJacobian_weak(dmom_w_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3748 ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_w[i]) +
3749 ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_w[i]) +
3750 ck.NumericalDiffusionJacobian(q_numDiff_w_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3757 for (
int i=0;i<nDOF_test_element;i++)
3759 register int eN_i = eN*nDOF_test_element+i;
3760 for (
int j=0;j<nDOF_trial_element;j++)
3762 register int eN_i_j = eN_i*nDOF_trial_element+j;
3763 globalJacobian.data()[csrRowIndeces_p_p.data()[eN_i] + csrColumnOffsets_p_p.data()[eN_i_j]] += elementJacobian_p_p[i][j];
3764 globalJacobian.data()[csrRowIndeces_p_u.data()[eN_i] + csrColumnOffsets_p_u.data()[eN_i_j]] += elementJacobian_p_u[i][j];
3765 globalJacobian.data()[csrRowIndeces_p_v.data()[eN_i] + csrColumnOffsets_p_v.data()[eN_i_j]] += elementJacobian_p_v[i][j];
3766 globalJacobian.data()[csrRowIndeces_p_w.data()[eN_i] + csrColumnOffsets_p_w.data()[eN_i_j]] += elementJacobian_p_w[i][j];
3768 globalJacobian.data()[csrRowIndeces_u_p.data()[eN_i] + csrColumnOffsets_u_p.data()[eN_i_j]] += elementJacobian_u_p[i][j];
3769 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
3770 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_u_v.data()[eN_i_j]] += elementJacobian_u_v[i][j];
3771 globalJacobian.data()[csrRowIndeces_u_w.data()[eN_i] + csrColumnOffsets_u_w.data()[eN_i_j]] += elementJacobian_u_w[i][j];
3773 globalJacobian.data()[csrRowIndeces_v_p.data()[eN_i] + csrColumnOffsets_v_p.data()[eN_i_j]] += elementJacobian_v_p[i][j];
3774 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_v_u.data()[eN_i_j]] += elementJacobian_v_u[i][j];
3775 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_v_v.data()[eN_i_j]] += elementJacobian_v_v[i][j];
3776 globalJacobian.data()[csrRowIndeces_v_w.data()[eN_i] + csrColumnOffsets_v_w.data()[eN_i_j]] += elementJacobian_v_w[i][j];
3778 globalJacobian.data()[csrRowIndeces_w_p.data()[eN_i] + csrColumnOffsets_w_p.data()[eN_i_j]] += elementJacobian_w_p[i][j];
3779 globalJacobian.data()[csrRowIndeces_w_u.data()[eN_i] + csrColumnOffsets_w_u.data()[eN_i_j]] += elementJacobian_w_u[i][j];
3780 globalJacobian.data()[csrRowIndeces_w_v.data()[eN_i] + csrColumnOffsets_w_v.data()[eN_i_j]] += elementJacobian_w_v[i][j];
3781 globalJacobian.data()[csrRowIndeces_w_w.data()[eN_i] + csrColumnOffsets_w_w.data()[eN_i_j]] += elementJacobian_w_w[i][j];
3788 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3790 register int ebN = exteriorElementBoundariesArray.data()[ebNE],
3791 eN = elementBoundaryElementsArray.data()[ebN*2+0],
3792 eN_nDOF_trial_element = eN*nDOF_trial_element,
3793 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
3794 register double eps_rho,eps_mu;
3795 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3797 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3798 ebNE_kb_nSpace = ebNE_kb*nSpace,
3799 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3800 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3802 register double p_ext=0.0,
3811 dmom_u_acc_u_ext=0.0,
3813 dmom_v_acc_v_ext=0.0,
3815 dmom_w_acc_w_ext=0.0,
3816 mass_adv_ext[nSpace],
3817 dmass_adv_u_ext[nSpace],
3818 dmass_adv_v_ext[nSpace],
3819 dmass_adv_w_ext[nSpace],
3820 mom_u_adv_ext[nSpace],
3821 dmom_u_adv_u_ext[nSpace],
3822 dmom_u_adv_v_ext[nSpace],
3823 dmom_u_adv_w_ext[nSpace],
3824 mom_v_adv_ext[nSpace],
3825 dmom_v_adv_u_ext[nSpace],
3826 dmom_v_adv_v_ext[nSpace],
3827 dmom_v_adv_w_ext[nSpace],
3828 mom_w_adv_ext[nSpace],
3829 dmom_w_adv_u_ext[nSpace],
3830 dmom_w_adv_v_ext[nSpace],
3831 dmom_w_adv_w_ext[nSpace],
3832 mom_uu_diff_ten_ext[nSpace],
3833 mom_vv_diff_ten_ext[nSpace],
3834 mom_ww_diff_ten_ext[nSpace],
3835 mom_uv_diff_ten_ext[1],
3836 mom_uw_diff_ten_ext[1],
3837 mom_vu_diff_ten_ext[1],
3838 mom_vw_diff_ten_ext[1],
3839 mom_wu_diff_ten_ext[1],
3840 mom_wv_diff_ten_ext[1],
3841 mom_u_source_ext=0.0,
3842 mom_v_source_ext=0.0,
3843 mom_w_source_ext=0.0,
3845 dmom_u_ham_grad_p_ext[nSpace],
3847 dmom_v_ham_grad_p_ext[nSpace],
3849 dmom_w_ham_grad_p_ext[nSpace],
3850 dmom_u_adv_p_ext[nSpace],
3851 dmom_v_adv_p_ext[nSpace],
3852 dmom_w_adv_p_ext[nSpace],
3853 dflux_mass_u_ext=0.0,
3854 dflux_mass_v_ext=0.0,
3855 dflux_mass_w_ext=0.0,
3856 dflux_mom_u_adv_p_ext=0.0,
3857 dflux_mom_u_adv_u_ext=0.0,
3858 dflux_mom_u_adv_v_ext=0.0,
3859 dflux_mom_u_adv_w_ext=0.0,
3860 dflux_mom_v_adv_p_ext=0.0,
3861 dflux_mom_v_adv_u_ext=0.0,
3862 dflux_mom_v_adv_v_ext=0.0,
3863 dflux_mom_v_adv_w_ext=0.0,
3864 dflux_mom_w_adv_p_ext=0.0,
3865 dflux_mom_w_adv_u_ext=0.0,
3866 dflux_mom_w_adv_v_ext=0.0,
3867 dflux_mom_w_adv_w_ext=0.0,
3872 bc_mom_u_acc_ext=0.0,
3873 bc_dmom_u_acc_u_ext=0.0,
3874 bc_mom_v_acc_ext=0.0,
3875 bc_dmom_v_acc_v_ext=0.0,
3876 bc_mom_w_acc_ext=0.0,
3877 bc_dmom_w_acc_w_ext=0.0,
3878 bc_mass_adv_ext[nSpace],
3879 bc_dmass_adv_u_ext[nSpace],
3880 bc_dmass_adv_v_ext[nSpace],
3881 bc_dmass_adv_w_ext[nSpace],
3882 bc_mom_u_adv_ext[nSpace],
3883 bc_dmom_u_adv_u_ext[nSpace],
3884 bc_dmom_u_adv_v_ext[nSpace],
3885 bc_dmom_u_adv_w_ext[nSpace],
3886 bc_mom_v_adv_ext[nSpace],
3887 bc_dmom_v_adv_u_ext[nSpace],
3888 bc_dmom_v_adv_v_ext[nSpace],
3889 bc_dmom_v_adv_w_ext[nSpace],
3890 bc_mom_w_adv_ext[nSpace],
3891 bc_dmom_w_adv_u_ext[nSpace],
3892 bc_dmom_w_adv_v_ext[nSpace],
3893 bc_dmom_w_adv_w_ext[nSpace],
3894 bc_mom_uu_diff_ten_ext[nSpace],
3895 bc_mom_vv_diff_ten_ext[nSpace],
3896 bc_mom_ww_diff_ten_ext[nSpace],
3897 bc_mom_uv_diff_ten_ext[1],
3898 bc_mom_uw_diff_ten_ext[1],
3899 bc_mom_vu_diff_ten_ext[1],
3900 bc_mom_vw_diff_ten_ext[1],
3901 bc_mom_wu_diff_ten_ext[1],
3902 bc_mom_wv_diff_ten_ext[1],
3903 bc_mom_u_source_ext=0.0,
3904 bc_mom_v_source_ext=0.0,
3905 bc_mom_w_source_ext=0.0,
3906 bc_mom_u_ham_ext=0.0,
3907 bc_dmom_u_ham_grad_p_ext[nSpace],
3908 bc_mom_v_ham_ext=0.0,
3909 bc_dmom_v_ham_grad_p_ext[nSpace],
3910 bc_mom_w_ham_ext=0.0,
3911 bc_dmom_w_ham_grad_p_ext[nSpace],
3912 fluxJacobian_p_p[nDOF_trial_element],
3913 fluxJacobian_p_u[nDOF_trial_element],
3914 fluxJacobian_p_v[nDOF_trial_element],
3915 fluxJacobian_p_w[nDOF_trial_element],
3916 fluxJacobian_u_p[nDOF_trial_element],
3917 fluxJacobian_u_u[nDOF_trial_element],
3918 fluxJacobian_u_v[nDOF_trial_element],
3919 fluxJacobian_u_w[nDOF_trial_element],
3920 fluxJacobian_v_p[nDOF_trial_element],
3921 fluxJacobian_v_u[nDOF_trial_element],
3922 fluxJacobian_v_v[nDOF_trial_element],
3923 fluxJacobian_v_w[nDOF_trial_element],
3924 fluxJacobian_w_p[nDOF_trial_element],
3925 fluxJacobian_w_u[nDOF_trial_element],
3926 fluxJacobian_w_v[nDOF_trial_element],
3927 fluxJacobian_w_w[nDOF_trial_element],
3928 jac_ext[nSpace*nSpace],
3930 jacInv_ext[nSpace*nSpace],
3931 boundaryJac[nSpace*(nSpace-1)],
3932 metricTensor[(nSpace-1)*(nSpace-1)],
3933 metricTensorDetSqrt,
3934 p_grad_trial_trace[nDOF_trial_element*nSpace],
3935 vel_grad_trial_trace[nDOF_trial_element*nSpace],
3937 p_test_dS[nDOF_test_element],
3938 vel_test_dS[nDOF_test_element],
3940 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3941 vel_grad_test_dS[nDOF_trial_element*nSpace],
3945 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
3946 ck.calculateMapping_elementBoundary(eN,
3952 mesh_trial_trace_ref.data(),
3953 mesh_grad_trial_trace_ref.data(),
3954 boundaryJac_ref.data(),
3960 metricTensorDetSqrt,
3964 ck.calculateMappingVelocity_elementBoundary(eN,
3968 mesh_velocity_dof.data(),
3970 mesh_trial_trace_ref.data(),
3971 xt_ext,yt_ext,zt_ext,
3978 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
3979 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3980 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
3982 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3983 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3987 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
3988 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3990 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);
3991 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);
3992 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);
3993 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);
3994 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
3995 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3996 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3997 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
3999 for (
int j=0;j<nDOF_trial_element;j++)
4001 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4002 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4003 for (
int I=0;I<nSpace;I++)
4004 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
4009 bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
4010 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
4011 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
4012 bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext;
4014 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
4025 elementDiameter.data()[eN],
4026 smagorinskyConstant,
4027 turbulenceClosureModel,
4030 ebqe_vf_ext.data()[ebNE_kb],
4031 ebqe_phi_ext.data()[ebNE_kb],
4032 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
4033 ebqe_kappa_phi_ext.data()[ebNE_kb],
4067 mom_uu_diff_ten_ext,
4068 mom_vv_diff_ten_ext,
4069 mom_ww_diff_ten_ext,
4070 mom_uv_diff_ten_ext,
4071 mom_uw_diff_ten_ext,
4072 mom_vu_diff_ten_ext,
4073 mom_vw_diff_ten_ext,
4074 mom_wu_diff_ten_ext,
4075 mom_wv_diff_ten_ext,
4080 dmom_u_ham_grad_p_ext,
4082 dmom_v_ham_grad_p_ext,
4084 dmom_w_ham_grad_p_ext,
4085 ebqe_dragBeam1.data()[ebNE_kb],
4086 ebqe_dragBeam2.data()[ebNE_kb],
4087 ebqe_dragBeam3.data()[ebNE_kb]);
4095 elementDiameter.data()[eN],
4096 smagorinskyConstant,
4097 turbulenceClosureModel,
4100 bc_ebqe_vf_ext.data()[ebNE_kb],
4101 bc_ebqe_phi_ext.data()[ebNE_kb],
4102 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
4103 ebqe_kappa_phi_ext.data()[ebNE_kb],
4116 bc_dmom_u_acc_u_ext,
4118 bc_dmom_v_acc_v_ext,
4120 bc_dmom_w_acc_w_ext,
4126 bc_dmom_u_adv_u_ext,
4127 bc_dmom_u_adv_v_ext,
4128 bc_dmom_u_adv_w_ext,
4130 bc_dmom_v_adv_u_ext,
4131 bc_dmom_v_adv_v_ext,
4132 bc_dmom_v_adv_w_ext,
4134 bc_dmom_w_adv_u_ext,
4135 bc_dmom_w_adv_v_ext,
4136 bc_dmom_w_adv_w_ext,
4137 bc_mom_uu_diff_ten_ext,
4138 bc_mom_vv_diff_ten_ext,
4139 bc_mom_ww_diff_ten_ext,
4140 bc_mom_uv_diff_ten_ext,
4141 bc_mom_uw_diff_ten_ext,
4142 bc_mom_vu_diff_ten_ext,
4143 bc_mom_vw_diff_ten_ext,
4144 bc_mom_wu_diff_ten_ext,
4145 bc_mom_wv_diff_ten_ext,
4146 bc_mom_u_source_ext,
4147 bc_mom_v_source_ext,
4148 bc_mom_w_source_ext,
4150 bc_dmom_u_ham_grad_p_ext,
4152 bc_dmom_v_ham_grad_p_ext,
4154 bc_dmom_w_ham_grad_p_ext,
4155 ebqe_dragBeam1.data()[ebNE_kb],
4156 ebqe_dragBeam2.data()[ebNE_kb],
4157 ebqe_dragBeam3.data()[ebNE_kb]);
4159 if (turbulenceClosureModel >= 3)
4161 const double turb_var_grad_0_dummy[3] = {0.,0.,0.};
4162 const double c_mu = 0.09;
4171 ebqe_vf_ext.data()[ebNE_kb],
4172 ebqe_phi_ext.data()[ebNE_kb],
4175 ebqe_turb_var_0.data()[ebNE_kb],
4176 ebqe_turb_var_1.data()[ebNE_kb],
4177 turb_var_grad_0_dummy,
4178 mom_uu_diff_ten_ext,
4179 mom_vv_diff_ten_ext,
4180 mom_ww_diff_ten_ext,
4181 mom_uv_diff_ten_ext,
4182 mom_uw_diff_ten_ext,
4183 mom_vu_diff_ten_ext,
4184 mom_vw_diff_ten_ext,
4185 mom_wu_diff_ten_ext,
4186 mom_wv_diff_ten_ext,
4199 ebqe_vf_ext.data()[ebNE_kb],
4200 ebqe_phi_ext.data()[ebNE_kb],
4203 ebqe_turb_var_0.data()[ebNE_kb],
4204 ebqe_turb_var_1.data()[ebNE_kb],
4205 turb_var_grad_0_dummy,
4206 bc_mom_uu_diff_ten_ext,
4207 bc_mom_vv_diff_ten_ext,
4208 bc_mom_ww_diff_ten_ext,
4209 bc_mom_uv_diff_ten_ext,
4210 bc_mom_uw_diff_ten_ext,
4211 bc_mom_vu_diff_ten_ext,
4212 bc_mom_vw_diff_ten_ext,
4213 bc_mom_wu_diff_ten_ext,
4214 bc_mom_wv_diff_ten_ext,
4215 bc_mom_u_source_ext,
4216 bc_mom_v_source_ext,
4217 bc_mom_w_source_ext);
4222 mass_adv_ext[0] -= MOVING_DOMAIN*xt_ext;
4223 mass_adv_ext[1] -= MOVING_DOMAIN*yt_ext;
4224 mass_adv_ext[2] -= MOVING_DOMAIN*zt_ext;
4226 mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext;
4227 mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext;
4228 mom_u_adv_ext[2] -= MOVING_DOMAIN*mom_u_acc_ext*zt_ext;
4229 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
4230 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
4231 dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext;
4233 mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext;
4234 mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext;
4235 mom_v_adv_ext[2] -= MOVING_DOMAIN*mom_v_acc_ext*zt_ext;
4236 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
4237 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
4238 dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext;
4240 mom_w_adv_ext[0] -= MOVING_DOMAIN*mom_w_acc_ext*xt_ext;
4241 mom_w_adv_ext[1] -= MOVING_DOMAIN*mom_w_acc_ext*yt_ext;
4242 mom_w_adv_ext[2] -= MOVING_DOMAIN*mom_w_acc_ext*zt_ext;
4243 dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext;
4244 dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext;
4245 dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext;
4248 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext;
4249 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext;
4250 bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*bc_mom_u_acc_ext*zt_ext;
4252 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext;
4253 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext;
4254 bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*bc_mom_v_acc_ext*zt_ext;
4256 bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*bc_mom_w_acc_ext*xt_ext;
4257 bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*bc_mom_w_acc_ext*yt_ext;
4258 bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*bc_mom_w_acc_ext*zt_ext;
4263 isDOFBoundary_u.data()[ebNE_kb],
4264 isDOFBoundary_v.data()[ebNE_kb],
4265 isDOFBoundary_w.data()[ebNE_kb],
4266 isAdvectiveFluxBoundary_p.data()[ebNE_kb],
4267 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
4268 isAdvectiveFluxBoundary_v.data()[ebNE_kb],
4269 isAdvectiveFluxBoundary_w.data()[ebNE_kb],
4270 dmom_u_ham_grad_p_ext[0],
4277 ebqe_bc_flux_mass_ext.data()[ebNE_kb],
4278 ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
4279 ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
4280 ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
4304 dflux_mom_u_adv_p_ext,
4305 dflux_mom_u_adv_u_ext,
4306 dflux_mom_u_adv_v_ext,
4307 dflux_mom_u_adv_w_ext,
4308 dflux_mom_v_adv_p_ext,
4309 dflux_mom_v_adv_u_ext,
4310 dflux_mom_v_adv_v_ext,
4311 dflux_mom_v_adv_w_ext,
4312 dflux_mom_w_adv_p_ext,
4313 dflux_mom_w_adv_u_ext,
4314 dflux_mom_w_adv_v_ext,
4315 dflux_mom_w_adv_w_ext);
4319 ck.calculateGScale(G,normal,h_penalty);
4320 penalty = useMetrics*C_b*h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
4321 for (
int j=0;j<nDOF_trial_element;j++)
4323 register int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
4324 fluxJacobian_p_p[j]=0.0;
4325 fluxJacobian_p_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
4326 fluxJacobian_p_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
4327 fluxJacobian_p_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
4329 fluxJacobian_u_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
4330 fluxJacobian_u_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4332 ebqe_phi_ext.data()[ebNE_kb],
4333 sdInfo_u_u_rowptr.data(),
4334 sdInfo_u_u_colind.data(),
4335 isDOFBoundary_u.data()[ebNE_kb],
4336 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4338 mom_uu_diff_ten_ext,
4339 vel_trial_trace_ref.data()[ebN_local_kb_j],
4340 &vel_grad_trial_trace[j_nSpace],
4342 fluxJacobian_u_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4344 ebqe_phi_ext.data()[ebNE_kb],
4345 sdInfo_u_v_rowptr.data(),
4346 sdInfo_u_v_colind.data(),
4347 isDOFBoundary_v.data()[ebNE_kb],
4348 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4350 mom_uv_diff_ten_ext,
4351 vel_trial_trace_ref.data()[ebN_local_kb_j],
4352 &vel_grad_trial_trace[j_nSpace],
4354 fluxJacobian_u_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4356 ebqe_phi_ext.data()[ebNE_kb],
4357 sdInfo_u_w_rowptr.data(),
4358 sdInfo_u_w_colind.data(),
4359 isDOFBoundary_w.data()[ebNE_kb],
4360 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4362 mom_uw_diff_ten_ext,
4363 vel_trial_trace_ref.data()[ebN_local_kb_j],
4364 &vel_grad_trial_trace[j_nSpace],
4367 fluxJacobian_v_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
4368 fluxJacobian_v_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4370 ebqe_phi_ext.data()[ebNE_kb],
4371 sdInfo_v_u_rowptr.data(),
4372 sdInfo_v_u_colind.data(),
4373 isDOFBoundary_u.data()[ebNE_kb],
4374 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4376 mom_vu_diff_ten_ext,
4377 vel_trial_trace_ref.data()[ebN_local_kb_j],
4378 &vel_grad_trial_trace[j_nSpace],
4380 fluxJacobian_v_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4382 ebqe_phi_ext.data()[ebNE_kb],
4383 sdInfo_v_v_rowptr.data(),
4384 sdInfo_v_v_colind.data(),
4385 isDOFBoundary_v.data()[ebNE_kb],
4386 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4388 mom_vv_diff_ten_ext,
4389 vel_trial_trace_ref.data()[ebN_local_kb_j],
4390 &vel_grad_trial_trace[j_nSpace],
4392 fluxJacobian_v_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4394 ebqe_phi_ext.data()[ebNE_kb],
4395 sdInfo_v_w_rowptr.data(),
4396 sdInfo_v_w_colind.data(),
4397 isDOFBoundary_w.data()[ebNE_kb],
4398 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4400 mom_vw_diff_ten_ext,
4401 vel_trial_trace_ref.data()[ebN_local_kb_j],
4402 &vel_grad_trial_trace[j_nSpace],
4405 fluxJacobian_w_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
4406 fluxJacobian_w_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4408 ebqe_phi_ext.data()[ebNE_kb],
4409 sdInfo_w_u_rowptr.data(),
4410 sdInfo_w_u_colind.data(),
4411 isDOFBoundary_u.data()[ebNE_kb],
4412 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4414 mom_wu_diff_ten_ext,
4415 vel_trial_trace_ref.data()[ebN_local_kb_j],
4416 &vel_grad_trial_trace[j_nSpace],
4418 fluxJacobian_w_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4420 ebqe_phi_ext.data()[ebNE_kb],
4421 sdInfo_w_v_rowptr.data(),
4422 sdInfo_w_v_colind.data(),
4423 isDOFBoundary_v.data()[ebNE_kb],
4424 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4426 mom_wv_diff_ten_ext,
4427 vel_trial_trace_ref.data()[ebN_local_kb_j],
4428 &vel_grad_trial_trace[j_nSpace],
4430 fluxJacobian_w_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4432 ebqe_phi_ext.data()[ebNE_kb],
4433 sdInfo_w_w_rowptr.data(),
4434 sdInfo_w_w_colind.data(),
4435 isDOFBoundary_w.data()[ebNE_kb],
4436 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4438 mom_ww_diff_ten_ext,
4439 vel_trial_trace_ref.data()[ebN_local_kb_j],
4440 &vel_grad_trial_trace[j_nSpace],
4446 for (
int i=0;i<nDOF_test_element;i++)
4448 register int eN_i = eN*nDOF_test_element+i;
4449 for (
int j=0;j<nDOF_trial_element;j++)
4451 register int ebN_i_j = ebN*4*
nDOF_test_X_trial_element + i*nDOF_trial_element + j,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
4453 globalJacobian.data()[csrRowIndeces_p_p.data()[eN_i] + csrColumnOffsets_eb_p_p.data()[ebN_i_j]] += fluxJacobian_p_p[j]*p_test_dS[i];
4454 globalJacobian.data()[csrRowIndeces_p_u.data()[eN_i] + csrColumnOffsets_eb_p_u.data()[ebN_i_j]] += fluxJacobian_p_u[j]*p_test_dS[i];
4455 globalJacobian.data()[csrRowIndeces_p_v.data()[eN_i] + csrColumnOffsets_eb_p_v.data()[ebN_i_j]] += fluxJacobian_p_v[j]*p_test_dS[i];
4456 globalJacobian.data()[csrRowIndeces_p_w.data()[eN_i] + csrColumnOffsets_eb_p_w.data()[ebN_i_j]] += fluxJacobian_p_w[j]*p_test_dS[i];
4458 globalJacobian.data()[csrRowIndeces_u_p.data()[eN_i] + csrColumnOffsets_eb_u_p.data()[ebN_i_j]] += fluxJacobian_u_p[j]*vel_test_dS[i];
4459 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*vel_test_dS[i]+
4460 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4461 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4463 vel_trial_trace_ref.data()[ebN_local_kb_j],
4465 sdInfo_u_u_rowptr.data(),
4466 sdInfo_u_u_colind.data(),
4467 mom_uu_diff_ten_ext,
4468 &vel_grad_test_dS[i*nSpace]);
4469 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_eb_u_v.data()[ebN_i_j]] += fluxJacobian_u_v[j]*vel_test_dS[i]+
4470 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4471 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4473 vel_trial_trace_ref.data()[ebN_local_kb_j],
4475 sdInfo_u_v_rowptr.data(),
4476 sdInfo_u_v_colind.data(),
4477 mom_uv_diff_ten_ext,
4478 &vel_grad_test_dS[i*nSpace]);
4479 globalJacobian.data()[csrRowIndeces_u_w.data()[eN_i] + csrColumnOffsets_eb_u_w.data()[ebN_i_j]] += fluxJacobian_u_w[j]*vel_test_dS[i]+
4480 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb],
4481 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4483 vel_trial_trace_ref.data()[ebN_local_kb_j],
4485 sdInfo_u_w_rowptr.data(),
4486 sdInfo_u_w_colind.data(),
4487 mom_uw_diff_ten_ext,
4488 &vel_grad_test_dS[i*nSpace]);
4490 globalJacobian.data()[csrRowIndeces_v_p.data()[eN_i] + csrColumnOffsets_eb_v_p.data()[ebN_i_j]] += fluxJacobian_v_p[j]*vel_test_dS[i];
4491 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_eb_v_u.data()[ebN_i_j]] += fluxJacobian_v_u[j]*vel_test_dS[i]+
4492 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4493 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4495 vel_trial_trace_ref.data()[ebN_local_kb_j],
4497 sdInfo_v_u_rowptr.data(),
4498 sdInfo_v_u_colind.data(),
4499 mom_vu_diff_ten_ext,
4500 &vel_grad_test_dS[i*nSpace]);
4501 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_eb_v_v.data()[ebN_i_j]] += fluxJacobian_v_v[j]*vel_test_dS[i]+
4502 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4503 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4505 vel_trial_trace_ref.data()[ebN_local_kb_j],
4507 sdInfo_v_v_rowptr.data(),
4508 sdInfo_v_v_colind.data(),
4509 mom_vv_diff_ten_ext,
4510 &vel_grad_test_dS[i*nSpace]);
4511 globalJacobian.data()[csrRowIndeces_v_w.data()[eN_i] + csrColumnOffsets_eb_v_w.data()[ebN_i_j]] += fluxJacobian_v_w[j]*vel_test_dS[i]+
4512 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb],
4513 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4515 vel_trial_trace_ref.data()[ebN_local_kb_j],
4517 sdInfo_v_w_rowptr.data(),
4518 sdInfo_v_w_colind.data(),
4519 mom_vw_diff_ten_ext,
4520 &vel_grad_test_dS[i*nSpace]);
4522 globalJacobian.data()[csrRowIndeces_w_p.data()[eN_i] + csrColumnOffsets_eb_w_p.data()[ebN_i_j]] += fluxJacobian_w_p[j]*vel_test_dS[i];
4523 globalJacobian.data()[csrRowIndeces_w_u.data()[eN_i] + csrColumnOffsets_eb_w_u.data()[ebN_i_j]] += fluxJacobian_w_u[j]*vel_test_dS[i]+
4524 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4525 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4527 vel_trial_trace_ref.data()[ebN_local_kb_j],
4529 sdInfo_w_u_rowptr.data(),
4530 sdInfo_w_u_colind.data(),
4531 mom_wu_diff_ten_ext,
4532 &vel_grad_test_dS[i*nSpace]);
4533 globalJacobian.data()[csrRowIndeces_w_v.data()[eN_i] + csrColumnOffsets_eb_w_v.data()[ebN_i_j]] += fluxJacobian_w_v[j]*vel_test_dS[i]+
4534 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4535 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4537 vel_trial_trace_ref.data()[ebN_local_kb_j],
4539 sdInfo_w_v_rowptr.data(),
4540 sdInfo_w_v_colind.data(),
4541 mom_wv_diff_ten_ext,
4542 &vel_grad_test_dS[i*nSpace]);
4543 globalJacobian.data()[csrRowIndeces_w_w.data()[eN_i] + csrColumnOffsets_eb_w_w.data()[ebN_i_j]] += fluxJacobian_w_w[j]*vel_test_dS[i]+
4544 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb],
4545 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4547 vel_trial_trace_ref.data()[ebN_local_kb_j],
4549 sdInfo_w_w_rowptr.data(),
4550 sdInfo_w_w_colind.data(),
4551 mom_ww_diff_ten_ext,
4552 &vel_grad_test_dS[i*nSpace]);
4561 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4562 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4563 int nInteriorElementBoundaries_global = args.
scalar<
int>(
"nInteriorElementBoundaries_global");
4564 xt::pyarray<int>& interiorElementBoundariesArray = args.
array<
int>(
"interiorElementBoundariesArray");
4565 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4566 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4567 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4568 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4569 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4570 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4571 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4572 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4573 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4574 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4575 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4576 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4577 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4578 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
4579 xt::pyarray<double>& velocityAverage = args.
array<
double>(
"velocityAverage");
4580 int permutations[nQuadraturePoints_elementBoundary];
4581 double xArray_left[nQuadraturePoints_elementBoundary*3],
4582 xArray_right[nQuadraturePoints_elementBoundary*3];
4583 for (
int i=0;i<nQuadraturePoints_elementBoundary;i++)
4585 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
4587 register int ebN = exteriorElementBoundariesArray.data()[ebNE];
4588 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4590 register int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace,
4591 ebNE_kb_nSpace = ebNE*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace;
4592 velocityAverage.data()[ebN_kb_nSpace+0]=ebqe_velocity.data()[ebNE_kb_nSpace+0];
4593 velocityAverage.data()[ebN_kb_nSpace+1]=ebqe_velocity.data()[ebNE_kb_nSpace+1];
4594 velocityAverage.data()[ebN_kb_nSpace+2]=ebqe_velocity.data()[ebNE_kb_nSpace+2];
4597 for (
int ebNI = 0; ebNI < nInteriorElementBoundaries_global; ebNI++)
4599 register int ebN = interiorElementBoundariesArray.data()[ebNI],
4600 left_eN_global = elementBoundaryElementsArray.data()[ebN*2+0],
4601 left_ebN_element = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
4602 right_eN_global = elementBoundaryElementsArray.data()[ebN*2+1],
4603 right_ebN_element = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+1],
4604 left_eN_nDOF_trial_element = left_eN_global*nDOF_trial_element,
4605 right_eN_nDOF_trial_element = right_eN_global*nDOF_trial_element;
4606 double jac[nSpace*nSpace],
4608 jacInv[nSpace*nSpace],
4609 boundaryJac[nSpace*(nSpace-1)],
4610 metricTensor[(nSpace-1)*(nSpace-1)],
4611 metricTensorDetSqrt,
4615 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4617 ck.calculateMapping_elementBoundary(left_eN_global,
4620 left_ebN_element*nQuadraturePoints_elementBoundary+kb,
4623 mesh_trial_trace_ref.data(),
4624 mesh_grad_trial_trace_ref.data(),
4625 boundaryJac_ref.data(),
4631 metricTensorDetSqrt,
4635 xArray_left[kb*3+0] = x;
4636 xArray_left[kb*3+1] = y;
4637 xArray_left[kb*3+2] =
z;
4638 ck.calculateMapping_elementBoundary(right_eN_global,
4641 right_ebN_element*nQuadraturePoints_elementBoundary+kb,
4644 mesh_trial_trace_ref.data(),
4645 mesh_grad_trial_trace_ref.data(),
4646 boundaryJac_ref.data(),
4652 metricTensorDetSqrt,
4656 xArray_right[kb*3+0] = x;
4657 xArray_right[kb*3+1] = y;
4658 xArray_right[kb*3+2] =
z;
4660 for (
int kb_left=0;kb_left<nQuadraturePoints_elementBoundary;kb_left++)
4662 double errorNormMin = 1.0;
4663 for (
int kb_right=0;kb_right<nQuadraturePoints_elementBoundary;kb_right++)
4665 double errorNorm=0.0;
4666 for (
int I=0;I<nSpace;I++)
4668 errorNorm += fabs(xArray_left[kb_left*3+I]
4670 xArray_right[kb_right*3+I]);
4672 if (errorNorm < errorNormMin)
4674 permutations[kb_right] = kb_left;
4675 errorNormMin = errorNorm;
4679 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4681 register int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace;
4682 register double u_left=0.0,
4688 register int left_kb = kb,
4689 right_kb = permutations[kb],
4690 left_ebN_element_kb_nDOF_test_element=(left_ebN_element*nQuadraturePoints_elementBoundary+left_kb)*nDOF_test_element,
4691 right_ebN_element_kb_nDOF_test_element=(right_ebN_element*nQuadraturePoints_elementBoundary+right_kb)*nDOF_test_element;
4695 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);
4696 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);
4697 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);
4699 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);
4700 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);
4701 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);
4703 velocityAverage.data()[ebN_kb_nSpace+0]=0.5*(u_left + u_right);
4704 velocityAverage.data()[ebN_kb_nSpace+1]=0.5*(v_left + v_right);
4705 velocityAverage.data()[ebN_kb_nSpace+2]=0.5*(w_left + w_right);
4712 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
4713 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
4714 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4715 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4716 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
4717 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
4718 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
4719 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
4720 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
4721 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
4722 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
4723 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
4724 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
4725 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4726 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4727 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
4728 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
4729 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
4730 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
4731 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
4732 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4733 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
4734 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
4735 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
4736 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4737 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4738 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
4739 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
4740 double hFactor = args.
scalar<
double>(
"hFactor");
4741 int nElements_global = args.
scalar<
int>(
"nElements_global");
4742 double useRBLES = args.
scalar<
double>(
"useRBLES");
4743 double useMetrics = args.
scalar<
double>(
"useMetrics");
4744 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
4745 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
4746 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
4747 double sigma = args.
scalar<
double>(
"sigma");
4752 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
4753 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
4754 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
4755 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
4756 double C_dc = args.
scalar<
double>(
"C_dc");
4757 double C_b = args.
scalar<
double>(
"C_b");
4758 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
4759 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4760 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
4761 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4762 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4763 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4764 xt::pyarray<double>& g = args.
array<
double>(
"g");
4765 xt::pyarray<double>& rho_init = args.
array<
double>(
"rho_init");
4766 const double useVF = args.
scalar<
double>(
"useVF");
4767 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
4768 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
4769 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
4770 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
4771 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
4772 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
4773 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
4774 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
4775 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
4776 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
4777 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
4778 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
4779 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
4780 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
4781 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
4782 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
4783 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
4784 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
4785 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
4786 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
4787 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
4788 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
4789 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
4790 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
4791 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
4792 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
4793 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
4794 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
4795 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
4796 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
4797 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
4798 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
4799 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
4800 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
4801 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
4802 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
4803 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
4804 int offset_p = args.
scalar<
int>(
"offset_p");
4805 int offset_u = args.
scalar<
int>(
"offset_u");
4806 int offset_v = args.
scalar<
int>(
"offset_v");
4807 int offset_w = args.
scalar<
int>(
"offset_w");
4808 int stride_p = args.
scalar<
int>(
"stride_p");
4809 int stride_u = args.
scalar<
int>(
"stride_u");
4810 int stride_v = args.
scalar<
int>(
"stride_v");
4811 int stride_w = args.
scalar<
int>(
"stride_w");
4812 xt::pyarray<double>& cg = args.
array<
double>(
"cg");
4813 xt::pyarray<double>& force = args.
array<
double>(
"force");
4814 xt::pyarray<double>& moment = args.
array<
double>(
"moment");
4815 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4816 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4817 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4818 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4819 xt::pyarray<int>& forceExtractionFaces = args.
array<
int>(
"forceExtractionFaces");
4820 int nForceExtractionFaces = args.
scalar<
int>(
"nForceExtractionFaces");
4821 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
4822 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
4823 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
4824 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
4825 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
4826 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
4827 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
4828 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
4829 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
4830 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
4831 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
4832 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
4833 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
4834 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
4835 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
4836 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
4837 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
4838 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
4839 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
4840 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
4841 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
4842 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
4843 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
4844 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
4845 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
4846 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
4847 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
4848 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
4849 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
4850 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
4851 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
4859 double tmp1[nSpace], tmp2[nSpace*nSpace];
4861 for (
int fef = 0; fef < nForceExtractionFaces; fef++)
4863 int ebNE = forceExtractionFaces.data()[fef];
4866 register int ebN = exteriorElementBoundariesArray.data()[ebNE],
4867 eN = elementBoundaryElementsArray.data()[ebN*2+0],
4868 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
4869 eN_nDOF_trial_element = eN*nDOF_trial_element;
4870 double eps_rho,eps_mu;
4874 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4877 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
4878 ebNE_kb_nSpace = ebNE_kb*nSpace,
4879 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
4880 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
4881 register double p_ext=0.0,
4889 jac_ext[nSpace*nSpace],
4891 jacInv_ext[nSpace*nSpace],
4892 boundaryJac[nSpace*(nSpace-1)],
4893 metricTensor[(nSpace-1)*(nSpace-1)],
4894 metricTensorDetSqrt,
4895 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],vel_grad_test_dS[nDOF_test_element*nSpace],
4896 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
4897 normal[3],x_ext,y_ext,z_ext,
4898 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,gi[nSpace], unormal,gnormal,uneg,
4899 H_rho,d_rho, H_mu,d_mu, rho, mu,nu;
4901 ck.calculateMapping_elementBoundary(eN,
4907 mesh_trial_trace_ref.data(),
4908 mesh_grad_trial_trace_ref.data(),
4909 boundaryJac_ref.data(),
4915 metricTensorDetSqrt,
4919 dS = metricTensorDetSqrt*dS_ref.data()[kb];
4922 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
4923 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
4925 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
4926 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
4928 H_rho = (1.0-useVF)*
RANS2P_IB::smoothedHeaviside(eps_rho,ebqe_phi_ext.data()[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,ebqe_vf_ext.data()[ebNE_kb]));
4930 H_mu = (1.0-useVF)*
RANS2P_IB::smoothedHeaviside(eps_mu,ebqe_phi_ext.data()[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,ebqe_vf_ext.data()[ebNE_kb]));
4936 ck.gradTrialFromRef(& p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
4937 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
4939 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);
4940 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);
4941 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);
4942 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);
4943 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
4944 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
4945 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
4946 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
4948 for (
int j=0;j<nDOF_trial_element;j++)
4950 int j_nSpace = j*nSpace;
4951 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4952 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4954 vel_grad_test_dS[j_nSpace+0] = vel_grad_trial_trace[j_nSpace+0]*dS;
4955 vel_grad_test_dS[j_nSpace+1] = vel_grad_trial_trace[j_nSpace+1]*dS;
4956 vel_grad_test_dS[j_nSpace+2] = vel_grad_trial_trace[j_nSpace+2]*dS;
4959 ebqe_velocity.data()[ebNE_kb_nSpace + 0 ] = u_ext;
4960 ebqe_velocity.data()[ebNE_kb_nSpace + 1 ] = v_ext;
4961 ebqe_velocity.data()[ebNE_kb_nSpace + 2 ] = w_ext;
4964 double norm_S,h_e=elementDiameter.data()[eN],t_nu_0,t_nu_1;
4967 switch (turbulenceClosureModel)
4971 norm_S = sqrt(2.0*(grad_u_ext[0]*grad_u_ext[0] + grad_v_ext[1]*grad_v_ext[1] + grad_w_ext[2]*grad_w_ext[2] +
4972 0.5*(grad_u_ext[1]+grad_v_ext[0])*(grad_u_ext[1]+grad_v_ext[0]) +
4973 0.5*(grad_u_ext[2]+grad_w_ext[0])*(grad_u_ext[2]+grad_w_ext[0]) +
4974 0.5*(grad_v_ext[2]+grad_w_ext[1])*(grad_v_ext[2]+grad_w_ext[1])));
4975 t_nu_0 += smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
4976 t_nu_1 += smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
4980 double re_0,cs_0,re_1,cs_1;
4981 norm_S = sqrt(2.0*(grad_u_ext[0]*grad_u_ext[0] + grad_v_ext[1]*grad_v_ext[1] + grad_w_ext[2]*grad_w_ext[2] +
4982 0.5*(grad_u_ext[1]+grad_v_ext[0])*(grad_u_ext[1]+grad_v_ext[0]) +
4983 0.5*(grad_u_ext[2]+grad_w_ext[0])*(grad_u_ext[2]+grad_w_ext[0]) +
4984 0.5*(grad_v_ext[2]+grad_w_ext[1])*(grad_v_ext[2]+grad_w_ext[1])));
4985 re_0 = h_e*h_e*norm_S/
nu_0;
4986 cs_0=0.027*pow(10.0,-3.23*pow(re_0,-0.92));
4987 t_nu_0 += cs_0*h_e*h_e*norm_S;
4988 re_1 = h_e*h_e*norm_S/
nu_1;
4989 cs_1=0.027*pow(10.0,-3.23*pow(re_1,-0.92));
4990 t_nu_1 += cs_1*h_e*h_e*norm_S;
4995 nu = t_nu_0*(1.0-H_mu)+t_nu_1*H_mu;
4996 mu =
rho_0*t_nu_0*(1.0-H_mu)+
rho_1*t_nu_1*H_mu;
5006 gi[0] = ebqe_bc_u_ext.data()[ebNE_kb];
5007 gi[1] = ebqe_bc_v_ext.data()[ebNE_kb];
5008 gi[2] = ebqe_bc_w_ext.data()[ebNE_kb];
5010 unormal = normal[0]*u_ext + normal[1]*v_ext + normal[2]*w_ext;
5011 gnormal = normal[0]*gi[0] + normal[1]*gi[1] + normal[2]*gi[2];
5012 uneg = 0.5*(unormal - fabs(unormal));
5017 ck.calculateGScale(G,normal,h_penalty);
5018 penalty = C_b*mu/h_penalty;
5023 tmp1[0] = -mu*(grad_u_ext[0]*normal[0] + grad_u_ext[0]*normal[0]+
5024 grad_u_ext[1]*normal[1] + grad_v_ext[0]*normal[1]+
5025 grad_u_ext[2]*normal[2] + grad_w_ext[0]*normal[2])
5026 + penalty*(u_ext-gi[0]) + p_ext*normal[0] - uneg*rho*(u_ext-gi[0]);
5028 tmp1[1] = -mu*(grad_v_ext[0]*normal[0] + grad_u_ext[1]*normal[0]+
5029 grad_v_ext[1]*normal[1] + grad_v_ext[1]*normal[1]+
5030 grad_v_ext[2]*normal[2] + grad_w_ext[1]*normal[2])
5031 + penalty*(v_ext-gi[1]) + p_ext*normal[1] - uneg*rho*(v_ext-gi[1]);
5033 tmp1[2] = -mu*(grad_w_ext[0]*normal[0] + grad_u_ext[2]*normal[0]+
5034 grad_w_ext[1]*normal[1] + grad_v_ext[2]*normal[1]+
5035 grad_w_ext[2]*normal[2] + grad_w_ext[2]*normal[2])
5036 + penalty*(w_ext-gi[2]) + p_ext*normal[2] - uneg*rho*(w_ext-gi[2]);
5038 tmp2[0*nSpace + 0] = 0.0;
5039 tmp2[0*nSpace + 1] = mu*( (v_ext - gi[1])*normal[0] + (u_ext - gi[0])*normal[1] );
5040 tmp2[0*nSpace + 2] = mu*( (w_ext - gi[2])*normal[0] + (u_ext - gi[0])*normal[2] );
5042 tmp2[1*nSpace + 0] = mu*( (u_ext - gi[0])*normal[1] + (v_ext - gi[1])*normal[0] );
5043 tmp2[1*nSpace + 1] = 0.0;
5044 tmp2[1*nSpace + 2] = mu*( (w_ext - gi[2])*normal[1] + (v_ext - gi[1])*normal[2] );
5046 tmp2[2*nSpace + 0] = mu*( (u_ext - gi[0])*normal[2] + (w_ext - gi[2])*normal[0] );
5047 tmp2[2*nSpace + 1] = mu*( (v_ext - gi[1])*normal[2] + (w_ext - gi[2])*normal[1] );
5048 tmp2[2*nSpace + 2] = 0.0;
5053 force.data()[0] += tmp1[0]*dS;
5054 force.data()[1] += tmp1[1]*dS;
5055 force.data()[2] += tmp1[2]*dS;
5060 moment.data()[0] += ( (y_ext-cg.data()[1])*tmp1[2] - (z_ext-cg.data()[2])*tmp1[1] + tmp2[2*nSpace + 1] - tmp2[1*nSpace + 2] )*dS;
5061 moment.data()[1] += ( (z_ext-cg.data()[2])*tmp1[0] - (x_ext-cg.data()[0])*tmp1[2] + tmp2[0*nSpace + 2] - tmp2[2*nSpace + 0] )*dS;
5062 moment.data()[2] += ( (x_ext-cg.data()[0])*tmp1[1] - (y_ext-cg.data()[1])*tmp1[0] + tmp2[1*nSpace + 0] - tmp2[0*nSpace + 1] )*dS;
5072 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
5073 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
5074 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
5075 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
5076 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
5077 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
5078 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
5079 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
5080 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
5081 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
5082 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
5083 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
5084 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
5085 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
5086 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
5087 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
5088 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
5089 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
5090 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
5091 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
5092 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
5093 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
5094 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
5095 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
5096 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
5097 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
5098 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
5099 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
5100 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
5101 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
5102 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
5103 double hFactor = args.
scalar<
double>(
"hFactor");
5104 int nElements_global = args.
scalar<
int>(
"nElements_global");
5105 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
5106 double useRBLES = args.
scalar<
double>(
"useRBLES");
5107 double useMetrics = args.
scalar<
double>(
"useMetrics");
5108 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
5109 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
5110 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
5111 double sigma = args.
scalar<
double>(
"sigma");
5116 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
5117 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
5118 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
5119 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
5120 double C_dc = args.
scalar<
double>(
"C_dc");
5121 double C_b = args.
scalar<
double>(
"C_b");
5122 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
5123 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
5124 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
5125 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
5126 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
5127 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
5128 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
5129 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
5130 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
5131 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
5132 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
5133 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
5134 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
5135 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
5136 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
5137 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
5138 xt::pyarray<double>& g = args.
array<
double>(
"g");
5139 const double useVF = args.
scalar<
double>(
"useVF");
5140 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
5141 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
5142 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
5143 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
5144 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
5145 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
5146 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
5147 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
5148 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
5149 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
5150 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
5151 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
5152 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
5153 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
5154 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
5155 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
5156 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
5157 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
5158 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
5159 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
5160 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
5161 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
5162 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
5163 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
5164 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
5165 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
5166 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
5167 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
5168 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
5169 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
5170 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
5171 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
5172 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
5173 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
5174 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
5175 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
5176 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
5177 int offset_p = args.
scalar<
int>(
"offset_p");
5178 int offset_u = args.
scalar<
int>(
"offset_u");
5179 int offset_v = args.
scalar<
int>(
"offset_v");
5180 int offset_w = args.
scalar<
int>(
"offset_w");
5181 int stride_p = args.
scalar<
int>(
"stride_p");
5182 int stride_u = args.
scalar<
int>(
"stride_u");
5183 int stride_v = args.
scalar<
int>(
"stride_v");
5184 int stride_w = args.
scalar<
int>(
"stride_w");
5185 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
5186 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
5187 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
5188 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
5189 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
5190 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
5191 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
5192 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
5193 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
5194 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
5195 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
5196 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
5197 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
5198 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
5199 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
5200 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
5201 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
5202 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
5203 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
5204 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
5205 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
5206 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
5207 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
5208 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
5209 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
5210 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
5211 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
5212 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
5213 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
5214 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
5215 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
5216 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
5217 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
5218 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
5219 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
5220 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
5221 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
5222 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
5223 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
5224 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
5225 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
5226 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
5227 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
5228 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
5229 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
5230 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
5231 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
5232 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
5233 xt::pyarray<double>& q_dragBeam1 = args.
array<
double>(
"q_dragBeam1");
5234 xt::pyarray<double>& q_dragBeam2 = args.
array<
double>(
"q_dragBeam2");
5235 xt::pyarray<double>& q_dragBeam3 = args.
array<
double>(
"q_dragBeam3");
5236 xt::pyarray<double>& ebqe_dragBeam1 = args.
array<
double>(
"ebqe_dragBeam1");
5237 xt::pyarray<double>& ebqe_dragBeam2 = args.
array<
double>(
"ebqe_dragBeam2");
5238 xt::pyarray<double>& ebqe_dragBeam3 = args.
array<
double>(
"ebqe_dragBeam3");
5239 int nBeams = args.
scalar<
int>(
"nBeams");
5240 int nBeamElements = args.
scalar<
int>(
"nBeamElements");
5241 int beam_quadOrder = args.
scalar<
int>(
"beam_quadOrder");
5242 double beam_Cd = args.
scalar<
double>(
"beam_Cd");
5243 xt::pyarray<double>& beamRadius = args.
array<
double>(
"beamRadius");
5244 xt::pyarray<double>& xq = args.
array<
double>(
"xq");
5245 xt::pyarray<double>& yq = args.
array<
double>(
"yq");
5246 xt::pyarray<double>& zq = args.
array<
double>(
"zq");
5247 xt::pyarray<double>& Beam_h = args.
array<
double>(
"Beam_h");
5248 xt::pyarray<double>& dV_beam = args.
array<
double>(
"dV_beam");
5249 xt::pyarray<double>& q1 = args.
array<
double>(
"q1");
5250 xt::pyarray<double>& q2 = args.
array<
double>(
"q2");
5251 xt::pyarray<double>& q3 = args.
array<
double>(
"q3");
5252 xt::pyarray<double>& vel_avg = args.
array<
double>(
"vel_avg");
5253 xt::pyarray<double>& netBeamDrag = args.
array<
double>(
"netBeamDrag");
5257 double globalConservationError=0.0;
5258 for(
int eN=0;eN<nElements_global;eN++)
5260 register double eps_rho,eps_mu;
5264 for(
int k=0;k<nQuadraturePoints_element;k++)
5267 register int eN_k = eN*nQuadraturePoints_element+k,
5268 eN_k_nSpace = eN_k*nSpace,
5269 eN_nDOF_trial_element = eN*nDOF_trial_element;
5270 register double p=0.0,
u=0.0,
v=0.0,
w=0.0,
5271 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
5279 dmass_adv_u[nSpace],
5280 dmass_adv_v[nSpace],
5281 dmass_adv_w[nSpace],
5283 dmom_u_adv_u[nSpace],
5284 dmom_u_adv_v[nSpace],
5285 dmom_u_adv_w[nSpace],
5287 dmom_v_adv_u[nSpace],
5288 dmom_v_adv_v[nSpace],
5289 dmom_v_adv_w[nSpace],
5291 dmom_w_adv_u[nSpace],
5292 dmom_w_adv_v[nSpace],
5293 dmom_w_adv_w[nSpace],
5294 mom_uu_diff_ten[nSpace],
5295 mom_vv_diff_ten[nSpace],
5296 mom_ww_diff_ten[nSpace],
5307 dmom_u_ham_grad_p[nSpace],
5309 dmom_v_ham_grad_p[nSpace],
5311 dmom_w_ham_grad_p[nSpace],
5322 Lstar_u_p[nDOF_test_element],
5323 Lstar_v_p[nDOF_test_element],
5324 Lstar_w_p[nDOF_test_element],
5325 Lstar_u_u[nDOF_test_element],
5326 Lstar_v_v[nDOF_test_element],
5327 Lstar_w_w[nDOF_test_element],
5328 Lstar_p_u[nDOF_test_element],
5329 Lstar_p_v[nDOF_test_element],
5330 Lstar_p_w[nDOF_test_element],
5335 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
5336 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
5339 jacInv[nSpace*nSpace],
5340 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
5341 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
5342 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
5348 dmom_u_source[nSpace],
5349 dmom_v_source[nSpace],
5350 dmom_w_source[nSpace],
5352 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
5354 ck.calculateMapping_element(eN,
5358 mesh_trial_ref.data(),
5359 mesh_grad_trial_ref.data(),
5364 ck.calculateH_element(eN,
5366 nodeDiametersArray.data(),
5368 mesh_trial_ref.data(),
5371 ck.calculateMappingVelocity_element(eN,
5373 mesh_velocity_dof.data(),
5375 mesh_trial_ref.data(),
5380 dV = fabs(jacDet)*dV_ref.data()[k];
5381 ck.calculateG(jacInv,G,G_dd_G,tr_G);
5384 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5385 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5388 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
5390 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
u);
5391 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
v);
5392 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
w);
5395 for (
int j=0;j<nDOF_trial_element;j++)
5397 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
5400 porosity = q_porosity.data()[eN_k];
5404 q_velocity.data()[eN_k_nSpace+0]=
u;
5405 q_velocity.data()[eN_k_nSpace+1]=
v;
5406 q_velocity.data()[eN_k_nSpace+2]=
w;
5407 q_x.data()[eN_k_nSpace+0]=x;
5408 q_x.data()[eN_k_nSpace+1]=y;
5409 q_x.data()[eN_k_nSpace+2]=
z;
5421 q_mom_u_acc_beta_bdf.data()[eN_k],
5427 q_mom_v_acc_beta_bdf.data()[eN_k],
5433 q_mom_w_acc_beta_bdf.data()[eN_k],
5459 q_dragBeam1.data()[eN_k],
5460 q_dragBeam2.data()[eN_k],
5461 q_dragBeam3.data()[eN_k],
5463 netBeamDrag.data());
5500 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
5502 register int ebN = exteriorElementBoundariesArray.data()[ebNE],
5503 eN = elementBoundaryElementsArray.data()[ebN*2+0],
5504 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
5505 eN_nDOF_trial_element = eN*nDOF_trial_element;
5506 register double eps_rho, eps_mu;
5507 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5509 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
5510 ebNE_kb_nSpace = ebNE_kb*nSpace,
5511 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5512 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5513 register double p_ext=0.0,
5522 dmom_u_acc_u_ext=0.0,
5524 dmom_v_acc_v_ext=0.0,
5526 dmom_w_acc_w_ext=0.0,
5527 mass_adv_ext[nSpace],
5528 dmass_adv_u_ext[nSpace],
5529 dmass_adv_v_ext[nSpace],
5530 dmass_adv_w_ext[nSpace],
5531 mom_u_adv_ext[nSpace],
5532 dmom_u_adv_u_ext[nSpace],
5533 dmom_u_adv_v_ext[nSpace],
5534 dmom_u_adv_w_ext[nSpace],
5535 mom_v_adv_ext[nSpace],
5536 dmom_v_adv_u_ext[nSpace],
5537 dmom_v_adv_v_ext[nSpace],
5538 dmom_v_adv_w_ext[nSpace],
5539 mom_w_adv_ext[nSpace],
5540 dmom_w_adv_u_ext[nSpace],
5541 dmom_w_adv_v_ext[nSpace],
5542 dmom_w_adv_w_ext[nSpace],
5543 mom_uu_diff_ten_ext[nSpace],
5544 mom_vv_diff_ten_ext[nSpace],
5545 mom_ww_diff_ten_ext[nSpace],
5546 mom_uv_diff_ten_ext[1],
5547 mom_uw_diff_ten_ext[1],
5548 mom_vu_diff_ten_ext[1],
5549 mom_vw_diff_ten_ext[1],
5550 mom_wu_diff_ten_ext[1],
5551 mom_wv_diff_ten_ext[1],
5552 mom_u_source_ext=0.0,
5553 mom_v_source_ext=0.0,
5554 mom_w_source_ext=0.0,
5556 dmom_u_ham_grad_p_ext[nSpace],
5558 dmom_v_ham_grad_p_ext[nSpace],
5560 dmom_w_ham_grad_p_ext[nSpace],
5561 dmom_u_adv_p_ext[nSpace],
5562 dmom_v_adv_p_ext[nSpace],
5563 dmom_w_adv_p_ext[nSpace],
5565 flux_mom_u_adv_ext=0.0,
5566 flux_mom_v_adv_ext=0.0,
5567 flux_mom_w_adv_ext=0.0,
5568 flux_mom_uu_diff_ext=0.0,
5569 flux_mom_uv_diff_ext=0.0,
5570 flux_mom_uw_diff_ext=0.0,
5571 flux_mom_vu_diff_ext=0.0,
5572 flux_mom_vv_diff_ext=0.0,
5573 flux_mom_vw_diff_ext=0.0,
5574 flux_mom_wu_diff_ext=0.0,
5575 flux_mom_wv_diff_ext=0.0,
5576 flux_mom_ww_diff_ext=0.0,
5581 bc_mom_u_acc_ext=0.0,
5582 bc_dmom_u_acc_u_ext=0.0,
5583 bc_mom_v_acc_ext=0.0,
5584 bc_dmom_v_acc_v_ext=0.0,
5585 bc_mom_w_acc_ext=0.0,
5586 bc_dmom_w_acc_w_ext=0.0,
5587 bc_mass_adv_ext[nSpace],
5588 bc_dmass_adv_u_ext[nSpace],
5589 bc_dmass_adv_v_ext[nSpace],
5590 bc_dmass_adv_w_ext[nSpace],
5591 bc_mom_u_adv_ext[nSpace],
5592 bc_dmom_u_adv_u_ext[nSpace],
5593 bc_dmom_u_adv_v_ext[nSpace],
5594 bc_dmom_u_adv_w_ext[nSpace],
5595 bc_mom_v_adv_ext[nSpace],
5596 bc_dmom_v_adv_u_ext[nSpace],
5597 bc_dmom_v_adv_v_ext[nSpace],
5598 bc_dmom_v_adv_w_ext[nSpace],
5599 bc_mom_w_adv_ext[nSpace],
5600 bc_dmom_w_adv_u_ext[nSpace],
5601 bc_dmom_w_adv_v_ext[nSpace],
5602 bc_dmom_w_adv_w_ext[nSpace],
5603 bc_mom_uu_diff_ten_ext[nSpace],
5604 bc_mom_vv_diff_ten_ext[nSpace],
5605 bc_mom_ww_diff_ten_ext[nSpace],
5606 bc_mom_uv_diff_ten_ext[1],
5607 bc_mom_uw_diff_ten_ext[1],
5608 bc_mom_vu_diff_ten_ext[1],
5609 bc_mom_vw_diff_ten_ext[1],
5610 bc_mom_wu_diff_ten_ext[1],
5611 bc_mom_wv_diff_ten_ext[1],
5612 bc_mom_u_source_ext=0.0,
5613 bc_mom_v_source_ext=0.0,
5614 bc_mom_w_source_ext=0.0,
5615 bc_mom_u_ham_ext=0.0,
5616 bc_dmom_u_ham_grad_p_ext[nSpace],
5617 bc_mom_v_ham_ext=0.0,
5618 bc_dmom_v_ham_grad_p_ext[nSpace],
5619 bc_mom_w_ham_ext=0.0,
5620 bc_dmom_w_ham_grad_p_ext[nSpace],
5621 jac_ext[nSpace*nSpace],
5623 jacInv_ext[nSpace*nSpace],
5624 boundaryJac[nSpace*(nSpace-1)],
5625 metricTensor[(nSpace-1)*(nSpace-1)],
5626 metricTensorDetSqrt,
5627 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
5628 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
5629 vel_grad_test_dS[nDOF_trial_element*nSpace],
5630 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
5634 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
5635 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, junk[nSpace];
5637 ck.calculateMapping_elementBoundary(eN,
5643 mesh_trial_trace_ref.data(),
5644 mesh_grad_trial_trace_ref.data(),
5645 boundaryJac_ref.data(),
5651 metricTensorDetSqrt,
5655 ck.calculateMappingVelocity_elementBoundary(eN,
5659 mesh_velocity_dof.data(),
5661 mesh_trial_trace_ref.data(),
5662 xt_ext,yt_ext,zt_ext,
5672 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
5675 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
5676 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
5678 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5679 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5683 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
5686 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);
5687 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);
5688 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);
5691 for (
int j=0;j<nDOF_trial_element;j++)
5693 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
5696 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
5697 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
5698 bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext;
5700 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
5718 phi.data()[ebNE_kb],
5721 ebqe_dragBeam1.data()[ebNE_kb],
5722 ebqe_dragBeam2.data()[ebNE_kb],
5723 ebqe_dragBeam3.data()[ebNE_kb],
5741 int nQuadraturePoints_elementIn,
5742 int nDOF_mesh_trial_elementIn,
5743 int nDOF_trial_elementIn,
5744 int nDOF_test_elementIn,
5745 int nQuadraturePoints_elementBoundaryIn,
5748 return proteus::chooseAndAllocateDiscretization<RANS2P_IB_base,RANS2P_IB,CompKernel>(nSpaceIn,
5749 nQuadraturePoints_elementIn,
5750 nDOF_mesh_trial_elementIn,
5751 nDOF_trial_elementIn,
5752 nDOF_test_elementIn,
5753 nQuadraturePoints_elementBoundaryIn,