diff --git a/.gitattributes b/.gitattributes index d3eac55dbb..d97c6d8906 100644 --- a/.gitattributes +++ b/.gitattributes @@ -18,3 +18,4 @@ proteus/tests/cylinder2D/sbm_method/comparison_files/T1_sbm_rans3p.h5 filter=lfs proteus/tests/HotStart_3P/comparison_files/T01P1_hotstart.h5 filter=lfs diff=lfs merge=lfs -text proteus/tests/HotStart_3P/comparison_files/T01P2_hotstart.h5 filter=lfs diff=lfs merge=lfs -text proteus/tests/cylinder2D/sbm_3Dmesh/comparison_files/T001_P1_sbm_3Dmesh.h5 filter=lfs diff=lfs merge=lfs -text +proteus/tests/cylinder2D/ibm_rans2p/comparison_files/T1_ibm_rans2p.h5 filter=lfs diff=lfs merge=lfs -text diff --git a/proteus/mprans/RANS2P.h b/proteus/mprans/RANS2P.h index 59bf66c210..ac91d2b8ed 100644 --- a/proteus/mprans/RANS2P.h +++ b/proteus/mprans/RANS2P.h @@ -2,12 +2,16 @@ #define RANS2P_H #include #include +#include #include "CompKernel.h" #include "ModelFactory.h" const double DM=0.0;//1-mesh conservation and divergence, 0 - weak div(v) only const double DM2=0.0;//1-point-wise mesh volume strong-residual, 0 - div(v) only const double DM3=1.0;//1-point-wise divergence, 0-point-wise rate of volume change #include "PyEmbeddedFunctions.h" + +#define USE_CYLINDER_AS_PARTICLE//just for debug + namespace proteus { class RANS2P_base @@ -18,7 +22,7 @@ namespace proteus double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABLIZATION, + double PRESSURE_PROJECTION_STABLIZATION, double* numerical_viscosity, //element double* mesh_trial_ref, @@ -97,7 +101,7 @@ namespace proteus double* w_dof, double* g, const double useVF, - double* q_rho, + double* q_rho, double* vf, double* phi, double* normal_phi, @@ -176,12 +180,30 @@ namespace proteus double* netForces_v, double* netMoments, double* velocityError, - double* velocityErrorNodal)=0; + double* velocityErrorNodal, + double* forcex, + double* forcey, + double* forcez, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + double *particle_netForces, + double *particle_netMoments, + double *particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant)=0; virtual void calculateJacobian(double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABLIZATION, + double PRESSURE_PROJECTION_STABLIZATION, //element double* mesh_trial_ref, double* mesh_grad_trial_ref, @@ -344,7 +366,19 @@ namespace proteus int* csrColumnOffsets_eb_w_v, int* csrColumnOffsets_eb_w_w, int* elementFlags, - int* boundaryFlags)=0; + int* boundaryFlags, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant)=0; virtual void calculateVelocityAverage(int nExteriorElementBoundaries_global, int* exteriorElementBoundariesArray, int nInteriorElementBoundaries_global, @@ -625,7 +659,10 @@ namespace proteus double& dmom_w_ham_u, double& dmom_w_ham_v, double& dmom_w_ham_w, - double& rho) + double& rho, + double forcex, + double forcey, + double forcez) { double nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n,nu_t0=0.0,nu_t1=0.0,nu_t; H_rho = (1.0-useVF)*smoothedHeaviside(eps_rho,phi) + useVF*fmin(1.0,fmax(0.0,vf)); @@ -990,7 +1027,373 @@ namespace proteus dmom_w_ham_v =0.0; dmom_w_ham_w =0.0; } + mom_u_source -= forcex; + mom_v_source -= forcey; + } + int get_distance_to_ball(int n_balls,const double* ball_center, const double* ball_radius, const double x, const double y, const double z, double& distance) + { + distance = 1e10; + int index = -1; + double d_ball_i; + for (int i=0; i 0.0) ? 0.0 : nu * penalty; + double C_vol = (phi_s > 0.0) ? 0.0 : (alpha + beta * rel_vel_norm); + + C = (D_s * C_surf + (1.0 - H_s) * C_vol); + // force_x = dV * D_s * (p * phi_s_normal[0] - porosity * mu * (phi_s_normal[0] * grad_u[0] + phi_s_normal[1] * grad_u[1]) + C_surf * (u - u_s) * rho) + + // dV * (1.0 - H_s) * C_vol * (u - u_s) * rho; + // force_y = dV * D_s * (p * phi_s_normal[1] - porosity * mu * (phi_s_normal[0] * grad_v[0] + phi_s_normal[1] * grad_v[1]) + C_surf * (v - v_s) * rho) + + // dV * (1.0 - H_s) * C_vol * (v - v_s) * rho; +// force_x = dV*D_s*(p*fluid_outward_normal[0] - porosity*mu*(fluid_outward_normal[0]*grad_u[0] + fluid_outward_normal[1]*grad_u[1]) + C_surf*rel_vel_norm*(u-u_s)*rho) + dV*(1.0 - H_s)*C_vol*(u-u_s)*rho; +// force_y = dV*D_s*(p*fluid_outward_normal[1] - porosity*mu*(fluid_outward_normal[0]*grad_v[0] + fluid_outward_normal[1]*grad_v[1]) + C_surf*rel_vel_norm*(v-v_s)*rho) + dV*(1.0 - H_s)*C_vol*(v-v_s)*rho; +// force_x = dV * D_s * (p * fluid_outward_normal[0] +// -mu * (fluid_outward_normal[0] * 2* grad_u[0] + fluid_outward_normal[1] * (grad_u[1]+grad_v[0])) +// ); +// force_y = dV * D_s * (p * fluid_outward_normal[1] +// -mu * (fluid_outward_normal[0] * (grad_u[1]+grad_v[0]) + fluid_outward_normal[1] * 2* grad_v[1]) +// ); + force_p_x = dV * D_s * p * fluid_outward_normal[0]; + force_stress_x = dV * D_s * (-mu) * (fluid_outward_normal[0] * 2* grad_u[0] + +fluid_outward_normal[1] * (grad_v[0]+grad_u[1]) + +fluid_outward_normal[2] * (grad_w[0]+grad_u[2])); + force_p_y = dV * D_s * p * fluid_outward_normal[1]; + force_stress_y = dV * D_s * (-mu) * (fluid_outward_normal[0] * (grad_u[1]+grad_v[0]) + +fluid_outward_normal[1] * 2* grad_v[1] + +fluid_outward_normal[2] * (grad_w[1]+grad_v[2])); + force_p_z = dV * D_s * p * fluid_outward_normal[2]; + force_stress_z = dV * D_s * (-mu) * (fluid_outward_normal[0] * (grad_u[2]+grad_w[0]) + +fluid_outward_normal[1] * (grad_v[2]+grad_w[1]) + +fluid_outward_normal[1] * 2* grad_w[2]); + + force_x = force_p_x + force_stress_x; + force_y = force_p_y + force_stress_y; + force_z = force_p_z + force_stress_z; + //always 3D for particle centroids + r_x = x - center[0]; + r_y = y - center[1]; + r_z = z - center[2]; + + if (element_owned) + { + particle_surfaceArea[i] += dV * D_s; + particle_netForces[i * 3 + 0] += force_x; + particle_netForces[i * 3 + 1] += force_y; + particle_netForces[i * 3 + 2] += force_z; + particle_netForces[(i+ nParticles)*3+0]+= force_stress_x; + particle_netForces[(i+2*nParticles)*3+0]+= force_p_x; + particle_netForces[(i+ nParticles)*3+1]+= force_stress_y; + particle_netForces[(i+2*nParticles)*3+1]+= force_p_y; + particle_netForces[(i+ nParticles)*3+2]+= force_stress_z; + particle_netForces[(i+2*nParticles)*3+2]+= force_p_z; + particle_netMoments[i*3+0] += (r_y*force_z - r_z*force_y); + particle_netMoments[i*3+1] += (r_z*force_x - r_x*force_z); + particle_netMoments[i*3+2] += (r_x*force_y - r_y*force_x); + } + + // These should be done inside to make sure the correct velocity of different particles are used + //(1) + mom_u_source += C * (u - u_s); + mom_v_source += C * (v - v_s); + mom_w_source += C * (w - w_s); + + dmom_u_source[0] += C; + dmom_v_source[1] += C; + dmom_w_source[2] += C; + + if (NONCONSERVATIVE_FORM > 0.0) + { + //(2) + 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]); + dmom_u_ham_grad_u[0] -= D_s * porosity * nu * fluid_outward_normal[0]; + dmom_u_ham_grad_u[1] -= D_s * porosity * nu * fluid_outward_normal[1]; + dmom_u_ham_grad_u[2] -= D_s * porosity * nu * fluid_outward_normal[2]; + + 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]); + dmom_v_ham_grad_v[0] -= D_s * porosity * nu * fluid_outward_normal[0]; + dmom_v_ham_grad_v[1] -= D_s * porosity * nu * fluid_outward_normal[1]; + dmom_v_ham_grad_v[2] -= D_s * porosity * nu * fluid_outward_normal[2]; + + mom_w_ham -= D_s * porosity * nu * (fluid_outward_normal[0] * grad_w[0] + fluid_outward_normal[1] * grad_w[1] + fluid_outward_normal[2]*grad_w[2]); + dmom_w_ham_grad_w[0] -= D_s * porosity * nu * fluid_outward_normal[0]; + dmom_w_ham_grad_w[1] -= D_s * porosity * nu * fluid_outward_normal[1]; + dmom_w_ham_grad_w[2] -= D_s * porosity * nu * fluid_outward_normal[2]; + + //(3) + mom_u_adv[0] += D_s * porosity * nu * fluid_outward_normal[0] * (u - u_s); + mom_u_adv[1] += D_s * porosity * nu * fluid_outward_normal[1] * (u - u_s); + mom_u_adv[2] += D_s * porosity * nu * fluid_outward_normal[2] * (u - u_s); + dmom_u_adv_u[0] += D_s * porosity * nu * fluid_outward_normal[0]; + dmom_u_adv_u[1] += D_s * porosity * nu * fluid_outward_normal[1]; + dmom_u_adv_u[2] += D_s * porosity * nu * fluid_outward_normal[2]; + + mom_v_adv[0] += D_s * porosity * nu * fluid_outward_normal[0] * (v - v_s); + mom_v_adv[1] += D_s * porosity * nu * fluid_outward_normal[1] * (v - v_s); + mom_v_adv[2] += D_s * porosity * nu * fluid_outward_normal[2] * (v - v_s); + dmom_v_adv_v[0] += D_s * porosity * nu * fluid_outward_normal[0]; + dmom_v_adv_v[1] += D_s * porosity * nu * fluid_outward_normal[1]; + dmom_v_adv_v[2] += D_s * porosity * nu * fluid_outward_normal[2]; + + mom_w_adv[0] += D_s * porosity * nu * fluid_outward_normal[0] * (w - w_s); + mom_w_adv[1] += D_s * porosity * nu * fluid_outward_normal[1] * (w - w_s); + mom_w_adv[2] += D_s * porosity * nu * fluid_outward_normal[2] * (w - w_s); + dmom_w_adv_w[0] += D_s * porosity * nu * fluid_outward_normal[0]; + dmom_w_adv_w[1] += D_s * porosity * nu * fluid_outward_normal[1]; + dmom_w_adv_w[2] += D_s * porosity * nu * fluid_outward_normal[2]; + } + else + { + //(2) + mom_u_ham -= D_s * porosity * nu/rho * (fluid_outward_normal[0] * grad_u[0] + fluid_outward_normal[1] * grad_u[1] + fluid_outward_normal[2]*grad_u[2]); + dmom_u_ham_grad_u[0] -= D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_u_ham_grad_u[1] -= D_s * porosity * nu/rho * fluid_outward_normal[1]; + dmom_u_ham_grad_u[2] -= D_s * porosity * nu/rho * fluid_outward_normal[2]; + + mom_v_ham -= D_s * porosity * nu/rho * (fluid_outward_normal[0] * grad_v[0] + fluid_outward_normal[1] * grad_v[1] + fluid_outward_normal[2]*grad_v[2]); + dmom_v_ham_grad_v[0] -= D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_v_ham_grad_v[1] -= D_s * porosity * nu/rho * fluid_outward_normal[1]; + dmom_v_ham_grad_v[2] -= D_s * porosity * nu/rho * fluid_outward_normal[2]; + + + mom_w_ham -= D_s * porosity * nu/rho * (fluid_outward_normal[0] * grad_w[0] + fluid_outward_normal[1] * grad_w[1] + fluid_outward_normal[2]*grad_w[2]); + dmom_w_ham_grad_w[0] -= D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_w_ham_grad_w[1] -= D_s * porosity * nu/rho * fluid_outward_normal[1]; + dmom_w_ham_grad_w[2] -= D_s * porosity * nu/rho * fluid_outward_normal[2]; + //(3) + mom_u_adv[0] += D_s * porosity * nu/rho * fluid_outward_normal[0] * (u - u_s); + mom_u_adv[1] += D_s * porosity * nu/rho * fluid_outward_normal[1] * (u - u_s); + mom_u_adv[2] += D_s * porosity * nu/rho * fluid_outward_normal[2] * (u - u_s); + dmom_u_adv_u[0] += D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_u_adv_u[1] += D_s * porosity * nu/rho * fluid_outward_normal[1]; + dmom_u_adv_u[2] += D_s * porosity * nu/rho * fluid_outward_normal[2]; + + mom_v_adv[0] += D_s * porosity * nu/rho * fluid_outward_normal[0] * (v - v_s); + mom_v_adv[1] += D_s * porosity * nu/rho * fluid_outward_normal[1] * (v - v_s); + mom_v_adv[2] += D_s * porosity * nu/rho * fluid_outward_normal[2] * (v - v_s); + dmom_v_adv_v[0] += D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_v_adv_v[1] += D_s * porosity * nu/rho * fluid_outward_normal[1]; + dmom_v_adv_v[2] += D_s * porosity * nu/rho * fluid_outward_normal[2]; + + mom_w_adv[0] += D_s * porosity * nu/rho * fluid_outward_normal[0] * (w - w_s); + mom_w_adv[1] += D_s * porosity * nu/rho * fluid_outward_normal[1] * (w - w_s); + mom_w_adv[2] += D_s * porosity * nu/rho * fluid_outward_normal[2] * (w - w_s); + dmom_w_adv_w[0] += D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_w_adv_w[1] += D_s * porosity * nu/rho * fluid_outward_normal[1]; + dmom_w_adv_w[2] += D_s * porosity * nu/rho * fluid_outward_normal[2]; + //(4) +// mom_u_ham += D_s * porosity * (fluid_outward_normal[0] * u + fluid_outward_normal[1] * v)*u; +// mom_v_ham += D_s * porosity * (fluid_outward_normal[0] * u + fluid_outward_normal[1] * v)*v; +// dmom_u_ham_u += D_s * porosity * fluid_outward_normal[0] * u * 2.0; +// dmom_u_ham_v += D_s * porosity * fluid_outward_normal[1] * u; +// dmom_v_ham_u += D_s * porosity * fluid_outward_normal[0] * v; +// dmom_v_ham_v += D_s * porosity * fluid_outward_normal[1] * v * 2.0; + + } + //(6) +// mass_ham += D_s * porosity * (fluid_outward_normal[0] * u + fluid_outward_normal[1] * v); +// dmass_ham_u += D_s * porosity * fluid_outward_normal[0]; +// dmass_ham_v += D_s * porosity * fluid_outward_normal[1]; + //(7) +// mass_ham += C_surf * D_s * porosity * (fluid_outward_normal[0] * (u - u_s) + fluid_outward_normal[1] * (v - v_s)); +// dmass_ham_u += C_surf * D_s * porosity * fluid_outward_normal[0]; +// dmass_ham_v += C_surf * D_s * porosity * fluid_outward_normal[1]; + } + } //VRANS specific inline void updateDarcyForchheimerTerms_Ergun(const double NONCONSERVATIVE_FORM, @@ -1942,41 +2345,41 @@ namespace proteus double PRESSURE_SGE, double VELOCITY_SGE, double PRESSURE_PROJECTION_STABLIZATION, - double* numerical_viscosity, + double *numerical_viscosity, //element - double* mesh_trial_ref, - double* mesh_grad_trial_ref, - double* mesh_dof, - double* mesh_velocity_dof, + double *mesh_trial_ref, + double *mesh_grad_trial_ref, + double *mesh_dof, + double *mesh_velocity_dof, double MOVING_DOMAIN, - int* mesh_l2g, - double* dV_ref, - double* p_trial_ref, - double* p_grad_trial_ref, - double* p_test_ref, - double* p_grad_test_ref, - double* vel_trial_ref, - double* vel_grad_trial_ref, - double* vel_test_ref, - double* vel_grad_test_ref, + int *mesh_l2g, + double *dV_ref, + double *p_trial_ref, + double *p_grad_trial_ref, + double *p_test_ref, + double *p_grad_test_ref, + double *vel_trial_ref, + double *vel_grad_trial_ref, + double *vel_test_ref, + double *vel_grad_test_ref, //element boundary - double* mesh_trial_trace_ref, - double* mesh_grad_trial_trace_ref, - double* dS_ref, - double* p_trial_trace_ref, - double* p_grad_trial_trace_ref, - double* p_test_trace_ref, - double* p_grad_test_trace_ref, - double* vel_trial_trace_ref, - double* vel_grad_trial_trace_ref, - double* vel_test_trace_ref, - double* vel_grad_test_trace_ref, - double* normal_ref, - double* boundaryJac_ref, + double *mesh_trial_trace_ref, + double *mesh_grad_trial_trace_ref, + double *dS_ref, + double *p_trial_trace_ref, + double *p_grad_trial_trace_ref, + double *p_test_trace_ref, + double *p_grad_test_trace_ref, + double *vel_trial_trace_ref, + double *vel_grad_trial_trace_ref, + double *vel_test_trace_ref, + double *vel_grad_test_trace_ref, + double *normal_ref, + double *boundaryJac_ref, //physics double eb_adjoint_sigma, - double* elementDiameter, - double* nodeDiametersArray, + double *elementDiameter, + double *nodeDiametersArray, double hFactor, int nElements_global, int nElementBoundaries_owned, @@ -1997,112 +2400,132 @@ namespace proteus double C_dc, double C_b, //VRANS - const double* eps_solid, - const double* phi_solid, - const double* q_velocity_solid, - const double* q_porosity, - const double* q_dragAlpha, - const double* q_dragBeta, - const double* q_mass_source, - const double* q_turb_var_0, - const double* q_turb_var_1, - const double* q_turb_var_grad_0, + const double *eps_solid, + const double *phi_solid, + const double *q_velocity_solid, + const double *q_porosity, + const double *q_dragAlpha, + const double *q_dragBeta, + const double *q_mass_source, + const double *q_turb_var_0, + const double *q_turb_var_1, + const double *q_turb_var_grad_0, const double LAG_LES, - double * q_eddy_viscosity, - double * q_eddy_viscosity_last, - double * ebqe_eddy_viscosity, - double * ebqe_eddy_viscosity_last, + double *q_eddy_viscosity, + double *q_eddy_viscosity_last, + double *ebqe_eddy_viscosity, + double *ebqe_eddy_viscosity_last, // - int* p_l2g, - int* vel_l2g, - double* p_dof, - double* u_dof, - double* v_dof, - double* w_dof, - double* g, + int *p_l2g, + int *vel_l2g, + double *p_dof, + double *u_dof, + double *v_dof, + double *w_dof, + double *g, const double useVF, - double* q_rho, - double* vf, - double* phi, - double* normal_phi, - double* kappa_phi, - double* q_mom_u_acc, - double* q_mom_v_acc, - double* q_mom_w_acc, - double* q_mass_adv, - double* q_mom_u_acc_beta_bdf, double* q_mom_v_acc_beta_bdf, double* q_mom_w_acc_beta_bdf, - double* q_dV, - double* q_dV_last, - double* q_velocity_sge, - double* q_cfl, - double* q_numDiff_u, double* q_numDiff_v, double* q_numDiff_w, - double* q_numDiff_u_last, double* q_numDiff_v_last, double* q_numDiff_w_last, - int* sdInfo_u_u_rowptr,int* sdInfo_u_u_colind, - int* sdInfo_u_v_rowptr,int* sdInfo_u_v_colind, - int* sdInfo_u_w_rowptr,int* sdInfo_u_w_colind, - int* sdInfo_v_v_rowptr,int* sdInfo_v_v_colind, - int* sdInfo_v_u_rowptr,int* sdInfo_v_u_colind, - int* sdInfo_v_w_rowptr,int* sdInfo_v_w_colind, - int* sdInfo_w_w_rowptr,int* sdInfo_w_w_colind, - int* sdInfo_w_u_rowptr,int* sdInfo_w_u_colind, - int* sdInfo_w_v_rowptr,int* sdInfo_w_v_colind, + double *q_rho, + double *vf, + double *phi, + double *normal_phi, + double *kappa_phi, + double *q_mom_u_acc, + double *q_mom_v_acc, + double *q_mom_w_acc, + double *q_mass_adv, + double *q_mom_u_acc_beta_bdf, double *q_mom_v_acc_beta_bdf, double *q_mom_w_acc_beta_bdf, + double *q_dV, + double *q_dV_last, + double *q_velocity_sge, + double *q_cfl, + double *q_numDiff_u, double *q_numDiff_v, double *q_numDiff_w, + double *q_numDiff_u_last, double *q_numDiff_v_last, double *q_numDiff_w_last, + int *sdInfo_u_u_rowptr, int *sdInfo_u_u_colind, + int *sdInfo_u_v_rowptr, int *sdInfo_u_v_colind, + int *sdInfo_u_w_rowptr, int *sdInfo_u_w_colind, + int *sdInfo_v_v_rowptr, int *sdInfo_v_v_colind, + int *sdInfo_v_u_rowptr, int *sdInfo_v_u_colind, + int *sdInfo_v_w_rowptr, int *sdInfo_v_w_colind, + int *sdInfo_w_w_rowptr, int *sdInfo_w_w_colind, + int *sdInfo_w_u_rowptr, int *sdInfo_w_u_colind, + int *sdInfo_w_v_rowptr, int *sdInfo_w_v_colind, int offset_p, int offset_u, int offset_v, int offset_w, int stride_p, int stride_u, int stride_v, int stride_w, - double* globalResidual, + double *globalResidual, int nExteriorElementBoundaries_global, - int* exteriorElementBoundariesArray, - int* elementBoundaryElementsArray, - int* elementBoundaryLocalElementBoundariesArray, - double* ebqe_vf_ext, - double* bc_ebqe_vf_ext, - double* ebqe_phi_ext, - double* bc_ebqe_phi_ext, - double* ebqe_normal_phi_ext, - double* ebqe_kappa_phi_ext, + int *exteriorElementBoundariesArray, + int *elementBoundaryElementsArray, + int *elementBoundaryLocalElementBoundariesArray, + double *ebqe_vf_ext, + double *bc_ebqe_vf_ext, + double *ebqe_phi_ext, + double *bc_ebqe_phi_ext, + double *ebqe_normal_phi_ext, + double *ebqe_kappa_phi_ext, //VRANS - const double* ebqe_porosity_ext, - const double* ebqe_turb_var_0, - const double* ebqe_turb_var_1, + const double *ebqe_porosity_ext, + const double *ebqe_turb_var_0, + const double *ebqe_turb_var_1, //VRANS end - int* isDOFBoundary_p, - int* isDOFBoundary_u, - int* isDOFBoundary_v, - int* isDOFBoundary_w, - int* isAdvectiveFluxBoundary_p, - int* isAdvectiveFluxBoundary_u, - int* isAdvectiveFluxBoundary_v, - int* isAdvectiveFluxBoundary_w, - int* isDiffusiveFluxBoundary_u, - int* isDiffusiveFluxBoundary_v, - int* isDiffusiveFluxBoundary_w, - double* ebqe_bc_p_ext, - double* ebqe_bc_flux_mass_ext, - double* ebqe_bc_flux_mom_u_adv_ext, - double* ebqe_bc_flux_mom_v_adv_ext, - double* ebqe_bc_flux_mom_w_adv_ext, - double* ebqe_bc_u_ext, - double* ebqe_bc_flux_u_diff_ext, - double* ebqe_penalty_ext, - double* ebqe_bc_v_ext, - double* ebqe_bc_flux_v_diff_ext, - double* ebqe_bc_w_ext, - double* ebqe_bc_flux_w_diff_ext, - double* q_x, - double* q_velocity, - double* ebqe_velocity, - double* flux, - double* elementResidual_p_save, - int* elementFlags, - int* boundaryFlags, - double* barycenters, - double* wettedAreas, - double* netForces_p, - double* netForces_v, - double* netMoments, - double* velocityError, - double* velocityErrorNodal) + int *isDOFBoundary_p, + int *isDOFBoundary_u, + int *isDOFBoundary_v, + int *isDOFBoundary_w, + int *isAdvectiveFluxBoundary_p, + int *isAdvectiveFluxBoundary_u, + int *isAdvectiveFluxBoundary_v, + int *isAdvectiveFluxBoundary_w, + int *isDiffusiveFluxBoundary_u, + int *isDiffusiveFluxBoundary_v, + int *isDiffusiveFluxBoundary_w, + double *ebqe_bc_p_ext, + double *ebqe_bc_flux_mass_ext, + double *ebqe_bc_flux_mom_u_adv_ext, + double *ebqe_bc_flux_mom_v_adv_ext, + double *ebqe_bc_flux_mom_w_adv_ext, + double *ebqe_bc_u_ext, + double *ebqe_bc_flux_u_diff_ext, + double *ebqe_penalty_ext, + double *ebqe_bc_v_ext, + double *ebqe_bc_flux_v_diff_ext, + double *ebqe_bc_w_ext, + double *ebqe_bc_flux_w_diff_ext, + double *q_x, + double *q_velocity, + double *ebqe_velocity, + double *flux, + double *elementResidual_p_save, + int *elementFlags, + int *boundaryFlags, + double *barycenters, + double *wettedAreas, + double *netForces_p, + double *netForces_v, + double *netMoments, + double *velocityError, + double *velocityErrorNodal, + double *forcex, + double *forcey, + double *forcez, + int use_ball_as_particle, + double *ball_center, + double *ball_radius, + double *ball_velocity, + double *ball_angular_velocity, + int nParticles, + double *particle_netForces, + double *particle_netMoments, + double *particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) { logEvent("Entered mprans 3D calculateResidual",6); + const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element); + // //loop over elements to compute volume integrals and load them into element and global residual // @@ -2153,6 +2576,10 @@ namespace proteus dmass_adv_u[nSpace], dmass_adv_v[nSpace], dmass_adv_w[nSpace], + mass_ham=0.0, + dmass_ham_u=0.0, + dmass_ham_v=0.0, + dmass_ham_w=0.0, mom_u_adv[nSpace], dmom_u_adv_u[nSpace], dmom_u_adv_v[nSpace], @@ -2404,7 +2831,10 @@ namespace proteus dmom_w_ham_u, dmom_w_ham_v, dmom_w_ham_w, - q_rho[eN_k]); + q_rho[eN_k], + forcex[eN_k], + forcey[eN_k], + forcez[eN_k]); //VRANS mass_source = q_mass_source[eN_k]; //todo: decide if these should be lagged or not? @@ -2441,7 +2871,84 @@ namespace proteus dmom_u_source, dmom_v_source, dmom_w_source); - + const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]); + if(nParticles > 0) + updateSolidParticleTerms(NONCONSERVATIVE_FORM, + eN < nElements_owned, + particle_nitsche, + dV, + nParticles, + nQuadraturePoints_global, +// &particle_signed_distances[eN_k], +// &particle_signed_distance_normals[eN_k_nSpace], +// particle_velocities, +// particle_centroids, + use_ball_as_particle, + ball_center, + ball_radius, + ball_velocity, + ball_angular_velocity, + porosity, + particle_penalty_constant/h_phi,//penalty, + particle_alpha, + particle_beta, + eps_rho, + eps_mu, + rho_0, + nu_0, + rho_1, + nu_1, + useVF, + vf[eN_k], + phi[eN_k], + x, + y, + z, + p, + u, + v, + w, + q_velocity_sge[eN_k_nSpace+0], + q_velocity_sge[eN_k_nSpace+1], + q_velocity_sge[eN_k_nSpace+1], + particle_eps, + grad_u, + grad_v, + grad_w, + mom_u_source, + mom_v_source, + mom_w_source, + dmom_u_source, + dmom_v_source, + dmom_w_source, + mom_u_adv, + mom_v_adv, + mom_w_adv, + dmom_u_adv_u, + dmom_v_adv_v, + dmom_w_adv_w, + mom_u_ham, + dmom_u_ham_grad_u, + dmom_u_ham_u, + dmom_u_ham_v, + dmom_u_ham_w, + mom_v_ham, + dmom_v_ham_grad_v, + dmom_v_ham_u, + dmom_v_ham_v, + dmom_v_ham_w, + mom_w_ham, + dmom_w_ham_grad_w, + dmom_w_ham_u, + dmom_w_ham_v, + dmom_w_ham_w, + mass_ham, + dmass_ham_u, + dmass_ham_v, + dmass_ham_w, + &particle_netForces[0], + &particle_netMoments[0], + &particle_surfaceArea[0]); //Turbulence closure model if (turbulenceClosureModel >= 3) { @@ -2721,8 +3228,9 @@ namespace proteus ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last[eN_k]/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace]); - elementResidual_p[i] += ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) + - DM*MOVING_DOMAIN*(ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) - + elementResidual_p[i] += ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) + + ck.Hamiltonian_weak(mass_ham, p_test_dV[i]) + + DM*MOVING_DOMAIN*(ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) - ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last[eN_k]/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])) + //VRANS @@ -3047,7 +3555,7 @@ namespace proteus //calculate the pde coefficients using the solution and the boundary values for the solution // double bc_eddy_viscosity_ext(0.); //not interested in saving boundary eddy viscosity for now - double rho; + double rho; evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -3131,7 +3639,10 @@ namespace proteus dmom_w_ham_u_ext, dmom_w_ham_v_ext, dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -3215,7 +3726,10 @@ namespace proteus bc_dmom_w_ham_u_ext, bc_dmom_w_ham_v_ext, bc_dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); //Turbulence closure model if (turbulenceClosureModel >= 3) @@ -3896,8 +4410,23 @@ namespace proteus int* csrColumnOffsets_eb_w_v, int* csrColumnOffsets_eb_w_w, int* elementFlags, - int* boundaryFlags) + int* boundaryFlags, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) { + const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element); + std::valarray particle_surfaceArea(nParticles), particle_netForces(nParticles*3*3), particle_netMoments(nParticles*3); + // //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian // @@ -3960,6 +4489,10 @@ namespace proteus dmass_adv_u[nSpace], dmass_adv_v[nSpace], dmass_adv_w[nSpace], + mass_ham=0.0, + dmass_ham_u=0.0, + dmass_ham_v=0.0, + dmass_ham_w=0.0, mom_u_adv[nSpace], dmom_u_adv_u[nSpace], dmom_u_adv_v[nSpace], @@ -4134,7 +4667,7 @@ namespace proteus //calculate pde coefficients and derivatives at quadrature points // double eddy_viscosity(0.);//not really interested in saving eddy_viscosity in jacobian - double rho; + double rho; evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -4218,7 +4751,10 @@ namespace proteus dmom_w_ham_u, dmom_w_ham_v, dmom_w_ham_w, - rho); + rho, + 0.0, + 0.0, + 0.0); //VRANS mass_source = q_mass_source[eN_k]; //todo: decide if these should be lagged or not @@ -4255,6 +4791,85 @@ namespace proteus dmom_u_source, dmom_v_source, dmom_w_source); + + const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]); + if(nParticles > 0) + updateSolidParticleTerms(NONCONSERVATIVE_FORM, + eN < nElements_owned, + particle_nitsche, + dV, + nParticles, + nQuadraturePoints_global, +// &particle_signed_distances[eN_k], +// &particle_signed_distance_normals[eN_k_nSpace], +// particle_velocities, +// particle_centroids, + use_ball_as_particle, + ball_center, + ball_radius, + ball_velocity, + ball_angular_velocity, + porosity, + particle_penalty_constant/h_phi,//penalty, + particle_alpha, + particle_beta, + eps_rho, + eps_mu, + rho_0, + nu_0, + rho_1, + nu_1, + useVF, + vf[eN_k], + phi[eN_k], + x, + y, + z, + p, + u, + v, + w, + q_velocity_sge[eN_k_nSpace+0], + q_velocity_sge[eN_k_nSpace+1], + q_velocity_sge[eN_k_nSpace+1], + particle_eps, + grad_u, + grad_v, + grad_w, + mom_u_source, + mom_v_source, + mom_w_source, + dmom_u_source, + dmom_v_source, + dmom_w_source, + mom_u_adv, + mom_v_adv, + mom_w_adv, + dmom_u_adv_u, + dmom_v_adv_v, + dmom_w_adv_w, + mom_u_ham, + dmom_u_ham_grad_u, + dmom_u_ham_u, + dmom_u_ham_v, + dmom_u_ham_w, + mom_v_ham, + dmom_v_ham_grad_v, + dmom_v_ham_u, + dmom_v_ham_v, + dmom_v_ham_w, + mom_w_ham, + dmom_w_ham_grad_w, + dmom_w_ham_u, + dmom_w_ham_v, + dmom_w_ham_w, + mass_ham, + dmass_ham_u, + dmass_ham_v, + dmass_ham_w, + &particle_netForces[0], + &particle_netMoments[0], + &particle_surfaceArea[0]); //Turbulence closure model if (turbulenceClosureModel >= 3) { @@ -4566,12 +5181,15 @@ namespace proteus ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_p[i]) + ck.SubgridErrorJacobian(dsubgridError_w_p[j],Lstar_w_p[i]); - elementJacobian_p_u[i][j] += ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + - ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_p[i]); - elementJacobian_p_v[i][j] += ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + - ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_p[i]); - elementJacobian_p_w[i][j] += ck.AdvectionJacobian_weak(dmass_adv_w,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + - ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_p[i]); + elementJacobian_p_u[i][j] += ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + + ck.MassJacobian_weak(dmass_ham_u,vel_trial_ref[k*nDOF_trial_element+j],p_test_dV[i]) + + ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_p[i]); + elementJacobian_p_v[i][j] += ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + + ck.MassJacobian_weak(dmass_ham_v,vel_trial_ref[k*nDOF_trial_element+j],p_test_dV[i]) + + ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_p[i]); + elementJacobian_p_w[i][j] += ck.AdvectionJacobian_weak(dmass_adv_w,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + + ck.MassJacobian_weak(dmass_ham_w,vel_trial_ref[k*nDOF_trial_element+j],p_test_dV[i]) + + ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_p[i]); elementJacobian_u_p[i][j] += ck.HamiltonianJacobian_weak(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) + MOMENTUM_SGE*VELOCITY_SGE*ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_u[i]); @@ -4953,7 +5571,7 @@ namespace proteus //calculate the internal and external trace of the pde coefficients // double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);//not interested in saving boundary eddy viscosity for now - double rho; + double rho; evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -5037,7 +5655,10 @@ namespace proteus dmom_w_ham_u_ext, dmom_w_ham_v_ext, dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -5121,7 +5742,10 @@ namespace proteus bc_dmom_w_ham_u_ext, bc_dmom_w_ham_v_ext, bc_dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); //Turbulence closure model if (turbulenceClosureModel >= 3) { diff --git a/proteus/mprans/RANS2P.py b/proteus/mprans/RANS2P.py index 1fade33eb7..6dd4d35cd2 100644 --- a/proteus/mprans/RANS2P.py +++ b/proteus/mprans/RANS2P.py @@ -208,7 +208,47 @@ def __init__(self, VELOCITY_SGE=1.0, PRESSURE_PROJECTION_STABILIZATION=0.0, phaseFunction=None, - LAG_LES=1.0): + LAG_LES=1.0, + use_ball_as_particle=1, + ball_center=None, + ball_radius=None, + ball_velocity=None, + ball_angular_velocity=None, + nParticles = 0, + particle_epsFact=3.0, + particle_alpha=1000.0, + particle_beta=1000.0, + particle_penalty_constant=1000.0, + particle_nitsche=1.0,): + self.use_ball_as_particle = use_ball_as_particle + self.nParticles = nParticles + self.particle_nitsche = particle_nitsche + self.particle_epsFact = particle_epsFact + self.particle_alpha = particle_alpha + self.particle_beta = particle_beta + self.particle_penalty_constant = particle_penalty_constant + self.particle_netForces = np.zeros((3*self.nParticles, 3), 'd')#####[total_force_1,total_force_2,...,stress_1,stress_2,...,pressure_1,pressure_2,...] + self.particle_netMoments = np.zeros((self.nParticles, 3), 'd') + self.particle_surfaceArea = np.zeros((self.nParticles,), 'd') + if ball_center is None: + self.ball_center = 1e10*numpy.ones((self.nParticles,3),'d') + else: + self.ball_center = ball_center + + if ball_radius is None: + self.ball_radius = 1e10*numpy.ones((self.nParticles,1),'d') + else: + self.ball_radius = ball_radius + + if ball_velocity is None: + self.ball_velocity = numpy.zeros((self.nParticles,3),'d') + else: + self.ball_velocity = ball_velocity + + if ball_angular_velocity is None: + self.ball_angular_velocity = numpy.zeros((self.nParticles,3),'d') + else: + self.ball_angular_velocity = ball_angular_velocity self.LAG_LES=LAG_LES self.phaseFunction=phaseFunction self.NONCONSERVATIVE_FORM=NONCONSERVATIVE_FORM @@ -505,6 +545,9 @@ def initializeMesh(self, mesh): self.forceHistory_p = open(os.path.join(proteus.Profiling.logDir, "forceHistory_p.txt"), "w") self.forceHistory_v = open(os.path.join(proteus.Profiling.logDir, "forceHistory_v.txt"), "w") self.momentHistory = open(os.path.join(proteus.Profiling.logDir, "momentHistory.txt"), "w") + if self.nParticles: + self.particle_forceHistory = open(os.path.join(proteus.Profiling.logDir, "particle_forceHistory.txt"), "w") + self.particle_momentHistory = open(os.path.join(proteus.Profiling.logDir, "particle_momentHistory.txt"), "w") self.comm = comm # initialize so it can run as single phase @@ -794,6 +837,11 @@ def postStep(self, t, firstStep=False): self.forceHistory_v.flush() self.momentHistory.write("%21.15e %21.16e %21.16e\n" % tuple(self.netMoments[-1, :])) self.momentHistory.flush() + if self.nParticles: + self.particle_forceHistory.write("%21.16e %21.16e %21.16e\n" % tuple(self.particle_netForces[0, :])) + self.particle_forceHistory.flush() + self.particle_momentHistory.write("%21.15e %21.16e %21.16e\n" % tuple(self.particle_netMoments[0, :])) + self.particle_momentHistory.flush() class LevelModel(proteus.Transport.OneLevelTransport): @@ -1409,6 +1457,14 @@ def __init__(self, self.velocityErrorNodal = self.u[0].dof.copy() logEvent('WARNING: The boundary fluxes at interpart boundaries are skipped if elementBoundaryMaterialType is 0 for RANS2P-based models. This means that DG methods are currently incompatible with RANS2P.') + + self.q[('force', 0)] = numpy.zeros( + (self.mesh.nElements_global, self.nQuadraturePoints_element), 'd') + self.q[('force', 1)] = numpy.zeros( + (self.mesh.nElements_global, self.nQuadraturePoints_element), 'd') + self.q[('force', 2)] = numpy.zeros( + (self.mesh.nElements_global, self.nQuadraturePoints_element), 'd') + def getResidual(self, u, r): """ Calculate the element residuals and add in to the global residual @@ -1451,6 +1507,9 @@ def getResidual(self, u, r): self.coefficients.netForces_p[:, :] = 0.0 self.coefficients.netForces_v[:, :] = 0.0 self.coefficients.netMoments[:, :] = 0.0 + self.coefficients.particle_netForces[:, :] = 0.0 + self.coefficients.particle_netMoments[:, :] = 0.0 + self.coefficients.particle_surfaceArea[:] = 0.0 if self.forceStrongConditions: for cj in range(len(self.dirichletConditionsForceDOF)): for dofN, g in self.dirichletConditionsForceDOF[cj].DOFBoundaryConditionsDict.iteritems(): @@ -1628,7 +1687,25 @@ def getResidual(self, u, r): self.coefficients.netForces_v, self.coefficients.netMoments, self.q['velocityError'], - self.velocityErrorNodal) + self.velocityErrorNodal, + self.q[('force', 0)], + self.q[('force', 1)], + self.q[('force', 2)], + self.coefficients.use_ball_as_particle, + self.coefficients.ball_center, + self.coefficients.ball_radius, + self.coefficients.ball_velocity, + self.coefficients.ball_angular_velocity, + self.coefficients.nParticles, + self.coefficients.particle_netForces, + self.coefficients.particle_netMoments, + self.coefficients.particle_surfaceArea, + self.mesh.nElements_owned, + self.coefficients.particle_nitsche, + self.coefficients.particle_epsFact, + self.coefficients.particle_alpha, + self.coefficients.particle_beta, + self.coefficients.particle_penalty_constant) from proteus.flcbdfWrappers import globalSum for i in range(self.coefficients.netForces_p.shape[0]): self.coefficients.wettedAreas[i] = globalSum(self.coefficients.wettedAreas[i]) @@ -1646,6 +1723,25 @@ def getResidual(self, u, r): # self.coefficients.netForces_p[i,I] = (125.0* math.pi**2 * 0.125*math.cos(self.timeIntegration.t*math.pi) + 125.0*9.81)/4.0 #if I==2: # self.coefficients.netMoments[i,I] = (4.05* math.pi**2 * (math.pi/4.0)*math.cos(self.timeIntegration.t*math.pi))/4.0 + for i in range(self.coefficients.nParticles): + for I in range(3): + self.coefficients.particle_netForces[i, I] = globalSum( + self.coefficients.particle_netForces[i, I]) + self.coefficients.particle_netForces[i+self.coefficients.nParticles, I] = globalSum( + self.coefficients.particle_netForces[i+self.coefficients.nParticles, I]) + self.coefficients.particle_netForces[i+2*self.coefficients.nParticles, I] = globalSum( + self.coefficients.particle_netForces[i+2*self.coefficients.nParticles, I]) + self.coefficients.particle_netMoments[i, I] = globalSum( + self.coefficients.particle_netMoments[i, I]) + self.coefficients.particle_surfaceArea[i] = globalSum( + self.coefficients.particle_surfaceArea[i]) + logEvent("particle i=" + `i`+ " force " + `self.coefficients.particle_netForces[i]`) + logEvent("particle i=" + `i`+ " moment " + `self.coefficients.particle_netMoments[i]`) + logEvent("particle i=" + `i`+ " surfaceArea " + `self.coefficients.particle_surfaceArea[i]`) + logEvent("particle i=" + `i`+ " stress force " + `self.coefficients.particle_netForces[i+self.coefficients.nParticles]`) + logEvent("particle i=" + `i`+ " pressure force " + `self.coefficients.particle_netForces[i+2*self.coefficients.nParticles]`) + + if self.forceStrongConditions: for cj in range(len(self.dirichletConditionsForceDOF)): for dofN, g in self.dirichletConditionsForceDOF[cj].DOFBoundaryConditionsDict.iteritems(): @@ -1865,7 +1961,19 @@ def getJacobian(self, jacobian): self.csrColumnOffsets_eb[(3, 2)], self.csrColumnOffsets_eb[(3, 3)], self.mesh.elementMaterialTypes, - self.mesh.elementBoundaryMaterialTypes) + self.mesh.elementBoundaryMaterialTypes, + self.coefficients.use_ball_as_particle, + self.coefficients.ball_center, + self.coefficients.ball_radius, + self.coefficients.ball_velocity, + self.coefficients.ball_angular_velocity, + self.coefficients.nParticles, + self.mesh.nElements_owned, + self.coefficients.particle_nitsche, + self.coefficients.particle_epsFact, + self.coefficients.particle_alpha, + self.coefficients.particle_beta, + self.coefficients.particle_penalty_constant) if not self.forceStrongConditions and max(numpy.linalg.norm(self.u[1].dof, numpy.inf), numpy.linalg.norm(self.u[2].dof, numpy.inf), numpy.linalg.norm(self.u[3].dof, numpy.inf)) < 1.0e-8: self.pp_hasConstantNullSpace = True diff --git a/proteus/mprans/RANS2P2D.h b/proteus/mprans/RANS2P2D.h index 9949a176da..96ff62dcbc 100644 --- a/proteus/mprans/RANS2P2D.h +++ b/proteus/mprans/RANS2P2D.h @@ -2,6 +2,7 @@ #define RANS2P2D_H #include #include +#include #include "CompKernel.h" #include "ModelFactory.h" #include "PyEmbeddedFunctions.h" @@ -18,7 +19,7 @@ namespace proteus double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABLIZATION, + double PRESSURE_PROJECTION_STABLIZATION, double* numerical_viscosity, //element double* mesh_trial_ref, @@ -97,7 +98,7 @@ namespace proteus double* w_dof, double* g, const double useVF, - double* q_rho, + double* q_rho, double* vf, double* phi, double* normal_phi, @@ -176,12 +177,30 @@ namespace proteus double* netForces_v, double* netMoments, double* velocityError, - double* velocityErrorNodal)=0; + double* velocityErrorNodal, + double* forcex, + double* forcey, + double* forcez, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + double *particle_netForces, + double *particle_netMoments, + double *particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant)=0; virtual void calculateJacobian(double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABILIZATION, + double PRESSURE_PROJECTION_STABILIZATION, //element double* mesh_trial_ref, double* mesh_grad_trial_ref, @@ -344,7 +363,19 @@ namespace proteus int* csrColumnOffsets_eb_w_v, int* csrColumnOffsets_eb_w_w, int* elementFlags, - int* boundaryFlags)=0; + int* boundaryFlags, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant)=0; virtual void calculateVelocityAverage(int nExteriorElementBoundaries_global, int* exteriorElementBoundariesArray, int nInteriorElementBoundaries_global, @@ -620,7 +651,10 @@ namespace proteus double& dmom_w_ham_u, double& dmom_w_ham_v, double& dmom_w_ham_w, - double& rho) + double& rho, + double forcex, + double forcey, + double forcez) { double nu,mu,H_rho,d_rho,H_mu,d_mu,norm_n,nu_t0=0.0,nu_t1=0.0,nu_t; H_rho = (1.0-useVF)*smoothedHeaviside(eps_rho,phi) + useVF*fmin(1.0,fmax(0.0,vf)); @@ -826,7 +860,291 @@ namespace proteus dmom_v_ham_u =0.0; dmom_v_ham_v =0.0; } + + mom_u_source -= forcex; + mom_v_source -= forcey; + } + + int get_distance_to_ball(int n_balls,const double* ball_center, const double* ball_radius, const double x, const double y, const double z, double& distance) + { + distance = 1e10; + int index = -1; + double d_ball_i; + for (int i=0; i 0.0) ? 0.0 : nu * penalty; + double C_vol = (phi_s > 0.0) ? 0.0 : (alpha + beta * rel_vel_norm); + + C = (D_s * C_surf + (1.0 - H_s) * C_vol); + // force_x = dV * D_s * (p * phi_s_normal[0] - porosity * mu * (phi_s_normal[0] * grad_u[0] + phi_s_normal[1] * grad_u[1]) + C_surf * (u - u_s) * rho) + + // dV * (1.0 - H_s) * C_vol * (u - u_s) * rho; + // force_y = dV * D_s * (p * phi_s_normal[1] - porosity * mu * (phi_s_normal[0] * grad_v[0] + phi_s_normal[1] * grad_v[1]) + C_surf * (v - v_s) * rho) + + // dV * (1.0 - H_s) * C_vol * (v - v_s) * rho; +// force_x = dV*D_s*(p*fluid_outward_normal[0] - porosity*mu*(fluid_outward_normal[0]*grad_u[0] + fluid_outward_normal[1]*grad_u[1]) + C_surf*rel_vel_norm*(u-u_s)*rho) + dV*(1.0 - H_s)*C_vol*(u-u_s)*rho; +// force_y = dV*D_s*(p*fluid_outward_normal[1] - porosity*mu*(fluid_outward_normal[0]*grad_v[0] + fluid_outward_normal[1]*grad_v[1]) + C_surf*rel_vel_norm*(v-v_s)*rho) + dV*(1.0 - H_s)*C_vol*(v-v_s)*rho; +// force_x = dV * D_s * (p * fluid_outward_normal[0] +// -mu * (fluid_outward_normal[0] * 2* grad_u[0] + fluid_outward_normal[1] * (grad_u[1]+grad_v[0])) +// ); +// force_y = dV * D_s * (p * fluid_outward_normal[1] +// -mu * (fluid_outward_normal[0] * (grad_u[1]+grad_v[0]) + fluid_outward_normal[1] * 2* grad_v[1]) +// ); + force_p_x = dV * D_s * p * fluid_outward_normal[0]; + force_stress_x = dV * D_s * (-mu) * (fluid_outward_normal[0] * 2* grad_u[0] + fluid_outward_normal[1] * (grad_u[1]+grad_v[0])); + force_p_y = dV * D_s * p * fluid_outward_normal[1]; + force_stress_y = dV * D_s * (-mu) * (fluid_outward_normal[0] * (grad_u[1]+grad_v[0]) + fluid_outward_normal[1] * 2* grad_v[1]); + force_x = force_p_x + force_stress_x; + force_y = force_p_y + force_stress_y; + //always 3D for particle centroids + r_x = x - center[0]; + r_y = y - center[1]; + + if (element_owned) + { + particle_surfaceArea[i] += dV * D_s; + particle_netForces[i * 3 + 0] += force_x; + particle_netForces[i * 3 + 1] += force_y; + particle_netForces[(i+ nParticles)*3+0]+= force_stress_x; + particle_netForces[(i+2*nParticles)*3+0]+= force_p_x; + particle_netForces[(i+ nParticles)*3+1]+= force_stress_y; + particle_netForces[(i+2*nParticles)*3+1]+= force_p_y; + particle_netMoments[i * 3 + 2] += (r_x * force_y - r_y * force_x); + } + + // These should be done inside to make sure the correct velocity of different particles are used + //(1) + mom_u_source += C * (u - u_s); + mom_v_source += C * (v - v_s); + + dmom_u_source[0] += C; + dmom_v_source[1] += C; + + if (NONCONSERVATIVE_FORM > 0.0) + { + //(2) + mom_u_ham -= D_s * porosity * nu * (fluid_outward_normal[0] * grad_u[0] + fluid_outward_normal[1] * grad_u[1]); + dmom_u_ham_grad_u[0] -= D_s * porosity * nu * fluid_outward_normal[0]; + dmom_u_ham_grad_u[1] -= D_s * porosity * nu * fluid_outward_normal[1]; + + mom_v_ham -= D_s * porosity * nu * (fluid_outward_normal[0] * grad_v[0] + fluid_outward_normal[1] * grad_v[1]); + dmom_v_ham_grad_v[0] -= D_s * porosity * nu * fluid_outward_normal[0]; + dmom_v_ham_grad_v[1] -= D_s * porosity * nu * fluid_outward_normal[1]; + //(3) + mom_u_adv[0] += D_s * porosity * nu * fluid_outward_normal[0] * (u - u_s); + mom_u_adv[1] += D_s * porosity * nu * fluid_outward_normal[1] * (u - u_s); + dmom_u_adv_u[0] += D_s * porosity * nu * fluid_outward_normal[0]; + dmom_u_adv_u[1] += D_s * porosity * nu * fluid_outward_normal[1]; + + mom_v_adv[0] += D_s * porosity * nu * fluid_outward_normal[0] * (v - v_s); + mom_v_adv[1] += D_s * porosity * nu * fluid_outward_normal[1] * (v - v_s); + dmom_v_adv_v[0] += D_s * porosity * nu * fluid_outward_normal[0]; + dmom_v_adv_v[1] += D_s * porosity * nu * fluid_outward_normal[1]; + } + else + { + //(2) + mom_u_ham -= D_s * porosity * nu/rho * (fluid_outward_normal[0] * grad_u[0] + fluid_outward_normal[1] * grad_u[1]); + dmom_u_ham_grad_u[0] -= D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_u_ham_grad_u[1] -= D_s * porosity * nu/rho * fluid_outward_normal[1]; + + mom_v_ham -= D_s * porosity * nu/rho * (fluid_outward_normal[0] * grad_v[0] + fluid_outward_normal[1] * grad_v[1]); + dmom_v_ham_grad_v[0] -= D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_v_ham_grad_v[1] -= D_s * porosity * nu/rho * fluid_outward_normal[1]; + //(3) + mom_u_adv[0] += D_s * porosity * nu/rho * fluid_outward_normal[0] * (u - u_s); + mom_u_adv[1] += D_s * porosity * nu/rho * fluid_outward_normal[1] * (u - u_s); + dmom_u_adv_u[0] += D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_u_adv_u[1] += D_s * porosity * nu/rho * fluid_outward_normal[1]; + + mom_v_adv[0] += D_s * porosity * nu/rho * fluid_outward_normal[0] * (v - v_s); + mom_v_adv[1] += D_s * porosity * nu/rho * fluid_outward_normal[1] * (v - v_s); + dmom_v_adv_v[0] += D_s * porosity * nu/rho * fluid_outward_normal[0]; + dmom_v_adv_v[1] += D_s * porosity * nu/rho * fluid_outward_normal[1]; + //(4) +// mom_u_ham += D_s * porosity * (fluid_outward_normal[0] * u + fluid_outward_normal[1] * v)*u; +// mom_v_ham += D_s * porosity * (fluid_outward_normal[0] * u + fluid_outward_normal[1] * v)*v; +// dmom_u_ham_u += D_s * porosity * fluid_outward_normal[0] * u * 2.0; +// dmom_u_ham_v += D_s * porosity * fluid_outward_normal[1] * u; +// dmom_v_ham_u += D_s * porosity * fluid_outward_normal[0] * v; +// dmom_v_ham_v += D_s * porosity * fluid_outward_normal[1] * v * 2.0; + + } + //(6) +// mass_ham += D_s * porosity * (fluid_outward_normal[0] * u + fluid_outward_normal[1] * v); +// dmass_ham_u += D_s * porosity * fluid_outward_normal[0]; +// dmass_ham_v += D_s * porosity * fluid_outward_normal[1]; + //(7) +// mass_ham += C_surf * D_s * porosity * (fluid_outward_normal[0] * (u - u_s) + fluid_outward_normal[1] * (v - v_s)); +// dmass_ham_u += C_surf * D_s * porosity * fluid_outward_normal[0]; +// dmass_ham_v += C_surf * D_s * porosity * fluid_outward_normal[1]; + } + } //VRANS specific inline void updateDarcyForchheimerTerms_Ergun(const double NONCONSERVATIVE_FORM, @@ -1563,7 +1881,7 @@ namespace proteus double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABILIZATION, + double PRESSURE_PROJECTION_STABILIZATION, double* numerical_viscosity, //element double* mesh_trial_ref, @@ -1643,7 +1961,7 @@ namespace proteus double* w_dof, double* g, const double useVF, - double* q_rho, + double* q_rho, double* vf, double* phi, double* normal_phi, @@ -1722,9 +2040,29 @@ namespace proteus double* netForces_v, double* netMoments, double* velocityError, - double* velocityErrorNodal) + double* velocityErrorNodal, + double* forcex, + double* forcey, + double* forcez, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + double *particle_netForces, + double *particle_netMoments, + double *particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) { logEvent("Entered mprans 2D calculateResidual",6); + + const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element); // //loop over elements to compute volume integrals and load them into element and global residual @@ -1777,6 +2115,10 @@ namespace proteus dmass_adv_u[nSpace], dmass_adv_v[nSpace], dmass_adv_w[nSpace], + mass_ham=0.0, + dmass_ham_u=0.0, + dmass_ham_v=0.0, + dmass_ham_w=0.0, mom_u_adv[nSpace], dmom_u_adv_u[nSpace], dmom_u_adv_v[nSpace], @@ -2027,7 +2369,10 @@ namespace proteus dmom_w_ham_u, dmom_w_ham_v, dmom_w_ham_w, - q_rho[eN_k]); + q_rho[eN_k], + forcex[eN_k], + forcey[eN_k], + forcez[eN_k]); //VRANS mass_source = q_mass_source[eN_k]; //todo: decide if these should be lagged or not? @@ -2065,6 +2410,84 @@ namespace proteus dmom_v_source, dmom_w_source); + const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]); + if(nParticles > 0) + updateSolidParticleTerms(NONCONSERVATIVE_FORM, + eN < nElements_owned, + particle_nitsche, + dV, + nParticles, + nQuadraturePoints_global, +// &particle_signed_distances[eN_k], +// &particle_signed_distance_normals[eN_k_nSpace], +// particle_velocities, +// particle_centroids, + use_ball_as_particle, + ball_center, + ball_radius, + ball_velocity, + ball_angular_velocity, + porosity, + particle_penalty_constant/h_phi,//penalty, + particle_alpha, + particle_beta, + eps_rho, + eps_mu, + rho_0, + nu_0, + rho_1, + nu_1, + useVF, + vf[eN_k], + phi[eN_k], + x, + y, + z, + p, + u, + v, + w, + q_velocity_sge[eN_k_nSpace+0], + q_velocity_sge[eN_k_nSpace+1], + q_velocity_sge[eN_k_nSpace+1], + particle_eps, + grad_u, + grad_v, + grad_w, + mom_u_source, + mom_v_source, + mom_w_source, + dmom_u_source, + dmom_v_source, + dmom_w_source, + mom_u_adv, + mom_v_adv, + mom_w_adv, + dmom_u_adv_u, + dmom_v_adv_v, + dmom_w_adv_w, + mom_u_ham, + dmom_u_ham_grad_u, + dmom_u_ham_u, + dmom_u_ham_v, + dmom_u_ham_w, + mom_v_ham, + dmom_v_ham_grad_v, + dmom_v_ham_u, + dmom_v_ham_v, + dmom_v_ham_w, + mom_w_ham, + dmom_w_ham_grad_w, + dmom_w_ham_u, + dmom_w_ham_v, + dmom_w_ham_w, + mass_ham, + dmass_ham_u, + dmass_ham_v, + dmass_ham_w, + &particle_netForces[0], + &particle_netMoments[0], + &particle_surfaceArea[0]); //Turbulence closure model if (turbulenceClosureModel >= 3) { @@ -2288,8 +2711,9 @@ namespace proteus ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last[eN_k]/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace]); - elementResidual_p[i] += ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) + - DM*MOVING_DOMAIN*(ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) - + elementResidual_p[i] += ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) + + ck.Hamiltonian_weak(mass_ham, p_test_dV[i]) + + DM*MOVING_DOMAIN*(ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) - ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last[eN_k]/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])) + //VRANS @@ -2599,7 +3023,7 @@ namespace proteus //calculate the pde coefficients using the solution and the boundary values for the solution // double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.); //not interested in saving boundary eddy viscosity for now - double rho; + double rho; evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -2683,7 +3107,10 @@ namespace proteus dmom_w_ham_u_ext, dmom_w_ham_v_ext, dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -2767,7 +3194,10 @@ namespace proteus bc_dmom_w_ham_u_ext, bc_dmom_w_ham_v_ext, bc_dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); //Turbulence closure model if (turbulenceClosureModel >= 3) @@ -3111,7 +3541,7 @@ namespace proteus double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABILIZATION, + double PRESSURE_PROJECTION_STABILIZATION, //element double* mesh_trial_ref, double* mesh_grad_trial_ref, @@ -3275,8 +3705,22 @@ namespace proteus int* csrColumnOffsets_eb_w_v, int* csrColumnOffsets_eb_w_w, int* elementFlags, - int* boundaryFlags) + int* boundaryFlags, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) { + const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element); + std::valarray particle_surfaceArea(nParticles), particle_netForces(nParticles*3*3), particle_netMoments(nParticles*3); // //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian // @@ -3339,6 +3783,10 @@ namespace proteus dmass_adv_u[nSpace], dmass_adv_v[nSpace], dmass_adv_w[nSpace], + mass_ham=0.0, + dmass_ham_u=0.0, + dmass_ham_v=0.0, + dmass_ham_w=0.0, mom_u_adv[nSpace], dmom_u_adv_u[nSpace], dmom_u_adv_v[nSpace], @@ -3510,7 +3958,7 @@ namespace proteus //calculate pde coefficients and derivatives at quadrature points // double eddy_viscosity(0.);//not really interested in saving eddy_viscosity in jacobian - double rho; + double rho; evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -3594,7 +4042,10 @@ namespace proteus dmom_w_ham_u, dmom_w_ham_v, dmom_w_ham_w, - rho); + rho, + 0.0, + 0.0, + 0.0); //VRANS mass_source = q_mass_source[eN_k]; //todo: decide if these should be lagged or not @@ -3631,6 +4082,85 @@ namespace proteus dmom_u_source, dmom_v_source, dmom_w_source); + + const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]); + if(nParticles > 0) + updateSolidParticleTerms(NONCONSERVATIVE_FORM, + eN < nElements_owned, + particle_nitsche, + dV, + nParticles, + nQuadraturePoints_global, +// &particle_signed_distances[eN_k], +// &particle_signed_distance_normals[eN_k_nSpace], +// particle_velocities, +// particle_centroids, + use_ball_as_particle, + ball_center, + ball_radius, + ball_velocity, + ball_angular_velocity, + porosity, + particle_penalty_constant/h_phi,//penalty, + particle_alpha, + particle_beta, + eps_rho, + eps_mu, + rho_0, + nu_0, + rho_1, + nu_1, + useVF, + vf[eN_k], + phi[eN_k], + x, + y, + z, + p, + u, + v, + w, + q_velocity_sge[eN_k_nSpace+0], + q_velocity_sge[eN_k_nSpace+1], + q_velocity_sge[eN_k_nSpace+1], + particle_eps, + grad_u, + grad_v, + grad_w, + mom_u_source, + mom_v_source, + mom_w_source, + dmom_u_source, + dmom_v_source, + dmom_w_source, + mom_u_adv, + mom_v_adv, + mom_w_adv, + dmom_u_adv_u, + dmom_v_adv_v, + dmom_w_adv_w, + mom_u_ham, + dmom_u_ham_grad_u, + dmom_u_ham_u, + dmom_u_ham_v, + dmom_u_ham_w, + mom_v_ham, + dmom_v_ham_grad_v, + dmom_v_ham_u, + dmom_v_ham_v, + dmom_v_ham_w, + mom_w_ham, + dmom_w_ham_grad_w, + dmom_w_ham_u, + dmom_w_ham_v, + dmom_w_ham_w, + mass_ham, + dmass_ham_u, + dmass_ham_v, + dmass_ham_w, + &particle_netForces[0], + &particle_netMoments[0], + &particle_surfaceArea[0]); //Turbulence closure model if (turbulenceClosureModel >= 3) { @@ -3880,12 +4410,17 @@ namespace proteus (1-PRESSURE_PROJECTION_STABILIZATION)*ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_p[i]) + PRESSURE_PROJECTION_STABILIZATION*ck.pressureProjection_weak(mom_uu_diff_ten[1], p_trial_ref[k*nDOF_trial_element+j], 1./3., p_test_ref[k*nDOF_test_element +i],dV); - elementJacobian_p_u[i][j] += ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + + elementJacobian_p_u[i][j] += ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + + ck.MassJacobian_weak(dmass_ham_u,vel_trial_ref[k*nDOF_trial_element+j],p_test_dV[i]) + + (1-PRESSURE_PROJECTION_STABILIZATION)*ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_p[i]); - elementJacobian_p_v[i][j] += ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + + elementJacobian_p_v[i][j] += ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) + + ck.MassJacobian_weak(dmass_ham_v,vel_trial_ref[k*nDOF_trial_element+j],p_test_dV[i]) + + (1-PRESSURE_PROJECTION_STABILIZATION)*ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_p[i]); - elementJacobian_u_p[i][j] += ck.HamiltonianJacobian_weak(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) + + elementJacobian_u_p[i][j] += ck.HamiltonianJacobian_weak(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) + + MOMENTUM_SGE*VELOCITY_SGE*ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_u[i]); elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) + @@ -4208,7 +4743,7 @@ namespace proteus //calculate the internal and external trace of the pde coefficients // double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);//not interested in saving boundary eddy viscosity for now - double rho; + double rho; evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -4292,7 +4827,10 @@ namespace proteus dmom_w_ham_u_ext, dmom_w_ham_v_ext, dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); evaluateCoefficients(NONCONSERVATIVE_FORM, eps_rho, eps_mu, @@ -4376,7 +4914,10 @@ namespace proteus bc_dmom_w_ham_u_ext, bc_dmom_w_ham_v_ext, bc_dmom_w_ham_w_ext, - rho); + rho, + 0.0, + 0.0, + 0.0); //Turbulence closure model if (turbulenceClosureModel >= 3) { diff --git a/proteus/mprans/cRANS2P.pyx b/proteus/mprans/cRANS2P.pyx index 7bcf7a32ec..1cd0eaf187 100644 --- a/proteus/mprans/cRANS2P.pyx +++ b/proteus/mprans/cRANS2P.pyx @@ -11,7 +11,7 @@ cdef extern from "RANS2P.h" namespace "proteus": double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABLIZATION, + double PRESSURE_PROJECTION_STABLIZATION, double * numerical_viscosity, double * mesh_trial_ref, double * mesh_grad_trial_ref, @@ -167,12 +167,30 @@ cdef extern from "RANS2P.h" namespace "proteus": double * netForces_v, double * netMoments, double * velocityError, - double * velocityErrorNodal) + double * velocityErrorNodal, + double* forcex, + double* forcey, + double* forcez, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + double *particle_netForces, + double *particle_netMoments, + double *particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) void calculateJacobian(double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABLIZATION, + double PRESSURE_PROJECTION_STABLIZATION, double * mesh_trial_ref, double * mesh_grad_trial_ref, double * mesh_dof, @@ -333,7 +351,19 @@ cdef extern from "RANS2P.h" namespace "proteus": int * csrColumnOffsets_eb_w_v, int * csrColumnOffsets_eb_w_w, int * elementFlags, - int * boundaryFlags) + int * boundaryFlags, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) void calculateVelocityAverage(int nExteriorElementBoundaries_global, int * exteriorElementBoundariesArray, int nInteriorElementBoundaries_global, @@ -489,7 +519,7 @@ cdef class cRANS2P_base: double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABLIZATION, + double PRESSURE_PROJECTION_STABLIZATION, numpy.ndarray numerical_viscosity, numpy.ndarray mesh_trial_ref, numpy.ndarray mesh_grad_trial_ref, @@ -645,12 +675,30 @@ cdef class cRANS2P_base: numpy.ndarray netForces_v, numpy.ndarray netMoments, numpy.ndarray velocityError, - numpy.ndarray velocityErrorNodal): + numpy.ndarray velocityErrorNodal, + numpy.ndarray forcex, + numpy.ndarray forcey, + numpy.ndarray forcez, + int use_ball_as_particle, + numpy.ndarray ball_center, + numpy.ndarray ball_radius, + numpy.ndarray ball_velocity, + numpy.ndarray ball_angular_velocity, + int nParticles, + numpy.ndarray particle_netForces, + numpy.ndarray particle_netMoments, + numpy.ndarray particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant): self.thisptr.calculateResidual(NONCONSERVATIVE_FORM, MOMENTUM_SGE, PRESSURE_SGE, VELOCITY_SGE, - PRESSURE_PROJECTION_STABLIZATION, + PRESSURE_PROJECTION_STABLIZATION, < double * > numerical_viscosity.data, < double * > mesh_trial_ref.data, < double * > mesh_grad_trial_ref.data, @@ -806,14 +854,32 @@ cdef class cRANS2P_base: < double * > netForces_v.data, < double * > netMoments.data, < double * > velocityError.data, - < double * > velocityErrorNodal.data) + < double * > velocityErrorNodal.data, + < double * > forcex.data, + < double * > forcey.data, + < double * > forcez.data, + use_ball_as_particle, + < double * > ball_center.data, + < double * > ball_radius.data, + < double * > ball_velocity.data, + < double * > ball_angular_velocity.data, + nParticles, + < double * > particle_netForces.data, + < double * > particle_netMoments.data, + < double * > particle_surfaceArea.data, + nElements_owned, + particle_nitsche, + particle_epsFact, + particle_alpha, + particle_beta, + particle_penalty_constant) def calculateJacobian(self, double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABLIZATION, + double PRESSURE_PROJECTION_STABLIZATION, numpy.ndarray mesh_trial_ref, numpy.ndarray mesh_grad_trial_ref, numpy.ndarray mesh_dof, @@ -974,14 +1040,26 @@ cdef class cRANS2P_base: numpy.ndarray csrColumnOffsets_eb_w_v, numpy.ndarray csrColumnOffsets_eb_w_w, numpy.ndarray elementFlags, - numpy.ndarray boundaryFlags): + numpy.ndarray boundaryFlags, + int use_ball_as_particle, + numpy.ndarray ball_center, + numpy.ndarray ball_radius, + numpy.ndarray ball_velocity, + numpy.ndarray ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant): cdef numpy.ndarray rowptr, colind, globalJacobian_a (rowptr, colind, globalJacobian_a) = globalJacobian.getCSRrepresentation() self.thisptr.calculateJacobian(NONCONSERVATIVE_FORM, MOMENTUM_SGE, PRESSURE_SGE, VELOCITY_SGE, - PRESSURE_PROJECTION_STABLIZATION, + PRESSURE_PROJECTION_STABLIZATION, < double * > mesh_trial_ref.data, < double * > mesh_grad_trial_ref.data, < double * > mesh_dof.data, @@ -1142,7 +1220,19 @@ cdef class cRANS2P_base: < int * > csrColumnOffsets_eb_w_v.data, < int * > csrColumnOffsets_eb_w_w.data, < int * > elementFlags.data, - < int * > boundaryFlags.data) + < int * > boundaryFlags.data, + use_ball_as_particle, + < double * > ball_center.data, + < double * > ball_radius.data, + < double * > ball_velocity.data, + < double * > ball_angular_velocity.data, + nParticles, + nElements_owned, + particle_nitsche, + particle_epsFact, + particle_alpha, + particle_beta, + particle_penalty_constant) def calculateVelocityAverage(self, int nExteriorElementBoundaries_global, diff --git a/proteus/mprans/cRANS2P2D.pyx b/proteus/mprans/cRANS2P2D.pyx index ef059e2dd3..e4b3c8c57e 100644 --- a/proteus/mprans/cRANS2P2D.pyx +++ b/proteus/mprans/cRANS2P2D.pyx @@ -11,7 +11,7 @@ cdef extern from "RANS2P2D.h" namespace "proteus": double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABILIZATION, + double PRESSURE_PROJECTION_STABILIZATION, double* numerical_viscosity, double* mesh_trial_ref, double* mesh_grad_trial_ref, @@ -167,12 +167,30 @@ cdef extern from "RANS2P2D.h" namespace "proteus": double * netForces_v, double * netMoments, double * velocityError, - double * velocityErrorNodal) + double * velocityErrorNodal, + double* forcex, + double* forcey, + double* forcez, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + double *particle_netForces, + double *particle_netMoments, + double *particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) void calculateJacobian(double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, - double PRESSURE_PROJECTION_STABILIZATION, + double PRESSURE_PROJECTION_STABILIZATION, double* mesh_trial_ref, double* mesh_grad_trial_ref, double* mesh_dof, @@ -333,7 +351,19 @@ cdef extern from "RANS2P2D.h" namespace "proteus": int * csrColumnOffsets_eb_w_v, int * csrColumnOffsets_eb_w_w, int * elementFlags, - int * boundaryFlags) + int * boundaryFlags, + int use_ball_as_particle, + double* ball_center, + double* ball_radius, + double* ball_velocity, + double* ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant) void calculateVelocityAverage(int nExteriorElementBoundaries_global, int * exteriorElementBoundariesArray, int nInteriorElementBoundaries_global, @@ -640,12 +670,30 @@ cdef class cRANS2P2D_base: numpy.ndarray netForces_v, numpy.ndarray netMoments, numpy.ndarray velocityError, - numpy.ndarray velocityErrorNodal): + numpy.ndarray velocityErrorNodal, + numpy.ndarray forcex, + numpy.ndarray forcey, + numpy.ndarray forcez, + int use_ball_as_particle, + numpy.ndarray ball_center, + numpy.ndarray ball_radius, + numpy.ndarray ball_velocity, + numpy.ndarray ball_angular_velocity, + int nParticles, + numpy.ndarray particle_netForces, + numpy.ndarray particle_netMoments, + numpy.ndarray particle_surfaceArea, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant): self.thisptr.calculateResidual(NONCONSERVATIVE_FORM, MOMENTUM_SGE, PRESSURE_SGE, VELOCITY_SGE, - PRESSURE_PROJECTION_STABILIZATION, + PRESSURE_PROJECTION_STABILIZATION, numerical_viscosity.data, mesh_trial_ref.data, mesh_grad_trial_ref.data, @@ -801,7 +849,25 @@ cdef class cRANS2P2D_base: < double * > netForces_v.data, < double * > netMoments.data, < double * > velocityError.data, - < double * > velocityErrorNodal.data) + < double * > velocityErrorNodal.data, + < double * > forcex.data, + < double * > forcey.data, + < double * > forcez.data, + use_ball_as_particle, + < double * > ball_center.data, + < double * > ball_radius.data, + < double * > ball_velocity.data, + < double * > ball_angular_velocity.data, + nParticles, + < double * > particle_netForces.data, + < double * > particle_netMoments.data, + < double * > particle_surfaceArea.data, + nElements_owned, + particle_nitsche, + particle_epsFact, + particle_alpha, + particle_beta, + particle_penalty_constant) def calculateJacobian(self, double NONCONSERVATIVE_FORM, @@ -969,14 +1035,26 @@ cdef class cRANS2P2D_base: numpy.ndarray csrColumnOffsets_eb_w_v, numpy.ndarray csrColumnOffsets_eb_w_w, numpy.ndarray elementFlags, - numpy.ndarray boundaryFlags): + numpy.ndarray boundaryFlags, + int use_ball_as_particle, + numpy.ndarray ball_center, + numpy.ndarray ball_radius, + numpy.ndarray ball_velocity, + numpy.ndarray ball_angular_velocity, + int nParticles, + int nElements_owned, + double particle_nitsche, + double particle_epsFact, + double particle_alpha, + double particle_beta, + double particle_penalty_constant): cdef numpy.ndarray rowptr, colind, globalJacobian_a (rowptr, colind, globalJacobian_a) = globalJacobian.getCSRrepresentation() self.thisptr.calculateJacobian(NONCONSERVATIVE_FORM, MOMENTUM_SGE, PRESSURE_SGE, VELOCITY_SGE, - PRESSURE_PROJECTION_STABILIZATION, + PRESSURE_PROJECTION_STABILIZATION, mesh_trial_ref.data, mesh_grad_trial_ref.data, mesh_dof.data, @@ -1137,7 +1215,19 @@ cdef class cRANS2P2D_base: < int * > csrColumnOffsets_eb_w_v.data, < int * > csrColumnOffsets_eb_w_w.data, < int * > elementFlags.data, - < int * > boundaryFlags.data) + < int * > boundaryFlags.data, + use_ball_as_particle, + < double * > ball_center.data, + < double * > ball_radius.data, + < double * > ball_velocity.data, + < double * > ball_angular_velocity.data, + nParticles, + nElements_owned, + particle_nitsche, + particle_epsFact, + particle_alpha, + particle_beta, + particle_penalty_constant) def calculateVelocityAverage(self, int nExteriorElementBoundaries_global, diff --git a/proteus/tests/cylinder2D/README.md b/proteus/tests/cylinder2D/README.md index f2be5eb01b..d44cae9568 100644 --- a/proteus/tests/cylinder2D/README.md +++ b/proteus/tests/cylinder2D/README.md @@ -1,8 +1,18 @@ # Flow around cylinder benchmark 2D-2 -The problem is to simulate the flow (density=1, viscosity=0.001) around the cylinder (center=(0.2, 0.2), radius=0.05) in the rectangular domain [0,2.2]x[0,0.41]. +The problem is to simulate the flow (density=1, viscosity=0.001) +around the cylinder (center=(0.2, 0.2), radius=0.05) in the +rectangular domain [0,2.2]x[0,0.41]. ## References -1. Schäfer, M., Turek, S., Durst, F., Krause, E., & Rannacher, R. (1996). Benchmark computations of laminar flow around a cylinder. In Flow simulation with high-performance computers II (pp. 547-566). Vieweg+ Teubner Verlag. https://doi.org/10.1007/978-3-322-89849-4_39 -2. John, V. (2002). Higher order finite element methods and multigrid solvers in a benchmark problem for the 3D Navier–Stokes equations. International Journal for Numerical Methods in Fluids, 40(6), 775-798. https://doi.org/10.1002/fld.377 +1. Schäfer, M., Turek, S., Durst, F., Krause, E., & Rannacher, +R. (1996). Benchmark computations of laminar flow around a +cylinder. In Flow simulation with high-performance computers II +(pp. 547-566). Vieweg+ Teubner +Verlag. https://doi.org/10.1007/978-3-322-89849-4_39 + +2. John, V. (2002). Higher order finite element methods and multigrid +solvers in a benchmark problem for the 3D Navier–Stokes +equations. International Journal for Numerical Methods in Fluids, +40(6), 775-798. https://doi.org/10.1002/fld.377 diff --git a/proteus/tests/cylinder2D/conforming_rans2p/README.md b/proteus/tests/cylinder2D/conforming_rans2p/README.md index bd17527540..6c5987aad9 100644 --- a/proteus/tests/cylinder2D/conforming_rans2p/README.md +++ b/proteus/tests/cylinder2D/conforming_rans2p/README.md @@ -1,19 +1,22 @@ # Flow around cylinder benchmark 2D-2 -The problelm is to simulate the flow (density=1, viscoity=0.001) around the cylinder (center=(0.2, 0.2), radius=0.05) in the rectangle domain [0,2.2]x[0,0.41]. Please refer to [1,2] for mote details. It seems none of finite element references cited in [1] can be found. - - - -The whole algorithm is to solve a saddle point problem, where velocity and pressure are solved together in RANS2P.py. - +This test simulates the flow (density=1, viscoity=0.001) around the +cylinder (center=(0.2, 0.2), radius=0.05) in the rectangular domain +[0,2.2]x[0,0.41]. Please refer to [1,2] for mote details. +The whole algorithm is to solve a saddle point problem, where velocity +and pressure are solved together in RANS2P.py. ## Preprocessing -The file cylinder2d.py is called first, where *he*, the maximum size of edge, is passed into *symmetric2D* to generate the mesh. To increase the accuracy of the cylinder, there are at least 40 points on the boundary of the cylinder. -The final time of the problem is *T=8*. The boundary flags are ['left', 'right', 'top', 'bottom', 'obstacle']. -For boundary condition, the 'left' side is of inflow type, the 'right' side is of outflow type, and other sides are enforced non-slip boundary condition. - +The file cylinder2d.py is called first, where *he*, the maximum size +of edge, is passed into *symmetric2D* to generate the mesh. To +increase the accuracy of the cylinder, there are at least 40 points on +the boundary of the cylinder. The final time of the problem is +*T=8*. The boundary flags are ['left', 'right', 'top', 'bottom', +'obstacle']. For boundary condition, the 'left' side is of inflow +type, the 'right' side is of outflow type, and other sides are +enforced non-slip boundary condition. ## Running @@ -23,14 +26,21 @@ This problem can be run using the following command parun -v -l 5 cylinder_so.py -C "he=0.01" ``` - ## Postprocessing -Several physics parameters can be used to evaluate the computation [1]. One is the drag/lift coefficient. The drag/lift coefficient is sensitive and gives us the direction to find the bug relating the penalty coefficient in Nitsches' term. It is computed in *plot-lift-drag-force-from-pv.ipynb* by using the file *forceHistory_p.txt* and *forceHistory_p.txt*. +Several physical parameters can be used to evaluate the computation +[1]. One is the drag/lift coefficient. It is computed in +*plot-lift-drag-force-from-pv.ipynb* by using the file +*forceHistory_p.txt* and *forceHistory_p.txt*. ## Reference -1. Turek, Schaefer; Benchmark computations of laminar flow around cylinder; in Flow Simulation with High-Performance Computers II, Notes on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996 -2. John; Higher order Finite element methods and multigrid solvers in a benchmark problem for the 3D Navier-Stokes equations; Int. J. Numer. Meth. Fluids 2002; 40: 775-798 DOI:10.1002/d.377 +1. Turek, Schaefer; Benchmark computations of laminar flow around +cylinder; in Flow Simulation with High-Performance Computers II, Notes +on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996 + +2. John; Higher order Finite element methods and multigrid solvers in +a benchmark problem for the 3D Navier-Stokes equations; +Int. J. Numer. Meth. Fluids 2002; 40: 775-798 DOI:10.1002/d.377 diff --git a/proteus/tests/cylinder2D/conforming_rans3p/README.md b/proteus/tests/cylinder2D/conforming_rans3p/README.md index c114b9c8ec..bb71f8d015 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/README.md +++ b/proteus/tests/cylinder2D/conforming_rans3p/README.md @@ -1,19 +1,25 @@ # Flow around cylinder benchmark 2D-2 -The problelm is to simulate the flow (density=1, viscoity=0.001) around the cylinder (center=(0.2, 0.2), radius=0.05) in the rectangle domain [0,2.2]x[0,0.41]. Please refer to [1,2] for mote details. It seems none of finite element references cited in [1] can be found. - - - -The whole algorithm is a projection scheme, where the momentum equation is solved using RANS3P.py, the pressure increment is computed in PreInc.py, and the pressure at the next step is updated in Pressure.py. +This test simulates the flow (density=1, viscoity=0.001) around the +cylinder (center=(0.2, 0.2), radius=0.05) in the rectangle domain +[0,2.2]x[0,0.41]. Please refer to [1,2] for mote details. It seems +none of finite element references cited in [1] can be found. +The whole algorithm is a projection scheme, where the momentum +equation is solved using RANS3P.py, the pressure increment is computed +in PreInc.py, and the pressure at the next step is updated in +Pressure.py. ## Preprocessing -The file cylinder.py is called first, where *he*, the maximum size of edge, is passed into *symmetric2D* to generate the mesh. To increase the accuracy of the cylinder, there are at least 40 points on the boundary of the cylinder. -The final time of the problem is *T=8*. The boundary flags are ['left', 'right', 'top', 'bottom', 'obstacle']. -For boundary condition, the 'left' side is of inflow type, the 'right' side is of outflow type, and other sides are enforced non-slip boundary condition. - - +The file cylinder.py is called first, where *he*, the maximum size of +edge, is passed into *symmetric2D* to generate the mesh. To increase +the accuracy of the cylinder, there are at least 40 points on the +boundary of the cylinder. The final time of the problem is *T=8*. The +boundary flags are ['left', 'right', 'top', 'bottom', 'obstacle']. +For boundary condition, the 'left' side is of inflow type, the 'right' +side is of outflow type, and other sides are enforced non-slip +boundary condition. ## Running @@ -25,11 +31,19 @@ This problem can be run using the following command ## Postprocessing -Several physics parameters can be used to evaluate the computation [1]. One is the drag/lift coefficient. The drag/lift coefficient is sensitive and gives us the direction to find the bug relating the penalty coefficient in Nitsches' term. It is computed in *plot-lift-drag-force-from-pv.ipynb* by using the file *forceHistory_p.txt* and *forceHistory_p.txt*. - - +Several physics parameters can be used to evaluate the computation +[1]. One is the drag/lift coefficient. The drag/lift coefficient is +sensitive and gives us the direction to find the bug relating the +penalty coefficient in Nitsches' term. It is computed in +*plot-lift-drag-force-from-pv.ipynb* by using the file +*forceHistory_p.txt* and *forceHistory_p.txt*. ## Reference -1. Turek, Schaefer; Benchmark computations of laminar flow around cylinder; in Flow Simulation with High-Performance Computers II, Notes on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996 -2. John; Higher order Finite element methods and multigrid solvers in a benchmark problem for the 3D Navier-Stokes equations; Int. J. Numer. Meth. Fluids 2002; 40: 775-798 DOI:10.1002/d.377 \ No newline at end of file +1. Turek, Schaefer; Benchmark computations of laminar flow around +cylinder; in Flow Simulation with High-Performance Computers II, Notes +on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996 + +2. John; Higher order Finite element methods and multigrid solvers in +a benchmark problem for the 3D Navier-Stokes equations; +Int. J. Numer. Meth. Fluids 2002; 40: 775-798 DOI:10.1002/d.377 \ No newline at end of file diff --git a/proteus/tests/cylinder2D/ibm_method/README.md b/proteus/tests/cylinder2D/ibm_method/README.md index fd44d3cf67..35a67509a3 100644 --- a/proteus/tests/cylinder2D/ibm_method/README.md +++ b/proteus/tests/cylinder2D/ibm_method/README.md @@ -1,17 +1,23 @@ # Flow around cylinder benchmark 2D-2 -The problelm is to simulate the flow (density=1, viscoity=0.001) around the cylinder (center=(0.2, 0.2), radius=0.05) in the rectangle domain [0,2.2]x[0,0.41]. Please refer to [1,2] for mote details. It seems none of finite element references cited in [1] can be found. - - - -The whole algorithm is a projection scheme, where the momentum equation is solved using RANS3P.py, the pressure increment is computed in PreInc.py, and the pressure at the next step is updated in Pressure.py. +This test simulates the flow (density=1, viscosity=0.001) around the +cylinder (center=(0.2, 0.2), radius=0.05) in the rectangle domain +[0,2.2]x[0,0.41]. Please refer to [1,2] for mote details. It seems +none of finite element references cited in [1] can be found. +The whole algorithm is a projection scheme, where the momentum +equation is solved using RANS3P.py, the pressure increment is computed +in PreInc.py, and the pressure at the next step is updated in +Pressure.py. ## Preprocessing -The file cylinder.py is called first, where `he`, depending on `Refinement`, is used to control the mesh generation. The final time of the problem is T=8. The boundary flags are ['left', 'right', 'top', 'bottom', 'obstacle']. For boundary condition, the 'left' side is of inflow type, the 'right' side is of outflow type, and other sides are enforced non-slip boundary condition. - - +The file cylinder.py is called first, where `he`, depending on +`Refinement`, is used to control the mesh generation. The final time +of the problem is T=8. The boundary flags are ['left', 'right', 'top', +'bottom', 'obstacle']. For boundary condition, the 'left' side is of +inflow type, the 'right' side is of outflow type, and other sides are +enforced non-slip boundary condition. ## Running @@ -20,14 +26,21 @@ This problem can be run using the following command parun -v -l 5 cylinder_so.py -C "Refinement=7" ``` - ## Postprocessing -Several physics parameters can be used to evaluate the computation [1]. One is the drag/lift coefficient. The drag/lift coefficient is sensitive and gives us the direction to find the bug relating the penalty coefficient in Nitsches' term. It is computed in *plot_lift_drag_force_from_particle.ipynb* by using the file *particle_forceHistory.txt*. - - +Several physics parameters can be used to evaluate the computation +[1]. One is the drag/lift coefficient. The drag/lift coefficient is +sensitive and gives us the direction to find the bug relating the +penalty coefficient in Nitsches' term. It is computed in +*plot_lift_drag_force_from_particle.ipynb* by using the file +*particle_forceHistory.txt*. ## Reference -1. Turek, Schaefer; Benchmark computations of laminar flow around cylinder; in Flow Simulation with High-Performance Computers II, Notes on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996 -2. John; Higher order Finite element methods and multigrid solvers in a benchmark problem for the 3D Navier-Stokes equations; Int. J. Numer. Meth. Fluids 2002; 40: 775-798 DOI:10.1002/d.377 \ No newline at end of file +1. Turek, Schaefer; Benchmark computations of laminar flow around +cylinder; in Flow Simulation with High-Performance Computers II, Notes +on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996 + +2. John; Higher order Finite element methods and multigrid solvers in +a benchmark problem for the 3D Navier-Stokes equations; +Int. J. Numer. Meth. Fluids 2002; 40: 775-798 DOI:10.1002/d.377 \ No newline at end of file diff --git a/proteus/tests/cylinder2D/ibm_rans2p/README.md b/proteus/tests/cylinder2D/ibm_rans2p/README.md new file mode 100644 index 0000000000..fb6495b882 --- /dev/null +++ b/proteus/tests/cylinder2D/ibm_rans2p/README.md @@ -0,0 +1,45 @@ +# Flow around cylinder benchmark 2D-2 + +This test simulates the flow (density=1, viscosity=0.001) around the +cylinder (center=(0.2, 0.2), radius=0.05) in the rectangular domain +[0,2.2]x[0,0.41]. Please refer to [1,2] for mote details. + + +## Preprocessing + +The file cylinder.py is called first, where `he`, depending on +`Refinement`, is used to control the mesh generation. The final time +of the problem is T=4. The boundary flags are ['left', 'right', 'top', +'bottom']. For boundary condition, the 'left' side is of inflow type, +the 'right' side is of outflow type, and other sides are enforced +no-slip boundary condition. The no-slip boundary condition over the +cylinder is implemented using the domain integral with the help of +delta function with support on the surface of the cylinder in the same +spirit as Nitsche method. + +## Running + +This problem can be run using the following command +```bash + parun -v -l 5 cylinder_so.py -C "Refinement=5" -D ibm_p1 +``` + + +## Postprocessing + +Several physical parameters can be used to evaluate the computation +[1]. One is the drag/lift coefficient. It is computed in +*plot_lift_drag_force_from_particle.ipynb* by using the file +*particle_forceHistory.txt*. + + + +## Reference + +1. Turek, Schaefer; Benchmark computations of laminar flow around +cylinder; in Flow Simulation with High-Performance Computers II, Notes +on Numerical Fluid Mechanics 52, 547-566, Vieweg 1996 + +2. John; Higher order Finite element methods and multigrid solvers in +a benchmark problem for the 3D Navier-Stokes equations; +Int. J. Numer. Meth. Fluids 2002; 40: 775-798 DOI:10.1002/d.377 diff --git a/proteus/tests/cylinder2D/ibm_rans2p/__init__.py b/proteus/tests/cylinder2D/ibm_rans2p/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/proteus/tests/cylinder2D/ibm_rans2p/cylinder.py b/proteus/tests/cylinder2D/ibm_rans2p/cylinder.py new file mode 100644 index 0000000000..f9d90c7868 --- /dev/null +++ b/proteus/tests/cylinder2D/ibm_rans2p/cylinder.py @@ -0,0 +1,431 @@ +from math import * +import proteus.MeshTools +from proteus import Domain +from proteus.default_n import * +from proteus.Profiling import logEvent +from proteus.Gauges import PointGauges, LineGauges, LineIntegralGauges + +from proteus import Context + +ct = Context.Options([ + ("T", 3.0, "Time interval [0, T]"), + ("onlySaveFinalSolution",False,"Only save the final solution"), + ("parallel",False,"Use parallel or not"), + ("dt_fixed",0.005,"fixed time step"), + ################################## + ("Refinement", 3, "refinement"), + ("StrongDirichlet", True,"weak or strong"), + ("use_sbm", 0,"use sbm instead of imb"), + ("spaceOrder", 1,"FE space for velocity"), + ("timeOrder", 2,"1=be or 2=vbdf"),#both works, but 2 gives better cd-cl + ("use_supg", 1,"Use supg or not"), + ("nonconservative", 1,"0=conservative"), + ("forceStrongDirichlet",False,"strong or weak"), + ################################## + ("use_ball_as_particle",1,"1 or 0 == use ball or particle"), + ("isHotStart",False,"Use hotStart or not"), +], mutable=True) + + +#=============================================================================== +# supg +#=============================================================================== +use_supg = ct.use_supg########## +#=============================================================================== +# use sbm or not +#=============================================================================== +USE_SBM=ct.use_sbm +#=============================================================================== +# Parameters +#=============================================================================== +# Discretization -- input options +#Refinement = 20#45min on a single core for spaceOrder=1, useHex=False +Refinement = ct.Refinement +sedimentDynamics=False +genMesh = True +movingDomain = False +applyRedistancing = True +useOldPETSc = False +useSuperlu = True + + +parallel = ct.parallel +if parallel: + usePETSc = True + useSuperlu=False +else: + usePETSc = False +parallelPartitioningType = MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 1 + +if ct.timeOrder==2: + timeDiscretization = 'vbdf'#vbdf'#'vbdf' # 'vbdf', 'be', 'flcbdf' +else: + timeDiscretization = 'be' + +spaceOrder = ct.spaceOrder +pspaceOrder = 1 +useHex = False +useRBLES = 0.0 +useMetrics = 1.0 +applyCorrection = True +useVF = 0.0 +useOnlyVF = False +useRANS = 0 # 0 -- None + # 1 -- K-Epsilon + # 2 -- K-Omega +openTop=True +fl_H = 0.41 +# Input checks +if spaceOrder not in [1, 2]: + print "INVALID: spaceOrder" + spaceOrder + sys.exit() + +if useRBLES not in [0.0, 1.0]: + print "INVALID: useRBLES" + useRBLES + sys.exit() + +if useMetrics not in [0.0, 1.0]: + print "INVALID: useMetrics" + sys.exit() + + +weak_bc_penalty_constant = 100.0 +nLevels = 1 +# parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element +parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.node +nLayersOfOverlapForParallel = 0 + +#=============================================================================== +# FE space +#=============================================================================== +# Discretization +nd = 2 + +if spaceOrder == 1: + hFactor = 1.0 + if useHex: + basis = C0_AffineLinearOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd, 2) + elementBoundaryQuadrature = CubeGaussQuadrature(nd - 1, 2) + else: + basis = C0_AffineLinearOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd, 3) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd - 1, 3) +elif spaceOrder == 2: + hFactor = 0.5 + if useHex: + basis = C0_AffineLagrangeOnCubeWithNodalBasis + elementQuadrature = CubeGaussQuadrature(nd, 4) + elementBoundaryQuadrature = CubeGaussQuadrature(nd - 1, 4) + else: + basis = C0_AffineQuadraticOnSimplexWithNodalBasis + elementQuadrature = SimplexGaussQuadrature(nd, 5) + elementBoundaryQuadrature = SimplexGaussQuadrature(nd - 1, 5) + +if pspaceOrder == 1: + if useHex: + pbasis = C0_AffineLinearOnCubeWithNodalBasis + else: + pbasis = C0_AffineLinearOnSimplexWithNodalBasis +elif pspaceOrder == 2: + if useHex: + pbasis = C0_AffineLagrangeOnCubeWithNodalBasis + else: + pbasis = C0_AffineQuadraticOnSimplexWithNodalBasis + + +#=============================================================================== +# PointGauges +#=============================================================================== +# new style PointGauges +# pointGauges = PointGauges(gauges = ((('u', 'v'), ((0.5, 0.5, 0), (1, 0.5, 0))), (('p',), ((0.5, 0.5, 0),))), +# activeTime=(0, 0.5), +# sampleRate=0, +# fileName='combined_gauge_0_0.5_sample_all.csv') +pressure_pointGauges = PointGauges(gauges = ((('p',), ((0.15, 0.2, 0.0),(0.25, 0.2, 0.0),),),), + fileName='front_back_pressure_sample_all.csv') +# lineGauges = LineGauges(gaugeEndpoints={'lineGauge_xtoH=0.825': ((0.495, 0.0, 0.0), (0.495, 1.8, 0.0))}, linePoints=20) +# #'lineGauge_x/H=1.653':((0.99,0.0,0.0),(0.99,1.8,0.0)) +# lineGauges_phi = LineGauges_phi(lineGauges.endpoints, linePoints=20) + + +#=============================================================================== +# Domain and mesh +#=============================================================================== +#L = (0.584,0.350) +L = (2.2, 0.41) +he = 1.0/2**Refinement +# he*=0.5 +# he*=0.5 +#he*=0.5 +#he*=0.5 +#he*=0.5 + +structured = False + +if useHex: + nnx = 4 * Refinement + 1 + nny = 2 * Refinement + 1 + hex = True + domain = Domain.RectangularDomain(L) +else: + boundaries = ['bottom', 'right', 'top', 'left', 'front', 'back'] + boundaryTags = dict([(key, i + 1) for (i, key) in enumerate(boundaries)]) + if structured: + nnx = 4 * Refinement + nny = 2 * Refinement + else: + vertices = [[0.0, 0.0], #0 + [L[0], 0.0], #1 + [L[0], L[1]], #2 + [0.0, L[1]], #3 + [0.2-0.16,L[1]*0.2], + [0.2-0.16,L[1]*0.8], + [0.2+0.3,L[1]*0.8], + [0.2+0.3,L[1]*0.2], + # the following are set for refining the mesh + [0.2-0.06,0.2-0.06], + [0.2-0.06,0.2+0.06], + [0.2+0.06,0.2+0.06], + [0.2+0.06,0.2-0.06]] + + + + vertexFlags = [boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['top'], + boundaryTags['top'], + # the interior vertices should be flaged to 0 + 0, 0, 0, 0, + 0, 0, 0, 0 ] + + segments = [[0, 1], + [1, 2], + [2, 3], + [3, 0], + #Interior segments + [4, 5], + [5, 6], + [6, 7], + [7,4], + [8,9], + [9,10], + [10,11], + [11,8]] + segmentFlags = [boundaryTags['bottom'], + boundaryTags['right'], + boundaryTags['top'], + boundaryTags['left'], + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0] + + regions = [[0.95*L[0], 0.2],[0.2-0.15,0.2],[0.2,0.2]] + regionFlags = [1,2,3] + regionConstraints=[0.5*he**2,0.5*(he/2.0)**2,0.5*(he/6.0)**2] + # for gaugeName,gaugeCoordinates in pointGauges.locations.iteritems(): + # vertices.append(gaugeCoordinates) + # vertexFlags.append(pointGauges.flags[gaugeName]) + + # for gaugeName, gaugeLines in lineGauges.linepoints.iteritems(): + # for gaugeCoordinates in gaugeLines: + # vertices.append(gaugeCoordinates) + # vertexFlags.append(lineGauges.flags[gaugeName]) + domain = Domain.PlanarStraightLineGraphDomain(vertices=vertices, + vertexFlags=vertexFlags, + segments=segments, + segmentFlags=segmentFlags, + regions=regions, + regionFlags=regionFlags, + regionConstraints=regionConstraints) + #go ahead and add a boundary tags member + domain.boundaryTags = boundaryTags + domain.writePoly("mesh") + domain.writePLY("mesh") + domain.writeAsymptote("mesh") + #triangleOptions = "VApq30Dena%8.8f" % ((he ** 2) / 2.0,) + triangleOptions = "VApq30Dena" + +logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions, domain.polyfile + ".poly")) +#=============================================================================== +# Time stepping +#=============================================================================== +T=ct.T +dt_fixed = ct.dt_fixed#0.03 +dt_init = 0.5*dt_fixed#min(0.1*dt_fixed,0.001) +if ct.isHotStart: + dt_init = dt_fixed +else: + dt_init = min(dt_init,0.5*dt_fixed) +runCFL=0.33 +nDTout = int(round(T/dt_fixed)) +if dt_init