在进行数值计算编程的过程中往往需要用到浮点数的计算,但浮点数的加减运算通常是会出现误差的。具体出现问题代码如下:
double begin = 0.0; //起始位置
double end = 20.0; //结束位置
int k = 2; //核函数支持半径的倍数
double initialDis = 0.1; //粒子间的初始距离(需要与粒子的控制体积保存一致性)
/*****************第一步初始化粒子**********/
for(double i=begin-(2*k-1)*initialDis/2; i<=end+(2*k)*initialDis/2; i=i+initialDis)
{
par = new LNode;//申请内存空间
(*par).data.changeCoordinate(i,0.0);//更新坐标
(*par).data.changeMass(mass); //更新质量
if(i<begin)
{
(*par).data.changeType(itghost); //更新粒子类型为itghost
(*par).data.changePressure(Pr);//更新压力(初始值)
(*par).data.changeVelocity(V0,0);//更新速度(初始值)
}
else if(i>=begin&&i<begin+initialDis)//inlet粒子
{
// cout<<"itoutlet<><><><><><>"<<i<<endl;
(*par).data.changeType(itinlet); //更新粒子类型为inlet
(*par).data.changePressure(Pr-i*Pf);//更新压力(初始值)
(*par).data.changeVelocity(V0,0);//更新速度(初始值)
}
else if(i<=end&&i>end-initialDis)//outlet粒子
{
(*par).data.changeType(itoutlet); //更新粒子类型为outlet
(*par).data.changePressure(Pr-i*Pf);//更新压力(初始值)
(*par).data.changeVelocity(V0,0);//更新速度(初始值)
}
else if(i>end)//outmirror粒子
{
(*par).data.changeType(itmirror); //更新粒子类型为itmirror
(*par).data.changePressure(Pr-(2*end-(40-i))*Pf);//更新压力
(*par).data.changeVelocity(0,0);//更新速度(初始值)
}
else
{
(*par).data.changeType(itfluid); //更新粒子类型为fluid
(*par).data.changePressure(Pr-i*Pf);//更新压力(初始值)
(*par).data.changeVelocity(V0,0);//更新速度(初始值)
}
(*par).next = (*totalPar).next;//头插法插入链表
(*totalPar).next = par;
}
如上述代码所示,在计算粒子的位置时需要进行浮点数的计算,但经过若干次加法计算,结果并不能和预期相同,inlet粒子的位置应为0.0,实际计算中,是通过0.1相加出来的,并不等于0.0,所以判断粒子为inlet粒子的条件根据具体的情况修改为(i>=begin&&i<begin+initialDis)。