📄 sph.cpp
字号:
txy[j] = txy[j] + mass[i]*hxy/rho[i];
tyy[i] = tyy[i] + mass[j]*hyy/rho[j];
tyy[j] = tyy[j] + mass[i]*hyy/rho[i];
}
else if (dim==3)
{
txx[i] = txx[i] + mass[j]*hxx/rho[j];
txx[j] = txx[j] + mass[i]*hxx/rho[i];
txy[i] = txy[i] + mass[j]*hxy/rho[j];
txy[j] = txy[j] + mass[i]*hxy/rho[i];txz[i] = txz[i] + mass[j]*hxz/rho[j];txz[j] = txz[j] + mass[i]*hxz/rho[i];
tyy[i] = tyy[i] + mass[j]*hyy/rho[j];tyy[j] = tyy[j] + mass[i]*hyy/rho[i];tyz[i] = tyz[i] + mass[j]*hyz/rho[j];
tyz[j] = tyz[j] + mass[i]*hyz/rho[i];tzz[i] = tzz[i] + mass[j]*hzz/rho[j];tzz[j] = tzz[j] + mass[i]*hzz/rho[i];
}
hvcc = 0;// Calculate SPH sum for vc,c = dvx1/dx + dvy/dy + dvz/dz:
for(d=1;d<=dim;d++)
hvcc = hvcc + dvx1[d]*dwdx[d][k];
vcc[i] = vcc[i] + mass[j]*hvcc/rho[j];
vcc[j] = vcc[j] + mass[i]*hvcc/rho[i];
}
}
for(i=1;i<=ntotal;i++)
{
if (visc)//c Viscous entropy Tds/dt = 1/2 eta/rho Tab Tab
{
if (dim==1)
{
tdsdt[i] = txx[i]*txx[i];
}
else if (dim==2)
{
tdsdt[i] = txx[i]*txx[i] + 2.0*txy[i]*txy[i]+ tyy[i]*tyy[i];
}
else if (dim==3)
{
tdsdt[i] = txx[i]*txx[i] + 2.0*txy[i]*txy[i]+ 2.0*txz[i]*txz[i]+ tyy[i]*tyy[i] + 2.0*tyz[i]*tyz[i]+ tzz[i]*tzz[i];
}
tdsdt[i] = 0.50*eta[i]/rho[i]*tdsdt[i];
}
//c Pressure from equation of state
if (abs(itype[i])==1)//粒子类型是在input中设定的
{
temp1=rho[i];temp2=u[i];
//temp3=p[i];temp4=c[i];
p_gas(temp1,temp2,temp3,temp4);
p[i]=temp3;c[i]=temp4;
}
else if (abs(itype[i])==2)
{
temp1=rho[i];
//temp3=p[i];temp4=c[i];
p_art_water(temp1,temp3,temp4);
p[i]=temp3;c[i]=temp4;
}
}
for(int k=1;k<=niac;k++) //c Calculate SPH sum for pressure force -p,a/rho and viscous force (eta Tab),b/rho and the internal energy change de/dt due to -p/rho vc,c
{
i = pair_i[k];
int j = pair_j[k];
he = 0.0;
rhoij = 1.0/(rho[i]*rho[j]); //c For SPH algorithm 1
if(pa_sph==1)
{
for(int d=1;d<=dim;d++)
{
h = -(p[i] + p[j])*dwdx[d][k];
he = he + (vx[d][j] - vx[d][i])*h; //c Pressure part
//c Viscous force
if (visc)
{
//c x-coordinate of acceleration
if (d==1)
{
h = h + (eta[i]*txx[i] + eta[j]*txx[j])*dwdx[1][k];
if (dim>=2)
{
h = h + (eta[i]*txy[i] + eta[j]*txy[j])*dwdx[2][k];
if (dim==3)
{
h = h + (eta[i]*txz[i] + eta[j]*txz[j])*dwdx[3][k];
}
}
}
//c y-coordinate of acceleration
else if (d==2)
{
h = h + (eta[i]*txy[i] + eta[j]*txy[j])*dwdx[1][k]+ (eta[i]*tyy[i] + eta[j]*tyy[j])*dwdx[2][k];
if (dim==3)
{
h = h + (eta[i]*tyz[i] + eta[j]*tyz[j])*dwdx[3][k];
}
}
//c z-coordinate of acceleration
else if (d==3)
{
h = h + (eta[i]*txz[i] + eta[j]*txz[j])*dwdx[1][k] + (eta[i]*tyz[i] + eta[j]*tyz[j])*dwdx[2][k]+(eta[i]*tzz[i] + eta[j]*tzz[j])*dwdx[3][k];
}
}
h = h*rhoij;
dvx1dt[d][i] = dvx1dt[d][i] + mass[j]*h;
dvx1dt[d][j] = dvx1dt[d][j] - mass[i]*h;
}
he = he*rhoij;
de1dt[i] = de1dt[i] + mass[j]*he;
de1dt[j] = de1dt[j] + mass[i]*he;
}
else if (pa_sph==2)//c For SPH algorithm 2
{
for(int d=1;d<=dim;d++)
{
h = -(p[i]/pow(rho[i],2) + p[j]/pow(rho[j],2))*dwdx[d][k];
he = he + (vx[d][j] - vx[d][i])*h;
//c Viscous force
if (visc)
{
//c x-coordinate of acceleration
if (d==1)
{
h = h + (eta[i]*txx[i]/pow(rho[i],2) +eta[j]*txx[j]/pow(rho[j],2))*dwdx[1][k];
if (dim>=2)
{
h = h + (eta[i]*txy[i]/pow(rho[i],2) +eta[j]*txy[j]/pow(rho[j],2))*dwdx[2][k];
if (dim==3)
{
h = h + (eta[i]*txz[i]/pow(rho[i],2) +eta[j]*txz[j]/pow(rho[j],2))*dwdx[3][k];
}
}
}
//c y-coordinate of acceleration
else if (d==2)
{
h = h + (eta[i]*txy[i]/pow(rho[i],2) + eta[j]*txy[j]/pow(rho[j],2))*dwdx[1][k]+ (eta[i]*tyy[i]/pow(rho[i],2)+eta[j]*tyy[j]/pow(rho[j],2))*dwdx[2][k];
if (dim==3)
{
h = h + (eta[i]*tyz[i]/pow(rho[i],2)+ eta[j]*tyz[j]/pow(rho[j],2))*dwdx[3][k];
}
}
//c z-coordinate of acceleration
else if (d==3)
{
h = h + (eta[i]*txz[i]/pow(rho[i],2) +eta[j]*txz[j]/pow(rho[j],2))*dwdx[1][k]+ (eta[i]*tyz[i]/pow(rho[i],2) +eta[j]*tyz[j]/pow(rho[j],2))*dwdx[2][k] + (eta[i]*tzz[i]/pow(rho[i],2) + eta[j]*tzz[j]/pow(rho[j],2))*dwdx[3][k];
}
}
dvx1dt[d][i] = dvx1dt[d][i] + mass[j]*h;dvx1dt[d][j] = dvx1dt[d][j] - mass[i]*h;
}
de1dt[i] = de1dt[i] + mass[j]*he;de1dt[j] = de1dt[j] + mass[i]*he;
}
}
for(i=1;i<=ntotal;i++)
de1dt[i] = tdsdt[i] + 0.5*de1dt[i]; // Change of specific internal energy de/dt = T ds/dt - p/rho vc,c:
//------------------------------------------------------------------------------------------------------------------------
//for( i=1;i<=ntotal;i++)
//{
// cout<<"int_force "<<"t[i]="<<t[i]<<"c[i]="<<c[i]<<"p[i]="<<p[i]<<"dvx1dt[1][i]="<<dvx1dt[1][i]<<"tdsdt[i]="<<tdsdt[i]<<"de1dt[i]="<<de1dt[i]<<endl;
//}
delete []dvx1;delete []txx;delete []tyy;delete []tzz;delete []txy;delete []txz;delete []tyz;delete []vcc;
}
void SPH::kernel(double r,double *dx,double hsml1,double &w1,double *(&dw1dx))
{//注意核函数的变量都 不是全局变量
// double dw;//未使用
double q, factor;
q = r/hsml1;
w1 = 0.0;
for(int d=1;d<=dim;d++)
dw1dx[d] = 0.0;
if (skf==1)
{
if (dim==1)
{
factor = 1.0/hsml1;
}
else if(dim==2)
{
factor = 15.0/(7.0*pi*hsml1*hsml1);
}
else if (dim==3)
{
factor = 3.0/(2.0*pi*hsml1*hsml1*hsml1);
}
else
{
cout<<"kernel=err"<<endl;
}
if ((q>=0)&&(q<=1.0))
{
w1 = factor * (2.0/3.0 - q*q + pow(q,3) / 2.0);
for(int d = 1;d<=dim;d++)
dw1dx[d] = factor * (-2.0+3.0/2.0*q)/pow(hsml1,2) * dx[d];
}
else if ((q>1.0)&&(q<=2.0))
{
w1 = factor * 1.0/6.0 * pow((2.0-q),3) ;
for(int d = 1; d<=dim;d++)
dw1dx[d] =-factor * 1.0/6.0 * 3.0*pow((2.0-q),2)/hsml1 * (dx[d]/r);
}
else
{
w1=0.0;
for(int d= 1;d<= dim;d++)
dw1dx[d] = 0.0;
}
}
else if (skf==2.0)
{
factor = 1.0 / (pow(hsml1,dim) * pow(pi,(dim/2.0)));
if((q>=0.0)&&(q<=3.0))
{
w1 = factor * exp(-q*q);
for(int d = 1;d<= dim;d++)
dw1dx[d] = w1 * ( -2.0* dx[d]/hsml1/hsml1);
}
else
{ w1 = 0.0;
for(int d = 1;d<= dim;d++)
dw1dx[d] = 0.0;
}
}
else if (skf==3.0)
{
if (dim==1.0)
{
factor = 1.0 / (120.0*hsml1);
}
else if (dim==2.0)
{
factor = 7.0 / (478.0*pi*hsml1*hsml1);
}
else if (dim==3.0)
{
factor = 1.0 / (120.0*pi*hsml1*hsml1*hsml1);
}
else
{
cout<<"kernel=Wrong dimension"<<endl;
}
if((q>=0.0)&&(q<=1.0))
{
w1 = factor * (pow((3.0-q),5) - 6.0*pow((2.0-q),5) + 15.0*pow((1.0-q),5));
for(int d= 1;d<=dim;d++)
dw1dx[d] = factor * ((-120.0 + 120.0*q - 50.0*pow(q,2))/ pow(hsml1,2) * dx[d] );
}
else if((q>1.0)&&(q<=2.0))
{
w1 = factor * ( pow((3.0-q),5) - 6.0*pow((2.0-q),5));
for(int d= 1;d<= dim;d++)
dw1dx[d] = factor * (-5.0*pow((3.0-q),4) + 30.0*pow((2.0-q),4))/ hsml1 * (dx[d]/r);
}
else if((q>2.0)&&(q<=3.0))
{
w1 = factor * pow((3.0-q),5);
for(int d= 1; d<=dim;d++)
dw1dx[d] = factor * (-5.0*pow((3.0-q),4)) / hsml1 * (dx[d]/r);
}
else
{
w1 = 0.0;
for(int d = 1;d<= dim;d++)
dw1dx[d] = 0.0;
}
}
}
void SPH::direct_find(int ntotal)
{
//这里ntotal应包含虚粒子,使用的是局部变量
// cout<<itimestep<<" "<<ntotal<<hsml[300]<<x[1][300]<<endl;
int sumiac, maxiac, miniac, noiac, maxp, minp, scale_k ;
double *dxiac, driac, r, mhsml, *tdwdx; //r有全局变量
dxiac=new double[dim+1];tdwdx=new double[dim+1];
double temp=0;
if (skf==1.0)
{
scale_k = 2;
}
else if (skf==2.0)
{
scale_k = 3;
}
else if (skf==3.0)
{
scale_k = 3;
}
for(int i=1;i<=ntotal;i++)
countiac[i] = 0;
niac = 0;
for( i=1;i<=ntotal-1;i++)
{
for(int j = i+1;j<= ntotal;j++)
{
dxiac[1] = x[1][i] - x[1][j];
driac = dxiac[1]*dxiac[1];
for(int d=2;d<=dim;d++)
{
dxiac[d] = x[d][i] - x[d][j];
driac = driac + dxiac[d]*dxiac[d];
}
mhsml = (hsml[i]+hsml[j])/2.0;
if (sqrt(driac)<(scale_k*mhsml))
{
if (niac<max_interaction)
{
//c Neighboring pair list, and totalinteraction number and
//c the interaction number for each particle
niac = niac + 1;
pair_i[niac] = i;
pair_j[niac] = j;
r = sqrt(driac);
countiac[i] = countiac[i] + 1;
countiac[j] = countiac[j] + 1;
//c Kernel and derivations of kernel
//temp=w[niac];//temp是空值???????
kernel(r,dxiac,mhsml,temp,tdwdx);//???????????,可能是程序出问题的地方
w[niac]=temp;//此结果仅在计算校正平均速度时用到，av_vel子程序
for(int d=1;d<=dim;d++)
dwdx[d][niac] = tdwdx[d];
// cout<<"dwdx[1]["<<niac<<"]"<<dwdx[1][niac]<<endl;
}
else
{
cout<<"direct_find=too many interactions"<<endl;
}
}
}
}
//c Statistics for the interaction
sumiac = 0;
maxiac = 0;
miniac = 1000;
noiac = 0;
for( i=1;i<=ntotal;i++)
{
sumiac = sumiac + countiac[i];
if (countiac[i]>maxiac)
{
maxiac = countiac[i];
maxp = i;
}
if (countiac[i]<miniac)
{
miniac = countiac[i];
minp = i;
}
if (countiac[i]==0)
{
noiac = noiac + 1;
}
}
//// if ((itimestep%print_step)==0)
//// {
//// if (int_stat)
//// {
//// cout<<"Statistics: interactions per particle:"<<endl;
//// cout<<"**** Particle:',maxp,"<<maxp<<" maximal interactions:"<<maxiac<<endl;
//// cout<<"**** Particle:',minp,"<<minp<<" minimal interactions:"<<miniac<<endl;
//// cout<<"**** Average :,"<<double(sumiac)/double(ntotal)<<endl;
//// cout<<"**** Total pairs :"<<niac<<endl;
//// cout<<"**** Particles with no interactions:"<<noiac<<endl;
//// }
//// }
//cout<<"direct_find1 "<<"niac="<<niac<<"pair_i="<<pair_i[niac]<<"pair_j="<<pair_j[niac]<<endl;
//cout<<"direct_find2 "<<"w[niac]="<<w[niac]<<"dwdx[1][niac]="<<dwdx[1][niac]<<"countiac[ntotal]="<<countiac[ntotal]<<endl;
delete []dxiac;delete []tdwdx;
}
void SPH::single_step()
{//single_step( u, s,rho,p, t, tdsdt, dx1, dvx, du, ds, drho,itype, av);
//cout<<u[300]<<" "<<s[300]<<" "<<rho[300]<<" "<<p[300]<<" "<<t[300]<<" "<<tdsdt[300]<<endl;
//cout<<dx1[300]<<" "<<dvx[300]<<" "<<du[300]<<" "<<ds[300]<<" "<<drho[300]<<" "<<av[300]<<endl;
//请核实这个子程序是否使用了*dx或**dx,c是表示温度还是声速？？这是程序一个明显的问题
//int nvirt, niac; //已定义全局变量
//int *pair_i,*pair_j;
int *ns;
// double *w,**dwdx
double **indvxdt,**exdvxdt;
double **ardvxdt,*avdudt,*ahdudt;
// double *c,*eta;
ns=new int[maxn+1];avdudt=new double[maxn+1];ahdudt=new double[maxn+1];
indvxdt=new double *[dim+1];for(int i=0;i<dim+1;i++)indvxdt[i]=new double[maxn+1];
exdvxdt=new double *[dim+1];for(i=0;i<dim+1;i++)exdvxdt[i]=new double[maxn+1];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -