//Conservation of mass, momentum and internal energy

//initial central high pressure/temperature

//velocity

let vx=[];

let vy=[];

let vz=[];

//momentum

let mx=[];

let my=[];

let mz=[];

//density

let rho=[];

//internal energy/temperature

let e=[];

//pressure

let p=[];

let phi=[];

//particle trace

let x=[];

let y=[];

let z=[];

//external force

let fx=[];

let fy=[];

let fz=[];

let N=40;

let dt;

let t=0;

let T=0;

let S=1;

let h;

let gamma=0.4;

​

​

function setup() {

  createCanvas(400, 400);

  for (var i=0;i<N;i++){

    x[i]=[];

    y[i]=[];

    z[i]=[];

    vx[i]=[];

    vy[i]=[];

    vz[i]=[];

    mx[i]=[];

    my[i]=[];

    mz[i]=[];

    rho[i]=[];

    e[i]=[];

    p[i]=[];

    phi[i]=[];

    fx[i]=[];

    fy[i]=[];

    fz[i]=[];

  for (var j=0;j<N;j++){

    x[i][j]=[];

    y[i][j]=[];

    z[i][j]=[];

    vx[i][j]=[];

    vy[i][j]=[];

    vz[i][j]=[];

    mx[i][j]=[];

    my[i][j]=[];

    mz[i][j]=[];

    fx[i][j]=[];

    fy[i][j]=[];

    fz[i][j]=[];

    rho[i][j]=[];

    e[i][j]=[];

    p[i][j]=[];

    phi[i][j]=[];

for (var k=0;k<N;k++){

    x[i][j][k]=10*i;

    y[i][j][k]=10*j;

    z[i][j][k]=10*k;

    vx[i][j][k]=0;

    vy[i][j][k]=0;

    vz[i][j][k]=0;

    mx[i][j][k]=0;

    my[i][j][k]=0;

    mz[i][j][k]=0;

    fx[i][j][k]=0;

    fy[i][j][k]=0;

    fz[i][j][k]=0;

    rho[i][j][k]=1;

    e[i][j][k]=1+pow(sin(PI*i/(N-1)),4)*pow(sin(PI*j/(N-1)),4)*pow(sin(PI*k/(N-1)),4);

    p[i][j][k]=gamma*e[i][j][k];

    phi[i][j][k]=0;

  }

  }

  }

  h=1/N;

  dt=0.05*h;

  for (var q=0;q<1000;q++){ 

for (var i=1;i<N-1;i++){

for (var j=1;j<N-1;j++){

for (var k=1;k<N-1;k++){

phi[i][j][k]=phi[i][j][k]+0.05*dt*((phi[i+1][j][k]+phi[i-1][j][k]-6*phi[i][j][k]+phi[i][j+1][k]+phi[i][j-1][k]+phi[i][j][k+1]+phi[i][j][k-1])/pow(h,2)+2*rho[i][j][k]); 

}

}

}

}

}

​

