14 #include "xtensor-python/pyarray.hpp"
17 return double((0.0 < val) - (val < 0.0));
34 #define CELL_BASED_EV_COEFF 1
35 #define POWER_SMOOTHNESS_INDICATOR 2
36 #define EPS_FOR_GAMMA_INDICATOR 1E-10
37 #define C_FOR_GAMMA_INDICATOR 0.25 // increase gamma to make the indicator more agressive (less dissipative)
38 #define USE_GAMMA_INDICATOR 0
39 #define ANISOTROPIC_DIFFUSION 0
41 static const double DM=0.0;
42 static const double DM2=0.0;
43 static const double DM3=1.0;
45 #define USE_CYLINDER_AS_PARTICLE//just for debug
49 template<
int nSpace,
int nP,
int nQ,
int nEBQ>
89 template<
class CompKernelType,
91 int nQuadraturePoints_element,
92 int nDOF_mesh_trial_element,
93 int nDOF_trial_element,
94 int nDOF_test_element,
95 int nQuadraturePoints_elementBoundary>
159 double mu_fr_limiter)
180 inline double Dot(
const double vec1[nSpace],
181 const double vec2[nSpace])
184 for (
int I=0; I<nSpace; I++)
185 dot += vec1[I]*vec2[I];
190 const double vel_grad[nSpace],
191 double vel_tgrad[nSpace])
193 double normal_dot_vel_grad =
Dot(normal,vel_grad);
194 for (
int I=0; I<nSpace; I++)
195 vel_tgrad[I] = vel_grad[I] - normal_dot_vel_grad*normal[I];
208 const double smagorinskyConstant,
209 const int turbulenceClosureModel,
210 const double g[nSpace],
214 const double n[nSpace],
215 const double distance_to_omega_solid,
217 const double porosity,
219 const double grad_p[nSpace],
220 const double grad_u[nSpace],
221 const double grad_v[nSpace],
222 const double grad_w[nSpace],
229 double& eddy_viscosity,
231 double& dmom_u_acc_u,
233 double& dmom_v_acc_v,
235 double& dmom_w_acc_w,
236 double mass_adv[nSpace],
237 double dmass_adv_u[nSpace],
238 double dmass_adv_v[nSpace],
239 double dmass_adv_w[nSpace],
240 double mom_u_adv[nSpace],
241 double dmom_u_adv_u[nSpace],
242 double dmom_u_adv_v[nSpace],
243 double dmom_u_adv_w[nSpace],
244 double mom_v_adv[nSpace],
245 double dmom_v_adv_u[nSpace],
246 double dmom_v_adv_v[nSpace],
247 double dmom_v_adv_w[nSpace],
248 double mom_w_adv[nSpace],
249 double dmom_w_adv_u[nSpace],
250 double dmom_w_adv_v[nSpace],
251 double dmom_w_adv_w[nSpace],
252 double mom_uu_diff_ten[nSpace],
253 double mom_vv_diff_ten[nSpace],
254 double mom_ww_diff_ten[nSpace],
255 double mom_uv_diff_ten[1],
256 double mom_uw_diff_ten[1],
257 double mom_vu_diff_ten[1],
258 double mom_vw_diff_ten[1],
259 double mom_wu_diff_ten[1],
260 double mom_wv_diff_ten[1],
261 double& mom_u_source,
262 double& mom_v_source,
263 double& mom_w_source,
265 double dmom_u_ham_grad_p[nSpace],
266 double dmom_u_ham_grad_u[nSpace],
268 double dmom_v_ham_grad_p[nSpace],
269 double dmom_v_ham_grad_v[nSpace],
271 double dmom_w_ham_grad_p[nSpace],
272 double dmom_w_ham_grad_w[nSpace],
275 int KILL_PRESSURE_TERM,
276 int MULTIPLY_EXTERNAL_FORCE_BY_DENSITY,
280 int MATERIAL_PARAMETERS_AS_FUNCTION,
281 double density_as_function,
282 double dynamic_viscosity_as_function,
284 double x,
double y,
double z,
285 int use_ball_as_particle,
288 double* ball_velocity,
289 double* ball_angular_velocity,
291 int INT_BY_PARTS_PRESSURE)
293 double rho,nu,mu,H_rho,ImH_rho,d_rho,H_mu,ImH_mu,d_mu,norm_n,nu_t0=0.0,nu_t1=0.0,nu_t;
294 H_rho = (1.0-useVF)*
gf.H(eps_rho,
phi) + useVF*fmin(1.0,fmax(0.0,vf));
295 ImH_rho = (1.0-useVF)*
gf.ImH(eps_rho,
phi) + useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
296 d_rho = (1.0-useVF)*
gf.D(eps_rho,
phi);
297 H_mu = (1.0-useVF)*
gf.H(eps_mu,
phi) + useVF*fmin(1.0,fmax(0.0,vf));
298 ImH_mu = (1.0-useVF)*
gf.ImH(eps_mu,
phi) + useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
299 d_mu = (1.0-useVF)*
gf.D(eps_mu,
phi);
302 switch (turbulenceClosureModel)
307 norm_S = sqrt(2.0*(grad_u[0]*grad_u[0] + grad_v[1]*grad_v[1] + grad_w[2]*grad_w[2] +
308 0.5*(grad_u[1]+grad_v[0])*(grad_u[1]+grad_v[0]) +
309 0.5*(grad_u[2]+grad_w[0])*(grad_u[2]+grad_w[0]) +
310 0.5*(grad_v[2]+grad_w[1])*(grad_v[2]+grad_w[1])));
311 nu_t0 = smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
312 nu_t1 = smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
316 double re_0,cs_0=0.0,re_1,cs_1=0.0;
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 re_0 = h_e*h_e*norm_S/
nu_0;
323 cs_0=0.027*pow(10.0,-3.23*pow(re_0,-0.92));
324 nu_t0 = cs_0*h_e*h_e*norm_S;
325 re_1 = h_e*h_e*norm_S/
nu_1;
327 cs_1=0.027*pow(10.0,-3.23*pow(re_1,-0.92));
328 nu_t1 = cs_1*h_e*h_e*norm_S;
332 if (MATERIAL_PARAMETERS_AS_FUNCTION==0)
335 nu_t= nu_t0*ImH_mu+nu_t1*H_mu;
342 rho = density_as_function;
344 mu = dynamic_viscosity_as_function;
351 eddy_viscosity = nu_t*rho;
355 double phi_s_effect = (distance_to_omega_solid > 0.0) ? 1.0 : 1e-10;
361 dmom_u_acc_u=phi_s_effect * rho*porosity;
365 dmom_v_acc_v=phi_s_effect * rho*porosity;
368 mom_w_acc=phi_s_effect*
w;
369 dmom_w_acc_w=phi_s_effect*rho*porosity;
372 mass_adv[0]=phi_s_effect * porosity*
u;
373 mass_adv[1]=phi_s_effect * porosity*
v;
374 mass_adv[2]=phi_s_effect * porosity*
w;
376 dmass_adv_u[0]=phi_s_effect * porosity;
381 dmass_adv_v[1]=phi_s_effect * porosity;
386 dmass_adv_w[2]=phi_s_effect * porosity;
441 mom_uu_diff_ten[0] = phi_s_effect*porosity*2.0*mu;
442 mom_uu_diff_ten[1] = phi_s_effect*porosity*mu;
443 mom_uu_diff_ten[2] = phi_s_effect*porosity*mu;
445 mom_uv_diff_ten[0]=phi_s_effect * porosity*mu;
447 mom_uw_diff_ten[0]=phi_s_effect*porosity*mu;
450 mom_vv_diff_ten[0] = phi_s_effect*porosity*mu;
451 mom_vv_diff_ten[1] = phi_s_effect*porosity*2.0*mu;
452 mom_vv_diff_ten[2] = phi_s_effect*porosity*mu;
454 mom_vu_diff_ten[0]=phi_s_effect * porosity*mu;
456 mom_vw_diff_ten[0]=phi_s_effect*porosity*mu;
459 mom_ww_diff_ten[0] = phi_s_effect*porosity*mu;
460 mom_ww_diff_ten[1] = phi_s_effect*porosity*mu;
461 mom_ww_diff_ten[2] = phi_s_effect*porosity*2.0*mu;
463 mom_wu_diff_ten[0]=phi_s_effect*porosity*mu;
465 mom_wv_diff_ten[0]=phi_s_effect*porosity*mu;
468 norm_n = sqrt(
n[0]*
n[0]+
n[1]*
n[1]+
n[2]*
n[2]);
470 mom_u_source = -phi_s_effect*porosity*rho*g[0];
471 mom_v_source = -phi_s_effect*porosity*rho*g[1];
472 mom_w_source = -phi_s_effect*porosity*rho*g[2];
475 mom_u_source -= (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcex;
476 mom_v_source -= (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcey;
477 mom_w_source -= (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcez;
480 double aux_pressure = (KILL_PRESSURE_TERM==1 ? 0. : 1.)*(INT_BY_PARTS_PRESSURE==1 ? 0. : 1.);
481 mom_u_ham = phi_s_effect * porosity*grad_p[0]*aux_pressure;
482 dmom_u_ham_grad_p[0]=phi_s_effect * porosity*aux_pressure;
483 dmom_u_ham_grad_p[1]=0.0;
484 dmom_u_ham_grad_p[2]=0.0;
487 mom_v_ham = phi_s_effect * porosity*grad_p[1]*aux_pressure;
488 dmom_v_ham_grad_p[0]=0.0;
489 dmom_v_ham_grad_p[1]=phi_s_effect * porosity*aux_pressure;
490 dmom_v_ham_grad_p[2]=0.0;
493 mom_w_ham = phi_s_effect*porosity*grad_p[2]*aux_pressure;
494 dmom_w_ham_grad_p[0]=0.0;
495 dmom_w_ham_grad_p[1]=0.0;
496 dmom_w_ham_grad_p[2]=phi_s_effect*porosity*aux_pressure;
499 mom_u_ham += phi_s_effect*rho*porosity*(uStar*grad_u[0]+vStar*grad_u[1]+wStar*grad_u[2]);
500 dmom_u_ham_grad_u[0]=phi_s_effect*rho*porosity*uStar;
501 dmom_u_ham_grad_u[1]=phi_s_effect*rho*porosity*vStar;
502 dmom_u_ham_grad_u[2]=phi_s_effect*rho*porosity*wStar;
505 mom_v_ham += phi_s_effect*rho*porosity*(uStar*grad_v[0]+vStar*grad_v[1]+wStar*grad_v[2]);
506 dmom_v_ham_grad_v[0]=phi_s_effect*rho*porosity*uStar;
507 dmom_v_ham_grad_v[1]=phi_s_effect*rho*porosity*vStar;
508 dmom_v_ham_grad_v[2]=phi_s_effect*rho*porosity*wStar;
511 mom_w_ham += phi_s_effect*rho*porosity*(uStar*grad_w[0]+vStar*grad_w[1]+wStar*grad_w[2]);
512 dmom_w_ham_grad_w[0]=phi_s_effect*rho*porosity*uStar;
513 dmom_w_ham_grad_w[1]=phi_s_effect*rho*porosity*vStar;
514 dmom_w_ham_grad_w[2]=phi_s_effect*rho*porosity*wStar;
525 const double eps_rho,
546 const double uStar_s,
547 const double vStar_s,
548 const double wStar_s,
549 double& mom_u_source,
550 double& mom_v_source,
551 double& mom_w_source,
552 double dmom_u_source[nSpace],
553 double dmom_v_source[nSpace],
554 double dmom_w_source[nSpace],
559 double rho, mu,nu,H_mu,ImH_mu,uc,duc_du,duc_dv,duc_dw,viscosity,H_s;
560 H_mu = (1.0-useVF)*
gf.H(eps_mu,
phi)+useVF*fmin(1.0,fmax(0.0,vf));
561 ImH_mu = (1.0-useVF)*
gf.ImH(eps_mu,
phi) + useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
567 duc_du =
u/(uc+1.0e-12);
568 duc_dv =
v/(uc+1.0e-12);
569 duc_dw =
w/(uc+1.0e-12);
570 double fluid_velocity[3]={uStar,vStar,wStar}, solid_velocity[3]={uStar_s,vStar_s,wStar_s};
577 mom_u_source += (1.0 - phi_s)*new_beta*( (
u - uStar_s) - nu_t*gradC_x/
closure.
sigmaC_ );
578 mom_v_source += (1.0 - phi_s)*new_beta*( (
v - vStar_s)- nu_t*gradC_y/
closure.
sigmaC_ );
579 mom_w_source += (1.0 - phi_s)*new_beta*( (
w - wStar_s)- nu_t*gradC_z/
closure.
sigmaC_ );;
581 dmom_u_source[0] = (1.0 - phi_s)*new_beta;
582 dmom_u_source[1] = 0.0;
583 dmom_u_source[2] = 0.0;
585 dmom_v_source[0] = 0.0;
586 dmom_v_source[1] = (1.0 - phi_s)*new_beta;
587 dmom_v_source[2] = 0.0;
589 dmom_w_source[0] = 0.0;
590 dmom_w_source[1] = 0.0;
591 dmom_w_source[2] = (1.0 - phi_s)*new_beta;
595 const double particle_nitsche,
597 const int nParticles,
599 double *particle_signed_distances,
600 double *particle_signed_distance_normals,
601 double *particle_velocities,
602 double *particle_centroids,
603 int use_ball_as_particle,
606 double* ball_velocity,
607 double* ball_angular_velocity,
608 const double porosity,
609 const double penalty,
612 const double eps_rho,
632 const double grad_u[nSpace],
633 const double grad_v[nSpace],
634 const double grad_w[nSpace],
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],
641 double mom_u_adv[nSpace],
642 double mom_v_adv[nSpace],
643 double mom_w_adv[nSpace],
644 double dmom_u_adv_u[nSpace],
645 double dmom_v_adv_v[nSpace],
646 double dmom_w_adv_w[nSpace],
648 double dmom_u_ham_grad_u[nSpace],
650 double dmom_v_ham_grad_v[nSpace],
652 double dmom_w_ham_grad_w[nSpace],
653 double *particle_netForces,
654 double *particle_netMoments,
655 double *particle_surfaceArea)
657 double C,rho, mu,nu,H_mu,ImH_mu,uc,duc_du,duc_dv,duc_dw,H_s,D_s,phi_s,u_s,v_s,w_s,force_x,force_y,force_z,r_x,r_y,r_z;
658 double phi_s_normal[3];
659 double fluid_outward_normal[3];
662 H_mu = (1.0-useVF)*
gf.H(eps_mu,
phi)+useVF*fmin(1.0,fmax(0.0,vf));
663 ImH_mu = (1.0-useVF)*
gf.ImH(eps_mu,
phi) + useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
668 for (
int i=0;i<nParticles;i++)
670 double* vel_pointer= &particle_velocities[i*sd_offset*nSpace];
671 if(use_ball_as_particle==1)
674 get_normal_to_ith_ball(nParticles,ball_center,ball_radius,i,x,y,
z,phi_s_normal[0],phi_s_normal[1],phi_s_normal[2]);
676 ball_velocity,ball_angular_velocity,
679 center[0] = ball_center[3*i + 0];
680 center[1] = ball_center[3*i + 1];
681 center[2] = ball_center[3*i + 2];
688 phi_s = particle_signed_distances[i * sd_offset];
689 phi_s_normal[0] = particle_signed_distance_normals[i * sd_offset * nSpace + 0];
690 phi_s_normal[1] = particle_signed_distance_normals[i * sd_offset * nSpace + 1];
691 phi_s_normal[2] = particle_signed_distance_normals[i * sd_offset * nSpace + 2];
692 vel[0] = particle_velocities[i * sd_offset * nSpace + 0];
693 vel[1] = particle_velocities[i * sd_offset * nSpace + 1];
694 vel[2] = particle_velocities[i * sd_offset * nSpace + 2];
695 u_s = vel_pointer[0];
696 v_s = vel_pointer[1];
697 w_s = vel_pointer[2];
698 center[0] = particle_centroids[3*i + 0];
699 center[1] = particle_centroids[3*i + 1];
700 center[2] = particle_centroids[3*i + 2];
702 fluid_outward_normal[0] = -phi_s_normal[0];
703 fluid_outward_normal[1] = -phi_s_normal[1];
704 fluid_outward_normal[2] = -phi_s_normal[2];
706 H_s =
gf_s.H(eps_s, phi_s);
707 D_s =
gf_s.D(eps_s, phi_s);
708 double rel_vel_norm=sqrt((uStar - u_s)*(uStar - u_s)+
709 (vStar - v_s)*(vStar - v_s)+
710 (wStar - w_s)*(wStar - w_s));
711 double C_surf = (phi_s > 0.0) ? 0.0 : nu*penalty;
712 double C_vol = (phi_s > 0.0) ? 0.0 : (alpha + beta*rel_vel_norm);
713 C = (D_s*C_surf +
gf_s.ImH(eps_s, phi_s) * C_vol);
714 force_x = dV*D_s*(p*fluid_outward_normal[0] - porosity*mu*(fluid_outward_normal[0]*grad_u[0] +
715 fluid_outward_normal[1]*grad_u[1] +
716 fluid_outward_normal[2]*grad_u[2]) +
717 C_surf*rel_vel_norm*(
u-u_s)*rho) + dV*
gf_s.ImH(eps_s, phi_s)*C_vol*(
u-u_s)*rho;
718 force_y = dV*D_s*(p*fluid_outward_normal[1] - porosity*mu*(fluid_outward_normal[0]*grad_v[0] +
719 fluid_outward_normal[1]*grad_v[1] +
720 fluid_outward_normal[2]*grad_v[2]) +
721 C_surf*rel_vel_norm*(
v-v_s)*rho) + dV*
gf_s.ImH(eps_s, phi_s)*C_vol*(
v-v_s)*rho;
722 force_z = dV*D_s*(p*fluid_outward_normal[2] - porosity*mu*(fluid_outward_normal[0]*grad_w[0] +
723 fluid_outward_normal[1]*grad_w[1] +
724 fluid_outward_normal[2]*grad_w[2]) +
725 C_surf*rel_vel_norm*(
v-v_s)*rho) + dV*
gf_s.ImH(eps_s, phi_s)*C_vol*(
v-v_s)*rho;
734 particle_netForces[i*3+0] += force_x;
735 particle_netForces[i*3+1] += force_y;
736 particle_netForces[i*3+2] += force_z;
737 particle_netMoments[i*3+0] += (r_y*force_z - r_z*force_y);
738 particle_netMoments[i*3+1] += (r_z*force_x - r_x*force_z);
739 particle_netMoments[i*3+2] += (r_x*force_y - r_y*force_x);
743 mom_u_source += C*(
u - u_s);
744 mom_v_source += C*(
v - v_s);
745 mom_w_source += C*(
w - w_s);
747 dmom_u_source[0] += C;
748 dmom_v_source[1] += C;
749 dmom_w_source[2] += C;
752 mom_u_ham -= D_s*porosity*nu*(fluid_outward_normal[0]*grad_u[0] + fluid_outward_normal[1]*grad_u[1] + fluid_outward_normal[2]*grad_u[2]);
753 dmom_u_ham_grad_u[0] -= D_s*porosity*nu*fluid_outward_normal[0];
754 dmom_u_ham_grad_u[1] -= D_s*porosity*nu*fluid_outward_normal[1];
755 dmom_u_ham_grad_u[2] -= D_s*porosity*nu*fluid_outward_normal[2];
757 mom_v_ham -= D_s*porosity*nu*(fluid_outward_normal[0]*grad_v[0] + fluid_outward_normal[1]*grad_v[1]+ fluid_outward_normal[2]*grad_v[2]);
758 dmom_v_ham_grad_v[0] -= D_s*porosity*nu*fluid_outward_normal[0];
759 dmom_v_ham_grad_v[1] -= D_s*porosity*nu*fluid_outward_normal[1];
760 dmom_v_ham_grad_v[2] -= D_s*porosity*nu*fluid_outward_normal[2];
762 mom_u_adv[0] += D_s*porosity*nu*fluid_outward_normal[0]*(
u - u_s);
763 mom_u_adv[1] += D_s*porosity*nu*fluid_outward_normal[1]*(
u - u_s);
764 mom_u_adv[2] += D_s*porosity*nu*fluid_outward_normal[2]*(
u - u_s);
765 dmom_u_adv_u[0] += D_s*porosity*nu*fluid_outward_normal[0];
766 dmom_u_adv_u[1] += D_s*porosity*nu*fluid_outward_normal[1];
767 dmom_u_adv_u[2] += D_s*porosity*nu*fluid_outward_normal[2];
769 mom_v_adv[0] += D_s*porosity*nu*fluid_outward_normal[0]*(
v - v_s);
770 mom_v_adv[1] += D_s*porosity*nu*fluid_outward_normal[1]*(
v - v_s);
771 mom_v_adv[2] += D_s*porosity*nu*fluid_outward_normal[2]*(
v - v_s);
772 dmom_v_adv_v[0] += D_s*porosity*nu*fluid_outward_normal[0];
773 dmom_v_adv_v[1] += D_s*porosity*nu*fluid_outward_normal[1];
774 dmom_v_adv_v[2] += D_s*porosity*nu*fluid_outward_normal[2];
776 mom_w_adv[0] += D_s*porosity*nu*fluid_outward_normal[0]*(
w - w_s);
777 mom_w_adv[1] += D_s*porosity*nu*fluid_outward_normal[1]*(
w - w_s);
778 mom_w_adv[2] += D_s*porosity*nu*fluid_outward_normal[2]*(
w - w_s);
779 dmom_w_adv_w[0] += D_s*porosity*nu*fluid_outward_normal[0];
780 dmom_w_adv_w[1] += D_s*porosity*nu*fluid_outward_normal[1];
781 dmom_w_adv_w[2] += D_s*porosity*nu*fluid_outward_normal[2];
787 const double& elementDiameter,
789 const double df[nSpace],
792 double h,density,nrm_df=0.0;
793 h = hFactor*elementDiameter;
795 for(
int I=0;I<nSpace;I++)
797 nrm_df = sqrt(nrm_df);
798 if (density > 1.0e-8)
799 cfl = nrm_df/(h*density);
806 const double eps_rho,
815 const double porosity,
816 const double eddy_visc_coef_0,
817 const double turb_var_0,
818 const double turb_var_1,
819 const double turb_grad_0[nSpace],
820 double &eddy_viscosity,
821 double mom_uu_diff_ten[nSpace],
822 double mom_vv_diff_ten[nSpace],
823 double mom_ww_diff_ten[nSpace],
824 double mom_uv_diff_ten[1],
825 double mom_uw_diff_ten[1],
826 double mom_vu_diff_ten[1],
827 double mom_vw_diff_ten[1],
828 double mom_wu_diff_ten[1],
829 double mom_wv_diff_ten[1],
830 double &mom_u_source,
831 double &mom_v_source,
832 double &mom_w_source)
840 assert (turbulenceClosureModel >=3);
841 double rho,nu,H_mu,ImH_mu,nu_t=0.0,nu_t_keps =0.0, nu_t_komega=0.0;
842 double isKEpsilon = 1.0;
843 if (turbulenceClosureModel == 4)
845 H_mu = (1.0-useVF)*
gf.H(eps_mu,
phi)+useVF*fmin(1.0,fmax(0.0,vf));
846 ImH_mu = (1.0-useVF)*
gf.ImH(eps_mu,
phi) + useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
850 const double twoThirds = 2.0/3.0;
const double div_zero = 1.0e-2*fmin(
nu_0,
nu_1);
851 mom_u_source += twoThirds*turb_grad_0[0];
852 mom_v_source += twoThirds*turb_grad_0[1];
853 mom_w_source += twoThirds*turb_grad_0[2];
857 nu_t_keps = eddy_visc_coef_0*turb_var_0*turb_var_0/(fabs(turb_var_1) + div_zero);
859 nu_t_komega = turb_var_0/(fabs(turb_var_1) + div_zero);
861 nu_t = isKEpsilon*nu_t_keps + (1.0-isKEpsilon)*nu_t_komega;
868 nu_t = fmax(nu_t,1.0e-4*nu);
870 nu_t = fmin(nu_t,1.0e6*nu);
872 eddy_viscosity = nu_t*rho;
874 mom_uu_diff_ten[0] += porosity*2.0*eddy_viscosity;
875 mom_uu_diff_ten[1] += porosity*eddy_viscosity;
876 mom_uu_diff_ten[2] += porosity*eddy_viscosity;
878 mom_uv_diff_ten[0]+=porosity*eddy_viscosity;
880 mom_uw_diff_ten[0]+=porosity*eddy_viscosity;
883 mom_vv_diff_ten[0] += porosity*eddy_viscosity;
884 mom_vv_diff_ten[1] += porosity*2.0*eddy_viscosity;
885 mom_vv_diff_ten[2] += porosity*eddy_viscosity;
887 mom_vu_diff_ten[0]+=porosity*eddy_viscosity;
889 mom_vw_diff_ten[0]+=porosity*eddy_viscosity;
892 mom_ww_diff_ten[0] += porosity*eddy_viscosity;
893 mom_ww_diff_ten[1] += porosity*eddy_viscosity;
894 mom_ww_diff_ten[2] += porosity*2.0*eddy_viscosity;
896 mom_wu_diff_ten[0]+=porosity*eddy_viscosity;
898 mom_wv_diff_ten[0]+=eddy_viscosity;
902 const double &elementDiameter,
905 const double df[nSpace],
912 double h, oneByAbsdt, density, viscosity, nrm_df;
913 h = hFactor * elementDiameter;
917 for (
int I = 0; I < nSpace; I++)
918 nrm_df +=
df[I] *
df[I];
919 nrm_df = sqrt(nrm_df);
920 if (density > 1.0e-8)
921 cfl = nrm_df/(h*density);
924 oneByAbsdt = fabs(dmt);
925 tau_v = 1.0/(4.0*viscosity/(h*h) + 2.0*nrm_df/h + oneByAbsdt);
926 tau_p = (4.0*viscosity + 2.0*nrm_df*h + oneByAbsdt*h*h)/pfac;
930 const double &Cd_sge,
931 const double G[nSpace * nSpace],
932 const double &G_dd_G,
935 const double Ai[nSpace],
943 for (
int I = 0; I < nSpace; I++)
944 for (
int J = 0; J < nSpace; J++)
945 v_d_Gv += Ai[I] * G[I * nSpace + J] * Ai[J];
946 tau_v = 1.0 / sqrt(Ct_sge * A0 * A0 + v_d_Gv + Cd_sge * Kij * Kij * G_dd_G + 1.0e-12);
947 tau_p = 1.0 / (pfac * tr_G * tau_v);
952 const double &pdeResidualP,
953 const double &pdeResidualU,
954 const double &pdeResidualV,
955 const double &pdeResidualW,
956 double &subgridErrorP,
957 double &subgridErrorU,
958 double &subgridErrorV,
959 double &subgridErrorW)
962 subgridErrorP = -tau_p * pdeResidualP;
964 subgridErrorU = -tau_v*pdeResidualU;
965 subgridErrorV = -tau_v*pdeResidualV;
966 subgridErrorW = -tau_v*pdeResidualW;
971 const double dpdeResidualP_du[nDOF_trial_element],
972 const double dpdeResidualP_dv[nDOF_trial_element],
973 const double dpdeResidualP_dw[nDOF_trial_element],
974 const double dpdeResidualU_dp[nDOF_trial_element],
975 const double dpdeResidualU_du[nDOF_trial_element],
976 const double dpdeResidualV_dp[nDOF_trial_element],
977 const double dpdeResidualV_dv[nDOF_trial_element],
978 const double dpdeResidualW_dp[nDOF_trial_element],
979 const double dpdeResidualW_dw[nDOF_trial_element],
980 double dsubgridErrorP_du[nDOF_trial_element],
981 double dsubgridErrorP_dv[nDOF_trial_element],
982 double dsubgridErrorP_dw[nDOF_trial_element],
983 double dsubgridErrorU_dp[nDOF_trial_element],
984 double dsubgridErrorU_du[nDOF_trial_element],
985 double dsubgridErrorV_dp[nDOF_trial_element],
986 double dsubgridErrorV_dv[nDOF_trial_element],
987 double dsubgridErrorW_dp[nDOF_trial_element],
988 double dsubgridErrorW_dw[nDOF_trial_element])
990 for (
int j = 0; j < nDOF_trial_element; j++)
993 dsubgridErrorP_du[j] = -tau_p*dpdeResidualP_du[j];
994 dsubgridErrorP_dv[j] = -tau_p*dpdeResidualP_dv[j];
995 dsubgridErrorP_dw[j] = -tau_p*dpdeResidualP_dw[j];
998 dsubgridErrorU_dp[j] = -tau_v * dpdeResidualU_dp[j];
999 dsubgridErrorU_du[j] = -tau_v * dpdeResidualU_du[j];
1001 dsubgridErrorV_dp[j] = -tau_v*dpdeResidualV_dp[j];
1002 dsubgridErrorV_dv[j] = -tau_v*dpdeResidualV_dv[j];
1004 dsubgridErrorW_dp[j] = -tau_v*dpdeResidualW_dp[j];
1005 dsubgridErrorW_dw[j] = -tau_v*dpdeResidualW_dw[j];
1011 const int& isDOFBoundary_u,
1012 const int& isDOFBoundary_v,
1013 const int& isDOFBoundary_w,
1014 const int& isFluxBoundary_p,
1015 const int& isFluxBoundary_u,
1016 const int& isFluxBoundary_v,
1017 const int& isFluxBoundary_w,
1018 const double& oneByRho,
1019 const double& bc_oneByRho,
1020 const double n[nSpace],
1021 const double& porosity,
1026 const double bc_f_mass[nSpace],
1027 const double bc_f_umom[nSpace],
1028 const double bc_f_vmom[nSpace],
1029 const double bc_f_wmom[nSpace],
1030 const double& bc_flux_mass,
1031 const double& bc_flux_umom,
1032 const double& bc_flux_vmom,
1033 const double& bc_flux_wmom,
1038 const double f_mass[nSpace],
1039 const double f_umom[nSpace],
1040 const double f_vmom[nSpace],
1041 const double f_wmom[nSpace],
1042 const double df_mass_du[nSpace],
1043 const double df_mass_dv[nSpace],
1044 const double df_mass_dw[nSpace],
1045 const double df_umom_dp[nSpace],
1046 const double df_umom_du[nSpace],
1047 const double df_umom_dv[nSpace],
1048 const double df_umom_dw[nSpace],
1049 const double df_vmom_dp[nSpace],
1050 const double df_vmom_du[nSpace],
1051 const double df_vmom_dv[nSpace],
1052 const double df_vmom_dw[nSpace],
1053 const double df_wmom_dp[nSpace],
1054 const double df_wmom_du[nSpace],
1055 const double df_wmom_dv[nSpace],
1056 const double df_wmom_dw[nSpace],
1061 double* velocity_star,
1064 double flowSpeedNormal;
1069 flowSpeedNormal=porosity*(
n[0]*velocity_star[0] +
1070 n[1]*velocity_star[1] +
1071 n[2]*velocity_star[2]);
1075 if (isDOFBoundary_u != 1)
1077 flux_mass +=
n[0]*f_mass[0];
1078 if (flowSpeedNormal < 0.0)
1080 flux_umom+=flowSpeedNormal*(0.0 -
u);
1085 flux_mass +=
n[0]*f_mass[0];
1086 if (flowSpeedNormal < 0.0)
1088 flux_umom+=flowSpeedNormal*(bc_u -
u);
1092 if (isDOFBoundary_v != 1)
1094 flux_mass+=
n[1]*f_mass[1];
1095 if (flowSpeedNormal < 0.0)
1097 flux_vmom+=flowSpeedNormal*(0.0 -
v);
1102 flux_mass+=
n[1]*f_mass[1];
1103 if (flowSpeedNormal < 0.0)
1105 flux_vmom+=flowSpeedNormal*(bc_v -
v);
1109 if (isDOFBoundary_w != 1)
1111 flux_mass+=
n[2]*f_mass[2];
1112 if (flowSpeedNormal < 0.0)
1114 flux_wmom+=flowSpeedNormal*(0.0 -
w);
1119 flux_mass +=
n[2]*f_mass[2];
1120 if (flowSpeedNormal < 0.0)
1122 flux_wmom+=flowSpeedNormal*(bc_w -
w);
1139 if (isFluxBoundary_u == 1)
1141 flux_umom = bc_flux_umom;
1144 if (isFluxBoundary_v == 1)
1146 flux_vmom = bc_flux_vmom;
1149 if (isFluxBoundary_w == 1)
1151 flux_wmom = bc_flux_wmom;
1158 const int& isDOFBoundary_u,
1159 const int& isDOFBoundary_v,
1160 const int& isDOFBoundary_w,
1161 const int& isFluxBoundary_p,
1162 const int& isFluxBoundary_u,
1163 const int& isFluxBoundary_v,
1164 const int& isFluxBoundary_w,
1165 const double& oneByRho,
1166 const double n[nSpace],
1167 const double& porosity,
1172 const double bc_f_mass[nSpace],
1173 const double bc_f_umom[nSpace],
1174 const double bc_f_vmom[nSpace],
1175 const double bc_f_wmom[nSpace],
1176 const double& bc_flux_mass,
1177 const double& bc_flux_umom,
1178 const double& bc_flux_vmom,
1179 const double& bc_flux_wmom,
1184 const double f_mass[nSpace],
1185 const double f_umom[nSpace],
1186 const double f_vmom[nSpace],
1187 const double f_wmom[nSpace],
1188 const double df_mass_du[nSpace],
1189 const double df_mass_dv[nSpace],
1190 const double df_mass_dw[nSpace],
1191 const double df_umom_dp[nSpace],
1192 const double df_umom_du[nSpace],
1193 const double df_umom_dv[nSpace],
1194 const double df_umom_dw[nSpace],
1195 const double df_vmom_dp[nSpace],
1196 const double df_vmom_du[nSpace],
1197 const double df_vmom_dv[nSpace],
1198 const double df_vmom_dw[nSpace],
1199 const double df_wmom_dp[nSpace],
1200 const double df_wmom_du[nSpace],
1201 const double df_wmom_dv[nSpace],
1202 const double df_wmom_dw[nSpace],
1203 double& dflux_mass_du,
1204 double& dflux_mass_dv,
1205 double& dflux_mass_dw,
1206 double& dflux_umom_dp,
1207 double& dflux_umom_du,
1208 double& dflux_umom_dv,
1209 double& dflux_umom_dw,
1210 double& dflux_vmom_dp,
1211 double& dflux_vmom_du,
1212 double& dflux_vmom_dv,
1213 double& dflux_vmom_dw,
1214 double& dflux_wmom_dp,
1215 double& dflux_wmom_du,
1216 double& dflux_wmom_dv,
1217 double& dflux_wmom_dw,
1218 double* velocity_star)
1220 double flowSpeedNormal;
1221 dflux_mass_du = 0.0;
1222 dflux_mass_dv = 0.0;
1223 dflux_mass_dw = 0.0;
1225 dflux_umom_dp = 0.0;
1226 dflux_umom_du = 0.0;
1227 dflux_umom_dv = 0.0;
1228 dflux_umom_dw = 0.0;
1230 dflux_vmom_dp = 0.0;
1231 dflux_vmom_du = 0.0;
1232 dflux_vmom_dv = 0.0;
1233 dflux_vmom_dw = 0.0;
1235 dflux_wmom_dp = 0.0;
1236 dflux_wmom_du = 0.0;
1237 dflux_wmom_dv = 0.0;
1238 dflux_wmom_dw = 0.0;
1239 flowSpeedNormal=porosity*(
n[0]*velocity_star[0] +
1240 n[1]*velocity_star[1] +
1241 n[2]*velocity_star[2]);
1242 if (isDOFBoundary_u != 1)
1244 dflux_mass_du +=
n[0]*df_mass_du[0];
1245 if (flowSpeedNormal < 0.0)
1246 dflux_umom_du -= flowSpeedNormal;
1250 dflux_mass_du +=
n[0]*df_mass_du[0];
1251 if (flowSpeedNormal < 0.0)
1252 dflux_umom_du -= flowSpeedNormal;
1254 if (isDOFBoundary_v != 1)
1256 dflux_mass_dv +=
n[1]*df_mass_dv[1];
1257 if (flowSpeedNormal < 0.0)
1258 dflux_vmom_dv -= flowSpeedNormal;
1262 dflux_mass_dv +=
n[1]*df_mass_dv[1];
1263 if (flowSpeedNormal < 0.0)
1264 dflux_vmom_dv -= flowSpeedNormal;
1266 if (isDOFBoundary_w != 1)
1268 dflux_mass_dw+=
n[2]*df_mass_dw[2];
1269 if (flowSpeedNormal < 0.0)
1270 dflux_wmom_dw -= flowSpeedNormal;
1274 dflux_mass_dw +=
n[2]*df_mass_dw[2];
1275 if (flowSpeedNormal < 0.0)
1276 dflux_wmom_dw -= flowSpeedNormal;
1290 if (isFluxBoundary_u == 1)
1292 dflux_umom_dp = 0.0;
1293 dflux_umom_du = 0.0;
1294 dflux_umom_dv = 0.0;
1295 dflux_umom_dw = 0.0;
1297 if (isFluxBoundary_v == 1)
1299 dflux_vmom_dp = 0.0;
1300 dflux_vmom_du = 0.0;
1301 dflux_vmom_dv = 0.0;
1302 dflux_vmom_dw = 0.0;
1304 if (isFluxBoundary_w == 1)
1306 dflux_wmom_dp = 0.0;
1307 dflux_wmom_du = 0.0;
1308 dflux_wmom_dv = 0.0;
1309 dflux_wmom_dw = 0.0;
1318 const int& isDOFBoundary,
1319 const int& isFluxBoundary,
1320 const double n[nSpace],
1323 const double& bc_flux,
1325 const double grad_potential[nSpace],
1327 const double& penalty,
1330 double diffusiveVelocityComponent_I,penaltyFlux,max_a;
1331 if(isFluxBoundary == 1)
1335 else if(isDOFBoundary == 1)
1339 for(
int I=0;I<nSpace;I++)
1341 diffusiveVelocityComponent_I=0.0;
1342 for(
int m=rowptr[I];m<rowptr[I+1];m++)
1344 diffusiveVelocityComponent_I -= a[m]*grad_potential[colind[m]];
1345 max_a = fmax(max_a,a[m]);
1347 flux+= diffusiveVelocityComponent_I*
n[I];
1349 penaltyFlux = max_a*penalty*(
u-bc_u);
1350 flux += penaltyFlux;
1356 std::cerr<<
"RANS3PF: warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0"<<std::endl;
1366 const int& isDOFBoundary,
1367 const int& isFluxBoundary,
1368 const double n[nSpace],
1371 const double grad_v[nSpace],
1372 const double& penalty)
1374 double dvel_I,tmp=0.0,max_a=0.0;
1375 if(isFluxBoundary==0 && isDOFBoundary==1)
1377 for(
int I=0;I<nSpace;I++)
1380 for(
int m=rowptr[I];m<rowptr[I+1];m++)
1382 dvel_I -= a[m]*grad_v[colind[m]];
1383 max_a = fmax(max_a,a[m]);
1387 tmp +=max_a*penalty*
v;
1399 res[0] = grad_u[0]*
n[0]+grad_u[1]*
n[1]+grad_u[2]*
n[2];
1400 res[1] = grad_v[0]*
n[0]+grad_v[1]*
n[1]+grad_v[2]*
n[2];
1401 res[2] = grad_w[0]*
n[0]+grad_w[1]*
n[1]+grad_w[2]*
n[2];;
1404 void get_stress_in_n(
const double *grad_u,
const double *grad_v,
const double *grad_w,
const double *
n,
double p,
double mu,
double f[3])
1406 f[0] = 2.0*grad_u[0]*
n[0]+(grad_u[1]+grad_v[0])*
n[1]+(grad_u[2]+grad_w[0])*
n[2];
1407 f[1] = (grad_v[0]+grad_u[1])*
n[0]+ 2*grad_v[1]*
n[1]+(grad_v[2]+grad_w[1])*
n[2];
1408 f[2] = (grad_w[0]+grad_u[2])*
n[0]+(grad_w[1]+grad_v[2])*
n[1]+ 2*grad_w[2]*
n[2];
1420 res[0] =
u[1]*
v[2]-
u[2]*
v[1];
1421 res[1] =
u[2]*
v[0]-
u[0]*
v[2];
1422 res[2] =
u[0]*
v[1]-
u[1]*
v[0];
1426 return u[0]*
v[0]+
u[1]*
v[1]+
u[2]*
v[2];
1428 int get_distance_to_ball(
int n_balls,
double* ball_center,
double* ball_radius,
double x,
double y,
double z,
double& distance)
1433 for (
int i=0; i<n_balls; ++i)
1435 d_ball_i = std::sqrt((ball_center[i*3+0]-x)*(ball_center[i*3+0]-x)
1436 +(ball_center[i*3+1]-y)*(ball_center[i*3+1]-y)
1438 +(ball_center[i*3+2]-
z)*(ball_center[i*3+2]-
z)
1441 if(d_ball_i<distance)
1443 distance = d_ball_i;
1451 double x,
double y,
double z,
1454 distance = std::sqrt((ball_center[I*3+0]-x)*(ball_center[I*3+0]-x)
1455 + (ball_center[I*3+1]-y)*(ball_center[I*3+1]-y)
1457 + (ball_center[I*3+2]-
z)*(ball_center[I*3+2]-
z)
1463 double x,
double y,
double z,
1464 double& nx,
double& ny,
double& nz)
1466 double distance = std::sqrt((ball_center[I*3+0]-x)*(ball_center[I*3+0]-x)
1467 + (ball_center[I*3+1]-y)*(ball_center[I*3+1]-y)
1469 + (ball_center[I*3+2]-
z)*(ball_center[I*3+2]-
z)
1472 nx = (x - ball_center[I*3+0])/(distance+1e-10);
1473 ny = (y - ball_center[I*3+1])/(distance+1e-10);
1474 #ifdef USE_CYLINDER_AS_PARTICLE
1477 nz = (
z - ball_center[I*3+2])/(distance+1e-10);
1481 double* ball_velocity,
double* ball_angular_velocity,
1483 double x,
double y,
double z,
1484 double&
vx,
double&
vy,
double& vz)
1486 #ifdef USE_CYLINDER_AS_PARTICLE
1487 double position[3]={x-ball_center[3*I + 0],y-ball_center[3*I + 1],0.0};
1489 double position[3]={x-ball_center[3*I + 0],y-ball_center[3*I + 1],
z-ball_center[3*I + 2]};
1491 double angular_cross_position[3];
1492 get_cross_product(&ball_angular_velocity[3*I + 0],position,angular_cross_position);
1493 vx = ball_velocity[3*I + 0] + angular_cross_position[0];
1494 vy = ball_velocity[3*I + 1] + angular_cross_position[1];
1495 vz = ball_velocity[3*I + 2] + angular_cross_position[2];
1501 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1502 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1503 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1504 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1505 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1506 double PSTAB = args.
scalar<
double>(
"PSTAB");
1507 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1508 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
1509 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1510 int nDOF_per_element_pressure = args.
scalar<
int>(
"nDOF_per_element_pressure");
1511 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
1512 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
1513 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
1514 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
1515 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
1516 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
1517 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
1518 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
1519 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
1520 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
1521 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
1522 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
1523 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
1524 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1525 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1526 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1527 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
1528 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
1529 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
1530 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
1531 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
1532 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
1533 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
1534 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
1535 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1536 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1537 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
1538 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1539 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1540 double hFactor = args.
scalar<
double>(
"hFactor");
1541 int nElements_global = args.
scalar<
int>(
"nElements_global");
1542 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
1543 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
1544 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
1545 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
1546 double useRBLES = args.
scalar<
double>(
"useRBLES");
1547 double useMetrics = args.
scalar<
double>(
"useMetrics");
1548 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1549 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
1550 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
1551 double sigma = args.
scalar<
double>(
"sigma");
1556 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
1557 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
1558 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
1559 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
1560 double C_dc = args.
scalar<
double>(
"C_dc");
1561 double C_b = args.
scalar<
double>(
"C_b");
1562 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
1563 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
1564 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
1565 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
1566 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
1567 xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
1568 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
1569 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
1570 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1571 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
1572 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
1573 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
1574 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
1575 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
1576 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
1577 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
1578 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
1579 xt::pyarray<double>& q_eddy_viscosity = args.
array<
double>(
"q_eddy_viscosity");
1580 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
1581 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
1582 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
1583 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1584 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
1585 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
1586 xt::pyarray<double>& u_dof_old = args.
array<
double>(
"u_dof_old");
1587 xt::pyarray<double>& v_dof_old = args.
array<
double>(
"v_dof_old");
1588 xt::pyarray<double>& w_dof_old = args.
array<
double>(
"w_dof_old");
1589 xt::pyarray<double>& u_dof_old_old = args.
array<
double>(
"u_dof_old_old");
1590 xt::pyarray<double>& v_dof_old_old = args.
array<
double>(
"v_dof_old_old");
1591 xt::pyarray<double>& w_dof_old_old = args.
array<
double>(
"w_dof_old_old");
1592 xt::pyarray<double>& uStar_dof = args.
array<
double>(
"uStar_dof");
1593 xt::pyarray<double>& vStar_dof = args.
array<
double>(
"vStar_dof");
1594 xt::pyarray<double>& wStar_dof = args.
array<
double>(
"wStar_dof");
1595 xt::pyarray<double>& g = args.
array<
double>(
"g");
1596 const double useVF = args.
scalar<
double>(
"useVF");
1597 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
1598 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
1599 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
1600 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
1601 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
1602 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
1603 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
1604 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
1605 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
1606 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
1607 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
1608 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
1609 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
1610 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
1611 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
1612 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
1613 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
1614 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1615 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
1616 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
1617 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1618 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
1619 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
1620 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
1621 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
1622 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
1623 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
1624 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
1625 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
1626 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
1627 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
1628 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
1629 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
1630 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
1631 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
1632 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
1633 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
1634 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
1635 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
1636 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
1637 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
1638 int offset_p = args.
scalar<
int>(
"offset_p");
1639 int offset_u = args.
scalar<
int>(
"offset_u");
1640 int offset_v = args.
scalar<
int>(
"offset_v");
1641 int offset_w = args.
scalar<
int>(
"offset_w");
1642 int stride_p = args.
scalar<
int>(
"stride_p");
1643 int stride_u = args.
scalar<
int>(
"stride_u");
1644 int stride_v = args.
scalar<
int>(
"stride_v");
1645 int stride_w = args.
scalar<
int>(
"stride_w");
1646 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1647 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1648 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1649 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
1650 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1651 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1652 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
1653 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
1654 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
1655 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
1656 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
1657 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
1658 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
1659 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
1660 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
1661 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
1662 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1663 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
1664 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
1665 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
1666 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
1667 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
1668 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
1669 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
1670 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
1671 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
1672 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
1673 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
1674 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
1675 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
1676 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
1677 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1678 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
1679 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1680 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
1681 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
1682 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
1683 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
1684 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
1685 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
1686 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
1687 xt::pyarray<double>& q_grad_u = args.
array<
double>(
"q_grad_u");
1688 xt::pyarray<double>& q_grad_v = args.
array<
double>(
"q_grad_v");
1689 xt::pyarray<double>& q_grad_w = args.
array<
double>(
"q_grad_w");
1690 xt::pyarray<double>& q_divU = args.
array<
double>(
"q_divU");
1691 xt::pyarray<double>& ebqe_grad_u = args.
array<
double>(
"ebqe_grad_u");
1692 xt::pyarray<double>& ebqe_grad_v = args.
array<
double>(
"ebqe_grad_v");
1693 xt::pyarray<double>& ebqe_grad_w = args.
array<
double>(
"ebqe_grad_w");
1694 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
1695 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
1696 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
1697 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
1698 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
1699 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
1700 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
1701 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
1702 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
1703 xt::pyarray<double>& q_rho = args.
array<
double>(
"q_rho");
1704 xt::pyarray<double>& ebqe_rho = args.
array<
double>(
"ebqe_rho");
1705 xt::pyarray<double>& q_nu = args.
array<
double>(
"q_nu");
1706 xt::pyarray<double>& ebqe_nu = args.
array<
double>(
"ebqe_nu");
1707 int nParticles = args.
scalar<
int>(
"nParticles");
1708 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
1709 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
1710 double particle_beta = args.
scalar<
double>(
"particle_beta");
1711 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
1712 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
1713 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
1714 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
1715 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
1716 xt::pyarray<double>& particle_netForces = args.
array<
double>(
"particle_netForces");
1717 xt::pyarray<double>& particle_netMoments = args.
array<
double>(
"particle_netMoments");
1718 xt::pyarray<double>& particle_surfaceArea = args.
array<
double>(
"particle_surfaceArea");
1719 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
1720 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
1721 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
1722 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
1723 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
1724 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
1725 xt::pyarray<double>& phisError = args.
array<
double>(
"phisError");
1726 xt::pyarray<double>& phisErrorNodal = args.
array<
double>(
"phisErrorNodal");
1727 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
1728 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
1730 double cE = args.
scalar<
double>(
"cE");
1731 int MULTIPLY_EXTERNAL_FORCE_BY_DENSITY = args.
scalar<
int>(
"MULTIPLY_EXTERNAL_FORCE_BY_DENSITY");
1732 xt::pyarray<double>& forcex = args.
array<
double>(
"forcex");
1733 xt::pyarray<double>& forcey = args.
array<
double>(
"forcey");
1734 xt::pyarray<double>& forcez = args.
array<
double>(
"forcez");
1735 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
1736 double dt = args.
scalar<
double>(
"dt");
1737 xt::pyarray<double>& quantDOFs = args.
array<
double>(
"quantDOFs");
1738 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
1739 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
1740 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
1741 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
1742 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
1743 double order_polynomial = args.
scalar<
double>(
"order_polynomial");
1744 xt::pyarray<double>& isActiveDOF = args.
array<
double>(
"isActiveDOF");
1745 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
1746 xt::pyarray<double>& ncDrag = args.
array<
double>(
"ncDrag");
1747 xt::pyarray<double>& betaDrag = args.
array<
double>(
"betaDrag");
1748 xt::pyarray<double>& vos_vel_nodes = args.
array<
double>(
"vos_vel_nodes");
1749 xt::pyarray<double>& entropyResidualPerNode = args.
array<
double>(
"entropyResidualPerNode");
1750 xt::pyarray<double>& laggedEntropyResidualPerNode = args.
array<
double>(
"laggedEntropyResidualPerNode");
1751 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
1752 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
1753 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
1754 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
1755 int NNZ_1D = args.
scalar<
int>(
"NNZ_1D");
1756 xt::pyarray<int>& csrRowIndeces_1D = args.
array<
int>(
"csrRowIndeces_1D");
1757 xt::pyarray<int>& csrColumnOffsets_1D = args.
array<
int>(
"csrColumnOffsets_1D");
1758 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
1759 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
1760 xt::pyarray<double>& isBoundary_1D = args.
array<
double>(
"isBoundary_1D");
1761 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
1762 gf.useExact=useExact;
1763 gf_s.useExact=useExact;
1764 register double element_uStar_He[nElements_global], element_vStar_He[nElements_global], element_wStar_He[nElements_global];
1768 den_hi.resize(numDOFs_1D,0.0);
1781 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1783 for (
int i=0; i<NNZ_1D; i++)
1785 uStar_dMatrix[i]=0.;
1786 vStar_dMatrix[i]=0.;
1787 wStar_dMatrix[i]=0.;
1791 for (
int i=0; i<numDOFs_1D; i++)
1796 entropyResidualPerNode[i]=0.;
1807 double mesh_volume_conservation=0.0,
1808 mesh_volume_conservation_weak=0.0,
1809 mesh_volume_conservation_err_max=0.0,
1810 mesh_volume_conservation_err_max_weak=0.0;
1811 double globalConservationError=0.0;
1812 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
1813 for(
int eN=0;eN<nElements_global;eN++)
1815 register double elementTransport[nDOF_test_element][nDOF_trial_element];
1816 register double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1818 register double elementResidual_p[nDOF_test_element],elementResidual_mesh[nDOF_test_element],
1819 elementResidual_u[nDOF_test_element],
1820 elementResidual_v[nDOF_test_element],
1821 mom_u_source_i[nDOF_test_element],
1822 mom_v_source_i[nDOF_test_element],
1823 mom_w_source_i[nDOF_test_element],
1824 betaDrag_i[nDOF_test_element],
1825 vos_i[nDOF_test_element],
1826 phisErrorElement[nDOF_test_element],
1827 elementResidual_w[nDOF_test_element],
1828 elementEntropyResidual[nDOF_test_element],
1830 double element_active=1.0;
1831 double mesh_volume_conservation_element=0.0,
1832 mesh_volume_conservation_element_weak=0.0;
1834 double linVisc_eN = 0, nlinVisc_eN_num = 0, nlinVisc_eN_den = 0;
1836 double det_hess_uStar_Ke=0.0, det_hess_vStar_Ke=0.0, det_hess_wStar_Ke=0.0, area_Ke=0.0;
1837 for (
int i=0;i<nDOF_test_element;i++)
1839 int eN_i = eN*nDOF_test_element+i;
1840 elementResidual_p_save[eN_i]=0.0;
1841 elementResidual_mesh[i]=0.0;
1842 elementResidual_p[i]=0.0;
1843 elementResidual_u[i]=0.0;
1844 elementResidual_v[i]=0.0;
1845 mom_u_source_i[i]=0.0;
1846 mom_v_source_i[i]=0.0;
1847 mom_w_source_i[i]=0.0;
1850 phisErrorElement[i]=0.0;
1851 elementResidual_w[i]=0.0;
1852 elementEntropyResidual[i]=0.0;
1853 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1855 for (
int j=0;j<nDOF_trial_element;j++)
1857 elementTransport[i][j]=0.0;
1858 elementTransposeTransport[i][j]=0.0;
1863 if(use_ball_as_particle==1)
1865 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1867 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1868 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1869 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1870 phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]]);
1879 double _distance[nDOF_mesh_trial_element]={0.0};
1881 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1883 if(use_ball_as_particle==1)
1886 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1887 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1888 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1893 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
1895 if ( _distance[I] >= 0)
1898 if (pos_counter == 3)
1902 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1904 if (_distance[I] < 0)
1907 assert(opp_node >=0);
1908 assert(opp_node <nDOF_mesh_trial_element);
1913 const int ebN = elementBoundariesArray[eN*nDOF_mesh_trial_element+opp_node];
1914 const int eN_oppo = (eN == elementBoundaryElementsArray[ebN*2+0])?elementBoundaryElementsArray[ebN*2+1]:elementBoundaryElementsArray[ebN*2+0];
1915 if((mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]<nNodes_owned
1916 || mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]<nNodes_owned
1917 || mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]<nNodes_owned
1924 if (eN == elementBoundaryElementsArray[ebN*2+0])
1936 double distance=1e10, distance_to_ith_particle;
1937 if(use_ball_as_particle==1)
1939 double middle_point_coord[3]={0.0};
1940 double middle_point_distance;
1941 middle_point_coord[0] = (mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]+0]
1942 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]+0]
1943 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]+0])/3.0;
1944 middle_point_coord[1] = (mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]+1]
1945 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]+1]
1946 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]+1])/3.0;
1947 middle_point_coord[2] = (mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]+2]
1948 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]+2]
1949 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]+2])/3.0;
1951 middle_point_coord[0],middle_point_coord[1],middle_point_coord[2],
1952 middle_point_distance);
1956 for (
int i=0;i<nParticles;++i)
1958 distance_to_ith_particle=particle_signed_distances[i*nElements_global*nQuadraturePoints_element
1959 +eN*nQuadraturePoints_element
1961 if (distance_to_ith_particle<distance)
1963 distance = distance_to_ith_particle;
1973 if(ebN<nElementBoundaries_owned)
1975 assert(eN_oppo==-1);
1979 else if (pos_counter == 4)
1982 for (
int i=0;i<nDOF_test_element;i++)
1984 isActiveDOF[offset_u+stride_u*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
1985 isActiveDOF[offset_v+stride_v*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
1986 isActiveDOF[offset_w+stride_w*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
1994 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
1995 for (
int j=0;j<nDOF_mesh_trial_element;j++)
1997 register int eN_j = eN*nDOF_mesh_trial_element+j;
1998 element_phi[j] = phi_dof[p_l2g[eN_j]];
1999 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
2001 double element_nodes[nDOF_mesh_trial_element*3];
2002 for (
int i=0;i<nDOF_mesh_trial_element;i++)
2004 register int eN_i=eN*nDOF_mesh_trial_element+i;
2005 for(
int I=0;I<3;I++)
2006 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
2008 gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),
false);
2009 gf.calculate(element_phi, element_nodes, x_ref.data(),
false);
2013 for(
int k=0;k<nQuadraturePoints_element;k++)
2018 register int eN_k = eN*nQuadraturePoints_element+k,
2019 eN_k_nSpace = eN_k*nSpace,
2021 eN_nDOF_trial_element = eN*nDOF_trial_element;
2022 register double p=0.0,
u=0.0,
v=0.0,
w=0.0,un=0.0,vn=0.0,wn=0.0,
2023 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
2032 dmass_adv_u[nSpace],
2033 dmass_adv_v[nSpace],
2034 dmass_adv_w[nSpace],
2036 dmom_u_adv_u[nSpace],
2037 dmom_u_adv_v[nSpace],
2038 dmom_u_adv_w[nSpace],
2040 dmom_v_adv_u[nSpace],
2041 dmom_v_adv_v[nSpace],
2042 dmom_v_adv_w[nSpace],
2044 dmom_w_adv_u[nSpace],
2045 dmom_w_adv_v[nSpace],
2046 dmom_w_adv_w[nSpace],
2047 mom_uu_diff_ten[nSpace],
2048 mom_vv_diff_ten[nSpace],
2049 mom_ww_diff_ten[nSpace],
2060 dmom_u_ham_grad_p[nSpace],
2061 dmom_u_ham_grad_u[nSpace],
2063 dmom_v_ham_grad_p[nSpace],
2064 dmom_v_ham_grad_v[nSpace],
2066 dmom_w_ham_grad_p[nSpace],
2067 dmom_w_ham_grad_w[nSpace],
2078 Lstar_u_p[nDOF_test_element],
2079 Lstar_v_p[nDOF_test_element],
2080 Lstar_w_p[nDOF_test_element],
2081 Lstar_u_u[nDOF_test_element],
2082 Lstar_v_v[nDOF_test_element],
2083 Lstar_w_w[nDOF_test_element],
2084 Lstar_p_u[nDOF_test_element],
2085 Lstar_p_v[nDOF_test_element],
2086 Lstar_p_w[nDOF_test_element],
2091 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
2092 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
2095 jacInv[nSpace*nSpace],
2096 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
2097 vel_hess_trial[nDOF_trial_element*
nSpace2],
2098 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
2099 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
2100 u_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2101 v_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2102 w_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2108 dmom_u_source[nSpace],
2109 dmom_v_source[nSpace],
2110 dmom_w_source[nSpace],
2114 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
2116 ck.calculateMapping_element(eN,
2120 mesh_trial_ref.data(),
2121 mesh_grad_trial_ref.data(),
2126 ck.calculateH_element(eN,
2128 nodeDiametersArray.data(),
2130 mesh_trial_ref.data(),
2132 ck.calculateMappingVelocity_element(eN,
2134 mesh_velocity_dof.data(),
2136 mesh_trial_ref.data(),
2141 dV = fabs(jacDet)*dV_ref[k];
2142 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2145 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2146 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2147 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2151 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
2152 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
2157 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
2158 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
2159 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
w);
2161 ck.valFromDOF(u_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],un);
2162 ck.valFromDOF(v_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],vn);
2163 ck.valFromDOF(w_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],wn);
2166 for (
int I=0;I<nSpace;I++)
2167 grad_p[I] = q_grad_p[eN_k_nSpace + I];
2168 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
2169 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
2170 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_w);
2171 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
2172 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
2173 ck.hessFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_w);
2174 ck.hessFromDOF(uStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_uStar);
2175 ck.hessFromDOF(vStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_vStar);
2176 ck.hessFromDOF(wStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_wStar);
2178 for (
int j=0;j<nDOF_trial_element;j++)
2181 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
2182 for (
int I=0;I<nSpace;I++)
2185 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
2186 if (ARTIFICIAL_VISCOSITY==4)
2189 u_times_vel_grad_test_dV[j*nSpace+I] =
2190 u*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_u[I];
2191 v_times_vel_grad_test_dV[j*nSpace+I] =
2192 v*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_v[I];
2193 w_times_vel_grad_test_dV[j*nSpace+I] =
2194 w*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_w[I];
2199 if (ARTIFICIAL_VISCOSITY==3)
2201 det_hess_uStar_Ke +=
2202 (hess_uStar[0]*(hess_uStar[4]*hess_uStar[8] - hess_uStar[7]*hess_uStar[5])
2203 -hess_uStar[1]*(hess_uStar[3]*hess_uStar[8] - hess_uStar[6]*hess_uStar[5])
2204 +hess_uStar[2]*(hess_uStar[3]*hess_uStar[7] - hess_uStar[6]*hess_uStar[4]))*dV;
2205 det_hess_vStar_Ke +=
2206 (hess_vStar[0]*(hess_vStar[4]*hess_vStar[8] - hess_vStar[7]*hess_vStar[5])
2207 -hess_vStar[1]*(hess_vStar[3]*hess_vStar[8] - hess_vStar[6]*hess_vStar[5])
2208 +hess_vStar[2]*(hess_vStar[3]*hess_vStar[7] - hess_vStar[6]*hess_vStar[4]))*dV;
2209 det_hess_wStar_Ke +=
2210 (hess_wStar[0]*(hess_wStar[4]*hess_wStar[8] - hess_wStar[7]*hess_wStar[5])
2211 -hess_wStar[1]*(hess_wStar[3]*hess_wStar[8] - hess_wStar[6]*hess_wStar[5])
2212 +hess_wStar[2]*(hess_wStar[3]*hess_wStar[7] - hess_wStar[6]*hess_wStar[4]))*dV;
2216 double div_mesh_velocity=0.0;
2217 int NDOF_MESH_TRIAL_ELEMENT=4;
2218 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
2220 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
2221 div_mesh_velocity +=
2222 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*nSpace+0] +
2223 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*nSpace+1] +
2224 mesh_velocity_dof[mesh_l2g[eN_j]*3+2]*vel_grad_trial[j*nSpace+2];
2226 mesh_volume_conservation_element += (alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity)*dV;
2227 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
2229 porosity = 1.0 - q_vos[eN_k];
2238 double distance_to_omega_solid = 1e10;
2239 if(use_ball_as_particle==1)
2243 distance_to_omega_solid);
2247 for (
int i = 0; i < nParticles; i++)
2249 double distance_to_i_th_solid = particle_signed_distances[i * nElements_global * nQuadraturePoints_element + eN_k];
2250 distance_to_omega_solid = (distance_to_i_th_solid < distance_to_omega_solid)?distance_to_i_th_solid:distance_to_omega_solid;
2253 phi_solid[eN_k] = distance_to_omega_solid;
2265 elementDiameter[eN],
2266 smagorinskyConstant,
2267 turbulenceClosureModel,
2272 &normal_phi[eN_k_nSpace],
2273 distance_to_omega_solid,
2286 q_velocity_sge[eN_k_nSpace+0],
2287 q_velocity_sge[eN_k_nSpace+1],
2288 q_velocity_sge[eN_k_nSpace+2],
2289 q_eddy_viscosity[eN_k],
2336 MULTIPLY_EXTERNAL_FORCE_BY_DENSITY,
2340 MATERIAL_PARAMETERS_AS_FUNCTION,
2341 density_as_function[eN_k],
2342 dynamic_viscosity_as_function[eN_k],
2345 use_ball_as_particle,
2348 ball_velocity.data(),
2349 ball_angular_velocity.data(),
2350 INT_BY_PARTS_PRESSURE);
2353 mass_source = q_mass_source[eN_k];
2354 for (
int I=0;I<nSpace;I++)
2356 dmom_u_source[I] = 0.0;
2357 dmom_v_source[I] = 0.0;
2358 dmom_w_source[I] = 0.0;
2369 q_eddy_viscosity[eN_k],
2376 q_velocity_sge[eN_k_nSpace+0],
2377 q_velocity_sge[eN_k_nSpace+1],
2378 q_velocity_sge[eN_k_nSpace+2],
2379 eps_solid[elementFlags[eN]],
2381 q_velocity_solid[eN_k_nSpace+0],
2382 q_velocity_solid[eN_k_nSpace+1],
2383 q_velocity_solid[eN_k_nSpace+2],
2384 q_velocityStar_solid[eN_k_nSpace+0],
2385 q_velocityStar_solid[eN_k_nSpace+1],
2386 q_velocityStar_solid[eN_k_nSpace+2],
2393 q_grad_vos[eN_k_nSpace+0],
2394 q_grad_vos[eN_k_nSpace+1],
2395 q_grad_vos[eN_k_nSpace+2]);
2396 double C_particles = 0.0;
2397 if (nParticles > 0 && USE_SBM==0)
2402 nQuadraturePoints_global,
2403 &particle_signed_distances[eN_k],
2404 &particle_signed_distance_normals[eN_k_nSpace],
2405 &particle_velocities[eN_k_nSpace],
2406 particle_centroids.data(),
2407 use_ball_as_particle,
2410 ball_velocity.data(),
2411 ball_angular_velocity.data(),
2413 particle_penalty_constant/h_phi,
2414 particle_alpha/h_phi,
2415 particle_beta/h_phi,
2432 q_velocity_sge[eN_k_nSpace + 0],
2433 q_velocity_sge[eN_k_nSpace + 1],
2434 q_velocity_sge[eN_k_nSpace + 2],
2457 particle_netForces.data(),
2458 particle_netMoments.data(),
2459 particle_surfaceArea.data());
2461 if (turbulenceClosureModel >= 3)
2463 const double c_mu = 0.09;
2478 &q_turb_var_grad_0[eN_k_nSpace],
2479 q_eddy_viscosity[eN_k],
2497 q_mom_u_acc[eN_k] = mom_u_acc;
2498 q_mom_v_acc[eN_k] = mom_v_acc;
2499 q_mom_w_acc[eN_k] = mom_w_acc;
2501 q_mass_adv[eN_k_nSpace+0] =
u;
2502 q_mass_adv[eN_k_nSpace+1] =
v;
2503 q_mass_adv[eN_k_nSpace+2] =
w;
2507 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
2508 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
2509 mom_u_adv[2] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*zt;
2510 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
2511 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
2512 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
2514 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
2515 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
2516 mom_v_adv[2] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*zt;
2517 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
2518 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
2519 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
2521 mom_w_adv[0] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*
xt;
2522 mom_w_adv[1] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*yt;
2523 mom_w_adv[2] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*zt;
2524 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
2525 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
2526 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
2531 if (q_dV_last[eN_k] <= -100)
2532 q_dV_last[eN_k] = dV;
2535 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2541 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2547 q_mom_w_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2553 mom_u_acc_t *= dmom_u_acc_u;
2554 mom_v_acc_t *= dmom_v_acc_v;
2555 mom_w_acc_t *= dmom_w_acc_w;
2563 ck.Mass_strong(-q_dvos_dt[eN_k]) +
2564 ck.Advection_strong(dmass_adv_u,grad_u) +
2565 ck.Advection_strong(dmass_adv_v,grad_v) +
2566 ck.Advection_strong(dmass_adv_w,grad_w) +
2567 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
2569 ck.Reaction_strong(mass_source);
2572 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
2573 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
2574 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
2577 ck.Mass_strong(mom_u_acc_t) +
2578 ck.Advection_strong(dmom_adv_sge,grad_u) +
2579 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
2580 ck.Reaction_strong(mom_u_source) -
2581 ck.Reaction_strong(
u*div_mesh_velocity);
2584 ck.Mass_strong(mom_v_acc_t) +
2585 ck.Advection_strong(dmom_adv_sge,grad_v) +
2586 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
2587 ck.Reaction_strong(mom_v_source) -
2588 ck.Reaction_strong(
v*div_mesh_velocity);
2591 ck.Mass_strong(mom_w_acc_t) +
2592 ck.Advection_strong(dmom_adv_sge,grad_w) +
2593 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
2594 ck.Reaction_strong(mom_w_source) -
2595 ck.Reaction_strong(
w*div_mesh_velocity);
2599 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
2601 elementDiameter[eN],
2606 dmom_u_ham_grad_p[0],
2616 dmom_u_ham_grad_p[0],
2621 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
2622 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
2635 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
2636 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
2637 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
2639 mom_u_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge[eN_k_nSpace+0]);
2640 mom_u_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge[eN_k_nSpace+1]);
2641 mom_u_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge[eN_k_nSpace+2]);
2644 for (
int i=0;i<nDOF_test_element;i++)
2646 register int i_nSpace = i*nSpace;
2651 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2652 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2653 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2654 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2655 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2656 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2659 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
2660 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
2661 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
2665 if (ARTIFICIAL_VISCOSITY==0 || ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2667 q_numDiff_u[eN_k] = 0;
2668 q_numDiff_v[eN_k] = 0;
2669 q_numDiff_w[eN_k] = 0;
2671 else if (ARTIFICIAL_VISCOSITY==1)
2673 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v + pdeResidual_w*pdeResidual_w);
2674 q_numDiff_u[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
2675 (1.0-useMetrics)*hFactor*hFactor*elementDiameter[eN]*elementDiameter[eN]);
2676 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
2677 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
2681 double rho = q_rho[eN_k];
2682 double mu = q_rho[eN_k]*q_nu[eN_k];
2684 double vel2 =
u*
u +
v*
v +
w*
w;
2688 porosity*rho*((
u-un)/dt + (
u*grad_u[0]+
v*grad_u[1]+
w*grad_u[2]) - g[0])
2689 + (KILL_PRESSURE_TERM == 1 ? 0 : 1.)*grad_p[0]
2690 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcex[eN_k]
2691 - mu*(hess_u[0] + hess_u[4] + hess_u[8])
2692 - mu*(hess_u[0] + hess_v[1] + hess_w[2]);
2694 porosity*rho*((
v-vn)/dt + (
u*grad_v[0]+
v*grad_v[1]+
w*grad_v[2]) - g[1])
2695 + (KILL_PRESSURE_TERM == 1 ? 0 : 1.)*grad_p[1]
2696 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcey[eN_k]
2697 - mu*(hess_v[0] + hess_v[4] + hess_v[8])
2698 - mu*(hess_u[1] + hess_v[4] + hess_w[5]);
2700 porosity*rho*((
w-wn)/dt + (
u*grad_w[0]+
v*grad_w[1]+
w*grad_w[2]) - g[2])
2701 + (KILL_PRESSURE_TERM == 1 ? 0 : 1.)*grad_p[2]
2702 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcez[eN_k]
2703 - mu*(hess_w[0] + hess_w[4] + hess_w[8])
2704 - mu*(hess_u[2] + hess_v[5] + hess_w[8]);
2707 double entRes_times_u = Res_in_x*
u + Res_in_y*
v + Res_in_z*
w;
2709 double hK = elementDiameter[eN]/order_polynomial;
2710 q_numDiff_u[eN_k] = fmin(
cMax*porosity*rho*hK*std::sqrt(vel2),
2711 cE*hK*hK*fabs(entRes_times_u)/(vel2+1E-10));
2712 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
2713 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
2717 linVisc_eN = fmax(porosity*rho*std::sqrt(vel2),linVisc_eN);
2718 nlinVisc_eN_num = fmax(fabs(entRes_times_u),nlinVisc_eN_num);
2719 nlinVisc_eN_den = fmax(vel2,nlinVisc_eN_den);
2730 q_velocity[eN_k_nSpace+0]=
u;
2731 q_velocity[eN_k_nSpace+1]=
v;
2732 q_velocity[eN_k_nSpace+2]=
w;
2733 for (
int I=0;I<nSpace;I++)
2735 q_grad_u[eN_k_nSpace+I] = grad_u[I];
2736 q_grad_v[eN_k_nSpace+I] = grad_v[I];
2737 q_grad_w[eN_k_nSpace+I] = grad_w[I];
2740 q_divU[eN_k] = q_grad_u[eN_k_nSpace+0] + q_grad_v[eN_k_nSpace+1] + q_grad_w[eN_k_nSpace+2];
2743 double unit_normal[nSpace];
2744 double norm_grad_phi = 0.;
2745 for (
int I=0;I<nSpace;I++)
2746 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
2747 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
2748 for (
int I=0;I<nSpace;I++)
2749 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
2753 v1[0]=1.-unit_normal[0]*unit_normal[0];
2754 v1[1]=-unit_normal[0]*unit_normal[1];
2755 v1[2]=-unit_normal[0]*unit_normal[2];
2758 v2[0]=-unit_normal[0]*unit_normal[1];
2759 v2[1]=1.-unit_normal[1]*unit_normal[1];
2760 v2[2]=-unit_normal[1]*unit_normal[2];
2763 v3[0]=-unit_normal[0]*unit_normal[2];
2764 v3[1]=-unit_normal[1]*unit_normal[2];
2765 v3[2]=1.-unit_normal[2]*unit_normal[2];
2766 double delta =
gf.D(eps_mu,
phi[eN_k]);
2767 register double vel_tgrad_test_i[nSpace],
2768 tgrad_u[nSpace], tgrad_v[nSpace], tgrad_w[nSpace];
2780 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2782 velStar[0] = q_velocity_sge[eN_k_nSpace+0];
2783 velStar[1] = q_velocity_sge[eN_k_nSpace+1];
2784 velStar[2] = q_velocity_sge[eN_k_nSpace+2];
2786 for(
int i=0;i<nDOF_test_element;i++)
2788 register int i_nSpace=i*nSpace;
2790 &vel_grad_trial[i_nSpace],
2792 phisErrorElement[i]+=std::abs(phisError[eN_k_nSpace+0])*p_test_dV[i];
2810 elementResidual_u[i] +=
2811 ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
2812 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
2813 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]) +
2814 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]) +
2815 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]) +
2816 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
2817 ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
2818 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+0] : 0.) +
2820 USE_SUPG*
ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
2821 ck.NumericalDiffusion(q_numDiff_u_last[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]) +
2823 ck.NumericalDiffusion(delta*sigma*dV,v1,vel_tgrad_test_i) +
2824 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_u,vel_tgrad_test_i);
2825 mom_u_source_i[i] +=
ck.Reaction_weak(mom_u_source,vel_test_dV[i]);
2826 betaDrag_i[i] +=
ck.Reaction_weak(dmom_u_source[0],
2828 vos_i[i] +=
ck.Reaction_weak(1.0-porosity,
2831 elementResidual_v[i] +=
2832 ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
2833 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
2834 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]) +
2835 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]) +
2836 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]) +
2837 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
2838 ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
2839 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+1] : 0.) +
2841 USE_SUPG*
ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
2842 ck.NumericalDiffusion(q_numDiff_v_last[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]) +
2844 ck.NumericalDiffusion(delta*sigma*dV,v2,vel_tgrad_test_i) +
2845 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_v,vel_tgrad_test_i);
2846 mom_v_source_i[i] +=
ck.Reaction_weak(mom_v_source,vel_test_dV[i]);
2848 elementResidual_w[i] +=
2849 ck.Mass_weak(mom_w_acc_t,vel_test_dV[i]) +
2850 ck.Advection_weak(mom_w_adv,&vel_grad_test_dV[i_nSpace]) +
2851 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]) +
2852 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]) +
2853 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]) +
2854 ck.Reaction_weak(mom_w_source,vel_test_dV[i]) +
2855 ck.Hamiltonian_weak(mom_w_ham,vel_test_dV[i]) +
2856 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+2] : 0.) +
2858 USE_SUPG*
ck.SubgridError(subgridError_w,Lstar_w_w[i]) +
2859 ck.NumericalDiffusion(q_numDiff_w_last[eN_k],grad_w,&vel_grad_test_dV[i_nSpace]) +
2861 ck.NumericalDiffusion(delta*sigma*dV,v3,vel_tgrad_test_i) +
2862 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_w,vel_tgrad_test_i);
2863 mom_w_source_i[i] +=
ck.Reaction_weak(mom_w_source,vel_test_dV[i]);
2865 if (ARTIFICIAL_VISCOSITY==4)
2869 elementEntropyResidual[i] +=
2871 ck.Mass_weak(mom_u_acc_t,
u*vel_test_dV[i]) +
2872 ck.Advection_weak(mom_u_adv,&u_times_vel_grad_test_dV[i_nSpace])+
2873 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),
2874 sdInfo_u_u_colind.data(),
2877 &u_times_vel_grad_test_dV[i_nSpace]) +
2878 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),
2879 sdInfo_u_v_colind.data(),
2882 &u_times_vel_grad_test_dV[i_nSpace]) +
2883 ck.Diffusion_weak(sdInfo_u_w_rowptr.data(),
2884 sdInfo_u_w_colind.data(),
2887 &u_times_vel_grad_test_dV[i_nSpace]) +
2888 ck.Reaction_weak(mom_u_source,
u*vel_test_dV[i]) +
2889 ck.Hamiltonian_weak(mom_u_ham,
u*vel_test_dV[i])
2891 ck.Mass_weak(mom_v_acc_t,
v*vel_test_dV[i]) +
2892 ck.Advection_weak(mom_v_adv,&v_times_vel_grad_test_dV[i_nSpace])+
2893 ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),
2894 sdInfo_v_u_colind.data(),
2897 &v_times_vel_grad_test_dV[i_nSpace])+
2898 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),
2899 sdInfo_v_v_colind.data(),
2902 &v_times_vel_grad_test_dV[i_nSpace])+
2903 ck.Diffusion_weak(sdInfo_v_w_rowptr.data(),
2904 sdInfo_v_w_colind.data(),
2907 &v_times_vel_grad_test_dV[i_nSpace]) +
2908 ck.Reaction_weak(mom_v_source,
v*vel_test_dV[i]) +
2909 ck.Hamiltonian_weak(mom_v_ham,
v*vel_test_dV[i])
2911 ck.Mass_weak(mom_w_acc_t,
w*vel_test_dV[i]) +
2912 ck.Advection_weak(mom_w_adv,&w_times_vel_grad_test_dV[i_nSpace])+
2913 ck.Diffusion_weak(sdInfo_w_u_rowptr.data(),
2914 sdInfo_w_u_colind.data(),
2917 &w_times_vel_grad_test_dV[i_nSpace]) +
2918 ck.Diffusion_weak(sdInfo_w_v_rowptr.data(),
2919 sdInfo_w_v_colind.data(),
2922 &w_times_vel_grad_test_dV[i_nSpace]) +
2923 ck.Diffusion_weak(sdInfo_w_w_rowptr.data(),
2924 sdInfo_w_w_colind.data(),
2927 &w_times_vel_grad_test_dV[i_nSpace]) +
2928 ck.Reaction_weak(mom_w_source,
w*vel_test_dV[i]) +
2929 ck.Hamiltonian_weak(mom_w_ham,
w*vel_test_dV[i]);
2931 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2933 for(
int j=0;j<nDOF_trial_element;j++)
2935 int j_nSpace = j*nSpace;
2936 int i_nSpace = i*nSpace;
2937 elementTransport[i][j] +=
2938 q_rho[eN_k]*porosity*
2939 ck.AdvectionJacobian_strong(velStar,
2940 &vel_grad_test_dV[j_nSpace])
2941 *vel_trial_ref[k*nDOF_trial_element+i];
2942 elementTransposeTransport[i][j] +=
2943 q_rho[eN_k]*porosity*
2944 ck.AdvectionJacobian_strong(velStar,
2945 &vel_grad_test_dV[i_nSpace])
2946 *vel_trial_ref[k*nDOF_trial_element+j];
2951 element_uStar_He[eN] = det_hess_uStar_Ke/area_Ke;
2952 element_vStar_He[eN] = det_hess_vStar_Ke/area_Ke;
2953 element_wStar_He[eN] = det_hess_wStar_Ke/area_Ke;
2958 double hK = elementDiameter[eN];
2959 double artVisc = fmin(
cMax*hK*linVisc_eN,
2960 cE*hK*hK*nlinVisc_eN_num/(nlinVisc_eN_den+1E-10));
2961 for(
int k=0;k<nQuadraturePoints_element;k++)
2963 register int eN_k = eN*nQuadraturePoints_element+k;
2964 q_numDiff_u[eN_k] = artVisc;
2965 q_numDiff_v[eN_k] = artVisc;
2966 q_numDiff_w[eN_k] = artVisc;
2972 for(
int i=0;i<nDOF_test_element;i++)
2974 register int eN_i=eN*nDOF_test_element+i;
2975 phisErrorNodal[vel_l2g[eN_i]]+= element_active*phisErrorElement[i];
2979 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=element_active*elementResidual_u[i];
2980 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=element_active*elementResidual_v[i];
2981 globalResidual[offset_w+stride_w*vel_l2g[eN_i]]+=element_active*elementResidual_w[i];
2982 ncDrag[offset_u+stride_u*vel_l2g[eN_i]]+=mom_u_source_i[i];
2983 ncDrag[offset_v+stride_v*vel_l2g[eN_i]]+=mom_v_source_i[i];
2984 ncDrag[offset_w+stride_w*vel_l2g[eN_i]]+=mom_w_source_i[i];
2985 betaDrag[vel_l2g[eN_i]] += betaDrag_i[i];
2986 vos_vel_nodes[vel_l2g[eN_i]] += vos_i[i];
2989 if (ARTIFICIAL_VISCOSITY==3)
2991 uStar_hi[vel_l2g[eN_i]] += element_uStar_He[eN];
2992 vStar_hi[vel_l2g[eN_i]] += element_vStar_He[eN];
2993 wStar_hi[vel_l2g[eN_i]] += element_wStar_He[eN];
2994 den_hi[vel_l2g[eN_i]] += 1;
2996 if (ARTIFICIAL_VISCOSITY==4)
2999 entropyResidualPerNode[vel_l2g[eN_i]] += elementEntropyResidual[i];
3001 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3003 for (
int j=0;j<nDOF_trial_element;j++)
3005 int eN_i_j = eN_i*nDOF_trial_element+j;
3007 + csrColumnOffsets_1D[eN_i_j]]
3008 += elementTransport[i][j];
3011 + csrColumnOffsets_1D[eN_i_j]]
3012 += elementTransposeTransport[i][j];
3023 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3026 for (
int i=0; i<numDOFs_1D; i++)
3028 if (ARTIFICIAL_VISCOSITY==4)
3031 double max_u2i = (std::pow(u_dof[i],2.) +
3032 std::pow(v_dof[i],2.) +
3033 std::pow(w_dof[i],2.));
3034 double min_u2i = max_u2i;
3035 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3037 int j = colind_1D[offset];
3038 double u2j = (std::pow(u_dof[j],2.) +
3039 std::pow(v_dof[j],2.) +
3040 std::pow(w_dof[j],2.));
3041 max_u2i = fmax(max_u2i,u2j);
3042 min_u2i = fmin(min_u2i,u2j);
3044 double normi = 0.5*(max_u2i + min_u2i) + 1E-10;
3045 entropyResidualPerNode[i] = fabs(entropyResidualPerNode[i])/normi;
3050 double uStari = uStar_dof[i];
3051 double vStari = vStar_dof[i];
3052 double wStari = wStar_dof[i];
3054 double u_beta_numerator = 0., u_beta_denominator = 0.;
3055 double v_beta_numerator = 0., v_beta_denominator = 0.;
3056 double w_beta_numerator = 0., w_beta_denominator = 0.;
3059 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3061 int j = colind_1D[offset];
3062 double uStarj = uStar_dof[j];
3063 double vStarj = vStar_dof[j];
3064 double wStarj = wStar_dof[j];
3067 u_beta_numerator += (uStarj - uStari);
3068 u_beta_denominator += fabs(uStarj - uStari);
3070 v_beta_numerator += (vStarj - vStari);
3071 v_beta_denominator += fabs(vStarj - vStari);
3073 w_beta_numerator += (wStarj - wStari);
3074 w_beta_denominator += fabs(wStarj - wStari);
3076 double u_beta = fabs(u_beta_numerator)/(u_beta_denominator+1E-10);
3077 double v_beta = fabs(v_beta_numerator)/(v_beta_denominator+1E-10);
3078 double w_beta = fabs(w_beta_numerator)/(w_beta_denominator+1E-10);
3100 if (ARTIFICIAL_VISCOSITY==3)
3102 for(
int eN=0;eN<nElements_global;eN++)
3104 double uStar_He = element_uStar_He[eN];
3105 double vStar_He = element_vStar_He[eN];
3106 double wStar_He = element_wStar_He[eN];
3107 for(
int i=0;i<nDOF_test_element;i++)
3109 register int eN_i=eN*nDOF_test_element+i;
3110 register int gi = vel_l2g[eN_i];
3119 if (ARTIFICIAL_VISCOSITY==3)
3121 for (
int i=0; i<numDOFs_1D; i++)
3127 if (isBoundary_1D[i] == 1)
3155 for (
int i=0; i<numDOFs_1D; i++)
3158 double uStar_dii = 0;
3159 double vStar_dii = 0;
3160 double wStar_dii = 0;
3161 double ui = u_dof[i];
3162 double vi = v_dof[i];
3163 double wi = w_dof[i];
3165 double ith_u_dissipative_term = 0;
3166 double ith_v_dissipative_term = 0;
3167 double ith_w_dissipative_term = 0;
3173 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3175 int j = colind_1D[offset];
3178 double uj = u_dof[j];
3179 double vj = v_dof[j];
3180 double wj = w_dof[j];
3186 if (ARTIFICIAL_VISCOSITY==4)
3188 double dEVij = fmax(laggedEntropyResidualPerNode[i],
3189 laggedEntropyResidualPerNode[j]);
3192 uStar_dMatrix[ij] = fmin(dLij,
cE*dEVij);
3193 vStar_dMatrix[i] = uStar_dMatrix[ij];
3194 wStar_dMatrix[i] = uStar_dMatrix[ij];
3205 uStar_dii -= uStar_dMatrix[ij];
3206 vStar_dii -= vStar_dMatrix[ij];
3207 wStar_dii -= wStar_dMatrix[ij];
3209 ith_u_dissipative_term += uStar_dMatrix[ij]*(uj-ui);
3210 ith_v_dissipative_term += vStar_dMatrix[ij]*(vj-vi);
3211 ith_w_dissipative_term += wStar_dMatrix[ij]*(wj-wi);
3220 uStar_dMatrix[ii] = uStar_dii;
3221 vStar_dMatrix[ii] = vStar_dii;
3222 wStar_dMatrix[ii] = wStar_dii;
3223 globalResidual[offset_u+stride_u*i] += -ith_u_dissipative_term;
3224 globalResidual[offset_v+stride_v*i] += -ith_v_dissipative_term;
3225 globalResidual[offset_w+stride_w*i] += -ith_w_dissipative_term;
3234 std::memset(particle_netForces.data(),0,nParticles*3*
sizeof(
double));
3235 std::memset(particle_netMoments.data(),0,nParticles*3*
sizeof(
double));
3241 eN_nDOF_trial_element = eN*nDOF_trial_element;
3242 register double elementResidual_mesh[nDOF_test_element],
3243 elementResidual_p[nDOF_test_element],
3244 elementResidual_u[nDOF_test_element],
3245 elementResidual_v[nDOF_test_element],
3246 elementResidual_w[nDOF_test_element],
3251 for (
int i=0;i<nDOF_test_element;i++)
3253 elementResidual_mesh[i]=0.0;
3254 elementResidual_p[i]=0.0;
3255 elementResidual_u[i]=0.0;
3256 elementResidual_v[i]=0.0;
3257 elementResidual_w[i]=0.0;
3259 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3261 register int ebN_kb = ebN*nQuadraturePoints_elementBoundary+kb,
3263 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3264 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3265 register double u_ext=0.0, v_ext=0.0, w_ext=0.0,
3266 bc_u_ext=0.0, bc_v_ext=0.0, bc_w_ext=0.0,
3267 grad_u_ext[nSpace], grad_v_ext[nSpace], grad_w_ext[nSpace],
3268 jac_ext[nSpace*nSpace],
3270 jacInv_ext[nSpace*nSpace],
3271 boundaryJac[nSpace*(nSpace-1)],
3272 metricTensor[(nSpace-1)*(nSpace-1)],
3273 metricTensorDetSqrt,
3275 p_test_dS[nDOF_test_element],p_grad_trial_trace[nDOF_trial_element*nSpace],
3276 vel_test_dS[nDOF_test_element],
3277 vel_grad_trial_trace[nDOF_trial_element*nSpace],
3278 vel_grad_test_dS[nDOF_trial_element*nSpace],
3280 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3281 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3282 force_x,force_y,force_z,
3283 force_p_x,force_p_y,force_p_z,
3284 force_v_x,force_v_y,force_v_z,
3287 ck.calculateMapping_elementBoundary(eN,
3293 mesh_trial_trace_ref.data(),
3294 mesh_grad_trial_trace_ref.data(),
3295 boundaryJac_ref.data(),
3301 metricTensorDetSqrt,
3305 ck.calculateMappingVelocity_elementBoundary(eN,
3309 mesh_velocity_dof.data(),
3311 mesh_trial_trace_ref.data(),
3312 xt_ext,yt_ext,zt_ext,
3317 dS = metricTensorDetSqrt*dS_ref[kb];
3319 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3322 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3324 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3325 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3326 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],w_ext);
3328 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3329 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3330 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
3331 for (
int j=0;j<nDOF_trial_element;j++)
3333 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3334 for (
int I=0;I<nSpace;I++)
3335 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3337 ck.calculateGScale(G,normal,h_penalty);
3342 double distance[3], P_normal[3], P_tangent[3];
3343 if(use_ball_as_particle==1)
3351 P_normal[0],P_normal[1],P_normal[2]);
3353 ball_velocity.data(),ball_angular_velocity.data(),
3355 x_ext-dist*P_normal[0],
3356 y_ext-dist*P_normal[1],
3357 z_ext-dist*P_normal[2],
3358 bc_u_ext,bc_v_ext,bc_w_ext);
3363 dist = ebq_global_phi_solid[ebN_kb];
3364 P_normal[0] = ebq_global_grad_phi_solid[ebN_kb*nSpace+0];
3365 P_normal[1] = ebq_global_grad_phi_solid[ebN_kb*nSpace+1];
3366 P_normal[2] = ebq_global_grad_phi_solid[ebN_kb*nSpace+2];
3367 bc_u_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+0];
3368 bc_v_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+1];
3369 bc_w_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+2];
3371 distance[0] = -P_normal[0]*dist;
3372 distance[1] = -P_normal[1]*dist;
3373 distance[2] = -P_normal[2]*dist;
3374 assert(h_penalty>0.0);
3375 if (h_penalty < std::abs(dist))
3376 h_penalty = std::abs(dist);
3380 for(
int i=0;i<3;++i)
3387 double C_adim =
C_sbm*visco/h_penalty;
3388 double beta_adim =
beta_sbm*h_penalty*visco;
3394 const double u_m_uD[3] = {u_ext - bc_u_ext,v_ext - bc_v_ext,w_ext - bc_w_ext};
3395 const double zero_vec[3]={0.,0.,0.};
3399 for (
int i=0;i<nDOF_test_element;i++)
3401 int eN_i = eN*nDOF_test_element+i;
3403 int GlobPos_u = offset_u+stride_u*vel_l2g[eN_i];
3404 int GlobPos_v = offset_v+stride_v*vel_l2g[eN_i];
3405 int GlobPos_w = offset_w+stride_w*vel_l2g[eN_i];
3406 const double phi_i = vel_test_dS[i];
3407 double *grad_phi_i = &vel_grad_test_dS[i*nSpace+0];
3411 globalResidual[GlobPos_u] += C_adim*phi_i*u_m_uD[0];
3412 globalResidual[GlobPos_v] += C_adim*phi_i*u_m_uD[1];
3413 globalResidual[GlobPos_w] += C_adim*phi_i*u_m_uD[2];
3417 globalResidual[GlobPos_u] -= visco * phi_i*res[0];
3418 globalResidual[GlobPos_v] -= visco * phi_i*res[1];
3419 globalResidual[GlobPos_w] -= visco * phi_i*res[2];
3430 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*u_m_uD[0];
3431 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*u_m_uD[1];
3432 globalResidual[GlobPos_w] += C_adim*grad_phi_i_dot_d*u_m_uD[2];
3435 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*grad_u_d[0];
3436 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*grad_u_d[1];
3437 globalResidual[GlobPos_w] += C_adim*grad_phi_i_dot_d*grad_u_d[2];
3440 globalResidual[GlobPos_u] += C_adim*phi_i*grad_u_d[0];
3441 globalResidual[GlobPos_v] += C_adim*phi_i*grad_u_d[1];
3442 globalResidual[GlobPos_w] += C_adim*phi_i*grad_u_d[2];
3463 for (
int i=0; i<nDOF_per_element_pressure;++i)
3465 p_ext += p_dof[p_l2g[eN*nDOF_per_element_pressure+i]]*p_trial_trace_ref[ebN_local_kb*nDOF_per_element_pressure+i];
3468 double force_quad_pt[3]={0.0,0.0,0.0},torque_quad_pt[3]={0.0,0.0,0.0},position_vector_to_mass_center[3];
3470 get_stress_in_n(grad_u_ext,grad_v_ext,grad_w_ext,P_normal,p_ext,visco,force_quad_pt);
3472 force_quad_pt[0] *= dS;
3473 force_quad_pt[1] *= dS;
3474 force_quad_pt[2] *= dS;
3475 if(use_ball_as_particle==1)
3488 if(ebN < nElementBoundaries_owned)
3507 gf_s.useExact=
false;
3508 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3510 register int ebN = exteriorElementBoundariesArray[ebNE],
3511 eN = elementBoundaryElementsArray[ebN*2+0],
3512 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
3513 eN_nDOF_trial_element = eN*nDOF_trial_element;
3514 register double elementResidual_mesh[nDOF_test_element],
3515 elementResidual_p[nDOF_test_element],
3516 elementResidual_u[nDOF_test_element],
3517 elementResidual_v[nDOF_test_element],
3518 elementResidual_w[nDOF_test_element],
3520 for (
int i=0;i<nDOF_test_element;i++)
3522 elementResidual_mesh[i]=0.0;
3523 elementResidual_p[i]=0.0;
3524 elementResidual_u[i]=0.0;
3525 elementResidual_v[i]=0.0;
3526 elementResidual_w[i]=0.0;
3528 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3530 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3531 ebNE_kb_nSpace = ebNE_kb*nSpace,
3532 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3533 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3534 register double p_ext=0.0,
3543 dmom_u_acc_u_ext=0.0,
3545 dmom_v_acc_v_ext=0.0,
3547 dmom_w_acc_w_ext=0.0,
3548 mass_adv_ext[nSpace],
3549 dmass_adv_u_ext[nSpace],
3550 dmass_adv_v_ext[nSpace],
3551 dmass_adv_w_ext[nSpace],
3552 mom_u_adv_ext[nSpace],
3553 dmom_u_adv_u_ext[nSpace],
3554 dmom_u_adv_v_ext[nSpace],
3555 dmom_u_adv_w_ext[nSpace],
3556 mom_v_adv_ext[nSpace],
3557 dmom_v_adv_u_ext[nSpace],
3558 dmom_v_adv_v_ext[nSpace],
3559 dmom_v_adv_w_ext[nSpace],
3560 mom_w_adv_ext[nSpace],
3561 dmom_w_adv_u_ext[nSpace],
3562 dmom_w_adv_v_ext[nSpace],
3563 dmom_w_adv_w_ext[nSpace],
3564 mom_uu_diff_ten_ext[nSpace],
3565 mom_vv_diff_ten_ext[nSpace],
3566 mom_ww_diff_ten_ext[nSpace],
3567 mom_uv_diff_ten_ext[1],
3568 mom_uw_diff_ten_ext[1],
3569 mom_vu_diff_ten_ext[1],
3570 mom_vw_diff_ten_ext[1],
3571 mom_wu_diff_ten_ext[1],
3572 mom_wv_diff_ten_ext[1],
3573 mom_u_source_ext=0.0,
3574 mom_v_source_ext=0.0,
3575 mom_w_source_ext=0.0,
3577 dmom_u_ham_grad_p_ext[nSpace],
3578 dmom_u_ham_grad_u_ext[nSpace],
3580 dmom_v_ham_grad_p_ext[nSpace],
3581 dmom_v_ham_grad_v_ext[nSpace],
3583 dmom_w_ham_grad_p_ext[nSpace],
3584 dmom_w_ham_grad_w_ext[nSpace],
3585 dmom_u_adv_p_ext[nSpace],
3586 dmom_v_adv_p_ext[nSpace],
3587 dmom_w_adv_p_ext[nSpace],
3589 flux_mom_u_adv_ext=0.0,
3590 flux_mom_v_adv_ext=0.0,
3591 flux_mom_w_adv_ext=0.0,
3592 flux_mom_uu_diff_ext=0.0,
3593 flux_mom_uv_diff_ext=0.0,
3594 flux_mom_uw_diff_ext=0.0,
3595 flux_mom_vu_diff_ext=0.0,
3596 flux_mom_vv_diff_ext=0.0,
3597 flux_mom_vw_diff_ext=0.0,
3598 flux_mom_wu_diff_ext=0.0,
3599 flux_mom_wv_diff_ext=0.0,
3600 flux_mom_ww_diff_ext=0.0,
3605 bc_mom_u_acc_ext=0.0,
3606 bc_dmom_u_acc_u_ext=0.0,
3607 bc_mom_v_acc_ext=0.0,
3608 bc_dmom_v_acc_v_ext=0.0,
3609 bc_mom_w_acc_ext=0.0,
3610 bc_dmom_w_acc_w_ext=0.0,
3611 bc_mass_adv_ext[nSpace],
3612 bc_dmass_adv_u_ext[nSpace],
3613 bc_dmass_adv_v_ext[nSpace],
3614 bc_dmass_adv_w_ext[nSpace],
3615 bc_mom_u_adv_ext[nSpace],
3616 bc_dmom_u_adv_u_ext[nSpace],
3617 bc_dmom_u_adv_v_ext[nSpace],
3618 bc_dmom_u_adv_w_ext[nSpace],
3619 bc_mom_v_adv_ext[nSpace],
3620 bc_dmom_v_adv_u_ext[nSpace],
3621 bc_dmom_v_adv_v_ext[nSpace],
3622 bc_dmom_v_adv_w_ext[nSpace],
3623 bc_mom_w_adv_ext[nSpace],
3624 bc_dmom_w_adv_u_ext[nSpace],
3625 bc_dmom_w_adv_v_ext[nSpace],
3626 bc_dmom_w_adv_w_ext[nSpace],
3627 bc_mom_uu_diff_ten_ext[nSpace],
3628 bc_mom_vv_diff_ten_ext[nSpace],
3629 bc_mom_ww_diff_ten_ext[nSpace],
3630 bc_mom_uv_diff_ten_ext[1],
3631 bc_mom_uw_diff_ten_ext[1],
3632 bc_mom_vu_diff_ten_ext[1],
3633 bc_mom_vw_diff_ten_ext[1],
3634 bc_mom_wu_diff_ten_ext[1],
3635 bc_mom_wv_diff_ten_ext[1],
3636 bc_mom_u_source_ext=0.0,
3637 bc_mom_v_source_ext=0.0,
3638 bc_mom_w_source_ext=0.0,
3639 bc_mom_u_ham_ext=0.0,
3640 bc_dmom_u_ham_grad_p_ext[nSpace],
3641 bc_dmom_u_ham_grad_u_ext[nSpace],
3642 bc_mom_v_ham_ext=0.0,
3643 bc_dmom_v_ham_grad_p_ext[nSpace],
3644 bc_dmom_v_ham_grad_v_ext[nSpace],
3645 bc_mom_w_ham_ext=0.0,
3646 bc_dmom_w_ham_grad_p_ext[nSpace],
3647 bc_dmom_w_ham_grad_w_ext[nSpace],
3648 jac_ext[nSpace*nSpace],
3650 jacInv_ext[nSpace*nSpace],
3651 boundaryJac[nSpace*(nSpace-1)],
3652 metricTensor[(nSpace-1)*(nSpace-1)],
3653 metricTensorDetSqrt,
3654 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
3655 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
3656 vel_grad_test_dS[nDOF_trial_element*nSpace],
3657 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3661 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3662 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;
3664 ck.calculateMapping_elementBoundary(eN,
3670 mesh_trial_trace_ref.data(),
3671 mesh_grad_trial_trace_ref.data(),
3672 boundaryJac_ref.data(),
3678 metricTensorDetSqrt,
3682 ck.calculateMappingVelocity_elementBoundary(eN,
3686 mesh_velocity_dof.data(),
3688 mesh_trial_trace_ref.data(),
3689 xt_ext,yt_ext,zt_ext,
3701 dS = metricTensorDetSqrt*dS_ref[kb];
3704 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3705 ck.calculateGScale(G,&ebqe_normal_phi_ext[ebNE_kb_nSpace],h_phi);
3707 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3708 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3709 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3714 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3718 p_ext = ebqe_p[ebNE_kb];
3719 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3720 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3721 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],w_ext);
3723 for (
int I=0;I<nSpace;I++)
3724 grad_p_ext[I] = ebqe_grad_p[ebNE_kb_nSpace + I];
3725 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3726 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3727 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
3729 for (
int j=0;j<nDOF_trial_element;j++)
3732 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3733 for (
int I=0;I<nSpace;I++)
3734 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3736 bc_p_ext = isDOFBoundary_p[ebNE_kb]*ebqe_bc_p_ext[ebNE_kb]+(1-isDOFBoundary_p[ebNE_kb])*p_ext;
3738 bc_u_ext = isDOFBoundary_u[ebNE_kb]*(ebqe_bc_u_ext[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u[ebNE_kb])*u_ext;
3739 bc_v_ext = isDOFBoundary_v[ebNE_kb]*(ebqe_bc_v_ext[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v[ebNE_kb])*v_ext;
3740 bc_w_ext = isDOFBoundary_w[ebNE_kb]*(ebqe_bc_w_ext[ebNE_kb] + MOVING_DOMAIN*zt_ext) + (1-isDOFBoundary_w[ebNE_kb])*w_ext;
3742 porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
3746 double distance_to_omega_solid = 1e10;
3747 if (use_ball_as_particle == 1)
3749 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(), x_ext, y_ext, z_ext, distance_to_omega_solid);
3753 for (
int i = 0; i < nParticles; i++)
3755 double distance_to_i_th_solid = ebq_global_phi_solid[i * nElementBoundaries_global * nQuadraturePoints_elementBoundary + ebNE_kb];
3756 distance_to_omega_solid = (distance_to_i_th_solid < distance_to_omega_solid)?distance_to_i_th_solid:distance_to_omega_solid;
3759 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);
3768 elementDiameter[eN],
3769 smagorinskyConstant,
3770 turbulenceClosureModel,
3773 ebqe_vf_ext[ebNE_kb],
3774 ebqe_phi_ext[ebNE_kb],
3775 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
3776 distance_to_omega_solid,
3777 ebqe_kappa_phi_ext[ebNE_kb],
3789 ebqe_velocity_star[ebNE_kb_nSpace+0],
3790 ebqe_velocity_star[ebNE_kb_nSpace+1],
3791 ebqe_velocity_star[ebNE_kb_nSpace+2],
3815 mom_uu_diff_ten_ext,
3816 mom_vv_diff_ten_ext,
3817 mom_ww_diff_ten_ext,
3818 mom_uv_diff_ten_ext,
3819 mom_uw_diff_ten_ext,
3820 mom_vu_diff_ten_ext,
3821 mom_vw_diff_ten_ext,
3822 mom_wu_diff_ten_ext,
3823 mom_wv_diff_ten_ext,
3828 dmom_u_ham_grad_p_ext,
3829 dmom_u_ham_grad_u_ext,
3831 dmom_v_ham_grad_p_ext,
3832 dmom_v_ham_grad_v_ext,
3834 dmom_w_ham_grad_p_ext,
3835 dmom_w_ham_grad_w_ext,
3843 MATERIAL_PARAMETERS_AS_FUNCTION,
3844 ebqe_density_as_function[ebNE_kb],
3845 ebqe_dynamic_viscosity_as_function[ebNE_kb],
3848 use_ball_as_particle,
3851 ball_velocity.data(),
3852 ball_angular_velocity.data(),
3853 INT_BY_PARTS_PRESSURE);
3862 elementDiameter[eN],
3863 smagorinskyConstant,
3864 turbulenceClosureModel,
3867 bc_ebqe_vf_ext[ebNE_kb],
3868 bc_ebqe_phi_ext[ebNE_kb],
3869 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
3870 distance_to_omega_solid,
3871 ebqe_kappa_phi_ext[ebNE_kb],
3883 ebqe_velocity_star[ebNE_kb_nSpace+0],
3884 ebqe_velocity_star[ebNE_kb_nSpace+1],
3885 ebqe_velocity_star[ebNE_kb_nSpace+2],
3886 bc_eddy_viscosity_ext,
3888 bc_dmom_u_acc_u_ext,
3890 bc_dmom_v_acc_v_ext,
3892 bc_dmom_w_acc_w_ext,
3898 bc_dmom_u_adv_u_ext,
3899 bc_dmom_u_adv_v_ext,
3900 bc_dmom_u_adv_w_ext,
3902 bc_dmom_v_adv_u_ext,
3903 bc_dmom_v_adv_v_ext,
3904 bc_dmom_v_adv_w_ext,
3906 bc_dmom_w_adv_u_ext,
3907 bc_dmom_w_adv_v_ext,
3908 bc_dmom_w_adv_w_ext,
3909 bc_mom_uu_diff_ten_ext,
3910 bc_mom_vv_diff_ten_ext,
3911 bc_mom_ww_diff_ten_ext,
3912 bc_mom_uv_diff_ten_ext,
3913 bc_mom_uw_diff_ten_ext,
3914 bc_mom_vu_diff_ten_ext,
3915 bc_mom_vw_diff_ten_ext,
3916 bc_mom_wu_diff_ten_ext,
3917 bc_mom_wv_diff_ten_ext,
3918 bc_mom_u_source_ext,
3919 bc_mom_v_source_ext,
3920 bc_mom_w_source_ext,
3922 bc_dmom_u_ham_grad_p_ext,
3923 bc_dmom_u_ham_grad_u_ext,
3925 bc_dmom_v_ham_grad_p_ext,
3926 bc_dmom_v_ham_grad_v_ext,
3928 bc_dmom_w_ham_grad_p_ext,
3929 bc_dmom_w_ham_grad_w_ext,
3937 MATERIAL_PARAMETERS_AS_FUNCTION,
3938 ebqe_density_as_function[ebNE_kb],
3939 ebqe_dynamic_viscosity_as_function[ebNE_kb],
3942 use_ball_as_particle,
3945 ball_velocity.data(),
3946 ball_angular_velocity.data(),
3947 INT_BY_PARTS_PRESSURE);
3950 if (turbulenceClosureModel >= 3)
3952 const double turb_var_grad_0_dummy[3] = {0.,0.,0.};
3953 const double c_mu = 0.09;
3962 ebqe_vf_ext[ebNE_kb],
3963 ebqe_phi_ext[ebNE_kb],
3966 ebqe_turb_var_0[ebNE_kb],
3967 ebqe_turb_var_1[ebNE_kb],
3968 turb_var_grad_0_dummy,
3970 mom_uu_diff_ten_ext,
3971 mom_vv_diff_ten_ext,
3972 mom_ww_diff_ten_ext,
3973 mom_uv_diff_ten_ext,
3974 mom_uw_diff_ten_ext,
3975 mom_vu_diff_ten_ext,
3976 mom_vw_diff_ten_ext,
3977 mom_wu_diff_ten_ext,
3978 mom_wv_diff_ten_ext,
3991 bc_ebqe_vf_ext[ebNE_kb],
3992 bc_ebqe_phi_ext[ebNE_kb],
3995 ebqe_turb_var_0[ebNE_kb],
3996 ebqe_turb_var_1[ebNE_kb],
3997 turb_var_grad_0_dummy,
3998 bc_eddy_viscosity_ext,
3999 bc_mom_uu_diff_ten_ext,
4000 bc_mom_vv_diff_ten_ext,
4001 bc_mom_ww_diff_ten_ext,
4002 bc_mom_uv_diff_ten_ext,
4003 bc_mom_uw_diff_ten_ext,
4004 bc_mom_vu_diff_ten_ext,
4005 bc_mom_vw_diff_ten_ext,
4006 bc_mom_wu_diff_ten_ext,
4007 bc_mom_wv_diff_ten_ext,
4008 bc_mom_u_source_ext,
4009 bc_mom_v_source_ext,
4010 bc_mom_w_source_ext);
4017 mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
4018 mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
4019 mom_u_adv_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*zt_ext;
4020 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
4021 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
4022 dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext;
4024 mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
4025 mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
4026 mom_v_adv_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*zt_ext;
4027 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
4028 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
4029 dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext;
4031 mom_w_adv_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*xt_ext;
4032 mom_w_adv_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*yt_ext;
4033 mom_w_adv_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*zt_ext;
4034 dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext;
4035 dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext;
4036 dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext;
4040 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
4041 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
4042 bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*zt_ext;
4044 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
4045 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
4046 bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*zt_ext;
4048 bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*xt_ext;
4049 bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*yt_ext;
4050 bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*zt_ext;
4055 ck.calculateGScale(G,normal,h_penalty);
4056 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext[ebNE_kb];
4058 isDOFBoundary_u[ebNE_kb],
4059 isDOFBoundary_v[ebNE_kb],
4060 isDOFBoundary_w[ebNE_kb],
4061 isAdvectiveFluxBoundary_p[ebNE_kb],
4062 isAdvectiveFluxBoundary_u[ebNE_kb],
4063 isAdvectiveFluxBoundary_v[ebNE_kb],
4064 isAdvectiveFluxBoundary_w[ebNE_kb],
4065 dmom_u_ham_grad_p_ext[0],
4066 bc_dmom_u_ham_grad_p_ext[0],
4068 porosity_ext*ebqe_rho[ebNE_kb],
4077 ebqe_bc_flux_mass_ext[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]+zt_ext*normal[2]),
4078 ebqe_bc_flux_mom_u_adv_ext[ebNE_kb],
4079 ebqe_bc_flux_mom_v_adv_ext[ebNE_kb],
4080 ebqe_bc_flux_mom_w_adv_ext[ebNE_kb],
4108 &ebqe_velocity_star[ebNE_kb_nSpace],
4109 &ebqe_velocity[ebNE_kb_nSpace]);
4111 for (
int I=0;I<nSpace;I++)
4113 ebqe_grad_u[ebNE_kb_nSpace+I] = grad_u_ext[I];
4114 ebqe_grad_v[ebNE_kb_nSpace+I] = grad_v_ext[I];
4115 ebqe_grad_w[ebNE_kb_nSpace+I] = grad_w_ext[I];
4118 ebqe_phi_ext[ebNE_kb],
4119 sdInfo_u_u_rowptr.data(),
4120 sdInfo_u_u_colind.data(),
4121 isDOFBoundary_u[ebNE_kb],
4122 isDiffusiveFluxBoundary_u[ebNE_kb],
4124 bc_mom_uu_diff_ten_ext,
4126 ebqe_bc_flux_u_diff_ext[ebNE_kb],
4127 mom_uu_diff_ten_ext,
4131 flux_mom_uu_diff_ext);
4133 ebqe_phi_ext[ebNE_kb],
4134 sdInfo_u_v_rowptr.data(),
4135 sdInfo_u_v_colind.data(),
4136 isDOFBoundary_v[ebNE_kb],
4137 isDiffusiveFluxBoundary_v[ebNE_kb],
4139 bc_mom_uv_diff_ten_ext,
4142 mom_uv_diff_ten_ext,
4146 flux_mom_uv_diff_ext);
4148 ebqe_phi_ext[ebNE_kb],
4149 sdInfo_u_w_rowptr.data(),
4150 sdInfo_u_w_colind.data(),
4151 isDOFBoundary_w[ebNE_kb],
4152 isDiffusiveFluxBoundary_w[ebNE_kb],
4154 bc_mom_uw_diff_ten_ext,
4157 mom_uw_diff_ten_ext,
4161 flux_mom_uw_diff_ext);
4163 ebqe_phi_ext[ebNE_kb],
4164 sdInfo_v_u_rowptr.data(),
4165 sdInfo_v_u_colind.data(),
4166 isDOFBoundary_u[ebNE_kb],
4167 isDiffusiveFluxBoundary_u[ebNE_kb],
4169 bc_mom_vu_diff_ten_ext,
4172 mom_vu_diff_ten_ext,
4176 flux_mom_vu_diff_ext);
4178 ebqe_phi_ext[ebNE_kb],
4179 sdInfo_v_v_rowptr.data(),
4180 sdInfo_v_v_colind.data(),
4181 isDOFBoundary_v[ebNE_kb],
4182 isDiffusiveFluxBoundary_v[ebNE_kb],
4184 bc_mom_vv_diff_ten_ext,
4186 ebqe_bc_flux_v_diff_ext[ebNE_kb],
4187 mom_vv_diff_ten_ext,
4191 flux_mom_vv_diff_ext);
4193 ebqe_phi_ext[ebNE_kb],
4194 sdInfo_v_w_rowptr.data(),
4195 sdInfo_v_w_colind.data(),
4196 isDOFBoundary_w[ebNE_kb],
4197 isDiffusiveFluxBoundary_w[ebNE_kb],
4199 bc_mom_vw_diff_ten_ext,
4202 mom_vw_diff_ten_ext,
4206 flux_mom_vw_diff_ext);
4208 ebqe_phi_ext[ebNE_kb],
4209 sdInfo_w_u_rowptr.data(),
4210 sdInfo_w_u_colind.data(),
4211 isDOFBoundary_u[ebNE_kb],
4212 isDiffusiveFluxBoundary_u[ebNE_kb],
4214 bc_mom_wu_diff_ten_ext,
4217 mom_wu_diff_ten_ext,
4221 flux_mom_wu_diff_ext);
4223 ebqe_phi_ext[ebNE_kb],
4224 sdInfo_w_v_rowptr.data(),
4225 sdInfo_w_v_colind.data(),
4226 isDOFBoundary_v[ebNE_kb],
4227 isDiffusiveFluxBoundary_v[ebNE_kb],
4229 bc_mom_wv_diff_ten_ext,
4232 mom_wv_diff_ten_ext,
4236 flux_mom_wv_diff_ext);
4238 ebqe_phi_ext[ebNE_kb],
4239 sdInfo_w_w_rowptr.data(),
4240 sdInfo_w_w_colind.data(),
4241 isDOFBoundary_w[ebNE_kb],
4242 isDiffusiveFluxBoundary_w[ebNE_kb],
4244 bc_mom_ww_diff_ten_ext,
4246 ebqe_bc_flux_w_diff_ext[ebNE_kb],
4247 mom_ww_diff_ten_ext,
4251 flux_mom_ww_diff_ext);
4252 flux[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
4260 if (ebN < nElementBoundaries_owned)
4262 force_v_x = (flux_mom_u_adv_ext + flux_mom_uu_diff_ext + flux_mom_uv_diff_ext + flux_mom_uw_diff_ext)/dmom_u_ham_grad_p_ext[0];
4263 force_v_y = (flux_mom_v_adv_ext + flux_mom_vu_diff_ext + flux_mom_vv_diff_ext + flux_mom_vw_diff_ext)/dmom_u_ham_grad_p_ext[0];
4264 force_v_z = (flux_mom_w_adv_ext + flux_mom_wu_diff_ext + flux_mom_wv_diff_ext + flux_mom_ww_diff_ext)/dmom_u_ham_grad_p_ext[0];
4266 force_p_x = p_ext*normal[0];
4267 force_p_y = p_ext*normal[1];
4268 force_p_z = p_ext*normal[2];
4270 force_x = force_p_x + force_v_x;
4271 force_y = force_p_y + force_v_y;
4272 force_z = force_p_z + force_v_z;
4274 r_x = x_ext - barycenters[3*boundaryFlags[ebN]+0];
4275 r_y = y_ext - barycenters[3*boundaryFlags[ebN]+1];
4276 r_z = z_ext - barycenters[3*boundaryFlags[ebN]+2];
4278 wettedAreas[boundaryFlags[ebN]] += dS*(1.0-ebqe_vf_ext[ebNE_kb]);
4280 netForces_p[3*boundaryFlags[ebN]+0] += force_p_x*dS;
4281 netForces_p[3*boundaryFlags[ebN]+1] += force_p_y*dS;
4282 netForces_p[3*boundaryFlags[ebN]+2] += force_p_z*dS;
4284 netForces_v[3*boundaryFlags[ebN]+0] += force_v_x*dS;
4285 netForces_v[3*boundaryFlags[ebN]+1] += force_v_y*dS;
4286 netForces_v[3*boundaryFlags[ebN]+2] += force_v_z*dS;
4288 netMoments[3*boundaryFlags[ebN]+0] += (r_y*force_z - r_z*force_y)*dS;
4289 netMoments[3*boundaryFlags[ebN]+1] += (r_z*force_x - r_x*force_z)*dS;
4290 netMoments[3*boundaryFlags[ebN]+2] += (r_x*force_y - r_y*force_x)*dS;
4295 for (
int i=0;i<nDOF_test_element;i++)
4297 elementResidual_mesh[i] -=
ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]+zt_ext*normal[2]),p_test_dS[i]);
4298 elementResidual_p[i] +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
4299 elementResidual_p[i] -=
DM*
ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]+zt_ext*normal[2]),p_test_dS[i]);
4300 globalConservationError +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
4301 elementResidual_u[i] +=
4302 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[0] : 0.) +
4303 ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
4304 ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
4305 ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
4306 ck.ExteriorElementBoundaryFlux(flux_mom_uw_diff_ext,vel_test_dS[i])+
4307 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4308 isDiffusiveFluxBoundary_u[ebNE_kb],
4313 sdInfo_u_u_rowptr.data(),
4314 sdInfo_u_u_colind.data(),
4315 mom_uu_diff_ten_ext,
4316 &vel_grad_test_dS[i*nSpace])+
4317 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4318 isDiffusiveFluxBoundary_u[ebNE_kb],
4323 sdInfo_u_v_rowptr.data(),
4324 sdInfo_u_v_colind.data(),
4325 mom_uv_diff_ten_ext,
4326 &vel_grad_test_dS[i*nSpace])+
4327 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w[ebNE_kb],
4328 isDiffusiveFluxBoundary_u[ebNE_kb],
4333 sdInfo_u_w_rowptr.data(),
4334 sdInfo_u_w_colind.data(),
4335 mom_uw_diff_ten_ext,
4336 &vel_grad_test_dS[i*nSpace]);
4337 elementResidual_v[i] +=
4338 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[1] : 0.) +
4339 ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
4340 ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
4341 ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
4342 ck.ExteriorElementBoundaryFlux(flux_mom_vw_diff_ext,vel_test_dS[i])+
4343 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4344 isDiffusiveFluxBoundary_v[ebNE_kb],
4349 sdInfo_v_u_rowptr.data(),
4350 sdInfo_v_u_colind.data(),
4351 mom_vu_diff_ten_ext,
4352 &vel_grad_test_dS[i*nSpace])+
4353 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4354 isDiffusiveFluxBoundary_v[ebNE_kb],
4359 sdInfo_v_v_rowptr.data(),
4360 sdInfo_v_v_colind.data(),
4361 mom_vv_diff_ten_ext,
4362 &vel_grad_test_dS[i*nSpace])+
4363 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w[ebNE_kb],
4364 isDiffusiveFluxBoundary_v[ebNE_kb],
4369 sdInfo_v_w_rowptr.data(),
4370 sdInfo_v_w_colind.data(),
4371 mom_vw_diff_ten_ext,
4372 &vel_grad_test_dS[i*nSpace]);
4374 elementResidual_w[i] +=
4375 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[2] : 0.) +
4376 ck.ExteriorElementBoundaryFlux(flux_mom_w_adv_ext,vel_test_dS[i]) +
4377 ck.ExteriorElementBoundaryFlux(flux_mom_wu_diff_ext,vel_test_dS[i])+
4378 ck.ExteriorElementBoundaryFlux(flux_mom_wv_diff_ext,vel_test_dS[i])+
4379 ck.ExteriorElementBoundaryFlux(flux_mom_ww_diff_ext,vel_test_dS[i])+
4380 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4381 isDiffusiveFluxBoundary_w[ebNE_kb],
4386 sdInfo_w_u_rowptr.data(),
4387 sdInfo_w_u_colind.data(),
4388 mom_wu_diff_ten_ext,
4389 &vel_grad_test_dS[i*nSpace])+
4390 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4391 isDiffusiveFluxBoundary_w[ebNE_kb],
4396 sdInfo_w_v_rowptr.data(),
4397 sdInfo_w_v_colind.data(),
4398 mom_wv_diff_ten_ext,
4399 &vel_grad_test_dS[i*nSpace])+
4400 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w[ebNE_kb],
4401 isDiffusiveFluxBoundary_w[ebNE_kb],
4406 sdInfo_w_w_rowptr.data(),
4407 sdInfo_w_w_colind.data(),
4408 mom_ww_diff_ten_ext,
4409 &vel_grad_test_dS[i*nSpace]);
4415 for (
int i=0;i<nDOF_test_element;i++)
4417 int eN_i = eN*nDOF_test_element+i;
4422 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=elementResidual_u[i];
4423 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=elementResidual_v[i];
4424 globalResidual[offset_w+stride_w*vel_l2g[eN_i]]+=elementResidual_w[i];
4427 gf.useExact=useExact;
4428 gf_s.useExact=useExact;
4439 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
4440 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
4441 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4442 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
4443 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
4444 double PSTAB = args.
scalar<
double>(
"PSTAB");
4445 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4446 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
4447 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
4448 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
4449 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
4450 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
4451 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
4452 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
4453 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
4454 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
4455 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
4456 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
4457 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
4458 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
4459 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
4460 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
4461 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4462 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4463 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
4464 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
4465 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
4466 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
4467 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
4468 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4469 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
4470 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
4471 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
4472 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4473 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4474 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
4475 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
4476 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
4477 double hFactor = args.
scalar<
double>(
"hFactor");
4478 int nElements_global = args.
scalar<
int>(
"nElements_global");
4479 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
4480 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
4481 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
4482 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
4483 double useRBLES = args.
scalar<
double>(
"useRBLES");
4484 double useMetrics = args.
scalar<
double>(
"useMetrics");
4485 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
4486 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
4487 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
4488 double sigma = args.
scalar<
double>(
"sigma");
4493 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
4494 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
4495 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
4496 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
4497 double C_dg = args.
scalar<
double>(
"C_dg");
4498 double C_b = args.
scalar<
double>(
"C_b");
4499 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
4500 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
4501 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
4502 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
4503 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
4504 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
4505 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
4506 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
4507 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
4508 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
4509 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
4510 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
4511 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
4512 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
4513 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
4514 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
4515 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
4516 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
4517 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4518 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
4519 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4520 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4521 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4522 xt::pyarray<double>& g = args.
array<
double>(
"g");
4523 const double useVF = args.
scalar<
double>(
"useVF");
4524 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
4525 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
4526 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
4527 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
4528 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
4529 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
4530 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
4531 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
4532 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
4533 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
4534 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
4535 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
4536 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
4537 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
4538 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
4539 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
4540 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
4541 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
4542 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
4543 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
4544 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
4545 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
4546 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
4547 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
4548 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
4549 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
4550 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
4551 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
4552 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
4553 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
4554 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
4555 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
4556 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
4557 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
4558 xt::pyarray<int>& csrRowIndeces_p_p = args.
array<
int>(
"csrRowIndeces_p_p");
4559 xt::pyarray<int>& csrColumnOffsets_p_p = args.
array<
int>(
"csrColumnOffsets_p_p");
4560 xt::pyarray<int>& csrRowIndeces_p_u = args.
array<
int>(
"csrRowIndeces_p_u");
4561 xt::pyarray<int>& csrColumnOffsets_p_u = args.
array<
int>(
"csrColumnOffsets_p_u");
4562 xt::pyarray<int>& csrRowIndeces_p_v = args.
array<
int>(
"csrRowIndeces_p_v");
4563 xt::pyarray<int>& csrColumnOffsets_p_v = args.
array<
int>(
"csrColumnOffsets_p_v");
4564 xt::pyarray<int>& csrRowIndeces_p_w = args.
array<
int>(
"csrRowIndeces_p_w");
4565 xt::pyarray<int>& csrColumnOffsets_p_w = args.
array<
int>(
"csrColumnOffsets_p_w");
4566 xt::pyarray<int>& csrRowIndeces_u_p = args.
array<
int>(
"csrRowIndeces_u_p");
4567 xt::pyarray<int>& csrColumnOffsets_u_p = args.
array<
int>(
"csrColumnOffsets_u_p");
4568 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
4569 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
4570 xt::pyarray<int>& csrRowIndeces_u_v = args.
array<
int>(
"csrRowIndeces_u_v");
4571 xt::pyarray<int>& csrColumnOffsets_u_v = args.
array<
int>(
"csrColumnOffsets_u_v");
4572 xt::pyarray<int>& csrRowIndeces_u_w = args.
array<
int>(
"csrRowIndeces_u_w");
4573 xt::pyarray<int>& csrColumnOffsets_u_w = args.
array<
int>(
"csrColumnOffsets_u_w");
4574 xt::pyarray<int>& csrRowIndeces_v_p = args.
array<
int>(
"csrRowIndeces_v_p");
4575 xt::pyarray<int>& csrColumnOffsets_v_p = args.
array<
int>(
"csrColumnOffsets_v_p");
4576 xt::pyarray<int>& csrRowIndeces_v_u = args.
array<
int>(
"csrRowIndeces_v_u");
4577 xt::pyarray<int>& csrColumnOffsets_v_u = args.
array<
int>(
"csrColumnOffsets_v_u");
4578 xt::pyarray<int>& csrRowIndeces_v_v = args.
array<
int>(
"csrRowIndeces_v_v");
4579 xt::pyarray<int>& csrColumnOffsets_v_v = args.
array<
int>(
"csrColumnOffsets_v_v");
4580 xt::pyarray<int>& csrRowIndeces_v_w = args.
array<
int>(
"csrRowIndeces_v_w");
4581 xt::pyarray<int>& csrColumnOffsets_v_w = args.
array<
int>(
"csrColumnOffsets_v_w");
4582 xt::pyarray<int>& csrRowIndeces_w_p = args.
array<
int>(
"csrRowIndeces_w_p");
4583 xt::pyarray<int>& csrColumnOffsets_w_p = args.
array<
int>(
"csrColumnOffsets_w_p");
4584 xt::pyarray<int>& csrRowIndeces_w_u = args.
array<
int>(
"csrRowIndeces_w_u");
4585 xt::pyarray<int>& csrColumnOffsets_w_u = args.
array<
int>(
"csrColumnOffsets_w_u");
4586 xt::pyarray<int>& csrRowIndeces_w_v = args.
array<
int>(
"csrRowIndeces_w_v");
4587 xt::pyarray<int>& csrColumnOffsets_w_v = args.
array<
int>(
"csrColumnOffsets_w_v");
4588 xt::pyarray<int>& csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
4589 xt::pyarray<int>& csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
4590 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
4591 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4592 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4593 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
4594 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4595 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4596 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
4597 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
4598 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
4599 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
4600 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
4601 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
4602 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
4603 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
4604 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
4605 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
4606 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
4607 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
4608 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
4609 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
4610 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
4611 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
4612 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
4613 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
4614 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
4615 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
4616 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
4617 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
4618 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
4619 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
4620 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
4621 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
4622 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
4623 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
4624 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
4625 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
4626 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
4627 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
4628 xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.
array<
int>(
"csrColumnOffsets_eb_p_p");
4629 xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.
array<
int>(
"csrColumnOffsets_eb_p_u");
4630 xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.
array<
int>(
"csrColumnOffsets_eb_p_v");
4631 xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.
array<
int>(
"csrColumnOffsets_eb_p_w");
4632 xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.
array<
int>(
"csrColumnOffsets_eb_u_p");
4633 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
4634 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.
array<
int>(
"csrColumnOffsets_eb_u_v");
4635 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.
array<
int>(
"csrColumnOffsets_eb_u_w");
4636 xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.
array<
int>(
"csrColumnOffsets_eb_v_p");
4637 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.
array<
int>(
"csrColumnOffsets_eb_v_u");
4638 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.
array<
int>(
"csrColumnOffsets_eb_v_v");
4639 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.
array<
int>(
"csrColumnOffsets_eb_v_w");
4640 xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.
array<
int>(
"csrColumnOffsets_eb_w_p");
4641 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.
array<
int>(
"csrColumnOffsets_eb_w_u");
4642 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.
array<
int>(
"csrColumnOffsets_eb_w_v");
4643 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.
array<
int>(
"csrColumnOffsets_eb_w_w");
4644 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
4645 int nParticles = args.
scalar<
int>(
"nParticles");
4646 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
4647 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
4648 double particle_beta = args.
scalar<
double>(
"particle_beta");
4649 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
4650 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
4651 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
4652 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
4653 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
4654 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
4655 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
4656 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
4657 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
4658 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
4659 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
4660 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
4661 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
4662 double dt = args.
scalar<
double>(
"dt");
4663 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
4664 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
4665 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
4666 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
4667 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
4668 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
4669 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
4670 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
4671 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
4672 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
4673 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
4674 int offset_u = args.
scalar<
int>(
"offset_u");
4675 int offset_v = args.
scalar<
int>(
"offset_v");
4676 int offset_w = args.
scalar<
int>(
"offset_w");
4677 int stride_u = args.
scalar<
int>(
"stride_u");
4678 int stride_v = args.
scalar<
int>(
"stride_v");
4679 int stride_w = args.
scalar<
int>(
"stride_w");
4680 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
4681 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
4682 xt::pyarray<int>& rowptr = args.
array<
int>(
"rowptr");
4683 xt::pyarray<int>& colind = args.
array<
int>(
"colind");
4684 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
4688 std::valarray<double> particle_surfaceArea(nParticles), particle_netForces(nParticles * 3), particle_netMoments(nParticles * 3);
4689 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
4690 for(
int eN=0;eN<nElements_global;eN++)
4692 register double eps_rho,eps_mu;
4693 double element_active=1.0;
4695 register double elementJacobian_p_p[nDOF_test_element][nDOF_trial_element],
4696 elementJacobian_p_u[nDOF_test_element][nDOF_trial_element],
4697 elementJacobian_p_v[nDOF_test_element][nDOF_trial_element],
4698 elementJacobian_p_w[nDOF_test_element][nDOF_trial_element],
4699 elementJacobian_u_p[nDOF_test_element][nDOF_trial_element],
4700 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
4701 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
4702 elementJacobian_u_w[nDOF_test_element][nDOF_trial_element],
4703 elementJacobian_v_p[nDOF_test_element][nDOF_trial_element],
4704 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
4705 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element],
4706 elementJacobian_v_w[nDOF_test_element][nDOF_trial_element],
4707 elementJacobian_w_p[nDOF_test_element][nDOF_trial_element],
4708 elementJacobian_w_u[nDOF_test_element][nDOF_trial_element],
4709 elementJacobian_w_v[nDOF_test_element][nDOF_trial_element],
4710 elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
4711 for (
int i=0;i<nDOF_test_element;i++)
4712 for (
int j=0;j<nDOF_trial_element;j++)
4714 elementJacobian_p_p[i][j]=0.0;
4715 elementJacobian_p_u[i][j]=0.0;
4716 elementJacobian_p_v[i][j]=0.0;
4717 elementJacobian_p_w[i][j]=0.0;
4718 elementJacobian_u_p[i][j]=0.0;
4719 elementJacobian_u_u[i][j]=0.0;
4720 elementJacobian_u_v[i][j]=0.0;
4721 elementJacobian_u_w[i][j]=0.0;
4722 elementJacobian_v_p[i][j]=0.0;
4723 elementJacobian_v_u[i][j]=0.0;
4724 elementJacobian_v_v[i][j]=0.0;
4725 elementJacobian_v_w[i][j]=0.0;
4726 elementJacobian_w_p[i][j]=0.0;
4727 elementJacobian_w_u[i][j]=0.0;
4728 elementJacobian_w_v[i][j]=0.0;
4729 elementJacobian_w_w[i][j]=0.0;
4736 double _distance[nDOF_mesh_trial_element]={0.0};
4738 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4740 if(use_ball_as_particle==1)
4743 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
4744 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
4745 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
4750 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
4752 if ( _distance[I] >= 0)
4755 if (pos_counter == 3)
4759 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4761 if (_distance[I] < 0)
4764 assert(opp_node >=0);
4765 assert(opp_node <nDOF_mesh_trial_element);
4767 else if (pos_counter == 4)
4776 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
4777 for (
int j=0;j<nDOF_mesh_trial_element;j++)
4779 register int eN_j = eN*nDOF_mesh_trial_element+j;
4780 element_phi[j] = phi_dof[p_l2g[eN_j]];
4781 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
4783 double element_nodes[nDOF_mesh_trial_element*3];
4784 for (
int i=0;i<nDOF_mesh_trial_element;i++)
4786 register int eN_i=eN*nDOF_mesh_trial_element+i;
4787 for(
int I=0;I<3;I++)
4788 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
4790 gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),
false);
4791 gf.calculate(element_phi, element_nodes, x_ref.data(),
false);
4792 for (
int k=0;k<nQuadraturePoints_element;k++)
4796 int eN_k = eN*nQuadraturePoints_element+k,
4797 eN_k_nSpace = eN_k*nSpace,
4798 eN_nDOF_trial_element = eN*nDOF_trial_element;
4801 register double p=0.0,
u=0.0,
v=0.0,
w=0.0,
4802 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
4811 dmass_adv_u[nSpace],
4812 dmass_adv_v[nSpace],
4813 dmass_adv_w[nSpace],
4815 dmom_u_adv_u[nSpace],
4816 dmom_u_adv_v[nSpace],
4817 dmom_u_adv_w[nSpace],
4819 dmom_v_adv_u[nSpace],
4820 dmom_v_adv_v[nSpace],
4821 dmom_v_adv_w[nSpace],
4823 dmom_w_adv_u[nSpace],
4824 dmom_w_adv_v[nSpace],
4825 dmom_w_adv_w[nSpace],
4826 mom_uu_diff_ten[nSpace],
4827 mom_vv_diff_ten[nSpace],
4828 mom_ww_diff_ten[nSpace],
4839 dmom_u_ham_grad_p[nSpace],
4840 dmom_u_ham_grad_u[nSpace],
4842 dmom_v_ham_grad_p[nSpace],
4843 dmom_v_ham_grad_v[nSpace],
4845 dmom_w_ham_grad_p[nSpace],
4846 dmom_w_ham_grad_w[nSpace],
4857 dpdeResidual_p_u[nDOF_trial_element],dpdeResidual_p_v[nDOF_trial_element],dpdeResidual_p_w[nDOF_trial_element],
4858 dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_trial_element],
4859 dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_trial_element],
4860 dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_trial_element],
4861 Lstar_u_p[nDOF_test_element],
4862 Lstar_v_p[nDOF_test_element],
4863 Lstar_w_p[nDOF_test_element],
4864 Lstar_u_u[nDOF_test_element],
4865 Lstar_v_v[nDOF_test_element],
4866 Lstar_w_w[nDOF_test_element],
4867 Lstar_p_u[nDOF_test_element],
4868 Lstar_p_v[nDOF_test_element],
4869 Lstar_p_w[nDOF_test_element],
4874 dsubgridError_p_u[nDOF_trial_element],
4875 dsubgridError_p_v[nDOF_trial_element],
4876 dsubgridError_p_w[nDOF_trial_element],
4877 dsubgridError_u_p[nDOF_trial_element],
4878 dsubgridError_u_u[nDOF_trial_element],
4879 dsubgridError_v_p[nDOF_trial_element],
4880 dsubgridError_v_v[nDOF_trial_element],
4881 dsubgridError_w_p[nDOF_trial_element],
4882 dsubgridError_w_w[nDOF_trial_element],
4883 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
4884 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
4887 jacInv[nSpace*nSpace],
4888 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
4889 vel_hess_trial[nDOF_trial_element*
nSpace2],
4891 p_test_dV[nDOF_test_element],vel_test_dV[nDOF_test_element],
4892 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
4897 dmom_u_source[nSpace],
4898 dmom_v_source[nSpace],
4899 dmom_w_source[nSpace],
4902 G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
4904 ck.calculateMapping_element(eN,
4908 mesh_trial_ref.data(),
4909 mesh_grad_trial_ref.data(),
4914 ck.calculateH_element(eN,
4916 nodeDiametersArray.data(),
4918 mesh_trial_ref.data(),
4920 ck.calculateMappingVelocity_element(eN,
4922 mesh_velocity_dof.data(),
4924 mesh_trial_ref.data(),
4929 dV = fabs(jacDet)*dV_ref[k];
4930 ck.calculateG(jacInv,G,G_dd_G,tr_G);
4933 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
4934 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
4935 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
4939 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
4940 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
4944 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
4945 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
4946 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
w);
4949 for (
int I=0;I<nSpace;I++)
4950 grad_p[I] = q_grad_p[eN_k_nSpace+I];
4951 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
4952 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
4953 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_w);
4954 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
4955 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
4956 ck.hessFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_w);
4958 for (
int j=0;j<nDOF_trial_element;j++)
4961 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
4962 for (
int I=0;I<nSpace;I++)
4965 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
4969 double div_mesh_velocity=0.0;
4970 int NDOF_MESH_TRIAL_ELEMENT=4;
4971 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
4973 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
4974 div_mesh_velocity +=
4975 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*3+0] +
4976 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*3+1] +
4977 mesh_velocity_dof[mesh_l2g[eN_j]*3+2]*vel_grad_trial[j*3+2];
4979 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
4982 porosity = 1.0 - q_vos[eN_k];
4987 double distance_to_omega_solid = phi_solid[eN_k];
4988 double eddy_viscosity(0.),rhoSave,nuSave;
4997 elementDiameter[eN],
4998 smagorinskyConstant,
4999 turbulenceClosureModel,
5004 &normal_phi[eN_k_nSpace],
5005 distance_to_omega_solid,
5018 q_velocity_sge[eN_k_nSpace+0],
5019 q_velocity_sge[eN_k_nSpace+1],
5020 q_velocity_sge[eN_k_nSpace+2],
5072 MATERIAL_PARAMETERS_AS_FUNCTION,
5073 density_as_function[eN_k],
5074 dynamic_viscosity_as_function[eN_k],
5077 use_ball_as_particle,
5080 ball_velocity.data(),
5081 ball_angular_velocity.data(),
5082 INT_BY_PARTS_PRESSURE);
5084 mass_source = q_mass_source[eN_k];
5085 for (
int I=0;I<nSpace;I++)
5087 dmom_u_source[I] = 0.0;
5088 dmom_v_source[I] = 0.0;
5089 dmom_w_source[I] = 0.0;
5110 q_velocity_sge[eN_k_nSpace+0],
5111 q_velocity_sge[eN_k_nSpace+1],
5112 q_velocity_sge[eN_k_nSpace+2],
5113 eps_solid[elementFlags[eN]],
5115 q_velocity_solid[eN_k_nSpace+0],
5116 q_velocity_solid[eN_k_nSpace+1],
5117 q_velocity_solid[eN_k_nSpace+2],
5118 q_velocityStar_solid[eN_k_nSpace+0],
5119 q_velocityStar_solid[eN_k_nSpace+1],
5120 q_velocityStar_solid[eN_k_nSpace+2],
5127 q_grad_vos[eN_k_nSpace+0],
5128 q_grad_vos[eN_k_nSpace+1],
5129 q_grad_vos[eN_k_nSpace+2]);
5130 double C_particles=0.0;
5131 if(nParticles > 0 && USE_SBM==0)
5136 nQuadraturePoints_global,
5137 &particle_signed_distances[eN_k],
5138 &particle_signed_distance_normals[eN_k_nSpace],
5139 &particle_velocities[eN_k_nSpace],
5140 particle_centroids.data(),
5141 use_ball_as_particle,
5144 ball_velocity.data(),
5145 ball_angular_velocity.data(),
5147 particle_penalty_constant/h_phi,
5148 particle_alpha/h_phi,
5149 particle_beta/h_phi,
5166 q_velocity_sge[eN_k_nSpace+0],
5167 q_velocity_sge[eN_k_nSpace+1],
5168 q_velocity_sge[eN_k_nSpace+2],
5191 &particle_netForces[0],
5192 &particle_netMoments[0],
5193 &particle_surfaceArea[0]);
5195 if (turbulenceClosureModel >= 3)
5197 const double c_mu = 0.09;
5212 &q_turb_var_grad_0[eN_k_nSpace],
5232 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
5233 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
5234 mom_u_adv[2] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*zt;
5235 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
5236 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
5237 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
5239 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
5240 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
5241 mom_v_adv[2] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*zt;
5242 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
5243 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
5244 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
5246 mom_w_adv[0] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*
xt;
5247 mom_w_adv[1] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*yt;
5248 mom_w_adv[2] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*zt;
5249 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
5250 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
5251 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
5256 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5262 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5268 q_mom_w_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5274 mom_u_acc_t *= dmom_u_acc_u;
5275 mom_v_acc_t *= dmom_v_acc_v;
5276 mom_w_acc_t *= dmom_w_acc_w;
5281 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
5282 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
5283 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
5288 ck.Mass_strong(-q_dvos_dt[eN_k]) +
5289 ck.Advection_strong(dmass_adv_u,grad_u) +
5290 ck.Advection_strong(dmass_adv_v,grad_v) +
5291 ck.Advection_strong(dmass_adv_w,grad_w) +
5292 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
5294 ck.Reaction_strong(mass_source);
5298 ck.Mass_strong(mom_u_acc_t) +
5299 ck.Advection_strong(dmom_adv_sge,grad_u) +
5300 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
5301 ck.Reaction_strong(mom_u_source) -
5302 ck.Reaction_strong(
u*div_mesh_velocity);
5305 ck.Mass_strong(mom_v_acc_t) +
5306 ck.Advection_strong(dmom_adv_sge,grad_v) +
5307 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
5308 ck.Reaction_strong(mom_v_source) -
5309 ck.Reaction_strong(
v*div_mesh_velocity);
5311 pdeResidual_w =
ck.Mass_strong(mom_w_acc_t) +
5312 ck.Advection_strong(dmom_adv_sge,grad_w) +
5313 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
5314 ck.Reaction_strong(mom_w_source) -
5315 ck.Reaction_strong(
w*div_mesh_velocity);
5318 for (
int j=0;j<nDOF_trial_element;j++)
5320 register int j_nSpace = j*nSpace;
5321 dpdeResidual_p_u[j]=
ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial[j_nSpace]);
5322 dpdeResidual_p_v[j]=
ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial[j_nSpace]);
5323 dpdeResidual_p_w[j]=
ck.AdvectionJacobian_strong(dmass_adv_w,&vel_grad_trial[j_nSpace]);
5325 dpdeResidual_u_p[j]=
ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace]);
5326 dpdeResidual_u_u[j]=
ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5327 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5328 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5330 dpdeResidual_v_p[j]=
ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace]);
5331 dpdeResidual_v_v[j]=
ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5332 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5333 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5335 dpdeResidual_w_p[j]=
ck.HamiltonianJacobian_strong(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace]);
5336 dpdeResidual_w_w[j]=
ck.MassJacobian_strong(dmom_w_acc_w_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5337 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5338 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5341 dpdeResidual_u_u[j]+=
ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial_ref[k*nDOF_trial_element+j]);
5342 dpdeResidual_v_v[j]+=
ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial_ref[k*nDOF_trial_element+j]);
5343 dpdeResidual_w_w[j]+=
ck.ReactionJacobian_strong(dmom_w_source[2],vel_trial_ref[k*nDOF_trial_element+j]);
5348 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
5350 elementDiameter[eN],
5355 dmom_u_ham_grad_p[0],
5365 dmom_u_ham_grad_p[0],
5371 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
5372 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
5405 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
5406 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
5407 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
5410 for (
int i=0;i<nDOF_test_element;i++)
5412 register int i_nSpace = i*nSpace;
5413 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
5414 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
5415 Lstar_w_p[i]=
ck.Advection_adjoint(dmass_adv_w,&p_grad_test_dV[i_nSpace]);
5416 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5417 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5418 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5419 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5420 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5421 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5423 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
5424 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
5425 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
5429 dmom_u_adv_u[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5430 dmom_u_adv_u[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5431 dmom_u_adv_u[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
5433 dmom_v_adv_v[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5434 dmom_v_adv_v[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5435 dmom_v_adv_v[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
5437 dmom_w_adv_w[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5438 dmom_w_adv_w[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5439 dmom_w_adv_w[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
5442 double unit_normal[nSpace];
5443 double norm_grad_phi = 0.;
5444 for (
int I=0;I<nSpace;I++)
5445 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
5446 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
5447 for (
int I=0;I<nSpace;I++)
5448 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
5449 double delta =
gf.D(eps_mu,
phi[eN_k]);
5450 register double vel_tgrad_test_i[nSpace], vel_tgrad_test_j[nSpace];
5454 for(
int i=0;i<nDOF_test_element;i++)
5456 register int i_nSpace = i*nSpace;
5458 &vel_grad_trial[i_nSpace],
5460 for(
int j=0;j<nDOF_trial_element;j++)
5462 register int j_nSpace = j*nSpace;
5464 &vel_grad_trial[j_nSpace],