From e0b37261e30c9c3ab00565cecd4054c243f33879 Mon Sep 17 00:00:00 2001 From: Chris Kees Date: Thu, 28 Mar 2019 09:08:32 -0500 Subject: [PATCH] Refactoring ibm/sbm/unfitted methods (#905) * Yy/debug pseudopenalty 2p (#904) * Cekees/cleanup ib (#906) * cleaned up some tests and removed some old code * moved notebooks to lfs * fix bug in particle_velocity and particle_signed_distance_normal sizing * added pressure integral to cut cell approach * use eN_k_3d instead of eN_k_nSpace for 2d particles (#915) * fix added mass global residual assembly; enable Pres.LevelModel * updated added mass comparison values * update some test files and turn off cut cell integration * disable MeshAdapt * turn off optimization for travis compile * reduce size of some tests * update optimization level to 1 --- .gitattributes | 5 + .travis.yml | 2 +- proteus/mprans/AddedMass.h | 20 +- proteus/mprans/Pres.py | 2 + proteus/mprans/PresInit.py | 18 - proteus/mprans/RANS2P.h | 336 ++- proteus/mprans/RANS2P.py | 155 +- proteus/mprans/RANS2P2D.h | 1851 ++++++++++------- proteus/mprans/RANS3PF.h | 91 +- proteus/mprans/RANS3PF.py | 31 +- proteus/mprans/RANS3PF2D.h | 487 ++--- proteus/mprans/cRANS2P.pyx | 114 +- proteus/mprans/cRANS2P2D.pyx | 114 +- proteus/mprans/cRANS3PF.pyx | 3 +- proteus/tests/AddedMass/test_addedmass2D.py | 4 +- proteus/tests/AddedMass/test_addedmass3D.py | 4 +- .../tests/CLSVOF/with_RANS3PF/pressure_p.py | 2 +- .../tests/ProjScheme_with_EV/pressure_p.py | 2 +- .../comparison_files/T1_rans2p.h5 | 4 +- .../conforming_rans2p/cylinder2d.py | 10 +- .../conforming_rans2p/cylinder_so.py | 1 - .../plot_lift_drag_from_pv_files.ipynb | 1746 +--------------- .../twp_navier_stokes_cylinder_2d_p.py | 9 +- .../comparison_files/T1P1.h5 | 4 +- .../cylinder2D/conforming_rans3p/cylinder.py | 8 +- .../conforming_rans3p/cylinder_so.py | 2 +- .../plot_lift_drag_from_pv_files.ipynb | 1746 +--------------- .../conforming_rans3p/pressureInitial_p.py | 2 + .../conforming_rans3p/pressureincrement_p.py | 10 +- .../conforming_rans3p/twp_navier_stokes_p.py | 16 +- .../ibm_method/comparison_files/T1_rans3p.h5 | 4 +- .../tests/cylinder2D/ibm_method/cylinder.py | 208 +- .../cylinder2D/ibm_method/cylinder_so.py | 2 +- .../plot_lift_drag_force_from_particle.ipynb | 1737 +--------------- .../tests/cylinder2D/ibm_method/pressure_p.py | 2 + .../ibm_method/twp_navier_stokes_p.py | 7 +- .../comparison_files/T1_ibm_rans2p.h5 | 4 +- .../tests/cylinder2D/ibm_rans2p/cylinder.py | 215 +- .../cylinder2D/ibm_rans2p/cylinder_so.py | 2 +- .../plot_lift_drag_force_from_particle.ipynb | 921 +------- .../ibm_rans2p/twp_navier_stokes_n.py | 6 +- .../ibm_rans2p/twp_navier_stokes_p.py | 54 +- .../comparison_files/T001_P1_sbm_3Dmesh.h5 | 4 +- .../tests/cylinder2D/sbm_3Dmesh/cylinder.py | 2 +- .../tests/cylinder2D/sbm_3Dmesh/pressure_p.py | 1 + .../comparison_files/T1_sbm_rans3p.h5 | 4 +- .../tests/cylinder2D/sbm_method/cylinder.py | 202 +- .../cylinder2D/sbm_method/cylinder_so.py | 2 +- .../plot_lift_drag_force_from_particle.ipynb | 1711 +-------------- .../sbm_method/pressureInitial_n.py | 5 +- .../sbm_method/pressureInitial_p.py | 5 +- .../tests/cylinder2D/sbm_method/pressure_n.py | 5 +- .../tests/cylinder2D/sbm_method/pressure_p.py | 6 +- .../sbm_method/pressureincrement_n.py | 5 +- .../sbm_method/pressureincrement_p.py | 7 +- .../sbm_method/twp_navier_stokes_n.py | 8 +- .../sbm_method/twp_navier_stokes_p.py | 5 +- .../elastoplastic_expected.h5 | 4 +- .../tests/griffiths_lane_6/re_gl_6_3d_p.py | 2 +- .../griffiths_lane_6/richards_expected.h5 | 4 +- .../tests/griffiths_lane_6/sm_gl_6_3d_p.py | 2 +- .../griffiths_lane_6/test_griffiths_lane6.py | 4 +- 62 files changed, 2589 insertions(+), 9360 deletions(-) diff --git a/.gitattributes b/.gitattributes index d97c6d8906..37c020e5d3 100644 --- a/.gitattributes +++ b/.gitattributes @@ -19,3 +19,8 @@ proteus/tests/HotStart_3P/comparison_files/T01P1_hotstart.h5 filter=lfs diff=lfs 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 +proteus/tests/cylinder2D/conforming_rans2p/plot_lift_drag_from_pv_files.ipynb filter=lfs diff=lfs merge=lfs -text +proteus/tests/cylinder2D/conforming_rans3p/plot_lift_drag_from_pv_files.ipynb filter=lfs diff=lfs merge=lfs -text +proteus/tests/cylinder2D/ibm_method/plot_lift_drag_force_from_particle.ipynb filter=lfs diff=lfs merge=lfs -text +proteus/tests/cylinder2D/ibm_rans2p/plot_lift_drag_force_from_particle.ipynb filter=lfs diff=lfs merge=lfs -text +proteus/tests/cylinder2D/sbm_method/plot_lift_drag_force_from_particle.ipynb filter=lfs diff=lfs merge=lfs -text diff --git a/.travis.yml b/.travis.yml index 897c91df9c..e90177f66e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -42,7 +42,7 @@ install: - export PATHSAVE=$PATH - export PATH=$PWD/linux2/bin:$PATH - export LD_LIBRARY_PATH=$PWD/linux2/lib:$LD_LIBRARY_PATH -- PROTEUS_OPT="-O2 -UNDEBUG -w" FC=gfortran CC=mpicc CXX=mpicxx make develop N=2 +- PROTEUS_OPT="-w -O1 -UNDEBUG" FC=gfortran CC=mpicc CXX=mpicxx make develop N=2 - export SSL_CERT_DIR=/etc/ssl/certs - pip install matplotlib diff --git a/proteus/mprans/AddedMass.h b/proteus/mprans/AddedMass.h index 0a7ecdc0f1..8d0a4433b0 100644 --- a/proteus/mprans/AddedMass.h +++ b/proteus/mprans/AddedMass.h @@ -511,16 +511,16 @@ namespace proteus Aij[36*eBMT+added_mass_i+6*4] += (rz*px-rx*pz)*dS; Aij[36*eBMT+added_mass_i+6*5] += (rx*py-ry*px)*dS; } - // - //update the element and global residual storage - // - for (int i=0;i0.0)?1:0; //calculate eddy viscosity switch (turbulenceClosureModel) { @@ -1119,21 +1163,40 @@ namespace proteus vz = ball_velocity[3*I + 2] + angular_cross_position[2]; } + void get_acceleration_to_ith_ball(int n_balls,const double* ball_center, const double* ball_radius, + const double* ball_velocity, const double* ball_angular_velocity, + const double* ball_center_acceleration, const double* ball_angular_acceleration, + int I, + const double x, const double y, const double z, + double& ax, double& ay, double& az) + { +#ifdef USE_CYLINDER_AS_PARTICLE + double position[3]={x-ball_center[3*I + 0],y-ball_center[3*I + 1],0.0}; +#else + double position[3]={x-ball_center[3*I + 0],y-ball_center[3*I + 1],z-ball_center[3*I + 2]}; +#endif + double angular_cross_position[3]; + get_cross_product(&ball_angular_velocity[3*I + 0],position,angular_cross_position); + //YY-incomplete + } inline void updateSolidParticleTerms(const double NONCONSERVATIVE_FORM, bool element_owned, const double particle_nitsche, const double dV, const int nParticles, const int sd_offset, -// double *particle_signed_distances, -// double *particle_signed_distance_normals, -// double *particle_velocities, -// double *particle_centroids, + double *particle_signed_distances, + double *particle_signed_distance_normals, + double *particle_velocities, + double *particle_centroids, const int use_ball_as_particle, const double* ball_center, const double* ball_radius, const double* ball_velocity, const double* ball_angular_velocity, + const double* ball_center_acceleration, + const double* ball_angular_acceleration, + const double* ball_density, const double porosity, //VRANS specific const double penalty, const double alpha, @@ -1194,7 +1257,8 @@ namespace proteus double &dmass_ham_w, double *particle_netForces, double *particle_netMoments, - double *particle_surfaceArea) + double *particle_surfaceArea, + const int use_pseudo_penalty) { double C, rho, mu, nu, H_mu, uc, duc_du, duc_dv, duc_dw, H_s, D_s, phi_s, u_s, v_s, w_s; double force_x, force_y, force_z, r_x, r_y, r_z, force_p_x, force_p_y, force_p_z, force_stress_x, force_stress_y, force_stress_z; @@ -1223,15 +1287,16 @@ namespace proteus } else { -// phi_s = particle_signed_distances[i * sd_offset]; -// phi_s_normal[0] = particle_signed_distance_normals[i * sd_offset * nSpace + 0]; -// phi_s_normal[1] = particle_signed_distance_normals[i * sd_offset * nSpace + 1]; -// vel[0] = particle_velocities[i * sd_offset * nSpace + 0]; -// vel[1] = particle_velocities[i * sd_offset * nSpace + 1]; -// center[0] = particle_centroids[3*i+0]; -// center[1] = particle_centroids[3*i+1]; - std::cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!YY: use_ball_as_particle should be 1"< 0) + phi_solid[eN_k] = distance_to_solids[eN_k]; // //calculate pde coefficients at quadrature points // @@ -2799,6 +2909,16 @@ namespace proteus kappa_phi[eN_k], //VRANS porosity, + phi_solid[eN_k],//distance to solid + p_old, + u_old, + v_old, + w_old, + grad_p_old, + grad_u_old, + grad_v_old, + grad_w_old, + use_pseudo_penalty, // p, grad_p, @@ -2912,15 +3032,18 @@ namespace proteus dV, nParticles, nQuadraturePoints_global, -// &particle_signed_distances[eN_k], -// &particle_signed_distance_normals[eN_k_nSpace], -// particle_velocities, -// particle_centroids, + &particle_signed_distances[eN_k], + &particle_signed_distance_normals[eN_k_3d], + &particle_velocities[eN_k_3d], + particle_centroids, use_ball_as_particle, ball_center, ball_radius, ball_velocity, ball_angular_velocity, + ball_center_acceleration, + ball_angular_acceleration, + ball_density, porosity, particle_penalty_constant/h_phi,//penalty, particle_alpha, @@ -2981,7 +3104,8 @@ namespace proteus dmass_ham_w, &particle_netForces[0], &particle_netMoments[0], - &particle_surfaceArea[0]); + &particle_surfaceArea[0], + use_pseudo_penalty); //Turbulence closure model if (turbulenceClosureModel >= 3) { @@ -3104,6 +3228,23 @@ namespace proteus dmom_w_acc_w, mom_w_acc_t, dmom_w_acc_w_t); + if(use_pseudo_penalty > 0 && phi_solid[eN_k]<0.0)//Do not have to change Jacobian + { + double distance,vx,vy,vz; + int index_ball = get_distance_to_ball(nParticles, ball_center, ball_radius,x,y,z,distance); + get_velocity_to_ith_ball(nParticles,ball_center,ball_radius,ball_velocity,ball_angular_velocity,index_ball,x,y,z,vx,vy,vz); + mom_u_acc_t = alphaBDF*(mom_u_acc - vx); + mom_v_acc_t = alphaBDF*(mom_v_acc - vy); + mom_w_acc_t = alphaBDF*(mom_w_acc - vz); + }else if(use_pseudo_penalty == -1 && phi_solid[eN_k]<0.0)//no derivative term inside the solid; Has to change Jacobian + { + mom_u_acc_t = 0.0; + mom_v_acc_t = 0.0; + mom_w_acc_t = 0.0; + dmom_u_acc_u= 0.0; + dmom_v_acc_v= 0.0; + dmom_w_acc_w= 0.0; + } if (NONCONSERVATIVE_FORM > 0.0) { mom_u_acc_t *= dmom_u_acc_u; @@ -3390,6 +3531,8 @@ namespace proteus grad_u_ext[nSpace], grad_v_ext[nSpace], grad_w_ext[nSpace], + p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0, + grad_p_old[nSpace],grad_u_old[nSpace],grad_v_old[nSpace],grad_w_old[nSpace], mom_u_acc_ext=0.0, dmom_u_acc_u_ext=0.0, mom_v_acc_ext=0.0, @@ -3586,10 +3729,18 @@ namespace proteus ck_v.valFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_ext); ck_v.valFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_ext); ck_v.valFromDOF(w_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],w_ext); + ck.valFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_trace_ref[ebN_local_kb*nDOF_test_element],p_old); + ck_v.valFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_old); + ck_v.valFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_old); + ck_v.valFromDOF(w_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],w_old); ck.gradFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext); ck_v.gradFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_ext); ck_v.gradFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_ext); ck_v.gradFromDOF(w_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_w_ext); + ck.gradFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_old); + ck_v.gradFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_old); + ck_v.gradFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_old); + ck_v.gradFromDOF(w_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_w_old); //precalculate test function products with integration weights for (int j=0;j particle_surfaceArea(nParticles), particle_netForces(nParticles*3*3), particle_netMoments(nParticles*3); @@ -4546,12 +4740,15 @@ namespace proteus { int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point eN_k_nSpace = eN_k*nSpace, + eN_k_3d = eN_k*3, eN_nDOF_trial_element = eN*nDOF_trial_element, //index to a vector at a quadrature point eN_nDOF_v_trial_element = eN*nDOF_v_trial_element; //index to a vector at a quadrature point //declare local storage register double p=0.0,u=0.0,v=0.0,w=0.0, grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace], + p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0, + grad_p_old[nSpace],grad_u_old[nSpace],grad_v_old[nSpace],grad_w_old[nSpace], mom_u_acc=0.0, dmom_u_acc_u=0.0, mom_v_acc=0.0, @@ -4704,11 +4901,19 @@ namespace proteus ck_v.valFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],u); ck_v.valFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],v); ck_v.valFromDOF(w_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],w); + ck.valFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_ref[k*nDOF_trial_element],p_old); + ck_v.valFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],u_old); + ck_v.valFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],v_old); + ck_v.valFromDOF(w_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],w_old); //get the solution gradients ck.gradFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial,grad_p); ck_v.gradFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_u); ck_v.gradFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_v); ck_v.gradFromDOF(w_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_w); + ck.gradFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial,grad_p_old); + ck_v.gradFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_u_old); + ck_v.gradFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_v_old); + ck_v.gradFromDOF(w_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_w_old); //precalculate test function products with integration weights for (int j=0;j= 3) { @@ -5058,6 +5277,15 @@ namespace proteus dmom_w_acc_w, mom_w_acc_t, dmom_w_acc_w_t); + if(use_pseudo_penalty == -1 && phi_solid[eN_k]<0.0)//no derivative term inside the solid; Has to change Jacobian + { + mom_u_acc_t = 0.0; + mom_v_acc_t = 0.0; + mom_w_acc_t = 0.0; + dmom_u_acc_u = 0.0; + dmom_v_acc_v = 0.0; + dmom_w_acc_w = 0.0; + } if (NONCONSERVATIVE_FORM > 0.0) { mom_u_acc_t *= dmom_u_acc_u; @@ -5470,6 +5698,8 @@ namespace proteus grad_u_ext[nSpace], grad_v_ext[nSpace], grad_w_ext[nSpace], + p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0, + grad_p_old[nSpace],grad_u_old[nSpace],grad_v_old[nSpace],grad_w_old[nSpace], mom_u_acc_ext=0.0, dmom_u_acc_u_ext=0.0, mom_v_acc_ext=0.0, @@ -5677,10 +5907,18 @@ namespace proteus ck_v.valFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_ext); ck_v.valFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_ext); ck_v.valFromDOF(w_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],w_ext); + ck.valFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_trace_ref[ebN_local_kb*nDOF_test_element],p_old); + ck_v.valFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_old); + ck_v.valFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_old); + ck_v.valFromDOF(w_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],w_old); ck.gradFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext); ck_v.gradFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_ext); ck_v.gradFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_ext); ck_v.gradFromDOF(w_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_w_ext); + ck.gradFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_old); + ck_v.gradFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_old); + ck_v.gradFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_old); + ck_v.gradFromDOF(w_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_w_old); //precalculate test function products with integration weights for (int j=0;j 0 and self.use_ball_as_particle == 0: + self.phi_s[:] = 1e10 + self.phisField[:] = 1e10 + self.ebq_global_phi_s[:] = 1e10 + self.ebq_global_grad_phi_s[:] = 1e10 + self.ebq_particle_velocity_s[:] = 1e10 + + for i in range(self.nParticles): + vel = lambda x: self.particle_velocityList[i](t, x) + sdf = lambda x: self.particle_sdfList[i](t, x) + + for j in range(self.mesh.nodeArray.shape[0]): + sdf_at_node, sdNormals = sdf(self.mesh.nodeArray[j, :]) + if (sdf_at_node < self.phi_s[j]): + self.phi_s[j] = sdf_at_node + for eN in range(self.model.q['x'].shape[0]): + for k in range(self.model.q['x'].shape[1]): + self.particle_signed_distances[i, eN, k], self.particle_signed_distance_normals[i, eN, k] = sdf(self.model.q['x'][eN, k]) + self.particle_velocities[i, eN, k] = vel(self.model.q['x'][eN, k]) + if (abs(self.particle_signed_distances[i, eN, k]) < abs(self.phisField[eN, k])): + self.phisField[eN, k] = self.particle_signed_distances[i, eN, k] + for ebN in range(self.model.ebq_global['x'].shape[0]): + for kb in range(self.model.ebq_global['x'].shape[1]): + sdf_ebN_kb,sdNormals = sdf(self.model.ebq_global['x'][ebN,kb]) + if ( sdf_ebN_kb < self.ebq_global_phi_s[ebN,kb]): + self.ebq_global_phi_s[ebN,kb]=sdf_ebN_kb + self.ebq_global_grad_phi_s[ebN,kb,:]=sdNormals + self.ebq_particle_velocity_s[ebN,kb,:] = vel(self.model.ebq_global['x'][ebN,kb]) # if self.comm.isMaster(): # print "wettedAreas" # print self.wettedAreas[:] @@ -684,7 +755,10 @@ def postStep(self, t, firstStep=False): `self.netForces_p[:,:]` + "\nForces_v\n" + `self.netForces_v[:,:]`) + self.timeHistory.write("%21.16e\n" % (t,)) + self.timeHistory.flush() self.wettedAreaHistory.write("%21.16e\n" % (self.wettedAreas[-1],)) + self.wettedAreaHistory.flush() self.forceHistory_p.write("%21.16e %21.16e %21.16e\n" % tuple(self.netForces_p[-1, :])) self.forceHistory_p.flush() self.forceHistory_v.write("%21.16e %21.16e %21.16e\n" % tuple(self.netForces_v[-1, :])) @@ -694,10 +768,13 @@ def postStep(self, t, firstStep=False): 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_pforceHistory.write("%21.16e %21.16e %21.16e\n" % tuple(self.particle_netForces[0+self.nParticles, :])) + self.particle_pforceHistory.flush() + self.particle_vforceHistory.write("%21.16e %21.16e %21.16e\n" % tuple(self.particle_netForces[0+2*self.nParticles, :])) + self.particle_vforceHistory.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): nCalls = 0 @@ -1323,6 +1400,13 @@ def __init__(self, self.q[('force', 2)] = numpy.zeros( (self.mesh.nElements_global, self.nQuadraturePoints_element), 'd') + if options is not None: + try: + self.chrono_model = options.chrono_model + except AttributeError: + logEvent('WARNING: did not find chrono model') + pass + def getResidual(self, u, r): """ Calculate the element residuals and add in to the global residual @@ -1458,6 +1542,10 @@ def getResidual(self, u, r): self.u[1].dof, self.u[2].dof, self.u[3].dof, + self.coefficients.p_old_dof, + self.coefficients.u_old_dof, + self.coefficients.v_old_dof, + self.coefficients.w_old_dof, self.coefficients.g, self.coefficients.useVF, self.q['rho'], @@ -1554,6 +1642,16 @@ def getResidual(self, u, r): self.coefficients.ball_radius, self.coefficients.ball_velocity, self.coefficients.ball_angular_velocity, + self.coefficients.ball_center_acceleration, + self.coefficients.ball_angular_acceleration, + self.coefficients.ball_density, + self.coefficients.particle_signed_distances, + self.coefficients.particle_signed_distance_normals, + self.coefficients.particle_velocities, + self.coefficients.particle_centroids, + self.coefficients.ebq_global_phi_s, + self.coefficients.ebq_global_grad_phi_s, + self.coefficients.ebq_particle_velocity_s, self.coefficients.nParticles, self.coefficients.particle_netForces, self.coefficients.particle_netMoments, @@ -1563,7 +1661,10 @@ def getResidual(self, u, r): self.coefficients.particle_epsFact, self.coefficients.particle_alpha, self.coefficients.particle_beta, - self.coefficients.particle_penalty_constant) + self.coefficients.particle_penalty_constant, + self.coefficients.phi_s, + self.coefficients.phisField, + self.coefficients.use_pseudo_penalty) for i in range(self.coefficients.netForces_p.shape[0]): self.coefficients.wettedAreas[i] = globalSum(self.coefficients.wettedAreas[i]) for I in range(3): @@ -1580,25 +1681,26 @@ 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 + + from mpi4py import MPI + comm = MPI.COMM_WORLD + + comm.Allreduce(self.coefficients.wettedAreas.copy(),self.coefficients.wettedAreas) + comm.Allreduce(self.coefficients.netForces_p.copy(),self.coefficients.netForces_p) + comm.Allreduce(self.coefficients.netForces_v.copy(),self.coefficients.netForces_v) + comm.Allreduce(self.coefficients.netMoments.copy(),self.coefficients.netMoments) + + comm.Allreduce(self.coefficients.particle_netForces.copy(),self.coefficients.particle_netForces) + comm.Allreduce(self.coefficients.particle_netMoments.copy(),self.coefficients.particle_netMoments) + comm.Allreduce(self.coefficients.particle_surfaceArea.copy(),self.coefficients.particle_surfaceArea) + 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=" + repr(i)+ " force " + repr(self.coefficients.particle_netForces[i])) logEvent("particle i=" + repr(i)+ " moment " + repr(self.coefficients.particle_netMoments[i])) logEvent("particle i=" + repr(i)+ " surfaceArea " + repr(self.coefficients.particle_surfaceArea[i])) logEvent("particle i=" + repr(i)+ " stress force " + repr(self.coefficients.particle_netForces[i+self.coefficients.nParticles])) logEvent("particle i=" + repr(i)+ " pressure force " + repr(self.coefficients.particle_netForces[i+2*self.coefficients.nParticles])) - if self.forceStrongConditions: for cj in range(len(self.dirichletConditionsForceDOF)): for dofN, g in list(self.dirichletConditionsForceDOF[cj].DOFBoundaryConditionsDict.items()): @@ -1722,6 +1824,10 @@ def getJacobian(self, jacobian): self.u[1].dof, self.u[2].dof, self.u[3].dof, + self.coefficients.p_old_dof, + self.coefficients.u_old_dof, + self.coefficients.v_old_dof, + self.coefficients.w_old_dof, self.coefficients.g, self.coefficients.useVF, self.coefficients.q_vf, @@ -1825,13 +1931,26 @@ def getJacobian(self, jacobian): self.coefficients.ball_radius, self.coefficients.ball_velocity, self.coefficients.ball_angular_velocity, + self.coefficients.ball_center_acceleration, + self.coefficients.ball_angular_acceleration, + self.coefficients.ball_density, + self.coefficients.particle_signed_distances, + self.coefficients.particle_signed_distance_normals, + self.coefficients.particle_velocities, + self.coefficients.particle_centroids, + self.coefficients.ebq_global_phi_s, + self.coefficients.ebq_global_grad_phi_s, + self.coefficients.ebq_particle_velocity_s, + self.coefficients.phi_s, + self.coefficients.phisField, 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) + self.coefficients.particle_penalty_constant, + self.coefficients.use_pseudo_penalty) 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 8f8b3cc115..c64c9b4c93 100644 --- a/proteus/mprans/RANS2P2D.h +++ b/proteus/mprans/RANS2P2D.h @@ -20,41 +20,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 MOVING_DOMAIN,//0 or 1 - 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, + double *mesh_trial_ref, + double *mesh_grad_trial_ref, + double *mesh_dof, + double *mesh_velocity_dof, + double MOVING_DOMAIN, //0 or 1 + 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, @@ -75,119 +75,133 @@ 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, + 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, - int* p_l2g, - int* vel_l2g, + double *q_eddy_viscosity, + double *q_eddy_viscosity_last, + double *ebqe_eddy_viscosity, + double *ebqe_eddy_viscosity_last, + int *p_l2g, + int *vel_l2g, int* rp_l2g, int* rvel_l2g, - double* p_dof, - double* u_dof, - double* v_dof, - double* w_dof, - double* g, + double *p_dof, + double *u_dof, + double *v_dof, + double *w_dof, + double *p_old_dof, + double *u_old_dof, + double *v_old_dof, + double *w_old_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, - 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 *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, + 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, + double* ball_center_acceleration, + double* ball_angular_acceleration, + double* ball_density, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, + double* ebq_global_phi_s, + double* ebq_global_grad_phi_s, + double* ebq_particle_velocity_s, int nParticles, double *particle_netForces, double *particle_netMoments, @@ -197,46 +211,49 @@ namespace proteus double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant)=0; + double particle_penalty_constant, + double* phi_solid_nodes, + double* distance_to_solids, + const int use_pseudo_penalty) = 0; virtual void calculateJacobian(double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, double VELOCITY_SGE, double PRESSURE_PROJECTION_STABILIZATION, //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, double useRBLES, @@ -256,128 +273,145 @@ namespace proteus double C_dg, 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_last, - 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, + double *q_eddy_viscosity_last, + double *ebqe_eddy_viscosity_last, + int *p_l2g, + int *vel_l2g, + double *p_dof, double *u_dof, double *v_dof, double *w_dof, + double *p_old_dof, + double *u_old_dof, + double *v_old_dof, + double *w_old_dof, + double *g, const double useVF, - double* vf, - double* phi, - double* normal_phi, - double* kappa_phi, - 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_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* csrRowIndeces_p_p,int* csrColumnOffsets_p_p, - int* csrRowIndeces_p_u,int* csrColumnOffsets_p_u, - int* csrRowIndeces_p_v,int* csrColumnOffsets_p_v, - int* csrRowIndeces_p_w,int* csrColumnOffsets_p_w, - int* csrRowIndeces_u_p,int* csrColumnOffsets_u_p, - int* csrRowIndeces_u_u,int* csrColumnOffsets_u_u, - int* csrRowIndeces_u_v,int* csrColumnOffsets_u_v, - int* csrRowIndeces_u_w,int* csrColumnOffsets_u_w, - int* csrRowIndeces_v_p,int* csrColumnOffsets_v_p, - int* csrRowIndeces_v_u,int* csrColumnOffsets_v_u, - int* csrRowIndeces_v_v,int* csrColumnOffsets_v_v, - int* csrRowIndeces_v_w,int* csrColumnOffsets_v_w, - int* csrRowIndeces_w_p,int* csrColumnOffsets_w_p, - int* csrRowIndeces_w_u,int* csrColumnOffsets_w_u, - int* csrRowIndeces_w_v,int* csrColumnOffsets_w_v, - int* csrRowIndeces_w_w,int* csrColumnOffsets_w_w, - double* globalJacobian, + double *vf, + double *phi, + double *normal_phi, + double *kappa_phi, + 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_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 *csrRowIndeces_p_p, int *csrColumnOffsets_p_p, + int *csrRowIndeces_p_u, int *csrColumnOffsets_p_u, + int *csrRowIndeces_p_v, int *csrColumnOffsets_p_v, + int *csrRowIndeces_p_w, int *csrColumnOffsets_p_w, + int *csrRowIndeces_u_p, int *csrColumnOffsets_u_p, + int *csrRowIndeces_u_u, int *csrColumnOffsets_u_u, + int *csrRowIndeces_u_v, int *csrColumnOffsets_u_v, + int *csrRowIndeces_u_w, int *csrColumnOffsets_u_w, + int *csrRowIndeces_v_p, int *csrColumnOffsets_v_p, + int *csrRowIndeces_v_u, int *csrColumnOffsets_v_u, + int *csrRowIndeces_v_v, int *csrColumnOffsets_v_v, + int *csrRowIndeces_v_w, int *csrColumnOffsets_v_w, + int *csrRowIndeces_w_p, int *csrColumnOffsets_w_p, + int *csrRowIndeces_w_u, int *csrColumnOffsets_w_u, + int *csrRowIndeces_w_v, int *csrColumnOffsets_w_v, + int *csrRowIndeces_w_w, int *csrColumnOffsets_w_w, + double *globalJacobian, 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, - int* csrColumnOffsets_eb_p_p, - int* csrColumnOffsets_eb_p_u, - int* csrColumnOffsets_eb_p_v, - int* csrColumnOffsets_eb_p_w, - int* csrColumnOffsets_eb_u_p, - int* csrColumnOffsets_eb_u_u, - int* csrColumnOffsets_eb_u_v, - int* csrColumnOffsets_eb_u_w, - int* csrColumnOffsets_eb_v_p, - int* csrColumnOffsets_eb_v_u, - int* csrColumnOffsets_eb_v_v, - int* csrColumnOffsets_eb_v_w, - int* csrColumnOffsets_eb_w_p, - int* csrColumnOffsets_eb_w_u, - int* csrColumnOffsets_eb_w_v, - int* csrColumnOffsets_eb_w_w, - int* elementFlags, - int* boundaryFlags, - int use_ball_as_particle, - double* ball_center, - double* ball_radius, - double* ball_velocity, - double* ball_angular_velocity, + 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, + int *csrColumnOffsets_eb_p_p, + int *csrColumnOffsets_eb_p_u, + int *csrColumnOffsets_eb_p_v, + int *csrColumnOffsets_eb_p_w, + int *csrColumnOffsets_eb_u_p, + int *csrColumnOffsets_eb_u_u, + int *csrColumnOffsets_eb_u_v, + int *csrColumnOffsets_eb_u_w, + int *csrColumnOffsets_eb_v_p, + int *csrColumnOffsets_eb_v_u, + int *csrColumnOffsets_eb_v_v, + int *csrColumnOffsets_eb_v_w, + int *csrColumnOffsets_eb_w_p, + int *csrColumnOffsets_eb_w_u, + int *csrColumnOffsets_eb_w_v, + int *csrColumnOffsets_eb_w_w, + int *elementFlags, + int *boundaryFlags, + int use_ball_as_particle, + double *ball_center, + double *ball_radius, + double *ball_velocity, + double *ball_angular_velocity, + double* ball_center_acceleration, + double* ball_angular_acceleration, + double* ball_density, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, + double* ebq_global_phi_s, + double* ebq_global_grad_phi_s, + double* ebq_particle_velocity_s, + double* phi_solid_nodes, + double* distance_to_solids, int nParticles, int nElements_owned, double particle_nitsche, double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant)=0; + double particle_penalty_constant, + const int use_pseudo_penalty) = 0; virtual void calculateVelocityAverage(int nExteriorElementBoundaries_global, int* exteriorElementBoundariesArray, int nInteriorElementBoundaries_global, @@ -550,7 +584,13 @@ namespace proteus else if (phi==0.0) H=0.5; else + { H = 0.5*(1.0 + phi/eps + sin(M_PI*phi/eps)/M_PI); + //phi /= eps; + //H = 0.5 + (9.0*phi-5.0*phi*phi*phi)/8.0; + //H = 0.5 + (45*phi-50*phi*phi*phi+21*phi*phi*phi*phi*phi)/32.0; + //H = 0.5 + (105*phi-175*phi*phi*phi+147*phi*phi*phi*phi*phi-45*phi*phi*phi*phi*phi*phi*phi)/64.0; + } return H; } @@ -583,8 +623,15 @@ namespace proteus else if (phi < -eps) d=0.0; else + { d = 0.5*(1.0 + cos(M_PI*phi/eps))/eps; + //phi /= eps; + //d = 3.0/4.0*(1.0-phi*phi); + //d = 15.0/32.0*(3.0-10.0*phi*phi+7.0*phi*phi*phi*phi); + //d = 315.0/512.0*(3.0-20.0*phi*phi+42.0*phi*phi*phi*phi-36.0*phi*phi*phi*phi*phi*phi+11.0*phi*phi*phi*phi*phi*phi*phi*phi); + } return d; + } inline @@ -606,6 +653,16 @@ namespace proteus const double n[nSpace], const double& kappa, const double porosity,//VRANS specific + const double phi_solid, + const double p_old, + const double u_old, + const double v_old, + const double w_old, + const double grad_p_old[nSpace], + const double grad_u_old[nSpace], + const double grad_v_old[nSpace], + const double grad_w_old[nSpace], + const int use_pseudo_penalty, const double& p, const double grad_p[nSpace], const double grad_u[nSpace], @@ -679,7 +736,7 @@ namespace proteus d_rho = (1.0-useVF)*smoothedDirac(eps_rho,phi); H_mu = (1.0-useVF)*smoothedHeaviside(eps_mu,phi) + useVF*fmin(1.0,fmax(0.0,vf)); d_mu = (1.0-useVF)*smoothedDirac(eps_mu,phi); - + const int in_fluid = (phi_solid>0.0)?1:0; //calculate eddy viscosity switch (turbulenceClosureModel) { @@ -773,7 +830,11 @@ namespace proteus norm_n = sqrt(n[0]*n[0]+n[1]*n[1]); mom_u_source = -porosity*rho*g[0];// - porosity*d_mu*sigma*kappa*n[0]; mom_v_source = -porosity*rho*g[1];// - porosity*d_mu*sigma*kappa*n[1]; - + if(use_pseudo_penalty>0) + { + mom_u_source = -in_fluid*porosity*rho*g[0]; + mom_v_source = -in_fluid*porosity*rho*g[1]; + } //u momentum Hamiltonian (pressure) mom_u_ham = porosity*grad_p[0]; @@ -783,23 +844,76 @@ namespace proteus //v momentum Hamiltonian (pressure) mom_v_ham = porosity*grad_p[1]; dmom_v_ham_grad_p[0]=0.0; - dmom_v_ham_grad_p[1]=porosity; + dmom_v_ham_grad_p[1] = porosity; + if (use_pseudo_penalty <= 0) + { + //u momentum Hamiltonian (advection) + mom_u_ham += rho * porosity * (u * grad_u[0] + v * grad_u[1]); + dmom_u_ham_grad_u[0] = rho * porosity * u; + dmom_u_ham_grad_u[1] = rho * porosity * v; + dmom_u_ham_u = rho * porosity * grad_u[0]; + dmom_u_ham_v = rho * porosity * grad_u[1]; + + //v momentum Hamiltonian (advection) + mom_v_ham += rho * porosity * (u * grad_v[0] + v * grad_v[1]); + dmom_v_ham_grad_v[0] = rho * porosity * u; + dmom_v_ham_grad_v[1] = rho * porosity * v; + dmom_v_ham_u = rho * porosity * grad_v[0]; + dmom_v_ham_v = rho * porosity * grad_v[1]; + } + else if (use_pseudo_penalty == 1) + { - //u momentum Hamiltonian (advection) - mom_u_ham += rho*porosity*(u*grad_u[0]+v*grad_u[1]); - dmom_u_ham_grad_u[0]=rho*porosity*u; - dmom_u_ham_grad_u[1]=rho*porosity*v; - dmom_u_ham_u =rho*porosity*grad_u[0]; - dmom_u_ham_v =rho*porosity*grad_u[1]; + //u momentum Hamiltonian (advection) + mom_u_ham += in_fluid * rho * porosity * (u * grad_u[0] + v * grad_u[1]); + dmom_u_ham_grad_u[0] = in_fluid * rho * porosity * u; + dmom_u_ham_grad_u[1] = in_fluid * rho * porosity * v; + dmom_u_ham_u = in_fluid * rho * porosity * grad_u[0]; + dmom_u_ham_v = in_fluid * rho * porosity * grad_u[1]; + + //v momentum Hamiltonian (advection) + mom_v_ham += in_fluid * rho * porosity * (u * grad_v[0] + v * grad_v[1]); + dmom_v_ham_grad_v[0] = in_fluid * rho * porosity * u; + dmom_v_ham_grad_v[1] = in_fluid * rho * porosity * v; + dmom_v_ham_u = in_fluid * rho * porosity * grad_v[0]; + dmom_v_ham_v = in_fluid * rho * porosity * grad_v[1]; + } + else if (use_pseudo_penalty == 2) + { - //v momentum Hamiltonian (advection) - mom_v_ham += rho*porosity*(u*grad_v[0]+v*grad_v[1]); - dmom_v_ham_grad_v[0]=rho*porosity*u; - dmom_v_ham_grad_v[1]=rho*porosity*v; - dmom_v_ham_u =rho*porosity*grad_v[0]; - dmom_v_ham_v =rho*porosity*grad_v[1]; + //u momentum Hamiltonian (advection) + mom_u_ham += in_fluid * rho * porosity * (u_old * grad_u_old[0] + v * grad_u_old[1]); + dmom_u_ham_grad_u[0] = 0.0; + dmom_u_ham_grad_u[1] = 0.0; + dmom_u_ham_u = 0.0; + dmom_u_ham_v = 0.0; + + //v momentum Hamiltonian (advection) + mom_v_ham += in_fluid * rho * porosity * (u_old * grad_v_old[0] + v_old * grad_v_old[1]); + dmom_v_ham_grad_v[0] = 0.0; + dmom_v_ham_grad_v[1] = 0.0; + dmom_v_ham_u = 0.0; + dmom_v_ham_v = 0.0; + } + else if (use_pseudo_penalty == 3) + { + + //u momentum Hamiltonian (advection) + mom_u_ham += 0.0; + dmom_u_ham_grad_u[0] = 0.0; + dmom_u_ham_grad_u[1] = 0.0; + dmom_u_ham_u = 0.0; + dmom_u_ham_v = 0.0; + + //v momentum Hamiltonian (advection) + mom_v_ham += 0.0; + dmom_v_ham_grad_v[0] = 0.0; + dmom_v_ham_grad_v[1] = 0.0; + dmom_v_ham_u = 0.0; + dmom_v_ham_v = 0.0; + } } - else + else { //u momentum accumulation mom_u_acc=porosity*u; @@ -839,6 +953,68 @@ namespace proteus dmom_v_adv_v[0]=porosity*u; dmom_v_adv_v[1]=2.0*porosity*v; + if(use_pseudo_penalty==1)///////////treat nonlinear term implicitly but only inside solid domain + { + //u momentum advective flux + mom_u_adv[0]=in_fluid*porosity*u*u; + mom_u_adv[1]=in_fluid*porosity*u*v; + + dmom_u_adv_u[0]=in_fluid*2.0*porosity*u; + dmom_u_adv_u[1]=in_fluid*porosity*v; + + dmom_u_adv_v[0]=0.0; + dmom_u_adv_v[1]=in_fluid*porosity*u; + + //v momentum advective_flux + mom_v_adv[0]=in_fluid*porosity*v*u; + mom_v_adv[1]=in_fluid*porosity*v*v; + + dmom_v_adv_u[0]=in_fluid*porosity*v; + dmom_v_adv_u[1]=0.0; + + dmom_v_adv_v[0]=in_fluid*porosity*u; + dmom_v_adv_v[1]=in_fluid*2.0*porosity*v; + }else if(use_pseudo_penalty==2){///////////treat nonlinear term explicitly + //u momentum advective flux + mom_u_adv[0]=in_fluid*porosity*u_old*u_old; + mom_u_adv[1]=in_fluid*porosity*u_old*v_old; + + dmom_u_adv_u[0]=0.0; + dmom_u_adv_u[1]=0.0; + + dmom_u_adv_v[0]=0.0; + dmom_u_adv_v[1]=0.0; + + //v momentum advective_flux + mom_v_adv[0]=in_fluid*porosity*v_old*u_old; + mom_v_adv[1]=in_fluid*porosity*v_old*v_old; + + dmom_v_adv_u[0]=0.0; + dmom_v_adv_u[1]=0.0; + + dmom_v_adv_v[0]=0.0; + dmom_v_adv_v[1]=0.0; + }else if(use_pseudo_penalty==3){///////////treat nonlinear term explicitly + //u momentum advective flux + mom_u_adv[0]=0.0; + mom_u_adv[1]=0.0; + + dmom_u_adv_u[0]=0.0; + dmom_u_adv_u[1]=0.0; + + dmom_u_adv_v[0]=0.0; + dmom_u_adv_v[1]=0.0; + + //v momentum advective_flux + mom_v_adv[0]=0.0; + mom_v_adv[1]=0.0; + + dmom_v_adv_u[0]=0.0; + dmom_v_adv_u[1]=0.0; + + dmom_v_adv_v[0]=0.0; + dmom_v_adv_v[1]=0.0; + } //u momentum diffusion tensor mom_uu_diff_ten[0] = 2.0*porosity*nu; mom_uu_diff_ten[1] = porosity*nu; @@ -853,9 +1029,13 @@ namespace proteus //momentum sources norm_n = sqrt(n[0]*n[0]+n[1]*n[1]);//+n[2]*n[2]); - mom_u_source = -porosity*g[0];// - porosity*d_mu*sigma*kappa*n[0]/(rho*(norm_n+1.0e-8)); - mom_v_source = -porosity*g[1];// - porosity*d_mu*sigma*kappa*n[1]/(rho*(norm_n+1.0e-8)); - + mom_u_source = -porosity*g[0]; + mom_v_source = -porosity*g[1]; + if(use_pseudo_penalty>0) + { + mom_u_source = -in_fluid*porosity*g[0]; + mom_v_source = -in_fluid*porosity*g[1]; + } //u momentum Hamiltonian (pressure) mom_u_ham = porosity*grad_p[0]/rho; dmom_u_ham_grad_p[0]=porosity/rho; @@ -930,21 +1110,36 @@ namespace proteus vx = ball_velocity[3*I + 0] - ball_angular_velocity[3*I + 2]*(y-ball_center[3*I + 1]); vy = ball_velocity[3*I + 1] + ball_angular_velocity[3*I + 2]*(x-ball_center[3*I + 0]); } + void get_acceleration_to_ith_ball(int n_balls,const double* ball_center, const double* ball_radius, + const double* ball_velocity, const double* ball_angular_velocity, + const double* ball_center_acceleration, const double* ball_angular_acceleration, + int I, + const double x, const double y, const double z, + double& ax, double& ay) + { + ax = ball_center_acceleration[3*I + 0] - ball_angular_acceleration[3*I + 2]*(y-ball_center[3*I + 1]) + - ball_angular_acceleration[3*I + 2]*ball_angular_acceleration[3*I + 2]*(x-ball_center[3*I + 0]); + ay = ball_center_acceleration[3*I + 1] + ball_angular_acceleration[3*I + 2]*(x-ball_center[3*I + 0]) + - ball_angular_acceleration[3*I + 2]*ball_angular_acceleration[3*I + 2]*(y-ball_center[3*I + 1]); + } inline void updateSolidParticleTerms(const double NONCONSERVATIVE_FORM, bool element_owned, const double particle_nitsche, const double dV, const int nParticles, const int sd_offset, -// double *particle_signed_distances, -// double *particle_signed_distance_normals, -// double *particle_velocities, -// double *particle_centroids, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, const int use_ball_as_particle, const double* ball_center, const double* ball_radius, const double* ball_velocity, const double* ball_angular_velocity, + const double* ball_center_acceleration, + const double* ball_angular_acceleration, + const double* ball_density, const double porosity, //VRANS specific const double penalty, const double alpha, @@ -1005,13 +1200,15 @@ namespace proteus double &dmass_ham_w, double *particle_netForces, double *particle_netMoments, - double *particle_surfaceArea) + double *particle_surfaceArea, + const int use_pseudo_penalty) { double C, rho, mu, nu, H_mu, uc, duc_du, duc_dv, duc_dw, H_s, D_s, phi_s, u_s, v_s, w_s; double force_x, force_y, r_x, r_y, force_p_x, force_p_y, force_stress_x, force_stress_y; double phi_s_normal[2]={0.0}; double fluid_outward_normal[2]; - double vel[2]; + double vel[2]={0.0}; + double acceleration[2]={0.0}; double center[2]; H_mu = (1.0 - useVF) * smoothedHeaviside(eps_mu, phi) + useVF * fmin(1.0, fmax(0.0, vf)); nu = nu_0 * (1.0 - H_mu) + nu_1 * H_mu; @@ -1019,9 +1216,9 @@ namespace proteus mu = rho_0 * nu_0 * (1.0 - H_mu) + rho_1 * nu_1 * H_mu; C = 0.0; for (int i = 0; i < nParticles; i++) - { + { if(use_ball_as_particle==1) - { + { get_distance_to_ith_ball(nParticles,ball_center,ball_radius,i,x,y,z,phi_s); get_normal_to_ith_ball(nParticles,ball_center,ball_radius,i,x,y,z,phi_s_normal[0],phi_s_normal[1]); get_velocity_to_ith_ball(nParticles,ball_center,ball_radius, @@ -1030,19 +1227,23 @@ namespace proteus vel[0],vel[1]); center[0] = ball_center[3*i+0]; center[1] = ball_center[3*i+1]; - } + get_acceleration_to_ith_ball(nParticles,ball_center,ball_radius, + ball_velocity,ball_angular_velocity, + ball_center_acceleration,ball_angular_acceleration, + i,x,y,z, + acceleration[0],acceleration[1]); + + } else - { -// phi_s = particle_signed_distances[i * sd_offset]; -// phi_s_normal[0] = particle_signed_distance_normals[i * sd_offset * nSpace + 0]; -// phi_s_normal[1] = particle_signed_distance_normals[i * sd_offset * nSpace + 1]; -// vel[0] = particle_velocities[i * sd_offset * nSpace + 0]; -// vel[1] = particle_velocities[i * sd_offset * nSpace + 1]; -// center[0] = particle_centroids[3*i+0]; -// center[1] = particle_centroids[3*i+1]; - std::cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!YY: use_ball_as_particle should be 1"< 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_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])) + +C_surf*(u-u_s)*rho + ); + 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]) + +C_surf*(v-v_s)*rho + ); 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; + 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])) + +C_surf*(u-u_s)*rho + ); + 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]) + +C_surf*(v-v_s)*rho + ); //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; @@ -1090,7 +1289,7 @@ namespace proteus 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) @@ -1099,65 +1298,57 @@ namespace proteus dmom_u_source[0] += C; dmom_v_source[1] += C; - + if(use_pseudo_penalty==-2)// not pseudo-penalty method == IBM + if (NONCONSERVATIVE_FORM > 0.0) + { + mom_u_source += (1.0-H_s)*(ball_density[i]-rho)*acceleration[0]; + mom_v_source += (1.0-H_s)*(ball_density[i]-rho)*acceleration[1]; + } + else + { + mom_u_source += (1.0-H_s)*(ball_density[i]-rho)*acceleration[0]/rho; + mom_v_source += (1.0-H_s)*(ball_density[i]-rho)*acceleration[1]/rho; + } 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]; + 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]; + + 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; - + 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]; + + 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]; } - //(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 @@ -1496,117 +1687,117 @@ namespace proteus flux_umom = 0.0; flux_vmom = 0.0; if (NONCONSERVATIVE_FORM > 0.0) - { - flowSpeedNormal=n[0]*df_vmom_dv[0]+n[1]*df_umom_du[1];//tricky, works for moving and fixed domains - flowSpeedNormal+=n[0]*dham_grad[0]+n[1]*dham_grad[1]; - } + { + flowSpeedNormal = n[0] * df_vmom_dv[0] + n[1] * df_umom_du[1]; //tricky, works for moving and fixed domains + flowSpeedNormal += n[0] * dham_grad[0] + n[1] * dham_grad[1]; + } else - flowSpeedNormal=n[0]*df_vmom_dv[0]+n[1]*df_umom_du[1];//tricky, works for moving and fixed domains + flowSpeedNormal = n[0] * df_vmom_dv[0] + n[1] * df_umom_du[1]; //tricky, works for moving and fixed domains if (isDOFBoundary_u != 1) + { + flux_mass += n[0] * f_mass[0]; + velocity[0] = f_mass[0]; + if (flowSpeedNormal >= 0.0) { - flux_mass += n[0]*f_mass[0]; - velocity[0] = f_mass[0]; - if (flowSpeedNormal >= 0.0) - { - flux_umom += n[0]*f_umom[0]; - flux_vmom += n[0]*f_vmom[0]; - } - else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - flux_umom+=(0.0-u)*flowSpeedNormal; - } - } + flux_umom += n[0] * f_umom[0]; + flux_vmom += n[0] * f_vmom[0]; + } + else + { + if (NONCONSERVATIVE_FORM > 0.0) + { + flux_umom += (0.0 - u) * flowSpeedNormal; + } } + } else + { + flux_mass += n[0] * f_mass[0]; + velocity[0] = f_mass[0]; + if (flowSpeedNormal >= 0.0) { - flux_mass += n[0]*f_mass[0]; - velocity[0] = f_mass[0]; - if (flowSpeedNormal >= 0.0) - { - flux_umom += n[0]*f_umom[0]; - flux_vmom += n[0]*f_vmom[0]; - } + flux_umom += n[0] * f_umom[0]; + flux_vmom += n[0] * f_vmom[0]; + } + else + { + if (NONCONSERVATIVE_FORM > 0.0) + { + flux_umom += (bc_u - u) * flowSpeedNormal; + } else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - flux_umom+=(bc_u-u)*flowSpeedNormal; - } - else - { - flux_umom+=n[0]*bc_f_umom[0]; - flux_vmom+=n[0]*bc_f_vmom[0]; - } - } + { + flux_umom += n[0] * bc_f_umom[0]; + flux_vmom += n[0] * bc_f_vmom[0]; + } } + } if (isDOFBoundary_v != 1) + { + flux_mass += n[1] * f_mass[1]; + velocity[1] = f_mass[1]; + if (flowSpeedNormal >= 0.0) { - flux_mass+=n[1]*f_mass[1]; - velocity[1] = f_mass[1]; - if (flowSpeedNormal >= 0.0) - { - flux_umom+=n[1]*f_umom[1]; - flux_vmom+=n[1]*f_vmom[1]; - } - else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - flux_vmom+=(0.0-v)*flowSpeedNormal; - } - } + flux_umom += n[1] * f_umom[1]; + flux_vmom += n[1] * f_vmom[1]; + } + else + { + if (NONCONSERVATIVE_FORM > 0.0) + { + flux_vmom += (0.0 - v) * flowSpeedNormal; + } } + } else + { + flux_mass += n[1] * f_mass[1]; + velocity[1] = f_mass[1]; + if (flowSpeedNormal >= 0.0) { - flux_mass+=n[1]*f_mass[1]; - velocity[1] = f_mass[1]; - if (flowSpeedNormal >= 0.0) - { - flux_umom+=n[1]*f_umom[1]; - flux_vmom+=n[1]*f_vmom[1]; - } - else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - flux_vmom+=(bc_v-v)*flowSpeedNormal; - } - else - { - flux_umom+=n[1]*bc_f_umom[1]; - flux_vmom+=n[1]*bc_f_vmom[1]; - } - } + flux_umom += n[1] * f_umom[1]; + flux_vmom += n[1] * f_vmom[1]; } - if (isDOFBoundary_p == 1) + else { if (NONCONSERVATIVE_FORM > 0.0) - { - flux_umom+= n[0]*(bc_p - p); - flux_vmom+= n[1]*(bc_p - p); - } + { + flux_vmom += (bc_v - v) * flowSpeedNormal; + } else - { - flux_umom+= n[0]*(bc_p*bc_oneByRho-p*oneByRho); - flux_vmom+= n[1]*(bc_p*bc_oneByRho-p*oneByRho); - } + { + flux_umom += n[1] * bc_f_umom[1]; + flux_vmom += n[1] * bc_f_vmom[1]; + } } - if (isFluxBoundary_p == 1) + } + if (isDOFBoundary_p == 1) + { + if (NONCONSERVATIVE_FORM > 0.0) { - velocity[0] += (bc_flux_mass - flux_mass)*n[0]; - velocity[1] += (bc_flux_mass - flux_mass)*n[1]; - flux_mass = bc_flux_mass; + flux_umom += n[0] * (bc_p - p); + flux_vmom += n[1] * (bc_p - p); } - if (isFluxBoundary_u == 1) + else { - flux_umom = bc_flux_umom; + flux_umom += n[0] * (bc_p * bc_oneByRho - p * oneByRho); + flux_vmom += n[1] * (bc_p * bc_oneByRho - p * oneByRho); } + } + if (isFluxBoundary_p == 1) + { + velocity[0] += (bc_flux_mass - flux_mass) * n[0]; + velocity[1] += (bc_flux_mass - flux_mass) * n[1]; + flux_mass = bc_flux_mass; + } + if (isFluxBoundary_u == 1) + { + flux_umom = bc_flux_umom; + } if (isFluxBoundary_v == 1) - { - flux_vmom = bc_flux_vmom; - } + { + flux_vmom = bc_flux_vmom; + } } inline @@ -1691,127 +1882,127 @@ namespace proteus flowSpeedNormal=n[0]*df_vmom_dv[0]+n[1]*df_umom_du[1];//tricky, works for moving and fixed domains flowSpeedNormal+=NONCONSERVATIVE_FORM*(n[0]*dham_grad[0]+n[1]*dham_grad[1]); if (isDOFBoundary_u != 1) + { + dflux_mass_du += n[0] * df_mass_du[0]; + if (flowSpeedNormal >= 0.0) { - dflux_mass_du += n[0]*df_mass_du[0]; - if (flowSpeedNormal >= 0.0) - { - dflux_umom_du += n[0]*df_umom_du[0]; - dflux_umom_dv += n[0]*df_umom_dv[0]; + dflux_umom_du += n[0] * df_umom_du[0]; + dflux_umom_dv += n[0] * df_umom_dv[0]; - dflux_vmom_du += n[0]*df_vmom_du[0]; - dflux_vmom_dv += n[0]*df_vmom_dv[0]; - } - else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - dflux_umom_du+=( dmom_u_acc_u*n[0]*(0.0 - u) - flowSpeedNormal ) ; - dflux_umom_dv+= dmom_u_acc_u * (0.0 - u) * n[1]; - } - } + dflux_vmom_du += n[0] * df_vmom_du[0]; + dflux_vmom_dv += n[0] * df_vmom_dv[0]; } + else + { + if (NONCONSERVATIVE_FORM > 0.0) + { + dflux_umom_du += (dmom_u_acc_u * n[0] * (0.0 - u) - flowSpeedNormal); + dflux_umom_dv += dmom_u_acc_u * (0.0 - u) * n[1]; + } + } + } else + { + //cek still upwind the advection for Dirichlet? + dflux_mass_du += n[0] * df_mass_du[0]; + if (flowSpeedNormal >= 0.0) { - //cek still upwind the advection for Dirichlet? - dflux_mass_du += n[0]*df_mass_du[0]; - if (flowSpeedNormal >= 0.0) - { - dflux_umom_du += n[0]*df_umom_du[0]; - dflux_umom_dv += n[0]*df_umom_dv[0]; + dflux_umom_du += n[0] * df_umom_du[0]; + dflux_umom_dv += n[0] * df_umom_dv[0]; - dflux_vmom_du += n[0]*df_vmom_du[0]; - dflux_vmom_dv += n[0]*df_vmom_dv[0]; - } + dflux_vmom_du += n[0] * df_vmom_du[0]; + dflux_vmom_dv += n[0] * df_vmom_dv[0]; + } + else + { + if (NONCONSERVATIVE_FORM > 0.0) + { + dflux_umom_du += (dmom_u_acc_u * n[0] * (bc_u - u) - flowSpeedNormal); + dflux_umom_dv += dmom_u_acc_u * (bc_u - u) * n[1]; + } else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - dflux_umom_du+=( dmom_u_acc_u*n[0]*(bc_u-u) - flowSpeedNormal ) ; - dflux_umom_dv+= dmom_u_acc_u * (bc_u - u) * n[1]; - } - else - { - if (isDOFBoundary_v != 1) - dflux_vmom_dv += n[0]*df_vmom_dv[0]; - } - } + { + if (isDOFBoundary_v != 1) + dflux_vmom_dv += n[0] * df_vmom_dv[0]; + } } + } if (isDOFBoundary_v != 1) + { + dflux_mass_dv += n[1] * df_mass_dv[1]; + if (flowSpeedNormal >= 0.0) { - dflux_mass_dv += n[1]*df_mass_dv[1]; - if (flowSpeedNormal >= 0.0) - { - dflux_umom_du += n[1]*df_umom_du[1]; - dflux_umom_dv += n[1]*df_umom_dv[1]; + dflux_umom_du += n[1] * df_umom_du[1]; + dflux_umom_dv += n[1] * df_umom_dv[1]; - dflux_vmom_du += n[1]*df_vmom_du[1]; - dflux_vmom_dv += n[1]*df_vmom_dv[1]; - } - else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - dflux_vmom_du += dmom_u_acc_u * n[0] * (0.0 - v); - dflux_vmom_dv += ( dmom_u_acc_u * n[1] * (0.0 - v) - flowSpeedNormal) ; - } - } + dflux_vmom_du += n[1] * df_vmom_du[1]; + dflux_vmom_dv += n[1] * df_vmom_dv[1]; + } + else + { + if (NONCONSERVATIVE_FORM > 0.0) + { + dflux_vmom_du += dmom_u_acc_u * n[0] * (0.0 - v); + dflux_vmom_dv += (dmom_u_acc_u * n[1] * (0.0 - v) - flowSpeedNormal); + } } + } else + { + //cek still upwind the advection for Dirichlet? + dflux_mass_dv += n[1] * df_mass_dv[1]; + if (flowSpeedNormal >= 0.0) { - //cek still upwind the advection for Dirichlet? - dflux_mass_dv += n[1]*df_mass_dv[1]; - if (flowSpeedNormal >= 0.0) - { - dflux_umom_du += n[1]*df_umom_du[1]; - dflux_umom_dv += n[1]*df_umom_dv[1]; + dflux_umom_du += n[1] * df_umom_du[1]; + dflux_umom_dv += n[1] * df_umom_dv[1]; - dflux_vmom_du += n[1]*df_vmom_du[1]; - dflux_vmom_dv += n[1]*df_vmom_dv[1]; - } - else - { - if (NONCONSERVATIVE_FORM > 0.0) - { - dflux_vmom_du += dmom_u_acc_u * n[0] * (bc_v - v); - dflux_vmom_dv += ( dmom_u_acc_u * n[1] * (bc_v - v) - flowSpeedNormal) ; - } - else - { - if (isDOFBoundary_u != 1) - dflux_umom_du += n[1]*df_umom_du[1]; - } - } + dflux_vmom_du += n[1] * df_vmom_du[1]; + dflux_vmom_dv += n[1] * df_vmom_dv[1]; } - if (isDOFBoundary_p == 1) + else { if (NONCONSERVATIVE_FORM > 0.0) - { - dflux_umom_dp= -n[0]; - dflux_vmom_dp= -n[1]; - } + { + dflux_vmom_du += dmom_u_acc_u * n[0] * (bc_v - v); + dflux_vmom_dv += (dmom_u_acc_u * n[1] * (bc_v - v) - flowSpeedNormal); + } else - { - dflux_umom_dp= -n[0]*oneByRho; - dflux_vmom_dp= -n[1]*oneByRho; - } + { + if (isDOFBoundary_u != 1) + dflux_umom_du += n[1] * df_umom_du[1]; + } } - if (isFluxBoundary_p == 1) + } + if (isDOFBoundary_p == 1) + { + if (NONCONSERVATIVE_FORM > 0.0) { - dflux_mass_du = 0.0; - dflux_mass_dv = 0.0; + dflux_umom_dp = -n[0]; + dflux_vmom_dp = -n[1]; } - if (isFluxBoundary_u == 1) + else { - dflux_umom_dp = 0.0; - dflux_umom_du = 0.0; - dflux_umom_dv = 0.0; + dflux_umom_dp = -n[0] * oneByRho; + dflux_vmom_dp = -n[1] * oneByRho; } + } + if (isFluxBoundary_p == 1) + { + dflux_mass_du = 0.0; + dflux_mass_dv = 0.0; + } + if (isFluxBoundary_u == 1) + { + dflux_umom_dp = 0.0; + dflux_umom_du = 0.0; + dflux_umom_dv = 0.0; + } if (isFluxBoundary_v == 1) - { - dflux_vmom_dp = 0.0; - dflux_vmom_du = 0.0; - dflux_vmom_dv = 0.0; - } + { + dflux_vmom_dp = 0.0; + dflux_vmom_du = 0.0; + dflux_vmom_dv = 0.0; + } } inline @@ -1959,7 +2150,7 @@ namespace proteus double C_b, //VRANS const double* eps_solid, - const double* phi_solid, + double* phi_solid, const double* q_velocity_solid, const double* q_porosity, const double* q_dragAlpha, @@ -1982,6 +2173,10 @@ namespace proteus double* u_dof, double* v_dof, double* w_dof, + double* p_old_dof, + double* u_old_dof, + double* v_old_dof, + double* w_old_dof, double* g, const double useVF, double* q_rho, @@ -2072,6 +2267,16 @@ namespace proteus double* ball_radius, double* ball_velocity, double* ball_angular_velocity, + double* ball_center_acceleration, + double* ball_angular_acceleration, + double* ball_density, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, + double* ebq_global_phi_s, + double* ebq_global_grad_phi_s, + double* ebq_particle_velocity_s, int nParticles, double *particle_netForces, double *particle_netMoments, @@ -2081,7 +2286,10 @@ namespace proteus double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant) + double particle_penalty_constant, + double* phi_solid_nodes, + double* distance_to_solids, + const int use_pseudo_penalty) { logEvent("Entered mprans calculateResidual",6); @@ -2119,6 +2327,20 @@ namespace proteus elementResidual_v[i]=0.0; velocityErrorElement[i]=0.0; }//i + //Use for plotting result + if(use_ball_as_particle==1) + { + for (int I=0;I 0) + phi_solid[eN_k] = distance_to_solids[eN_k]; // //calculate pde coefficients at quadrature points // @@ -2338,6 +2580,16 @@ namespace proteus kappa_phi[eN_k], //VRANS porosity, + phi_solid[eN_k],//distance to solid + p_old, + u_old, + v_old, + w_old, + grad_p_old, + grad_u_old, + grad_v_old, + grad_w_old, + use_pseudo_penalty, // p, grad_p, @@ -2445,82 +2697,86 @@ namespace proteus 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]); + updateSolidParticleTerms(NONCONSERVATIVE_FORM, + eN < nElements_owned, + particle_nitsche, + dV, + nParticles, + nQuadraturePoints_global, + &particle_signed_distances[eN_k], + &particle_signed_distance_normals[eN_k_3d], + &particle_velocities[eN_k_3d], + particle_centroids, + use_ball_as_particle, + ball_center, + ball_radius, + ball_velocity, + ball_angular_velocity, + ball_center_acceleration, + ball_angular_acceleration, + ball_density, + 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], + use_pseudo_penalty); //Turbulence closure model if (turbulenceClosureModel >= 3) { @@ -2613,6 +2869,20 @@ namespace proteus dmom_v_acc_v, mom_v_acc_t, dmom_v_acc_v_t); + if(use_pseudo_penalty > 0 && phi_solid[eN_k]<0.0)//Do not have to change Jacobian + { + double distance,vx,vy; + int index_ball = get_distance_to_ball(nParticles, ball_center, ball_radius,x,y,z,distance); + get_velocity_to_ith_ball(nParticles,ball_center,ball_radius,ball_velocity,ball_angular_velocity,index_ball,x,y,z,vx,vy); + mom_u_acc_t = alphaBDF*(mom_u_acc - vx); + mom_v_acc_t = alphaBDF*(mom_v_acc - vy); + }else if(use_pseudo_penalty == -1 && phi_solid[eN_k]<0.0)//no derivative term inside the solid; Has to change Jacobian + { + mom_u_acc_t = 0.0; + mom_v_acc_t = 0.0; + dmom_u_acc_u= 0.0; + dmom_v_acc_v= 0.0; + } if (NONCONSERVATIVE_FORM > 0.0) { mom_u_acc_t *= dmom_u_acc_u; @@ -2861,6 +3131,8 @@ namespace proteus grad_u_ext[nSpace], grad_v_ext[nSpace], grad_w_ext[nSpace], + p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0, + grad_p_old[nSpace],grad_u_old[nSpace],grad_v_old[nSpace],grad_w_old[nSpace], mom_u_acc_ext=0.0, dmom_u_acc_u_ext=0.0, mom_v_acc_ext=0.0, @@ -3056,9 +3328,15 @@ namespace proteus ck.valFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_trace_ref[ebN_local_kb*nDOF_test_element],p_ext); ck_v.valFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_ext); ck_v.valFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_ext); + ck.valFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_trace_ref[ebN_local_kb*nDOF_test_element],p_old); + ck_v.valFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_old); + ck_v.valFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_old); ck.gradFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext); ck_v.gradFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_ext); ck_v.gradFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_ext); + ck.gradFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_old); + ck_v.gradFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_old); + ck_v.gradFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_old); //precalculate test function products with integration weights for (int j=0;j particle_surfaceArea(nParticles), particle_netForces(nParticles*3*3), particle_netMoments(nParticles*3); @@ -3840,12 +4160,15 @@ namespace proteus { int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point eN_k_nSpace = eN_k*nSpace, + eN_k_3d = eN_k*3, eN_nDOF_trial_element = eN*nDOF_trial_element, //index to a vector at a quadrature point eN_nDOF_v_trial_element = eN*nDOF_v_trial_element; //index to a vector at a quadrature point //declare local storage register double p=0.0,u=0.0,v=0.0,w=0.0, grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace], + p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0, + grad_p_old[nSpace],grad_u_old[nSpace],grad_v_old[nSpace],grad_w_old[nSpace], mom_u_acc=0.0, dmom_u_acc_u=0.0, mom_v_acc=0.0, @@ -3997,10 +4320,16 @@ namespace proteus ck.valFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_ref[k*nDOF_trial_element],p); ck_v.valFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],u); ck_v.valFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],v); + ck.valFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_ref[k*nDOF_trial_element],p_old); + ck_v.valFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],u_old); + ck_v.valFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_ref[k*nDOF_v_trial_element],v_old); //get the solution gradients ck.gradFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial,grad_p); ck_v.gradFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_u); ck_v.gradFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_v); + ck.gradFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial,grad_p_old); + ck_v.gradFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_u_old); + ck_v.gradFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial,grad_v_old); //precalculate test function products with integration weights for (int j=0;j 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]); + updateSolidParticleTerms(NONCONSERVATIVE_FORM, + eN < nElements_owned, + particle_nitsche, + dV, + nParticles, + nQuadraturePoints_global, + &particle_signed_distances[eN_k], + &particle_signed_distance_normals[eN_k_3d], + &particle_velocities[eN_k_3d], + particle_centroids, + use_ball_as_particle, + ball_center, + ball_radius, + ball_velocity, + ball_angular_velocity, + ball_center_acceleration, + ball_angular_acceleration, + ball_density, + 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], + use_pseudo_penalty); //Turbulence closure model if (turbulenceClosureModel >= 3) { @@ -4320,6 +4663,13 @@ namespace proteus dmom_v_acc_v, mom_v_acc_t, dmom_v_acc_v_t); + if(use_pseudo_penalty == -1 && phi_solid[eN_k]<0.0)//no derivative term inside the solid; Has to change Jacobian + { + mom_u_acc_t = 0.0; + mom_v_acc_t = 0.0; + dmom_u_acc_u = 0.0; + dmom_v_acc_v = 0.0; + } if (NONCONSERVATIVE_FORM > 0.0) { mom_u_acc_t *= dmom_u_acc_u; @@ -4645,6 +4995,8 @@ namespace proteus grad_u_ext[nSpace], grad_v_ext[nSpace], grad_w_ext[nSpace], + p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0, + grad_p_old[nSpace],grad_u_old[nSpace],grad_v_old[nSpace],grad_w_old[nSpace], mom_u_acc_ext=0.0, dmom_u_acc_u_ext=0.0, mom_v_acc_ext=0.0, @@ -4851,10 +5203,16 @@ namespace proteus ck.valFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_trace_ref[ebN_local_kb*nDOF_test_element],p_ext); ck_v.valFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_ext); ck_v.valFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_ext); + ck.valFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],&p_trial_trace_ref[ebN_local_kb*nDOF_test_element],p_old); + ck_v.valFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],u_old); + ck_v.valFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_v_test_element],v_old); ck.gradFromDOF(p_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext); ck_v.gradFromDOF(u_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_ext); ck_v.gradFromDOF(v_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_ext); - //precalculate test function products with integration weights + ck.gradFromDOF(p_old_dof,&p_l2g[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_old); + ck_v.gradFromDOF(u_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_old); + ck_v.gradFromDOF(v_old_dof,&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_old); + //precalculate test function products with integration weights for (int j=0;j surrogate_boundaries, surrogate_boundary_elements, surrogate_boundary_particle; + double C_sbm, beta_sbm; cppHsuSedStress<3> closure; const int nDOF_test_X_trial_element, nSpace2=9; @@ -610,7 +612,9 @@ namespace proteus 0.05, 1.00), nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element), - ck() + ck(), + C_sbm(2.0), + beta_sbm(0.0) {/* std::cout<<"Constructing cppRANS3PF surrogate_boundaries, surrogate_boundary_elements, surrogate_boundary_particle; for(int eN=0;eN particle_surfaceArea(nParticles), particle_netForces(nParticles * 3), particle_netMoments(nParticles * 3); const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element); - std::vector surrogate_boundaries, surrogate_boundary_elements,surrogate_boundary_particle; for(int eN=0;eN=0); assert(opp_node 0.0); if (h_penalty < std::abs(dist)) h_penalty = std::abs(dist); + //hack: this won't work for two-phase flow, need mixture viscosity double visco = nu_0*rho_0; - double Csb=10; - double C_adim = Csb*visco/h_penalty; - double beta = 0.0; - double beta_adim = beta*h_penalty*visco; + double C_adim = C_sbm*visco/h_penalty; + double beta_adim = beta_sbm*h_penalty*visco; double res[3]; diff --git a/proteus/mprans/RANS3PF.py b/proteus/mprans/RANS3PF.py index 539ed00f96..4ab599d0f6 100644 --- a/proteus/mprans/RANS3PF.py +++ b/proteus/mprans/RANS3PF.py @@ -472,13 +472,13 @@ def attachModels(self, modelList): self.particle_surfaceArea = np.zeros((self.nParticles,), 'd') self.particle_centroids = np.zeros((self.nParticles, 3), 'd') self.particle_signed_distances = np.zeros((self.nParticles,) + self.model.q[('u', 0)].shape, 'd') - self.particle_signed_distance_normals = np.zeros((self.nParticles,) + self.model.q[('velocity', 0)].shape, 'd') - self.particle_velocities = np.zeros((self.nParticles,) + self.model.q[('velocity', 0)].shape, 'd') + self.particle_signed_distance_normals = np.zeros((self.nParticles,) + self.model.q[('velocity', 0)].shape[:-1]+(3,), 'd') + self.particle_velocities = np.zeros((self.nParticles,) + self.model.q[('velocity', 0)].shape[:-1]+(3,), 'd') self.phisField = np.ones(self.model.q[('u', 0)].shape, 'd') * 1e10 self.ebq_global_phi_s = numpy.ones_like(self.model.ebq_global[('totalFlux',0)]) * 1.e10 - self.ebq_global_grad_phi_s = numpy.ones_like(self.model.ebq_global[('velocityAverage',0)]) * 1.e10 - self.ebq_particle_velocity_s = numpy.ones_like(self.model.ebq_global[('velocityAverage',0)]) * 1.e10 + self.ebq_global_grad_phi_s = numpy.ones(self.model.ebq_global[('velocityAverage',0)].shape[:-1]+(3,),'d') * 1.e10 + self.ebq_particle_velocity_s = numpy.ones(self.model.ebq_global[('velocityAverage',0)].shape[:-1]+(3,),'d') * 1.e10 # This is making a special case for granular material simulations # if the user inputs a list of position/velocities then the sdf are calculated based on the "spherical" particles @@ -742,12 +742,16 @@ def initializeMesh(self, mesh): self.barycenters = numpy.zeros((nBoundariesMax, 3), 'd') comm = Comm.get() if comm.isMaster(): + self.timeHistory = open(os.path.join(proteus.Profiling.logDir, "timeHistory.txt"), "w") self.wettedAreaHistory = open(os.path.join(proteus.Profiling.logDir, "wettedAreaHistory.txt"), "w") 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_surfaceAreaHistory = open(os.path.join(proteus.Profiling.logDir, "particle_surfaceAreaHistory.txt"), "w") self.particle_forceHistory = open(os.path.join(proteus.Profiling.logDir, "particle_forceHistory.txt"), "w") + self.particle_pforceHistory = open(os.path.join(proteus.Profiling.logDir, "particle_pforceHistory.txt"), "w") + self.particle_vforceHistory = open(os.path.join(proteus.Profiling.logDir, "particle_vforceHistory.txt"), "w") self.particle_momentHistory = open(os.path.join(proteus.Profiling.logDir, "particle_momentHistory.txt"), "w") # initialize so it can run as single phase @@ -1072,7 +1076,10 @@ def postStep(self, t, firstStep=False): if self.model.isActiveDOF[dof] < 0.5: self.model.u[ci].dof[ci_g_dof] = vel_at_xyz[ci] if self.model.comm.isMaster(): + self.timeHistory.write("%21.16e\n" % (t,)) + self.timeHistory.flush() self.wettedAreaHistory.write("%21.16e\n" % (self.wettedAreas[-1],)) + self.wettedAreaHistory.flush() self.forceHistory_p.write("%21.16e %21.16e %21.16e\n" % tuple(self.netForces_p[-1, :])) self.forceHistory_p.flush() self.forceHistory_v.write("%21.16e %21.16e %21.16e\n" % tuple(self.netForces_v[-1, :])) @@ -1080,8 +1087,14 @@ def postStep(self, t, firstStep=False): self.momentHistory.write("%21.15e %21.16e %21.16e\n" % tuple(self.netMoments[-1, :])) self.momentHistory.flush() if self.nParticles: + self.particle_surfaceAreaHistory.write("%21.16e\n" % tuple(self.particle_surfaceArea[:])) + self.particle_surfaceAreaHistory.flush() self.particle_forceHistory.write("%21.16e %21.16e %21.16e\n" % tuple(self.particle_netForces[0, :])) self.particle_forceHistory.flush() + self.particle_pforceHistory.write("%21.16e %21.16e %21.16e\n" % tuple(self.particle_netForces[0+self.nParticles, :])) + self.particle_pforceHistory.flush() + self.particle_vforceHistory.write("%21.16e %21.16e %21.16e\n" % tuple(self.particle_netForces[0+2*self.nParticles, :])) + self.particle_vforceHistory.flush() self.particle_momentHistory.write("%21.15e %21.16e %21.16e\n" % tuple(self.particle_netMoments[0, :])) self.particle_momentHistory.flush() @@ -2143,8 +2156,8 @@ def getSparsityPatternForComponents(self): # fill vector rowptr_scalar for i in range(1, self.rowptr_1D.size): self.rowptr_1D[i]=(self.rowptr_1D[i - 1] + - old_div((self.rowptr[nSpace * (i - 1) + 1] - - self.rowptr[nSpace * (i - 1)]), nSpace)) + (self.rowptr[nSpace * (i - 1) + 1] - + self.rowptr[nSpace * (i - 1)])//nSpace) # fill vector colind_cMatrix ith_row = 0 for i in range(len(self.rowptr)-1): # 0 to total num of DOFs (i.e. num of rows of jacobian) @@ -2153,7 +2166,7 @@ def getSparsityPatternForComponents(self): offset_1D = list(range(self.rowptr_1D[ith_row], self.rowptr_1D[ith_row + 1])) if (j % nSpace == 0): - self.colind_1D[offset_1D[old_div(j, nSpace)]] = old_div(self.colind[offset], nSpace) + self.colind_1D[offset_1D[j//nSpace]] = self.colind[offset]//nSpace ith_row += 1 def getResidual(self, u, r): @@ -2355,7 +2368,7 @@ def getResidual(self, u, r): self.q['eddy_viscosity'], self.pressureModel.u[0].femSpace.dofMap.l2g, self.u[0].femSpace.dofMap.l2g, - self.pressureModel.u[0].dof, + self.pressureModel.p_sharp_dof, self.u[0].dof, self.u[1].dof, self.u[2].dof, @@ -2526,7 +2539,7 @@ def getResidual(self, u, r): self.numDOFs_1D, self.nnz_1D, self.csrRowIndeces[(0, 0)] // self.nSpace_global // self.nSpace_global, - old_div(self.csrColumnOffsets[(0, 0)], self.nSpace_global), + self.csrColumnOffsets[(0, 0)]//self.nSpace_global, self.rowptr_1D, self.colind_1D, self.isBoundary_1D, diff --git a/proteus/mprans/RANS3PF2D.h b/proteus/mprans/RANS3PF2D.h index 2e47e9b79f..938d84a372 100644 --- a/proteus/mprans/RANS3PF2D.h +++ b/proteus/mprans/RANS3PF2D.h @@ -11,7 +11,10 @@ #include "SedClosure.h" #define DRAG_FAC 1.0 #define TURB_FORCE_FAC 0.0 -#define CUT_CELL_INTEGRATION 0.0 +#define CUT_CELL_INTEGRATION 0 +double sgn(double val) { + return double((0.0 < val) - (val < 0.0)); +} ////////////////////// // ***** TODO ***** // ////////////////////// @@ -591,10 +594,9 @@ namespace proteus class cppRANS3PF2D : public cppRANS3PF2D_base { public: -// std::vector surrogate_boundaries, surrogate_boundary_elements, surrogate_boundary_particle; - const double C_sbm;//penalty constant for sbm - const double beta_sbm;//tangent penalty constant for sbm - + std::vector surrogate_boundaries, surrogate_boundary_elements, surrogate_boundary_particle; + std::valarray TransportMatrix, TransposeTransportMatrix, psi; + double C_sbm, beta_sbm; cppHsuSedStress<2> closure; const int nDOF_test_X_trial_element, nSpace2; @@ -617,9 +619,9 @@ namespace proteus 5.0, M_PI/6., 0.05, 1.00), nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element), - C_sbm(1000), - beta_sbm(0.0), - ck() + ck(), + C_sbm(10.0), + beta_sbm(0.0) {/* std::cout<<"Constructing cppRANS3PF2D surrogate_boundaries, surrogate_boundary_elements, surrogate_boundary_particle; - surrogate_boundaries.clear(); - surrogate_boundary_elements.clear(); - surrogate_boundary_particle.clear(); //std::set active_velocity_dof; for(int eN=0;eN ls_nodes; for (int I=0;I eps)//level set does NOT lie on edge (intersects line SOMEWHERE) + if (fabs(delta_phi) > eps)//level sets are not parallel to edge //need tolerance selection guidance { theta = -_distance[I]/delta_phi;//zero level set is at theta*xIp1+(1-theta)*xI - if (theta > 1.0-eps || theta < eps)//zero level does NOT intersect between nodes + if (theta > 1.0-eps || theta < eps)//zero level does NOT intersect between nodes; it may got through a node { - theta = 0.5;//just put the subelement node at midpoint + if (theta > 1.0-eps && theta <= 1.0)// + { + ls_nodes.push_back((I+1)%3); + //todo, fix connectivity for this case--can't use 4T + assert(false); + } + else if (theta > 0.0 && theta < eps)// + { + ls_nodes.push_back(I); + assert(false); + } + else + theta = 0.5;//just put the subelement node at midpoint } else - boundaryNodes[3+I]=1; + { + boundaryNodes[3+I]=1; + ls_nodes.push_back(3+I); + } } else //level set lies on edge { @@ -2516,6 +2539,8 @@ namespace proteus boundaryNodes[I]=1; boundaryNodes[3+I]=1; boundaryNodes[(I+1)%3]=1; + ls_nodes.push_back(I); + ls_nodes.push_back((I+1)%3); } } assert(theta <= 1.0); @@ -2523,6 +2548,13 @@ namespace proteus GI[3*3 + I*3 + (I+1)%3] = theta; GI[3*3 + I*3 + (I+2)%3] = 0.0; } + if (ls_nodes.size() != 2) + { + std::cout<<"level set nodes not 2 "< 1.0-1.0e-4 && dot_test < 1.0 + 1.0e-4); + cut_cell_boundary_length += DS; + p_force_x += sub_p_dof[L]*nx*0.5*DS + sub_p_dof[R]*nx*0.5*DS; + p_force_y += sub_p_dof[L]*ny*0.5*DS + sub_p_dof[R]*ny*0.5*DS; //TODO for P2 //1. Now define the Lagrange nodes for P2 on the submesh X //2. Define and evaluate the P2 trial functions for the parent element at the new submesh P2 nodes. X @@ -2598,123 +2656,123 @@ namespace proteus } for (int I=0; I<6; I++) { - std::cout<0) { - // - //detect cut cells - // - double _distance[nDOF_mesh_trial_element]={0.0}; - int pos_counter=0; - for (int I=0;I= 0) - pos_counter++; + if ( _distance[I] >= 0) + pos_counter++; } - if (pos_counter == 2) + if (pos_counter == 2) { - element_active=0.0; - int opp_node=-1; - for (int I=0;I=0); - assert(opp_node =0); + assert(opp_node 0) - std::cout< 0) + std::cout< globalResidual[GlobPos_u] += beta_adim*grad_u_t[0]*grad_phi_i_dot_t; globalResidual[GlobPos_v] += beta_adim*grad_u_t[1]*grad_phi_i_dot_t; + Fx += beta_adim*grad_u_t[0]*grad_phi_i_dot_t; + Fy += beta_adim*grad_u_t[1]*grad_phi_i_dot_t; }//i @@ -4097,21 +4173,11 @@ namespace proteus } double nx = P_normal[0]; //YY: normal direction outward of the solid. double ny = P_normal[1]; - - double S_xx = 2*visco*grad_u_ext[0]; - double S_xy = visco*(grad_u_ext[1] + grad_v_ext[0]); // sym tensor -> S_yx = S_xy - double S_yy = 2*visco*grad_v_ext[1]; - Fx -= p_ext*nx*dS; - Fx += (S_xx*nx + S_xy*ny)*dS; - // Fx += dS*(C_adim*(u_ext - bc_u_ext) - // - visco * (normal[0]*2*grad_u_ext[0] + normal[1]*(grad_u_ext[1]+grad_v_ext[0])) - // + C_adim*dd1); Fy -= p_ext*ny*dS; - Fy += (S_xy*nx + S_yy*ny)*dS; - // Fy += dS*(C_adim*(v_ext - bc_v_ext) - // - visco * (normal[0]*(grad_u_ext[1]+grad_v_ext[0]) + normal[1]*2*grad_v_ext[1]) - // + C_adim*dd2); + Fxp -= p_ext*nx*dS; + Fyp -= p_ext*ny*dS; + surfaceArea += dS; if(use_ball_as_particle==1) { r_x = x_ext - ball_center[surrogate_boundary_particle[ebN_s] * 3 + 0]; @@ -4125,13 +4191,17 @@ namespace proteus Mz += r_x*Fy-r_y*Fx; }//kb if(USE_SBM==1 - && ebN < nElementBoundaries_owned)//avoid double counting - { + && ebN < nElementBoundaries_owned)//avoid double counting + { + particle_surfaceArea[surrogate_boundary_particle[ebN_s]] += surfaceArea; particle_netForces[3*surrogate_boundary_particle[ebN_s]+0] += Fx; particle_netForces[3*surrogate_boundary_particle[ebN_s]+1] += Fy; + particle_netForces[3*( nParticles+surrogate_boundary_particle[ebN_s])+0] += Fxp; + particle_netForces[3*(2*nParticles+surrogate_boundary_particle[ebN_s])+0] += (Fx-Fxp); + particle_netForces[3*( nParticles+surrogate_boundary_particle[ebN_s])+1] += Fyp; + particle_netForces[3*(2*nParticles+surrogate_boundary_particle[ebN_s])+1] += (Fy-Fyp); particle_netMoments[3*surrogate_boundary_particle[ebN_s]+2]+= Mz; - } - + } }//ebN_s //std::cout<<" sbm force over surrogate boundary is: "< particle_surfaceArea(nParticles), particle_netForces(nParticles*3*3), particle_netMoments(nParticles*3); const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element); - std::vector surrogate_boundaries,surrogate_boundary_elements,surrogate_boundary_particle; //std::set active_velocity_dof; for(int eN=0;eN= 0) - pos_counter++; - } + pos_counter++; + } if (pos_counter == 2) { element_active=0.0; @@ -5367,65 +5434,6 @@ namespace proteus } assert(opp_node >=0); assert(opp_node 0.0); if (h_penalty < std::abs(dist)) h_penalty = std::abs(dist); + //hack: this won't work for two-phase flow, need mixture viscosity double visco = nu_0*rho_0; double C_adim = C_sbm*visco/h_penalty; double beta_adim = beta_sbm*visco/h_penalty; diff --git a/proteus/mprans/cRANS2P.pyx b/proteus/mprans/cRANS2P.pyx index 7df58a5207..83e9555444 100644 --- a/proteus/mprans/cRANS2P.pyx +++ b/proteus/mprans/cRANS2P.pyx @@ -88,6 +88,10 @@ cdef extern from "RANS2P.h" namespace "proteus": double * u_dof, double * v_dof, double * w_dof, + double * p_old_dof, + double * u_old_dof, + double * v_old_dof, + double * w_old_dof, double * g, double useVF, double * rho, @@ -178,6 +182,16 @@ cdef extern from "RANS2P.h" namespace "proteus": double* ball_radius, double* ball_velocity, double* ball_angular_velocity, + double* ball_center_acceleration, + double* ball_angular_acceleration, + double* ball_density, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, + double* ebq_global_phi_s, + double* ebq_global_grad_phi_s, + double* ebq_particle_velocity_s, int nParticles, double *particle_netForces, double *particle_netMoments, @@ -187,7 +201,10 @@ cdef extern from "RANS2P.h" namespace "proteus": double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant) + double particle_penalty_constant, + double *phi_solid_nodes, + double* distance_to_solids, + int use_pseudo_penalty) void calculateJacobian(double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, @@ -260,6 +277,10 @@ cdef extern from "RANS2P.h" namespace "proteus": int * p_l2g, int * vel_l2g, double * p_dof, double * u_dof, double * v_dof, double * w_dof, + double * p_old_dof, + double * u_old_dof, + double * v_old_dof, + double * w_old_dof, double * g, double useVF, double * vf, @@ -359,13 +380,26 @@ cdef extern from "RANS2P.h" namespace "proteus": double* ball_radius, double* ball_velocity, double* ball_angular_velocity, + double* ball_center_acceleration, + double* ball_angular_acceleration, + double* ball_density, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, + double* ebq_global_phi_s, + double* ebq_global_grad_phi_s, + double* ebq_particle_velocity_s, + double* phi_solid_nodes, + double* distance_to_solids, int nParticles, int nElements_owned, double particle_nitsche, double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant) + double particle_penalty_constant, + int use_pseudo_penalty) void calculateVelocityAverage(int nExteriorElementBoundaries_global, int * exteriorElementBoundariesArray, int nInteriorElementBoundaries_global, @@ -604,6 +638,10 @@ cdef class cRANS2P_base: numpy.ndarray u_dof, numpy.ndarray v_dof, numpy.ndarray w_dof, + numpy.ndarray p_old_dof, + numpy.ndarray u_old_dof, + numpy.ndarray v_old_dof, + numpy.ndarray w_old_dof, numpy.ndarray g, double useVF, numpy.ndarray rho, @@ -694,6 +732,16 @@ cdef class cRANS2P_base: numpy.ndarray ball_radius, numpy.ndarray ball_velocity, numpy.ndarray ball_angular_velocity, + numpy.ndarray ball_center_acceleration, + numpy.ndarray ball_angular_acceleration, + numpy.ndarray ball_density, + numpy.ndarray particle_signed_distances, + numpy.ndarray particle_signed_distance_normals, + numpy.ndarray particle_velocities, + numpy.ndarray particle_centroids, + numpy.ndarray ebq_global_phi_s, + numpy.ndarray ebq_global_grad_phi_s, + numpy.ndarray ebq_particle_velocity_s, int nParticles, numpy.ndarray particle_netForces, numpy.ndarray particle_netMoments, @@ -703,7 +751,10 @@ cdef class cRANS2P_base: double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant): + double particle_penalty_constant, + numpy.ndarray phi_solid_nodes, + numpy.ndarray distance_to_solids, + int use_pseudo_penalty): self.thisptr.calculateResidual(NONCONSERVATIVE_FORM, MOMENTUM_SGE, PRESSURE_SGE, @@ -785,6 +836,10 @@ cdef class cRANS2P_base: < double * > u_dof.data, < double * > v_dof.data, < double * > w_dof.data, + < double * > p_old_dof.data, + < double * > u_old_dof.data, + < double * > v_old_dof.data, + < double * > w_old_dof.data, < double * > g.data, useVF, < double * > rho.data, @@ -875,6 +930,16 @@ cdef class cRANS2P_base: < double * > ball_radius.data, < double * > ball_velocity.data, < double * > ball_angular_velocity.data, + < double * > ball_center_acceleration.data, + < double * > ball_angular_acceleration.data, + < double * > ball_density.data, + < double * > particle_signed_distances.data, + < double * > particle_signed_distance_normals.data, + < double * > particle_velocities.data, + < double * > particle_centroids.data, + < double * > ebq_global_phi_s.data, + < double * > ebq_global_grad_phi_s.data, + < double * > ebq_particle_velocity_s.data, nParticles, < double * > particle_netForces.data, < double * > particle_netMoments.data, @@ -884,7 +949,10 @@ cdef class cRANS2P_base: particle_epsFact, particle_alpha, particle_beta, - particle_penalty_constant) + particle_penalty_constant, + < double * > phi_solid_nodes.data, + < double * > distance_to_solids.data, + use_pseudo_penalty) def calculateJacobian(self, double NONCONSERVATIVE_FORM, @@ -959,6 +1027,10 @@ cdef class cRANS2P_base: numpy.ndarray p_l2g, numpy.ndarray vel_l2g, numpy.ndarray p_dof, numpy.ndarray u_dof, numpy.ndarray v_dof, numpy.ndarray w_dof, + numpy.ndarray p_old_dof, + numpy.ndarray u_old_dof, + numpy.ndarray v_old_dof, + numpy.ndarray w_old_dof, numpy.ndarray g, double useVF, numpy.ndarray vf, @@ -1058,13 +1130,26 @@ cdef class cRANS2P_base: numpy.ndarray ball_radius, numpy.ndarray ball_velocity, numpy.ndarray ball_angular_velocity, + numpy.ndarray ball_center_acceleration, + numpy.ndarray ball_angular_acceleration, + numpy.ndarray ball_density, + numpy.ndarray particle_signed_distances, + numpy.ndarray particle_signed_distance_normals, + numpy.ndarray particle_velocities, + numpy.ndarray particle_centroids, + numpy.ndarray ebq_global_phi_s, + numpy.ndarray ebq_global_grad_phi_s, + numpy.ndarray ebq_particle_velocity_s, + numpy.ndarray phi_solid_nodes, + numpy.ndarray distance_to_solids, int nParticles, int nElements_owned, double particle_nitsche, double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant): + double particle_penalty_constant, + int use_pseudo_penalty): cdef numpy.ndarray rowptr, colind, globalJacobian_a (rowptr, colind, globalJacobian_a) = globalJacobian.getCSRrepresentation() self.thisptr.calculateJacobian(NONCONSERVATIVE_FORM, @@ -1139,6 +1224,10 @@ cdef class cRANS2P_base: < int * > p_l2g.data, < int * > vel_l2g.data, < double * > p_dof.data, < double * > u_dof.data, < double * > v_dof.data, < double * > w_dof.data, + < double * > p_old_dof.data, + < double * > u_old_dof.data, + < double * > v_old_dof.data, + < double * > w_old_dof.data, < double * > g.data, useVF, < double * > vf.data, @@ -1238,13 +1327,26 @@ cdef class cRANS2P_base: < double * > ball_radius.data, < double * > ball_velocity.data, < double * > ball_angular_velocity.data, + < double * > ball_center_acceleration.data, + < double * > ball_angular_acceleration.data, + < double * > ball_density.data, + < double * > particle_signed_distances.data, + < double * > particle_signed_distance_normals.data, + < double * > particle_velocities.data, + < double * > particle_centroids.data, + < double * > ebq_global_phi_s.data, + < double * > ebq_global_grad_phi_s.data, + < double * > ebq_particle_velocity_s.data, + < double * > phi_solid_nodes.data, + < double * > distance_to_solids.data, nParticles, nElements_owned, particle_nitsche, particle_epsFact, particle_alpha, particle_beta, - particle_penalty_constant) + particle_penalty_constant, + use_pseudo_penalty) def calculateVelocityAverage(self, int nExteriorElementBoundaries_global, diff --git a/proteus/mprans/cRANS2P2D.pyx b/proteus/mprans/cRANS2P2D.pyx index 2374caf8d1..a691e2feaa 100644 --- a/proteus/mprans/cRANS2P2D.pyx +++ b/proteus/mprans/cRANS2P2D.pyx @@ -88,6 +88,10 @@ cdef extern from "RANS2P2D.h" namespace "proteus": double * u_dof, double * v_dof, double * w_dof, + double * p_old_dof, + double * u_old_dof, + double * v_old_dof, + double * w_old_dof, double * g, double useVF, double * rho, @@ -178,6 +182,16 @@ cdef extern from "RANS2P2D.h" namespace "proteus": double* ball_radius, double* ball_velocity, double* ball_angular_velocity, + double* ball_center_acceleration, + double* ball_angular_acceleration, + double* ball_density, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, + double* ebq_global_phi_s, + double* ebq_global_grad_phi_s, + double* ebq_particle_velocity_s, int nParticles, double *particle_netForces, double *particle_netMoments, @@ -187,7 +201,10 @@ cdef extern from "RANS2P2D.h" namespace "proteus": double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant) + double particle_penalty_constant, + double* phi_solid_nodes, + double* distance_to_solids, + int use_pseudo_penalty) void calculateJacobian(double NONCONSERVATIVE_FORM, double MOMENTUM_SGE, double PRESSURE_SGE, @@ -260,6 +277,10 @@ cdef extern from "RANS2P2D.h" namespace "proteus": int * p_l2g, int * vel_l2g, double * p_dof, double * u_dof, double * v_dof, double * w_dof, + double * p_old_dof, + double * u_old_dof, + double * v_old_dof, + double * w_old_dof, double * g, double useVF, double * vf, @@ -359,13 +380,26 @@ cdef extern from "RANS2P2D.h" namespace "proteus": double* ball_radius, double* ball_velocity, double* ball_angular_velocity, + double* ball_center_acceleration, + double* ball_angular_acceleration, + double* ball_density, + double* particle_signed_distances, + double* particle_signed_distance_normals, + double* particle_velocities, + double* particle_centroids, + double* ebq_global_phi_s, + double* ebq_global_grad_phi_s, + double* ebq_particle_velocity_s, + double* phi_solid_nodes, + double* distance_to_solids, int nParticles, int nElements_owned, double particle_nitsche, double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant) + double particle_penalty_constant, + int use_pseudo_penalty) void calculateVelocityAverage(int nExteriorElementBoundaries_global, int * exteriorElementBoundariesArray, int nInteriorElementBoundaries_global, @@ -604,6 +638,10 @@ cdef class cRANS2P2D_base: numpy.ndarray u_dof, numpy.ndarray v_dof, numpy.ndarray w_dof, + numpy.ndarray p_old_dof, + numpy.ndarray u_old_dof, + numpy.ndarray v_old_dof, + numpy.ndarray w_old_dof, numpy.ndarray g, double useVF, numpy.ndarray rho, @@ -694,6 +732,16 @@ cdef class cRANS2P2D_base: numpy.ndarray ball_radius, numpy.ndarray ball_velocity, numpy.ndarray ball_angular_velocity, + numpy.ndarray ball_center_acceleration, + numpy.ndarray ball_angular_acceleration, + numpy.ndarray ball_density, + numpy.ndarray particle_signed_distances, + numpy.ndarray particle_signed_distance_normals, + numpy.ndarray particle_velocities, + numpy.ndarray particle_centroids, + numpy.ndarray ebq_global_phi_s, + numpy.ndarray ebq_global_grad_phi_s, + numpy.ndarray ebq_particle_velocity_s, int nParticles, numpy.ndarray particle_netForces, numpy.ndarray particle_netMoments, @@ -703,7 +751,10 @@ cdef class cRANS2P2D_base: double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant): + double particle_penalty_constant, + numpy.ndarray phi_solid_nodes, + numpy.ndarray distance_to_solids, + int use_pseudo_penalty): self.thisptr.calculateResidual(NONCONSERVATIVE_FORM, MOMENTUM_SGE, PRESSURE_SGE, @@ -785,6 +836,10 @@ cdef class cRANS2P2D_base: < double * > u_dof.data, < double * > v_dof.data, < double * > w_dof.data, + < double * > p_old_dof.data, + < double * > u_old_dof.data, + < double * > v_old_dof.data, + < double * > w_old_dof.data, < double * > g.data, useVF, < double * > rho.data, @@ -875,6 +930,16 @@ cdef class cRANS2P2D_base: < double * > ball_radius.data, < double * > ball_velocity.data, < double * > ball_angular_velocity.data, + < double * > ball_center_acceleration.data, + < double * > ball_angular_acceleration.data, + < double * > ball_density.data, + < double * > particle_signed_distances.data, + < double * > particle_signed_distance_normals.data, + < double * > particle_velocities.data, + < double * > particle_centroids.data, + < double * > ebq_global_phi_s.data, + < double * > ebq_global_grad_phi_s.data, + < double * > ebq_particle_velocity_s.data, nParticles, < double * > particle_netForces.data, < double * > particle_netMoments.data, @@ -884,7 +949,10 @@ cdef class cRANS2P2D_base: particle_epsFact, particle_alpha, particle_beta, - particle_penalty_constant) + particle_penalty_constant, + < double * > phi_solid_nodes.data, + < double * > distance_to_solids.data, + use_pseudo_penalty) def calculateJacobian(self, double NONCONSERVATIVE_FORM, @@ -959,6 +1027,10 @@ cdef class cRANS2P2D_base: numpy.ndarray p_l2g, numpy.ndarray vel_l2g, numpy.ndarray p_dof, numpy.ndarray u_dof, numpy.ndarray v_dof, numpy.ndarray w_dof, + numpy.ndarray p_old_dof, + numpy.ndarray u_old_dof, + numpy.ndarray v_old_dof, + numpy.ndarray w_old_dof, numpy.ndarray g, double useVF, numpy.ndarray vf, @@ -1058,13 +1130,26 @@ cdef class cRANS2P2D_base: numpy.ndarray ball_radius, numpy.ndarray ball_velocity, numpy.ndarray ball_angular_velocity, + numpy.ndarray ball_center_acceleration, + numpy.ndarray ball_angular_acceleration, + numpy.ndarray ball_density, + numpy.ndarray particle_signed_distances, + numpy.ndarray particle_signed_distance_normals, + numpy.ndarray particle_velocities, + numpy.ndarray particle_centroids, + numpy.ndarray ebq_global_phi_s, + numpy.ndarray ebq_global_grad_phi_s, + numpy.ndarray ebq_particle_velocity_s, + numpy.ndarray phi_solid_nodes, + numpy.ndarray distance_to_solids, int nParticles, int nElements_owned, double particle_nitsche, double particle_epsFact, double particle_alpha, double particle_beta, - double particle_penalty_constant): + double particle_penalty_constant, + int use_pseudo_penalty): cdef numpy.ndarray rowptr, colind, globalJacobian_a (rowptr, colind, globalJacobian_a) = globalJacobian.getCSRrepresentation() self.thisptr.calculateJacobian(NONCONSERVATIVE_FORM, @@ -1139,6 +1224,10 @@ cdef class cRANS2P2D_base: < int * > p_l2g.data, < int * > vel_l2g.data, < double * > p_dof.data, < double * > u_dof.data, < double * > v_dof.data, < double * > w_dof.data, + < double * > p_old_dof.data, + < double * > u_old_dof.data, + < double * > v_old_dof.data, + < double * > w_old_dof.data, < double * > g.data, useVF, < double * > vf.data, @@ -1238,13 +1327,26 @@ cdef class cRANS2P2D_base: < double * > ball_radius.data, < double * > ball_velocity.data, < double * > ball_angular_velocity.data, + < double * > ball_center_acceleration.data, + < double * > ball_angular_acceleration.data, + < double * > ball_density.data, + < double * > particle_signed_distances.data, + < double * > particle_signed_distance_normals.data, + < double * > particle_velocities.data, + < double * > particle_centroids.data, + < double * > ebq_global_phi_s.data, + < double * > ebq_global_grad_phi_s.data, + < double * > ebq_particle_velocity_s.data, + < double * > phi_solid_nodes.data, + < double * > distance_to_solids.data, nParticles, nElements_owned, particle_nitsche, particle_epsFact, particle_alpha, particle_beta, - particle_penalty_constant) + particle_penalty_constant, + use_pseudo_penalty) def calculateVelocityAverage(self, int nExteriorElementBoundaries_global, diff --git a/proteus/mprans/cRANS3PF.pyx b/proteus/mprans/cRANS3PF.pyx index 71ccfec812..f98b8d2a9f 100644 --- a/proteus/mprans/cRANS3PF.pyx +++ b/proteus/mprans/cRANS3PF.pyx @@ -2600,7 +2600,8 @@ cdef class RANS3PF2D: wStar_dMatrix.data, numDOFs_1D, NNZ_1D, - csrRowIndeces_1D.data,csrColumnOffsets_1D.data, + csrRowIndeces_1D.data, + csrColumnOffsets_1D.data, rowptr_1D.data, colind_1D.data, isBoundary_1D.data, diff --git a/proteus/tests/AddedMass/test_addedmass2D.py b/proteus/tests/AddedMass/test_addedmass2D.py index 9f81ffaec3..9b80c4653c 100644 --- a/proteus/tests/AddedMass/test_addedmass2D.py +++ b/proteus/tests/AddedMass/test_addedmass2D.py @@ -85,8 +85,8 @@ def test_AddedMass_2D(self): ns = NumericalSolution.NS_base(so,pList,nList,so.sList,opts) ns.calculateSolution('addedmass2D') Aij = ns.so.ct.body.Aij - npt.assert_almost_equal(Aij[0,0], 500.1800561050536, decimal=5) - npt.assert_almost_equal(Aij[1,1], 1299.5160407630476, decimal=4) + npt.assert_almost_equal(Aij[0,0], 231.9854077642408, decimal=5) + npt.assert_almost_equal(Aij[1,1], 600.4898230969727, decimal=4) self.teardown_method(self) if __name__ == "__main__": diff --git a/proteus/tests/AddedMass/test_addedmass3D.py b/proteus/tests/AddedMass/test_addedmass3D.py index 788f2d0f8e..769faf7b26 100644 --- a/proteus/tests/AddedMass/test_addedmass3D.py +++ b/proteus/tests/AddedMass/test_addedmass3D.py @@ -85,8 +85,8 @@ def test_AddedMass_3D(self): ns = NumericalSolution.NS_base(so,pList,nList,so.sList,opts) ns.calculateSolution('addedmass3D') Aij = ns.so.ct.body.Aij - npt.assert_almost_equal(Aij[0,0], 284.27506488228266,decimal=5) - npt.assert_almost_equal(Aij[1,1], 281.1899393312541,decimal=5) + npt.assert_almost_equal(Aij[0,0], 92.23859079369586,decimal=5) + npt.assert_almost_equal(Aij[1,1], 94.3816831753334,decimal=5) self.teardown_method(self) if __name__ == "__main__": diff --git a/proteus/tests/CLSVOF/with_RANS3PF/pressure_p.py b/proteus/tests/CLSVOF/with_RANS3PF/pressure_p.py index 3963a2838a..7f002b95ce 100644 --- a/proteus/tests/CLSVOF/with_RANS3PF/pressure_p.py +++ b/proteus/tests/CLSVOF/with_RANS3PF/pressure_p.py @@ -10,7 +10,7 @@ from proteus.mprans import Pres name = "pressure" -#LevelModelType = Pres.LevelModel +LevelModelType = Pres.LevelModel coefficients=Pres.Coefficients(modelIndex=PRESSURE_model, fluidModelIndex=V_model, pressureIncrementModelIndex=PINC_model, diff --git a/proteus/tests/ProjScheme_with_EV/pressure_p.py b/proteus/tests/ProjScheme_with_EV/pressure_p.py index a2ed6898ef..7b62e95e1d 100644 --- a/proteus/tests/ProjScheme_with_EV/pressure_p.py +++ b/proteus/tests/ProjScheme_with_EV/pressure_p.py @@ -7,7 +7,7 @@ from proteus.mprans import Pres name = "pressure" -#LevelModelType = Pres.LevelModel +LevelModelType = Pres.LevelModel coefficients=Pres.Coefficients(modelIndex=PRESSURE_model, fluidModelIndex=V_model, pressureIncrementModelIndex=PINC_model, diff --git a/proteus/tests/cylinder2D/conforming_rans2p/comparison_files/T1_rans2p.h5 b/proteus/tests/cylinder2D/conforming_rans2p/comparison_files/T1_rans2p.h5 index 4d18f66d95..2615b4a354 100644 --- a/proteus/tests/cylinder2D/conforming_rans2p/comparison_files/T1_rans2p.h5 +++ b/proteus/tests/cylinder2D/conforming_rans2p/comparison_files/T1_rans2p.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:87c15c6d3254e63710f1628b6fd1f5effb03676562c70208115b7e808601ab88 -size 3088398 +oid sha256:3d1d90e75ecb932e10c54673030488d25057befdd34561aeb7aae240bc482ea4 +size 5035712 diff --git a/proteus/tests/cylinder2D/conforming_rans2p/cylinder2d.py b/proteus/tests/cylinder2D/conforming_rans2p/cylinder2d.py index 7d603cbc03..918b3a09b0 100644 --- a/proteus/tests/cylinder2D/conforming_rans2p/cylinder2d.py +++ b/proteus/tests/cylinder2D/conforming_rans2p/cylinder2d.py @@ -15,7 +15,7 @@ opts = Context.Options([ ("T", 4.0, "Time interval [0, T]"), - ("he",0.02, "maximum size of edges"), + ("he",0.01, "maximum size of edges"), ("backwardEuler",False,"use backward Euler or not"), ("onlySaveFinalSolution",False,"Only save the final solution") ], mutable=True) @@ -24,9 +24,7 @@ nd = 2 spaceOrder=1 -Refinement=1 useHex=False -points_on_grain = 21 DX = opts.he usePETSc = False#True @@ -86,7 +84,7 @@ T= opts.T runCFL = 0.9 dt_fixed = 0.005 -dt_init = 0.0025 +dt_init = 0.005 nDTout = int(old_div(T,dt_fixed)) dt_init = min(dt_init,0.5*dt_fixed) tnList = [0.0,dt_init]+[i*dt_fixed for i in range(1,nDTout+1)] @@ -116,9 +114,7 @@ # Gravity g = [0.0,0.0] -#triangleOptions="pAq30.0Dena%f" % (.5*DX**2) #% (0.5*(DX)**2,) -triangleOptions="pAq30.0Dena" -print(triangleOptions) +triangleOptions="pAq30ena"#D=Delaunay gives bad results for this composite meshing approach genMesh=True domain.writePLY('cylinder2D') domain.writePoly('cylinder2D') diff --git a/proteus/tests/cylinder2D/conforming_rans2p/cylinder_so.py b/proteus/tests/cylinder2D/conforming_rans2p/cylinder_so.py index 62d94d87db..8b93d90e5a 100644 --- a/proteus/tests/cylinder2D/conforming_rans2p/cylinder_so.py +++ b/proteus/tests/cylinder2D/conforming_rans2p/cylinder_so.py @@ -31,5 +31,4 @@ def __init__(self,modelList,system=defaultSystem,stepExact=True): # # tnList = [0.0,cylinder.dt_init]+[i*cylinder.dt_fixed for i in range(1,cylinder.nDTout+1)] tnList = cylinder.tnList -info = open("TimeList.txt","w") #archiveFlag = ArchiveFlags.EVERY_SEQUENCE_STEP diff --git a/proteus/tests/cylinder2D/conforming_rans2p/plot_lift_drag_from_pv_files.ipynb b/proteus/tests/cylinder2D/conforming_rans2p/plot_lift_drag_from_pv_files.ipynb index 296d9cfe5d..9187aa24af 100644 --- a/proteus/tests/cylinder2D/conforming_rans2p/plot_lift_drag_from_pv_files.ipynb +++ b/proteus/tests/cylinder2D/conforming_rans2p/plot_lift_drag_from_pv_files.ipynb @@ -1,1743 +1,3 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib notebook\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "directory=\"./\"" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "p_force = np.genfromtxt(directory+\"forceHistory_p.txt\", dtype=float)\n", - "v_force = np.genfromtxt(directory+\"forceHistory_v.txt\", dtype=float)" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "force=p_force+v_force" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "(801, 3)\n" - ] - } - ], - "source": [ - "print(force.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "N = force.shape[0]" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "scrolled": false - }, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure()\n", - "p1, = plt.plot(range(N), Cd, 'b', label='Cd')\n", - "p2, = plt.plot(range(N), Cl, 'y', label='Cl')\n", - "plt.axis([0, N, -1, 4])\n", - "plt.legend(handles=[p1,p2])\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "3.198053929138693" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "Cd[-100:].sum()/100.0" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 2", - "language": "python", - "name": "python2" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.11" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} +version https://git-lfs.github.com/spec/v1 +oid sha256:54992159a86fbdf6c4246ee0d62ebbe9845d2e426931c9f352e819e4ea98cb95 +size 575894 diff --git a/proteus/tests/cylinder2D/conforming_rans2p/twp_navier_stokes_cylinder_2d_p.py b/proteus/tests/cylinder2D/conforming_rans2p/twp_navier_stokes_cylinder_2d_p.py index 9df5d2f75e..b2f6d20a4d 100644 --- a/proteus/tests/cylinder2D/conforming_rans2p/twp_navier_stokes_cylinder_2d_p.py +++ b/proteus/tests/cylinder2D/conforming_rans2p/twp_navier_stokes_cylinder_2d_p.py @@ -35,11 +35,12 @@ def vel(x,t): U = Um*x[1]*(fl_H-x[1])/(old_div(fl_H,2.0))**2 -# if t < 2.0: -# return t*U/2.0 -# else: -# return U + if t < 2.0: + return t*U/2.0 + else: + return U return U + def getDBC_p(x,flag): if flag == boundaryTags['right']: return lambda x,t: 0.0 diff --git a/proteus/tests/cylinder2D/conforming_rans3p/comparison_files/T1P1.h5 b/proteus/tests/cylinder2D/conforming_rans3p/comparison_files/T1P1.h5 index 63d99553ee..7816360965 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/comparison_files/T1P1.h5 +++ b/proteus/tests/cylinder2D/conforming_rans3p/comparison_files/T1P1.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:46f3f62eb68a831f53b345cf250c3959533e95ace37b752db7715325d370edb1 -size 1469915 +oid sha256:6585b1ee70c978923c5ba859c36e18ee0181106c57024ef4299dec419f654446 +size 541846 diff --git a/proteus/tests/cylinder2D/conforming_rans3p/cylinder.py b/proteus/tests/cylinder2D/conforming_rans3p/cylinder.py index 969a698f3e..d12501cc62 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/cylinder.py +++ b/proteus/tests/cylinder2D/conforming_rans3p/cylinder.py @@ -88,7 +88,7 @@ pbasis = C0_AffineQuadraticOnSimplexWithNodalBasis -weak_bc_penalty_constant = 100.0 +weak_bc_penalty_constant = 10.0 nLevels = 1 #parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.node @@ -119,13 +119,13 @@ #domain.writeAsymptote("mesh") #triangleOptions = "VApq30Dena%8.8f" % ((he ** 2) / 2.0,) #triangleOptions = "VApq30Dena" -triangleOptions= "VApq30Dena" +triangleOptions= "Apq30ena" logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions, domain.polyfile + ".poly")) # Time stepping T=opts.T dt_fixed = 0.005#0.03 -dt_init = 0.0025#min(0.1*dt_fixed,0.001) +dt_init = 0.005#min(0.1*dt_fixed,0.001) runCFL=0.33 nDTout = int(round(old_div(T,dt_fixed))) tnList = [0.0,dt_init]+[i*dt_fixed for i in range(1,nDTout+1)] @@ -138,7 +138,7 @@ ns_forceStrongDirichlet = False ns_sed_forceStrongDirichlet = False if useMetrics: - ns_shockCapturingFactor = 0.5 + ns_shockCapturingFactor = 0.0 ns_lag_shockCapturing = True ns_lag_subgridError = True ls_shockCapturingFactor = 0.5 diff --git a/proteus/tests/cylinder2D/conforming_rans3p/cylinder_so.py b/proteus/tests/cylinder2D/conforming_rans3p/cylinder_so.py index 85ced6483e..7766acfc52 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/cylinder_so.py +++ b/proteus/tests/cylinder2D/conforming_rans3p/cylinder_so.py @@ -88,7 +88,7 @@ def __init__(self,modelList,system=defaultSystem,stepExact=True): needEBQ_GLOBAL = False needEBQ = False -# systemStepExact=False +systemStepExact=False tnList = cylinder.tnList diff --git a/proteus/tests/cylinder2D/conforming_rans3p/plot_lift_drag_from_pv_files.ipynb b/proteus/tests/cylinder2D/conforming_rans3p/plot_lift_drag_from_pv_files.ipynb index 009f028071..cd41641937 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/plot_lift_drag_from_pv_files.ipynb +++ b/proteus/tests/cylinder2D/conforming_rans3p/plot_lift_drag_from_pv_files.ipynb @@ -1,1743 +1,3 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib notebook\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "directory=\"./\"" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "p_force = np.genfromtxt(directory+\"forceHistory_p.txt\", dtype=float)\n", - "v_force = np.genfromtxt(directory+\"forceHistory_v.txt\", dtype=float)" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "force=p_force+v_force" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "(801, 3)\n" - ] - } - ], - "source": [ - "print(force.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "N = force.shape[0]" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "scrolled": false - }, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure()\n", - "p1, = plt.plot(range(N), Cd, 'b', label='Cd')\n", - "p2, = plt.plot(range(N), Cl, 'y', label='Cl')\n", - "plt.axis([0, N, -1, 4])\n", - "plt.legend(handles=[p1,p2])\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "3.2646816691441867" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "Cd[-100:].sum()/100.0" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 2", - "language": "python", - "name": "python2" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.11" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} +version https://git-lfs.github.com/spec/v1 +oid sha256:a08bf0edf6d47f277cb9ebc9dbe4888643c787fe3c32d8235f1bca65d4fce14a +size 57932 diff --git a/proteus/tests/cylinder2D/conforming_rans3p/pressureInitial_p.py b/proteus/tests/cylinder2D/conforming_rans3p/pressureInitial_p.py index f57e9a355b..09e7e17c1d 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/pressureInitial_p.py +++ b/proteus/tests/cylinder2D/conforming_rans3p/pressureInitial_p.py @@ -19,6 +19,8 @@ fluidModelIndex=V_model, pressureModelIndex=PRESSURE_model) +#LevelModelType = PresInit.LevelModel + #pressure increment should be zero on any pressure dirichlet boundaries def getDBC_pInit(x,flag): if flag == boundaryTags['right']: diff --git a/proteus/tests/cylinder2D/conforming_rans3p/pressureincrement_p.py b/proteus/tests/cylinder2D/conforming_rans3p/pressureincrement_p.py index eae82d9747..49d1094da2 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/pressureincrement_p.py +++ b/proteus/tests/cylinder2D/conforming_rans3p/pressureincrement_p.py @@ -38,9 +38,17 @@ def getDBC_phi(x,flag): from math import cos,pi #the advectiveFlux should be zero on any no-flow boundaries +def flux(x,t): + U = Um*x[1]*(fl_H-x[1])/(old_div(fl_H,2.0))**2 + if t < 2.0: + return -t*U/2.0 + else: + return -U + return -U + def getAdvectiveFlux_qt(x,flag): if flag == boundaryTags['left']: - return lambda x,t: -Um*x[1]*(fl_H-x[1])/(old_div(fl_H,2.0))**2 + return flux elif flag == boundaryTags['right']: return None else: diff --git a/proteus/tests/cylinder2D/conforming_rans3p/twp_navier_stokes_p.py b/proteus/tests/cylinder2D/conforming_rans3p/twp_navier_stokes_p.py index 086b459573..f3314338e9 100644 --- a/proteus/tests/cylinder2D/conforming_rans3p/twp_navier_stokes_p.py +++ b/proteus/tests/cylinder2D/conforming_rans3p/twp_navier_stokes_p.py @@ -52,16 +52,24 @@ turbulenceClosureModel=ns_closure, movingDomain=movingDomain, dragAlpha=dragAlpha, - PSTAB=0.0) - + PSTAB=0.0, + USE_SUPG=0.0, + ARTIFICIAL_VISCOSITY=1) +def vel(x,t): + U = Um*x[1]*(fl_H-x[1])/(old_div(fl_H,2.0))**2 + if t < 2.0: + return t*U/2.0 + else: + return U + return U def getDBC_u(x,flag): if flag in[ boundaryTags['left']]: - return lambda x,t: Um*x[1]*(fl_H-x[1])/(old_div(fl_H,2.0))**2 + return vel elif flag in [boundaryTags['front'],boundaryTags['back'],boundaryTags['top'],boundaryTags['bottom'],boundaryTags['obstacle']]: return lambda x,t: 0.0 elif ns_forceStrongDirichlet==False and flag == boundaryTags['right']: - return lambda x,t: Um*x[1]*(fl_H-x[1])/(old_div(fl_H,2.0))**2 + return vel def getDBC_v(x,flag): if flag in[boundaryTags['left']]: diff --git a/proteus/tests/cylinder2D/ibm_method/comparison_files/T1_rans3p.h5 b/proteus/tests/cylinder2D/ibm_method/comparison_files/T1_rans3p.h5 index 48829dbc6c..3ed18a69be 100644 --- a/proteus/tests/cylinder2D/ibm_method/comparison_files/T1_rans3p.h5 +++ b/proteus/tests/cylinder2D/ibm_method/comparison_files/T1_rans3p.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:94cc187470745a4584eba949533b19b8fdfbef2dad197de6cf0359a19eda075f -size 2886949 +oid sha256:afd45ee2e07c1d2abec744ef4e2aeacb35c557dc4f586a13bc0209b1d1c91526 +size 1279644 diff --git a/proteus/tests/cylinder2D/ibm_method/cylinder.py b/proteus/tests/cylinder2D/ibm_method/cylinder.py index c2b4801cc2..0b2a8732e1 100644 --- a/proteus/tests/cylinder2D/ibm_method/cylinder.py +++ b/proteus/tests/cylinder2D/ibm_method/cylinder.py @@ -8,20 +8,18 @@ from proteus.default_n import * from proteus.Profiling import logEvent - from proteus import Context ct = Context.Options([ ("T", 4.0, "Time interval [0, T]"), - ("Refinement",4, "refinement"), + ("he",0.04, "maximum size of edges"), ("onlySaveFinalSolution",False,"Only save the final solution"), ("vspaceOrder",2,"FE space for velocity"), - ("pspaceOrder",1,"FE space for pressure") + ("pspaceOrder",1,"FE space for pressure"), + ("use_regions",False, "use refinement regions") ], mutable=True) # Discretization -- input options -#Refinement = 20#45min on a single core for spaceOrder=1, useHex=False -Refinement = ct.Refinement sedimentDynamics=False genMesh = True movingDomain = False @@ -93,12 +91,7 @@ # Domain and mesh #L = (0.584,0.350) L = (2.2, 0.41) -he = old_div(L[0],float(4*Refinement-1)) -he*=0.5 -he*=0.5 -#he*=0.5 -#he*=0.5 -#he*=0.5 +he = ct.he weak_bc_penalty_constant = 100.0 nLevels = 1 #parallelPartitioningType = proteus.MeshTools.MeshParallelPartitioningTypes.element @@ -119,98 +112,136 @@ # lineGauges_phi = LineGauges_phi(lineGauges.endpoints, linePoints=20) if useHex: - nnx = 4 * Refinement + 1 - nny = 2 * Refinement + 1 + nnx = ceil(L[0]/ct.he) + nny = ceil(L[1]/ct.he) 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 + nnx = ceil(L[0]/ct.he) + nny = ceil(L[1]/ct.he) 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 + if ct.use_regions: + 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]] + [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*(old_div(he,2.0))**2,0.5*(old_div(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) + 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*(old_div(he,2.0))**2,0.5*(old_div(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) + else: + vertices = [[0.0, 0.0], #0 + [L[0], 0.0], #1 + [L[0], L[1]], #2 + [0.0, L[1]]] #3 + vertexFlags = [boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['top'], + boundaryTags['top']] + segments = [[0, 1], + [1, 2], + [2, 3], + [3, 0]] + segmentFlags = [boundaryTags['bottom'], + boundaryTags['right'], + boundaryTags['top'], + boundaryTags['left']] + + regions = [[0.95*L[0], 0.2]] + regionFlags = [1] + regionConstraints=[0.5*he**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 = "VApq30ena%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 = 0.005#0.03 -dt_init = 0.0025#min(0.1*dt_fixed,0.001) +dt_init = 0.005#min(0.1*dt_fixed,0.001) runCFL=0.33 nDTout = int(round(old_div(T,dt_fixed))) tnList = [0.0,dt_init]+[i*dt_fixed for i in range(1,nDTout+1)] @@ -218,10 +249,10 @@ if ct.onlySaveFinalSolution == True: tnList = [0.0,dt_init,ct.T] # Numerical parameters -ns_forceStrongDirichlet = True +ns_forceStrongDirichlet = False ns_sed_forceStrongDirichlet = False if useMetrics: - ns_shockCapturingFactor = 0.5 + ns_shockCapturingFactor = 0.0 ns_lag_shockCapturing = True ns_lag_subgridError = True ls_shockCapturingFactor = 0.5 @@ -328,13 +359,10 @@ U = 1.5 # this is the inlet max velocity not the mean velocity def velRamp(t): - return U -# if t < 0.25: -# return U*(1.0+math.cos((t-0.25)*math.pi/0.25))/2.0 -# else: -# return U - - + if t < 2.0: + return t*U/2.0 + else: + return U def signedDistance(x): return x[1]-old_div(L[1],2) diff --git a/proteus/tests/cylinder2D/ibm_method/cylinder_so.py b/proteus/tests/cylinder2D/ibm_method/cylinder_so.py index 80cc5799e5..1eec5089e8 100644 --- a/proteus/tests/cylinder2D/ibm_method/cylinder_so.py +++ b/proteus/tests/cylinder2D/ibm_method/cylinder_so.py @@ -90,7 +90,7 @@ def __init__(self,modelList,system=defaultSystem,stepExact=True): # systemStepControllerType = Sequential_FixedStep #Sequential_FixedStep_Simple # uses time steps in so.tnList # dt_system_fixed = 0.01; -# systemStepExact=False; +systemStepExact=False needEBQ_GLOBAL = False diff --git a/proteus/tests/cylinder2D/ibm_method/plot_lift_drag_force_from_particle.ipynb b/proteus/tests/cylinder2D/ibm_method/plot_lift_drag_force_from_particle.ipynb index 850b16b4f0..6c96d02c80 100644 --- a/proteus/tests/cylinder2D/ibm_method/plot_lift_drag_force_from_particle.ipynb +++ b/proteus/tests/cylinder2D/ibm_method/plot_lift_drag_force_from_particle.ipynb @@ -1,1734 +1,3 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib notebook\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "directory=\"./\"" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "force =np.genfromtxt(directory+\"particle_forceHistory.txt\", dtype=float)" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "N = force.shape[0]" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure()\n", - "p1, = plt.plot(range(N), Cd, 'b', label='Cd')\n", - "p2, = plt.plot(range(N), Cl, 'y', label='Cl')\n", - "plt.axis([0, N, -1, 4])\n", - "plt.legend(handles=[p1,p2])\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "2.852776035773133" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "Cd[-100:].sum()/100.0" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "2.914198732616997" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "2.914198732616997" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 2", - "language": "python", - "name": "python2" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.11" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} +version https://git-lfs.github.com/spec/v1 +oid sha256:09131e29051cffc22a9aa2d37ef734addbae712e17f8d49a8d1719c5c3c5879c +size 310787 diff --git a/proteus/tests/cylinder2D/ibm_method/pressure_p.py b/proteus/tests/cylinder2D/ibm_method/pressure_p.py index cdd1cdbc8f..f66a6e4e64 100644 --- a/proteus/tests/cylinder2D/ibm_method/pressure_p.py +++ b/proteus/tests/cylinder2D/ibm_method/pressure_p.py @@ -11,6 +11,8 @@ name = "pressure" +LevelModelType = Pres.LevelModel + coefficients=Pres.Coefficients(modelIndex=PRESSURE_model, fluidModelIndex=V_model, pressureIncrementModelIndex=PINC_model, diff --git a/proteus/tests/cylinder2D/ibm_method/twp_navier_stokes_p.py b/proteus/tests/cylinder2D/ibm_method/twp_navier_stokes_p.py index c4eec31d2e..e3a4ab0a82 100644 --- a/proteus/tests/cylinder2D/ibm_method/twp_navier_stokes_p.py +++ b/proteus/tests/cylinder2D/ibm_method/twp_navier_stokes_p.py @@ -50,12 +50,13 @@ turbulenceClosureModel=ns_closure, movingDomain=movingDomain, dragAlpha=dragAlpha, + USE_SUPG=0.0, PSTAB=0.0, nParticles=1, particle_epsFact=1.5, - particle_alpha=1e6, - particle_beta=1e6, - particle_penalty_constant=100.0, + particle_alpha=1000.0, + particle_beta=1000.0, + particle_penalty_constant=10.0, particle_sdfList=[particle_sdf], particle_velocityList=[particle_vel], use_ball_as_particle=1, diff --git a/proteus/tests/cylinder2D/ibm_rans2p/comparison_files/T1_ibm_rans2p.h5 b/proteus/tests/cylinder2D/ibm_rans2p/comparison_files/T1_ibm_rans2p.h5 index 15bf7d8e74..3d0bcd8a16 100644 --- a/proteus/tests/cylinder2D/ibm_rans2p/comparison_files/T1_ibm_rans2p.h5 +++ b/proteus/tests/cylinder2D/ibm_rans2p/comparison_files/T1_ibm_rans2p.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:21a78c1d8677c97c0ef30e6e26d236b39a4a40a2abdcca2891d860419acd5a88 -size 119042 +oid sha256:96819de8f3656dfd5dc3f124e1d0d2507c39c304168165a8c5039a46e40f5c59 +size 1299402 diff --git a/proteus/tests/cylinder2D/ibm_rans2p/cylinder.py b/proteus/tests/cylinder2D/ibm_rans2p/cylinder.py index f9d90c7868..62ca0d2c3a 100644 --- a/proteus/tests/cylinder2D/ibm_rans2p/cylinder.py +++ b/proteus/tests/cylinder2D/ibm_rans2p/cylinder.py @@ -8,22 +8,23 @@ from proteus import Context ct = Context.Options([ - ("T", 3.0, "Time interval [0, T]"), + ("T", 4.0, "Time interval [0, T]"), + ("he",0.02, "maximum size of edges"), ("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"), + ("StrongDirichlet", False,"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"), + ("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"), + ("use_regions",False,"Use refinement regions"), ], mutable=True) @@ -39,8 +40,6 @@ # 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 @@ -155,101 +154,134 @@ #=============================================================================== #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 +he = ct.he structured = False if useHex: - nnx = 4 * Refinement + 1 - nny = 2 * Refinement + 1 + nnx = ceil(L[0]/ct.he) + nny = ceil(L[1]/ct.he) 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 + nnx = ceil(L[0]/ct.he) + nny = ceil(L[1]/ct.he) 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]] + if ct.use_regions: + 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) + 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*(old_div(he,2.0))**2,0.5*(old_div(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) + else: + vertices = [[0.0, 0.0], #0 + [L[0], 0.0], #1 + [L[0], L[1]], #2 + [0.0, L[1]]] #3 + vertexFlags = [boundaryTags['bottom'], + boundaryTags['bottom'], + boundaryTags['top'], + boundaryTags['top']] + segments = [[0, 1], + [1, 2], + [2, 3], + [3, 0]] + segmentFlags = [boundaryTags['bottom'], + boundaryTags['right'], + boundaryTags['top'], + boundaryTags['left']] + + regions = [[0.95*L[0], 0.2]] + regionFlags = [1] + regionConstraints=[0.5*he**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 = "VApq30ena%8.8f" % ((he ** 2) / 2.0,) triangleOptions = "VApq30Dena" logEvent("""Mesh generated using: tetgen -%s %s""" % (triangleOptions, domain.polyfile + ".poly")) @@ -258,11 +290,9 @@ #=============================================================================== T=ct.T dt_fixed = ct.dt_fixed#0.03 -dt_init = 0.5*dt_fixed#min(0.1*dt_fixed,0.001) +dt_init = ct.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');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp = $('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " event.shiftKey = false;\n", - " // Send a \"J\" for go to next cell\n", - " event.which = 74;\n", - " event.keyCode = 74;\n", - " manager.command_mode();\n", - " manager.handle_keydown(event);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.plot(range(N), force[:,0], 'b', range(N),force[:,1], 'y')\n", - "plt.axis([0, N, -1, 1])\n", - "plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "D = 0.1\n", - "U = 1.0\n", - "rho=1.0\n", - "Cd = force[:,0]*2/rho/D/U/U\n", - "Cl = force[:,1]*2/rho/D/U/U" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " fig.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('