function draw() {

  background(255,255,255,10);

  noStroke();

  t=t+dt;

  T=T+dt;

for (var q=0;q<1;q++){ 

for (var i=1;i<N-1;i++){

for (var j=1;j<N-1;j++){

for (var k=1;k<N-1;k++){

phi[i][j][k]=phi[i][j][k]+0.05*dt*((phi[i+1][j][k]+phi[i-1][j][k]-6*phi[i][j][k]+phi[i][j+1][k]+phi[i][j-1][k]+phi[i][j][k+1]+phi[i][j][k-1])/pow(h,2)+2*rho[i][j][k]); 

}

}

}

}

​

for (var i=1;i<N-1;i++){

for (var j=1;j<N-1;j++){

for (var k=1;k<N-1;k++){

vx[i][j][k]=mx[i][j][k]/rho[i][j][k];

vy[i][j][k]=my[i][j][k]/rho[i][j][k];

vz[i][j][k]=mz[i][j][k]/rho[i][j][k];

p[i][j][k]=gamma*e[i][j][k]+2*pow(rho[i][j][k],2);

​

mx[i][j][k]=mx[i][j][k]+dt*(vx[i+1][j][k]+vx[i-1][j][k]-6*vx[i][j][k]+vx[i][j+1][k]+vx[i][j-1][k]+vx[i][j][k+1]+vx[i][j][k-1])/h-0.5*dt*(vx[i+1][j][k]*mx[i+1][j][k]-vx[i-1][j][k]*mx[i-1][j][k])/h-0.5*dt*(vy[i][j+1][k]*mx[i][j+1][k]-vy[i][j-1][k]*mx[i][j-1][k])/h-0.5*dt*(vz[i][j][k+1]*mx[i][j][k+1]-vz[i][j][k-1]*mx[i][j][k-1])/h-dt*gamma*(e[i+1][j][k]-e[i-1][j][k])/h+dt*rho[i][j][k]*(phi[i+1][j][k]-phi[i-1][j][k])/h;;

my[i][j][k]=my[i][j][k]+dt*(vy[i+1][j][k]+vy[i-1][j][k]-4*vy[i][j][k]+vy[i][j+1][k]+vy[i][j-1][k])/h-0.5*dt*(vx[i+1][j][k]*my[i+1][j][k]-vx[i-1][j][k]*my[i-1][j][k])/h-0.5*dt*(vy[i][j+1][k]*my[i][j+1][k]-vy[i][j-1][k]*my[i][j-1][k])/h-0.5*dt*(vz[i][j][k+1]*my[i][j][k+1]-vz[i][j][k-1]*my[i][j][k-1])/h-dt*gamma*(e[i][j+1][k]-e[i][j-1][k])/h+dt*rho[i][j][k]*(phi[i][j+1][k]-phi[i][j-1][k])/h;;

mz[i][j][k]=mz[i][j][k]+dt*(vz[i+1][j][k]+vz[i-1][j][k]-4*vz[i][j][k]+vz[i][j+1][k]+vz[i][j-1][k])/h-0.5*dt*(vx[i+1][j][k]*mz[i+1][j][k]-vx[i-1][j][k]*mz[i-1][j][k])/h-0.5*dt*(vy[i][j+1][k]*mz[i][j+1][k]-vy[i][j-1][k]*mz[i][j-1][k])/h-0.5*dt*(vz[i][j][k+1]*mz[i][j][k+1]-vz[i][j][k-1]*mz[i][j][k-1])/h-dt*gamma*(e[i][j][k+1]-e[i][j][k-1])/h+dt*rho[i][j][k]*(phi[i][j][k+1]-phi[i][j][k-1])/h;;

e[i][j][k]=e[i][j][k]+2*dt*(e[i+1][j][k]+e[i-1][j][k]-6*e[i][j][k]+e[i][j+1][k]+e[i][j-1][k]+e[i][j][k+1]+e[i][j][k-1])/h;

-0.5*dt*(vx[i+1][j][k]*e[i+1][j][k]-vx[i-1][j][k]*e[i-1][j][k])/h;

-0.5*dt*(vy[i][j+1][k]*e[i][j+1][k]-vy[i][j-1][k]*e[i][j-1][k])/h;

-0.5*dt*(vz[i][j][k+1]*e[i][j][k+1]-vz[i][j][k-1]*e[i][j][k-1])/h

-dt*gamma*e[i][j][k]*(vx[i+1][j][k]-vx[i-1][j][k]+vy[i][j+1][k]-vy[i][j-1][k]+vz[i][j][k+1]-vz[i][j][k-1])/h+dt*phi[i][j][k]*(mx[i+1][j][k]-mx[i-1][j][k]+my[i][j+1][k]-my[i][j-1][k]+mz[i][j][k+1]-mz[i][j][k-1])/h;;

rho[i][j][k]=rho[i][j][k]+2*dt*(rho[i+1][j][k]+rho[i-1][j][k]-6*rho[i][j][k]+rho[i][j+1][k]+rho[i][j-1][k]+rho[i][j][k+1]+rho[i][j][k-1])/h

-dt*(mx[i+1][j][k]-mx[i-1][j][k])/h;

-dt*(my[i][j+1][k]-my[i][j-1][k])/h

-dt*(mz[i][j][k+1]-mz[i][j][k-1])/h;

}

}

}

if (t>0.1){

for (var i=0;i<N;i++){

for (var j=0;j<N;j++){

    x[i][j][20]=10*i;

    y[i][j][20]=10*j;

}

}

t=0;

background(255,255,255,255);

}

if (T>1 & S==0){

S=0;

for (var i=0;i<N;i++){

for (var j=0;j<N;j++){

for (var k=0;k<N;k++){

//    vx[i][j][k]=-vx[i][j][k];

//    vy[i][j][k]=-vy[i][j][k];

//    vz[i][j][k]=-vz[i][j][k];

    mx[i][j][k]=-mx[i][j][k];

    my[i][j][k]=-my[i][j][k];

    mz[i][j][k]=-mz[i][j][k];

}

}

}

}

​

fill(255,0,0,255);

for (var i=0;i<N;i++){

for (var j=0;j<N;j++){

  x[i][j][20]=x[i][j][20]+5000*vx[i][j][20]*dt;  

  y[i][j][20]=y[i][j][20]+5000*vy[i][j][20]*dt;

  ellipse(x[i][j][20],y[i][j][20],3);

  }

  }

for (var i=0;i<N;i++){

  fill(0,0,255,255);

  ellipse(10*i,300-100*rho[i][20][20],10);

  fill(0,255,0,255);

  ellipse(10*i,300-100*e[i][20][20],10);

  fill(0,255,255,255);

  ellipse(10*i,300-2000*phi[i][20][20],10);

}

​

  fill(0);

  text("Cosmology 3d: Gravitation (contraction) vs Fluid Pressure (expansion)",10,30);

  text("momentum (mx,my,mz), density rho, internal energy e, grav pot phi",10,70);

  text("drho/dt+d(mx)/dx+d(my)/dy+d(mz)/dz=0",10,110);

  text("dmx/dt+d(mx*vx)/dx+d(mx*vy)/dy+d(mx*vz)/dz+dp/dx-rho*dphi/dx = 0",10,130);

  text("dmy/dt+d(my*vx)/dx+d(my*vy)/dy+d(my*vz)/dz+dp/dy-rho*dphi/dy=0",10,150);

  text("dmz/dt+d(mz*vx)/dx+d(mz*vy)/dy+d(mz*vz)/dz+dp/dz- rho*dphi/dz=0",10,170);

  text("de/dt+d(e*vx)/dx+d(e*vy)/dy+d(e*vz)/dz+p*(dvx/dx+dvy/dy+dvz/dz)",10,190);

  text("-phi*(dmx/dx+dmy/dy+dmz/dz) = 0",10,210);

  text("-Laplacian phi = 0.5*rho",1,230);

  text("vx = mx/rho, vy = my/rho, vz = mz/rho p =gamma*e, 0 < gamma < 1.",10,250);

  text("Conservation of mass, momentum, internal energy:" ,10,50);

  text("Initial data: zero velocity, high central pressure/temperature",10,270);

  text("Test different initial and boundary conditions.",10,370);

  text("Blue: Density (line crosscut), Light-Blue: Gravitational Energy",10,310);

  text("Green: Internal Energy (line crosscut).",10,330);

  text("Red: Trace of fluid particles (plane crosscut updated).",10,350);

}