📄 sph.cpp
字号:
ardvxdt=new double *[dim+1];for(i=0;i<dim+1;i++)ardvxdt[i]=new double[maxn+1];
for(i=1;i<=ntotal;i++)
{
avdudt[i] = 0;
ahdudt[i] = 0;
for(int d=1;d<=dim;d++)
{
indvxdt[d][i] = 0;
ardvxdt[d][i] = 0;
exdvxdt[d][i] = 0;
}
}
//c--- Positions of virtual (boundary) particles:
nvirt = 0;
// c--- Interaction parameters, calculating neighboring particles and optimzing smoothing length
if (nnps==1)direct_find(ntotal+nvirt);
if (summation_density)sum_density(ntotal+nvirt);
//c--- Internal forces:
int_force(ntotal+nvirt,indvxdt,tdsdt,du); //这可能是程序出错的地方，如果indvxdt,du改变将导致dedt,dvxdt同时改变。
//c--- Artificial viscosity:
if (visc_artificial)art_visc(ntotal+nvirt,ardvxdt,avdudt);//其中c表示温度还是声速，这里以声速计算
//c--- Convert velocity, force, and energy to f and dfdt
for( i=1;i<=ntotal;i++)
{
// cout<<indvxdt[1][i]<<" "<<du[i]<<endl;
for(int d=1;d<=dim;d++)
{
dvx[d][i] = indvxdt[d][i] + exdvxdt[d][i] + ardvxdt[d][i];
}
du[i] = du[i] + avdudt[i] + ahdudt[i];
}
// for(i=1;i<=ntotal;i++)cout<<"single_step "<<"dvx[1][i]="<<dvx[1][i]<<"du[i]="<<du[i]<<endl;
// if ((itimestep%print_step).eq.0) then
// write(*,*)
// write(*,*) '**** Information for particle ****',
// & moni_particle
// write(*,101)'internal a ','artifical a=',
// & 'external a ','total a '
// write(*,100)indvxdt(1,moni_particle),ardvxdt(1,moni_particle),
// & exdvxdt(1,moni_particle),dvx(1,moni_particle)
// endif
//101 format(1x,4(2x,a12))
//100 format(1x,4(2x,e12.6))
// end
delete []ns; delete []avdudt; delete []ahdudt;
for(i=0;i<dim+1;i++)
delete []indvxdt[i];
delete []indvxdt;
for( i=0;i<dim+1;i++)
delete []exdvxdt[i];
delete []exdvxdt;
for(i=0;i<dim+1;i++)
delete []ardvxdt[i];
delete []ardvxdt;
}
void SPH::time_integration()
{
// int current_ts;
int nstart=0;//源程序中没有定义，这肯定是个问题，这里先定义为0；
double **x_min, **v_min,*u_min,*rho_min;
double **dx1;//这个数组在程序中并没有使用，为区别*dx，改为**dx1；
u_min=new double[maxn+1];rho_min=new double[maxn+1];
x_min=new double *[dim+1];for(int i=0;i<dim+1;i++)x_min[i]=new double[maxn+1];v_min=new double *[dim+1];for(i=0;i<dim+1;i++)v_min[i]=new double[maxn+1];
dx1=new double *[dim+1];for(i=0;i<dim+1;i++)dx1[i]=new double[maxn+1];
double time, temp_rho, temp_u;
time=0;//统计迭代次数
for( i = 1;i<=ntotal;i++){for(int d = 1;d<=dim;d++) av[d][i] = 0;}
for(int j = nstart+1; j<=nstart+maxtimestep;j++)
{
itimestep=itimestep+1;
// current_ts=current_ts+1;
//// if ((itimestep%print_step)==0)
//// {
//// cout<<"______________________________________________"<<endl;
//// cout<<" current number of time step ="<< itimestep<<" "<<"current time="<<double(time+dt)<<endl;
//// cout<<"______________________________________________"<<endl;
//// //c If not first time step, then update thermal energy, density and
//// //c velocity half a time step
//// }
if (itimestep!=1)
{
for(int i = 1;i<= ntotal;i++)
{
u_min[i] = u[i];
temp_u=0.0;
if (dim==1)
temp_u=-nsym*p[i]*vx[1][i]/x[1][i]/rho[i];
u[i] = u[i] + (dt/2.0)* (du[i]+temp_u);
if(u[i]<0)
u[i] = 0;
if (!summation_density)
{
rho_min[i] = rho[i];
temp_rho=0.0;
if (dim==1)
temp_rho=-nsym*rho[i]*vx[1][i]/x[1][i];
rho[i] = rho[i] +(dt/2.0)*( drho[i]+ temp_rho);
}
for(int d = 1;d<=dim;d++)
{
v_min[d][i] = vx[d][i];
vx[d][i] = vx[d][i] + (dt/2.0)*dvx[d][i];
}
}
}
//c--- Definition of variables out of the function vector:
single_step();
if (itimestep==1)
{
for(int i=1;i<=ntotal;i++)
{
temp_u=0;
if (dim==1)
temp_u=-nsym*p[i]*vx[1][i]/x[1][i]/rho[i];
u[i] = u[i] + (dt/2.0)*(du[i] + temp_u);
if(u[i]<0)
u[i] = 0.0;
if(!summation_density )
{
temp_rho=0;
if (dim==1)
temp_rho=-nsym*rho[i]*vx[1][i]/x[1][i];
rho[i] = rho[i] + (dt/2.0)* (drho[i]+temp_rho);
}
for(int d = 1;d<=dim;d++)
{
vx[d][i] = vx[d][i] + (dt/2.0) * dvx[d][i] + av[d][i];
x[d][i] = x[d][i] + dt * vx[d][i];
}
}
}
else
{
for(int i=1;i<=ntotal;i++)
{
temp_u=0.0;
if (dim==1)
temp_u=-nsym*p[i]*vx[1][i]/x[1][i]/rho[i];
u[i] = u_min[i] + dt*(du[i]+temp_u);
if(u[i]<0)
u[i] = 0;
if (!summation_density )
{
temp_rho=0;
if (dim==1)
temp_rho=-nsym*rho[i]*vx[1][i]/x[1][i];
rho[i] = rho_min[i] + dt*(drho[i]+temp_rho);
}
for(int d = 1;d<= dim;d++)
{
vx[d][i] = v_min[d][i] + dt * dvx[d][i] + av[d][i];
x[d][i] = x[d][i] + dt * vx[d][i];
}
}
}
//for(int i=1;i<=ntotal;i++)
//{
//cout<<"time_integration1 "<<"x[1][i]="<<x[1][i]<<"vx[1][i]="<<vx[1][i]<<"rho[i]="<<rho[i]<<"p[i]="<<p[i]<<endl;
//cout<<"time_integration2 "<<"u[i]="<<u[i]<<"c[i]="<<c[i]<<"e[i]="<<e[i]<<"hsml[i]="<<hsml[i]<<endl;
//}
// time = time + dt;//time统计迭代次数？？？
// if ((itimestep%save_step)==0)
// {
// output(x, vx, mass, rho, p, u, c, itype, hsml, ntotal);//??????????????????
// }
//// if ((itimestep%print_step)==0){cout<<"x(1,moni_particle)="<<x[1][moni_particle]<<" "<<"vx(1,moni_particle)="<<x[1][moni_particle]<<"dvx(1,moni_particle) dvx="<<dvx[1][moni_particle]<<endl<<endl;}
}
// nstart=current_ts;
delete []u_min;delete []rho_min;
for( i=0;i<dim+1;i++)delete []x_min[i];delete []x_min;
for( i=0;i<dim+1;i++)delete []v_min[i];delete []v_min;
for( i=0;i<dim+1;i++)delete []dx1[i];delete []dx1;
}
void SPH::input()
{
//此处省略了读入文件数据，下面是生成数据的程序
if(shocktube)
{
//subroutine shock_tube(x, vx, mass, rho, p, u, itype, hsml, ntotal)
//c----------------------------------------------------------------------
//c This subroutine is used to generate initial data for the
//c 1 d noh shock tube problem
//c x-- coordinates of particles [out]
//c vx-- velocities of particles [out]
//c mass-- mass of particles [out]
//c rho-- dnesities of particles [out]
//c p-- pressure of particles [out]
//c u-- internal energy of particles [out]
//c itype-- types of particles [out]
//c =1 ideal gas
//c hsml-- smoothing lengths of particles [out]
//c ntotal-- total particle number [out]
double space_x;
ntotal=400;
space_x=0.6/80;
for(int i=1;i<=ntotal;i++)
{
mass[i]=0.75/400.0;
hsml[i]=0.015;
itype[i]=1;
for(int d = 1;d<=dim;d++)
{
x[d][i] = 0.0;
vx[d][i] = 0.0;
}
}
for( i=1;i<=320;i++)
x[1][i]=-0.6+space_x/4.0*(i-1);
for( i=320+1;i<=ntotal;i++)
x[1][i]=0.0+space_x*(i-320);
for( i=1;i<=ntotal;i++)
{
if (x[1][i]<=0.00000001) //1.e-8???????????????????????
{
u[i]=2.5;
rho[i]=1.0;
p[i]=1.0;
}
if (x[1][i]>0.00000001)//1.e-8??????????????????????
{
u[i]=1.795;
rho[i]=0.25;
p[i]=0.1795;
}
}
}
if(shearcavity)
{
// subroutine shear_cavity(x, vx, mass, rho, p, u, itype, hsml, ntotal)
//c----------------------------------------------------------------------
//c This subroutine is used to generate initial data for the
//c 2 d shear driven cavity probem with Re = 1
//c x-- coordinates of particles [out]
//c vx-- velocities of particles [out]
//c mass-- mass of particles [out]
//c rho-- dnesities of particles [out]
//c p-- pressure of particles [out]
//c u-- internal energy of particles [out]
//c itype-- types of particles [out]
//c =2 water
//c h-- smoothing lengths of particles [out]
//c ntotal-- total particle number [out]
int m, n, mp, np, k;
double xl, yl, dx, dy;
//c Giving mass and smoothing length as well as other data.
m = 41;
n = 41;
mp = m-1;
np = n-1;
ntotal = mp * np;
xl = 0.001;// 1.e-3????????????
yl = 0.001;// 1.e-3????????????
dx = xl/mp;
dy = yl/np;
for(int i = 1;i<= mp;i++)
{
for(int j = 1;j<= np;j++)
{
k = j + (i-1)*np;x[1][k] = (i-1)*dx + dx/2.0;x[2][k] = (j-1)*dy + dy/2.0;
}
}
for(i = 1;i<=mp*np;i++)
{
vx[1][i] = 0;
vx[2][i] = 0;
rho[i] = 1000;
mass[i] = dx*dy*rho[i];
p[i]= 0;
u[i]=357.1;
itype[i] = 2;
hsml[i] = dx;
}
}
//c,s,e均没有赋初始值,临时赋零,有待调整;????????????????????????????????????????
// for(int i=1;i<=ntotal+1;i++){c[i]=0;s[i]=0;e[i]=0;}
time_integration();
output();
}
void SPH::output()
{
// subroutine output(x, vx, mass, rho, p, u, c, itype, hsml, ntotal)
//c----------------------------------------------------------------------
//c Subroutine for saving particle information to external disk file
//c x-- coordinates of particles [in]
//c vx-- velocities of particles [in]
//c mass-- mass of particles [in]
//c rho-- dnesities of particles [in]
//c p-- pressure of particles [in]
//c u-- internal energy of particles [in]
//c c-- sound velocity of particles [in]
//c itype-- types of particles [in]
//c hsml-- smoothing lengths of particles [in]
//c ntotal-- total particle number [in]
ofstream out1("f_xv.plt");//创建输出流并创建打开文件
for(int i=1;i<=ntotal;i++)//分行列存储
{
out1<<setw(10)<<setprecision(30)<<i<<" ";
out1<<setw(10)<<setprecision(30)<<x[1][i]<<" ";
// out1<<setw(10)<<setprecision(30)<<x[2][i]<<" ";
// out1<<setw(10)<<setprecision(30)<<x[3][i]<<" ";
out1<<setw(10)<<setprecision(30)<<vx[1][i]<<endl;
// out1<<setw(10)<<setprecision(30)<<vx[2][i]<<" ";
// out1<<setw(10)<<setprecision(30)<<vx[3][i]<<endl;
}
ofstream out2("f_state.plt");//创建输出流并创建打开文件
for( i=1;i<=ntotal;i++)//分行列存储
{
out2<<setw(10)<<setprecision(30)<<i<<" ";
out2<<setw(10)<<setprecision(30)<<x[1][i]<<" ";
out2<<setw(10)<<setprecision(30)<<mass[i]<<" ";
out2<<setw(10)<<setprecision(30)<<rho[i]<<" ";
out2<<setw(10)<<setprecision(30)<<p[i]<<" ";
out2<<setw(10)<<setprecision(30)<<u[i]<<endl;
}
ofstream out3("f_other.plt");//创建输出流并创建打开文件
for( i=1;i<=ntotal;i++)//分行列存储
{
out3<<setw(10)<<setprecision(30)<<i<<" ";
out3<<setw(10)<<setprecision(30)<<itype[i]<<" ";
out3<<setw(10)<<setprecision(30)<<hsml[i]<<endl;
}
}
void main()
{
// program SPH
//c----------------------------------------------------------------------
//c This is a three dimensional SPH code. the followings are the
//c basic parameters needed in this codeor calculated by this code
//c mass-- mass of particles [in]
//c ntotal-- total particle number ues [in]
//c dt--- Time step used in the time integration [in]
//c itype-- types of particles [in]
//c x-- coordinates of particles [in/out]
//c vx-- velocities of particles [in/out]
//c rho-- dnesities of particles [in/out]
//c p-- pressure of particles [in/out]
//c u-- internal energy of particles [in/out]
//c hsml-- smoothing lengths of particles [in/out]
//c c-- sound velocity of particles [out]
//c s-- entropy of particles [out]
//c e-- total energy of particles [out]
int dim1=1;
int maxn1=12000;
// double s1, s2,dt;//用于计算CPU的时间
// call time_print
// call time_elapsed(s1)
SPH test(dim1,maxn1);
if (test.shocktube) test.dt = 0.005;
if (test.shearcavity) test.dt = 0.00005;
cin>>test.maxtimestep;
test.input();
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -