primitive equations on square mesh is diverging (computational fluid dynamics)
up vote
0
down vote
favorite
I am trying to make a basic program for simulating the primitive equations, but something has gone wrong as it rapidly diverges.
I used this paper as a guide:
http://www.met.reading.ac.uk/~ross/Science/PrimEqs.html It is really good but does little to help me translate it into working code. I get the concepts(Divergence,derivatives,Material derivative(sort of),Pressure gradiant, and so on...), but it is wierd as my code keeps messing up all the time.
Both viscocity and diffusion are 0 so dont worry about them.(I have tested them separatly and both are stable)
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*DT;
Density[i][j]-= Velocity[i][j].x*DT;
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT;
Density[i][j]-= Velocity[i][j].y*DT;
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
I am pretty sure that there are no problem with the sub-functions but here they are just in case.
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
Can anyone see the problem with it?
Minimum code needed for failiure:
void update()
{
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
if(i<amountX-1)
{
Density[i+1][j]+= Velocity[i][j].x*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].x*DT*Density[i][j];
}
if(j<amountY-1)
{
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j];
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
Simple 1D psudocode:
for(int i =0;i<N-1;i++)
Gradient[i]=(Pressure[i]-Pressure[i+1])*dx;
for(int i =0;i<N;i++)
{
Accel = -Gradient[i]/Density[i];
Velocity[i]+=Accel*dt;
if(i<N-1)
{
Density[i+1]+=Velocity*dt;
Density[i]-=Velocity*dt;
}
Pressure[i]=Density[i];
}
Here is the whole thing for testing. Contains temperature,viscocity and diffusion. They can be disabled if needed. Requires Processing https://processing.org/
int amountX=40;
int amountY=20;
float Temp=new float[amountX][amountY];
float Density=new float[amountX][amountY];
float Pressure=new float[amountX][amountY];
PVector Velocity=new PVector[amountX][amountY];
float Grav=0.05;
float GasConst=1;
float HeatCapacity=5;
float Diffusion=0.5;
float Viscosity=0.4;
float DX=1;
float DY=1;
float DT=0.001;
int Warp = 1;
float VDiv;
float TempChange;
PVector PGrad;
void setup()
{
size(1000,500);
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Temp[i][j]=3;
Pressure[i][j]=exp(j*Grav/(GasConst*Temp[i][j]))*40;
Density[i][j]=Pressure[i][j]/(GasConst*Temp[i][j]);
}
}
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Velocity[i][j]=new PVector(0,0);
}
}
int Min=10;
int Max=15;
for(int i=Min;i<Max;i++)
{
for(int j=Min;j<Max;j++)
{
float x=(i-Min)/(float)(Max-Min);
float y=(j-Min)/(float)(Max-Min);
Temp[i][j]=lerp(3,6,exp(-(sq(2*x-1)+sq(2*y-1))*3));
}
}
for(int i =0;i<0;i++)
{
update();
}
}
void draw()
{
background(201);
for(int i = 0;i<Warp;i++)
{
update();
}
show();
}
void show()
{
float dx=width/(float)amountX;
float dy=height/(float)amountY;
float dv=10*dx;
noStroke();
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
fill(Temp[i][j]*10,0,Density[i][j]*10);
rect(dx*i,dy*j,dx,dy);
}
}
stroke(255);
for(int i=0;i<amountX;i+=1)
{
for(int j=0;j<amountY;j+=1)
{
line(dx*(i+0.5),dy*(j+0.5),dx*(i+0.5)+Velocity[i][j].x*dv,dy*(j+0.5)+Velocity[i][j].y*dv);
}
}
int MouseI=(int)(mouseX/dx);
int MouseJ=(int)(mouseY/dy);
fill(255);
if(MouseI>0&&MouseJ>0&&MouseI<amountX-1&&MouseJ<amountY-1)
{
text("Temp:"+Temp[MouseI][MouseJ],10,20);
text("Pres:"+Pressure[MouseI][MouseJ],10,40);
text("Dens:"+Density[MouseI][MouseJ],10,60);
}
}
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*Density[i][j];
Density[i][j]-= Velocity[i][j].x*Density[i][j];
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
java physics
asked Nov 17 at 21:59
Eriksonn
42
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
|
show 8 more comments
up vote
0
down vote
favorite
I am trying to make a basic program for simulating the primitive equations, but something has gone wrong as it rapidly diverges.
I used this paper as a guide:
http://www.met.reading.ac.uk/~ross/Science/PrimEqs.html It is really good but does little to help me translate it into working code. I get the concepts(Divergence,derivatives,Material derivative(sort of),Pressure gradiant, and so on...), but it is wierd as my code keeps messing up all the time.
Both viscocity and diffusion are 0 so dont worry about them.(I have tested them separatly and both are stable)
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*DT;
Density[i][j]-= Velocity[i][j].x*DT;
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT;
Density[i][j]-= Velocity[i][j].y*DT;
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
I am pretty sure that there are no problem with the sub-functions but here they are just in case.
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
Can anyone see the problem with it?
Minimum code needed for failiure:
void update()
{
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
if(i<amountX-1)
{
Density[i+1][j]+= Velocity[i][j].x*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].x*DT*Density[i][j];
}
if(j<amountY-1)
{
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j];
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
Simple 1D psudocode:
for(int i =0;i<N-1;i++)
Gradient[i]=(Pressure[i]-Pressure[i+1])*dx;
for(int i =0;i<N;i++)
{
Accel = -Gradient[i]/Density[i];
Velocity[i]+=Accel*dt;
if(i<N-1)
{
Density[i+1]+=Velocity*dt;
Density[i]-=Velocity*dt;
}
Pressure[i]=Density[i];
}
Here is the whole thing for testing. Contains temperature,viscocity and diffusion. They can be disabled if needed. Requires Processing https://processing.org/
int amountX=40;
int amountY=20;
float Temp=new float[amountX][amountY];
float Density=new float[amountX][amountY];
float Pressure=new float[amountX][amountY];
PVector Velocity=new PVector[amountX][amountY];
float Grav=0.05;
float GasConst=1;
float HeatCapacity=5;
float Diffusion=0.5;
float Viscosity=0.4;
float DX=1;
float DY=1;
float DT=0.001;
int Warp = 1;
float VDiv;
float TempChange;
PVector PGrad;
void setup()
{
size(1000,500);
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Temp[i][j]=3;
Pressure[i][j]=exp(j*Grav/(GasConst*Temp[i][j]))*40;
Density[i][j]=Pressure[i][j]/(GasConst*Temp[i][j]);
}
}
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Velocity[i][j]=new PVector(0,0);
}
}
int Min=10;
int Max=15;
for(int i=Min;i<Max;i++)
{
for(int j=Min;j<Max;j++)
{
float x=(i-Min)/(float)(Max-Min);
float y=(j-Min)/(float)(Max-Min);
Temp[i][j]=lerp(3,6,exp(-(sq(2*x-1)+sq(2*y-1))*3));
}
}
for(int i =0;i<0;i++)
{
update();
}
}
void draw()
{
background(201);
for(int i = 0;i<Warp;i++)
{
update();
}
show();
}
void show()
{
float dx=width/(float)amountX;
float dy=height/(float)amountY;
float dv=10*dx;
noStroke();
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
fill(Temp[i][j]*10,0,Density[i][j]*10);
rect(dx*i,dy*j,dx,dy);
}
}
stroke(255);
for(int i=0;i<amountX;i+=1)
{
for(int j=0;j<amountY;j+=1)
{
line(dx*(i+0.5),dy*(j+0.5),dx*(i+0.5)+Velocity[i][j].x*dv,dy*(j+0.5)+Velocity[i][j].y*dv);
}
}
int MouseI=(int)(mouseX/dx);
int MouseJ=(int)(mouseY/dy);
fill(255);
if(MouseI>0&&MouseJ>0&&MouseI<amountX-1&&MouseJ<amountY-1)
{
text("Temp:"+Temp[MouseI][MouseJ],10,20);
text("Pres:"+Pressure[MouseI][MouseJ],10,40);
text("Dens:"+Density[MouseI][MouseJ],10,60);
}
}
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*Density[i][j];
Density[i][j]-= Velocity[i][j].x*Density[i][j];
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
java physics
asked Nov 17 at 21:59
Eriksonn
42
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
Your code isn't behaving as expected, and so the first thing that you must do (and do before you come here) is to isolate the source of the bug -- meaning that you need to do intense debugging. Consider stepping through the program using your IDE's debugger, inspecting the state of key variables as you do so. This will give you valuable insights and quite possibly a decent solution. Even if you don't get a full answer, you'll be closer and can ask a better question. If you still need our help, you will need to create and post a Minimal, Complete, and Verifiable example, and tell the details of expected and observed behaviors
– Hovercraft Full Of Eels
Nov 17 at 22:02
@HovercraftFullOfEels To be fair: For many sorts of simulations, the usual debugging approaches are essentially useless. There often is no point in stepping through some simulation time step that involves complex math where you simply cannot say whether a certain value is "right" or "wrong" at a certain point (except forNaNor so). Even more important is a small MCVE, and often the good ol'printlns, where you can at least print messages like"Now I'm adding these 1000 elements [...] to these [...]"to have the chance to spot where one vector goes off the rails.
– Marco13
Nov 17 at 22:12
Can you be more specific than "my code keeps messing up all the time."? What is the behavior you expect vs the behavior being observed?
– John Camerin
Nov 17 at 22:13
@Marco13: I have to bow to your expertise since I don't do simulation programming, but something sure needs to be done to first isolate the error, somehow
– Hovercraft Full Of Eels
Nov 17 at 22:23
The problem manifests as spatial oscillations growing over time. What i mean about that is that almost immideatly there becomes a "checkerboard pattern" of high/low pressure that grows and after a few seconds turn Everything into nan values. By making the temperature constant it works a bit longer and produce a shockwave pattern that looks very promesing but after a while also fails. It feels like something is overshooting in order to produce these oscillations.
– Eriksonn
Nov 17 at 22:48
|
show 8 more comments
up vote
0
down vote
favorite
up vote
0
down vote
favorite
I am trying to make a basic program for simulating the primitive equations, but something has gone wrong as it rapidly diverges.
I used this paper as a guide:
http://www.met.reading.ac.uk/~ross/Science/PrimEqs.html It is really good but does little to help me translate it into working code. I get the concepts(Divergence,derivatives,Material derivative(sort of),Pressure gradiant, and so on...), but it is wierd as my code keeps messing up all the time.
Both viscocity and diffusion are 0 so dont worry about them.(I have tested them separatly and both are stable)
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*DT;
Density[i][j]-= Velocity[i][j].x*DT;
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT;
Density[i][j]-= Velocity[i][j].y*DT;
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
I am pretty sure that there are no problem with the sub-functions but here they are just in case.
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
Can anyone see the problem with it?
Minimum code needed for failiure:
void update()
{
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
if(i<amountX-1)
{
Density[i+1][j]+= Velocity[i][j].x*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].x*DT*Density[i][j];
}
if(j<amountY-1)
{
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j];
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
Simple 1D psudocode:
for(int i =0;i<N-1;i++)
Gradient[i]=(Pressure[i]-Pressure[i+1])*dx;
for(int i =0;i<N;i++)
{
Accel = -Gradient[i]/Density[i];
Velocity[i]+=Accel*dt;
if(i<N-1)
{
Density[i+1]+=Velocity*dt;
Density[i]-=Velocity*dt;
}
Pressure[i]=Density[i];
}
Here is the whole thing for testing. Contains temperature,viscocity and diffusion. They can be disabled if needed. Requires Processing https://processing.org/
int amountX=40;
int amountY=20;
float Temp=new float[amountX][amountY];
float Density=new float[amountX][amountY];
float Pressure=new float[amountX][amountY];
PVector Velocity=new PVector[amountX][amountY];
float Grav=0.05;
float GasConst=1;
float HeatCapacity=5;
float Diffusion=0.5;
float Viscosity=0.4;
float DX=1;
float DY=1;
float DT=0.001;
int Warp = 1;
float VDiv;
float TempChange;
PVector PGrad;
void setup()
{
size(1000,500);
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Temp[i][j]=3;
Pressure[i][j]=exp(j*Grav/(GasConst*Temp[i][j]))*40;
Density[i][j]=Pressure[i][j]/(GasConst*Temp[i][j]);
}
}
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Velocity[i][j]=new PVector(0,0);
}
}
int Min=10;
int Max=15;
for(int i=Min;i<Max;i++)
{
for(int j=Min;j<Max;j++)
{
float x=(i-Min)/(float)(Max-Min);
float y=(j-Min)/(float)(Max-Min);
Temp[i][j]=lerp(3,6,exp(-(sq(2*x-1)+sq(2*y-1))*3));
}
}
for(int i =0;i<0;i++)
{
update();
}
}
void draw()
{
background(201);
for(int i = 0;i<Warp;i++)
{
update();
}
show();
}
void show()
{
float dx=width/(float)amountX;
float dy=height/(float)amountY;
float dv=10*dx;
noStroke();
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
fill(Temp[i][j]*10,0,Density[i][j]*10);
rect(dx*i,dy*j,dx,dy);
}
}
stroke(255);
for(int i=0;i<amountX;i+=1)
{
for(int j=0;j<amountY;j+=1)
{
line(dx*(i+0.5),dy*(j+0.5),dx*(i+0.5)+Velocity[i][j].x*dv,dy*(j+0.5)+Velocity[i][j].y*dv);
}
}
int MouseI=(int)(mouseX/dx);
int MouseJ=(int)(mouseY/dy);
fill(255);
if(MouseI>0&&MouseJ>0&&MouseI<amountX-1&&MouseJ<amountY-1)
{
text("Temp:"+Temp[MouseI][MouseJ],10,20);
text("Pres:"+Pressure[MouseI][MouseJ],10,40);
text("Dens:"+Density[MouseI][MouseJ],10,60);
}
}
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*Density[i][j];
Density[i][j]-= Velocity[i][j].x*Density[i][j];
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
java physics
asked Nov 17 at 21:59
Eriksonn
42
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
I am trying to make a basic program for simulating the primitive equations, but something has gone wrong as it rapidly diverges.
I used this paper as a guide:
http://www.met.reading.ac.uk/~ross/Science/PrimEqs.html It is really good but does little to help me translate it into working code. I get the concepts(Divergence,derivatives,Material derivative(sort of),Pressure gradiant, and so on...), but it is wierd as my code keeps messing up all the time.
Both viscocity and diffusion are 0 so dont worry about them.(I have tested them separatly and both are stable)
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*DT;
Density[i][j]-= Velocity[i][j].x*DT;
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT;
Density[i][j]-= Velocity[i][j].y*DT;
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
I am pretty sure that there are no problem with the sub-functions but here they are just in case.
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
Can anyone see the problem with it?
Minimum code needed for failiure:
void update()
{
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
if(i<amountX-1)
{
Density[i+1][j]+= Velocity[i][j].x*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].x*DT*Density[i][j];
}
if(j<amountY-1)
{
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j];
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
Simple 1D psudocode:
for(int i =0;i<N-1;i++)
Gradient[i]=(Pressure[i]-Pressure[i+1])*dx;
for(int i =0;i<N;i++)
{
Accel = -Gradient[i]/Density[i];
Velocity[i]+=Accel*dt;
if(i<N-1)
{
Density[i+1]+=Velocity*dt;
Density[i]-=Velocity*dt;
}
Pressure[i]=Density[i];
}
Here is the whole thing for testing. Contains temperature,viscocity and diffusion. They can be disabled if needed. Requires Processing https://processing.org/
int amountX=40;
int amountY=20;
float Temp=new float[amountX][amountY];
float Density=new float[amountX][amountY];
float Pressure=new float[amountX][amountY];
PVector Velocity=new PVector[amountX][amountY];
float Grav=0.05;
float GasConst=1;
float HeatCapacity=5;
float Diffusion=0.5;
float Viscosity=0.4;
float DX=1;
float DY=1;
float DT=0.001;
int Warp = 1;
float VDiv;
float TempChange;
PVector PGrad;
void setup()
{
size(1000,500);
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Temp[i][j]=3;
Pressure[i][j]=exp(j*Grav/(GasConst*Temp[i][j]))*40;
Density[i][j]=Pressure[i][j]/(GasConst*Temp[i][j]);
}
}
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Velocity[i][j]=new PVector(0,0);
}
}
int Min=10;
int Max=15;
for(int i=Min;i<Max;i++)
{
for(int j=Min;j<Max;j++)
{
float x=(i-Min)/(float)(Max-Min);
float y=(j-Min)/(float)(Max-Min);
Temp[i][j]=lerp(3,6,exp(-(sq(2*x-1)+sq(2*y-1))*3));
}
}
for(int i =0;i<0;i++)
{
update();
}
}
void draw()
{
background(201);
for(int i = 0;i<Warp;i++)
{
update();
}
show();
}
void show()
{
float dx=width/(float)amountX;
float dy=height/(float)amountY;
float dv=10*dx;
noStroke();
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
fill(Temp[i][j]*10,0,Density[i][j]*10);
rect(dx*i,dy*j,dx,dy);
}
}
stroke(255);
for(int i=0;i<amountX;i+=1)
{
for(int j=0;j<amountY;j+=1)
{
line(dx*(i+0.5),dy*(j+0.5),dx*(i+0.5)+Velocity[i][j].x*dv,dy*(j+0.5)+Velocity[i][j].y*dv);
}
}
int MouseI=(int)(mouseX/dx);
int MouseJ=(int)(mouseY/dy);
fill(255);
if(MouseI>0&&MouseJ>0&&MouseI<amountX-1&&MouseJ<amountY-1)
{
text("Temp:"+Temp[MouseI][MouseJ],10,20);
text("Pres:"+Pressure[MouseI][MouseJ],10,40);
text("Dens:"+Density[MouseI][MouseJ],10,60);
}
}
void update()
{
VDiv = Divergence(Velocity);
PVector ViscosityVector = VectorLaplace(Velocity);
TempChange = Laplace(Temp);//diffusion operator
PGrad = Gradient(Pressure);//pressure force
for(int i=amountX-1;i>=0;i--)
{
for(int j=amountY-1;j>=0;j--)
{
TempChange[i][j]*=Diffusion;
TempChange[i][j]-=VDiv[i][j]*Pressure[i][j];// temp changes as volume changes
PVector Accel=new PVector(0,Grav).sub(PVector.div(PGrad[i][j],Density[i][j]));
Accel.add(PVector.mult(ViscosityVector[i][j],Viscosity));//viscocity
float TS=(Density[i][j]*HeatCapacity);
//the material derivative
Velocity[i][j].add(Accel.mult(DT));
Temp[i][j]+=DT*TempChange[i][j]/TS;
if(i<amountX-1)
{
Temp[i+1][1]+= Velocity[i][j].x*DT;
Temp[i][j]-= Velocity[i][j].x*DT;
Density[i+1][j]+= Velocity[i][j].x*Density[i][j];
Density[i][j]-= Velocity[i][j].x*Density[i][j];
}
if(j<amountY-1)
{
Temp[i][j+1]+= Velocity[i][j].y*DT;
Temp[i][j]-= Velocity[i][j].y*DT;
Density[i][j+1]+= Velocity[i][j].y*DT*Density[i][j];
Density[i][j]-= Velocity[i][j].y*DT*Density[i][j];
}
// ideal gas law
Pressure[i][j]=Density[i][j]*(GasConst*Temp[i][j]);
}
}
}
PVector Gradient(float Input)
{
PVector Out = new PVector[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
Out[i][j]=new PVector();
if(j<amountY-1)
Out[i][j].y=(-Input[i][j]+Input[i][j+1])*DY;
if(i<amountX-1)
Out[i][j].x=(-Input[i][j]+Input[i+1][j])*DX;
}
}
return Out;
}
PVector VectorLaplace(PVector Input)
{
float InputX=new float[amountX][amountY];
float InputY=new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
InputX[i][j]=Input[i][j].x;
InputY[i][j]=Input[i][j].y;
}
}
InputX = Laplace(InputX);
InputY = Laplace(InputY);
PVector Out = new PVector[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
Out[i][j]=new PVector(InputX[i][j],InputY[i][j]);
}
}
return Out;
}
float Laplace(float Input)
{
float Flow = new float[amountX][amountY];
for(int i =0;i<amountX;i++)
{
for(int j =0;j<amountY;j++)
{
if(i>0)
{
Flow[i][j]-=Input[i][j];
Flow[i-1][j]+=Input[i][j];
}
if(j>0)
{
Flow[i][j]-=Input[i][j];
Flow[i][j-1]+=Input[i][j];
}
if(i<amountX-1)
{
Flow[i][j]-=Input[i][j];
Flow[i+1][j]+=Input[i][j];
}
if(j<amountY-1)
{
Flow[i][j]-=Input[i][j];
Flow[i][j+1]+=Input[i][j];
}
}
}
return Flow;
}
float Divergence(PVector Input)
{
float Out = new float[amountX][amountY];
for(int i=0;i<amountX;i++)
{
for(int j=0;j<amountY;j++)
{
if(i>0)
{
Out[i][j]-=(Input[i-1][j].x)*DX;
}
if(j>0)
{
Out[i][j]-=(Input[i][j-1].y)*DY;
}
if(i<amountX-1)
{
Out[i][j]+=(Input[i][j].x)*DX;
}
if(j<amountY-1)
{
Out[i][j]+=(Input[i][j].y)*DY;
}
}
}
return Out;
}
java physics
java physics
asked Nov 17 at 21:59
Eriksonn
42
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
asked Nov 17 at 21:59
Eriksonn
42
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
edited Nov 18 at 0:42
asked Nov 17 at 21:59
Eriksonn
42
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
asked Nov 17 at 21:59
Eriksonn
42
asked Nov 17 at 21:59
Eriksonn
42
42
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
New contributor
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
Eriksonn is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
Your code isn't behaving as expected, and so the first thing that you must do (and do before you come here) is to isolate the source of the bug -- meaning that you need to do intense debugging. Consider stepping through the program using your IDE's debugger, inspecting the state of key variables as you do so. This will give you valuable insights and quite possibly a decent solution. Even if you don't get a full answer, you'll be closer and can ask a better question. If you still need our help, you will need to create and post a Minimal, Complete, and Verifiable example, and tell the details of expected and observed behaviors
– Hovercraft Full Of Eels
Nov 17 at 22:02
@HovercraftFullOfEels To be fair: For many sorts of simulations, the usual debugging approaches are essentially useless. There often is no point in stepping through some simulation time step that involves complex math where you simply cannot say whether a certain value is "right" or "wrong" at a certain point (except forNaNor so). Even more important is a small MCVE, and often the good ol'printlns, where you can at least print messages like"Now I'm adding these 1000 elements [...] to these [...]"to have the chance to spot where one vector goes off the rails.
– Marco13
Nov 17 at 22:12
Can you be more specific than "my code keeps messing up all the time."? What is the behavior you expect vs the behavior being observed?
– John Camerin
Nov 17 at 22:13
@Marco13: I have to bow to your expertise since I don't do simulation programming, but something sure needs to be done to first isolate the error, somehow
– Hovercraft Full Of Eels
Nov 17 at 22:23
The problem manifests as spatial oscillations growing over time. What i mean about that is that almost immideatly there becomes a "checkerboard pattern" of high/low pressure that grows and after a few seconds turn Everything into nan values. By making the temperature constant it works a bit longer and produce a shockwave pattern that looks very promesing but after a while also fails. It feels like something is overshooting in order to produce these oscillations.
– Eriksonn
Nov 17 at 22:48
|
show 8 more comments
Your code isn't behaving as expected, and so the first thing that you must do (and do before you come here) is to isolate the source of the bug -- meaning that you need to do intense debugging. Consider stepping through the program using your IDE's debugger, inspecting the state of key variables as you do so. This will give you valuable insights and quite possibly a decent solution. Even if you don't get a full answer, you'll be closer and can ask a better question. If you still need our help, you will need to create and post a Minimal, Complete, and Verifiable example, and tell the details of expected and observed behaviors
– Hovercraft Full Of Eels
Nov 17 at 22:02
@HovercraftFullOfEels To be fair: For many sorts of simulations, the usual debugging approaches are essentially useless. There often is no point in stepping through some simulation time step that involves complex math where you simply cannot say whether a certain value is "right" or "wrong" at a certain point (except forNaNor so). Even more important is a small MCVE, and often the good ol'printlns, where you can at least print messages like"Now I'm adding these 1000 elements [...] to these [...]"to have the chance to spot where one vector goes off the rails.
– Marco13
Nov 17 at 22:12
Can you be more specific than "my code keeps messing up all the time."? What is the behavior you expect vs the behavior being observed?
– John Camerin
Nov 17 at 22:13
@Marco13: I have to bow to your expertise since I don't do simulation programming, but something sure needs to be done to first isolate the error, somehow
– Hovercraft Full Of Eels
Nov 17 at 22:23
The problem manifests as spatial oscillations growing over time. What i mean about that is that almost immideatly there becomes a "checkerboard pattern" of high/low pressure that grows and after a few seconds turn Everything into nan values. By making the temperature constant it works a bit longer and produce a shockwave pattern that looks very promesing but after a while also fails. It feels like something is overshooting in order to produce these oscillations.
– Eriksonn
Nov 17 at 22:48
Your code isn't behaving as expected, and so the first thing that you must do (and do before you come here) is to isolate the source of the bug -- meaning that you need to do intense debugging. Consider stepping through the program using your IDE's debugger, inspecting the state of key variables as you do so. This will give you valuable insights and quite possibly a decent solution. Even if you don't get a full answer, you'll be closer and can ask a better question. If you still need our help, you will need to create and post a Minimal, Complete, and Verifiable example, and tell the details of expected and observed behaviors
– Hovercraft Full Of Eels
Nov 17 at 22:02
Your code isn't behaving as expected, and so the first thing that you must do (and do before you come here) is to isolate the source of the bug -- meaning that you need to do intense debugging. Consider stepping through the program using your IDE's debugger, inspecting the state of key variables as you do so. This will give you valuable insights and quite possibly a decent solution. Even if you don't get a full answer, you'll be closer and can ask a better question. If you still need our help, you will need to create and post a Minimal, Complete, and Verifiable example, and tell the details of expected and observed behaviors
– Hovercraft Full Of Eels
Nov 17 at 22:02
@HovercraftFullOfEels To be fair: For many sorts of simulations, the usual debugging approaches are essentially useless. There often is no point in stepping through some simulation time step that involves complex math where you simply cannot say whether a certain value is "right" or "wrong" at a certain point (except for
NaN or so). Even more important is a small MCVE, and often the good ol' printlns, where you can at least print messages like "Now I'm adding these 1000 elements [...] to these [...]" to have the chance to spot where one vector goes off the rails.– Marco13
Nov 17 at 22:12
@HovercraftFullOfEels To be fair: For many sorts of simulations, the usual debugging approaches are essentially useless. There often is no point in stepping through some simulation time step that involves complex math where you simply cannot say whether a certain value is "right" or "wrong" at a certain point (except for
NaN or so). Even more important is a small MCVE, and often the good ol' printlns, where you can at least print messages like "Now I'm adding these 1000 elements [...] to these [...]" to have the chance to spot where one vector goes off the rails.– Marco13
Nov 17 at 22:12
Can you be more specific than "my code keeps messing up all the time."? What is the behavior you expect vs the behavior being observed?
– John Camerin
Nov 17 at 22:13
Can you be more specific than "my code keeps messing up all the time."? What is the behavior you expect vs the behavior being observed?
– John Camerin
Nov 17 at 22:13
@Marco13: I have to bow to your expertise since I don't do simulation programming, but something sure needs to be done to first isolate the error, somehow
– Hovercraft Full Of Eels
Nov 17 at 22:23
@Marco13: I have to bow to your expertise since I don't do simulation programming, but something sure needs to be done to first isolate the error, somehow
– Hovercraft Full Of Eels
Nov 17 at 22:23
The problem manifests as spatial oscillations growing over time. What i mean about that is that almost immideatly there becomes a "checkerboard pattern" of high/low pressure that grows and after a few seconds turn Everything into nan values. By making the temperature constant it works a bit longer and produce a shockwave pattern that looks very promesing but after a while also fails. It feels like something is overshooting in order to produce these oscillations.
– Eriksonn
Nov 17 at 22:48
The problem manifests as spatial oscillations growing over time. What i mean about that is that almost immideatly there becomes a "checkerboard pattern" of high/low pressure that grows and after a few seconds turn Everything into nan values. By making the temperature constant it works a bit longer and produce a shockwave pattern that looks very promesing but after a while also fails. It feels like something is overshooting in order to produce these oscillations.
– Eriksonn
Nov 17 at 22:48
|
show 8 more comments
active
oldest
votes
active
oldest
votes
active
oldest
votes
active
oldest
votes
active
oldest
votes
Eriksonn is a new contributor. Be nice, and check out our Code of Conduct.
draft saved
draft discarded
Eriksonn is a new contributor. Be nice, and check out our Code of Conduct.
Eriksonn is a new contributor. Be nice, and check out our Code of Conduct.
Eriksonn is a new contributor. Be nice, and check out our Code of Conduct.
draft saved
draft discarded
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53355940%2fprimitive-equations-on-square-mesh-is-diverging-computational-fluid-dynamics%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Your code isn't behaving as expected, and so the first thing that you must do (and do before you come here) is to isolate the source of the bug -- meaning that you need to do intense debugging. Consider stepping through the program using your IDE's debugger, inspecting the state of key variables as you do so. This will give you valuable insights and quite possibly a decent solution. Even if you don't get a full answer, you'll be closer and can ask a better question. If you still need our help, you will need to create and post a Minimal, Complete, and Verifiable example, and tell the details of expected and observed behaviors
– Hovercraft Full Of Eels
Nov 17 at 22:02
@HovercraftFullOfEels To be fair: For many sorts of simulations, the usual debugging approaches are essentially useless. There often is no point in stepping through some simulation time step that involves complex math where you simply cannot say whether a certain value is "right" or "wrong" at a certain point (except for
NaNor so). Even more important is a small MCVE, and often the good ol'printlns, where you can at least print messages like"Now I'm adding these 1000 elements [...] to these [...]"to have the chance to spot where one vector goes off the rails.– Marco13
Nov 17 at 22:12
Can you be more specific than "my code keeps messing up all the time."? What is the behavior you expect vs the behavior being observed?
– John Camerin
Nov 17 at 22:13
@Marco13: I have to bow to your expertise since I don't do simulation programming, but something sure needs to be done to first isolate the error, somehow
– Hovercraft Full Of Eels
Nov 17 at 22:23
The problem manifests as spatial oscillations growing over time. What i mean about that is that almost immideatly there becomes a "checkerboard pattern" of high/low pressure that grows and after a few seconds turn Everything into nan values. By making the temperature constant it works a bit longer and produce a shockwave pattern that looks very promesing but after a while also fails. It feels like something is overshooting in order to produce these oscillations.
– Eriksonn
Nov 17 at 22:48