MathOnco / HAL

HAL (Hybrid Automata Library) is the one-stop shop for all of your hybrid modeling needs
MIT License
37 stars 12 forks source link

Bug in advection solver: Negative y speeds not treated correctly #5

Closed JulianoGianlupi closed 2 months ago

JulianoGianlupi commented 2 months ago

Hello, while playing around with Examples/_4PDEexample/PDEexample.java I noticed that the jet never went down. I confirmed that it should go down by modifying PDEexample.Step in with the following:

// double advectionX=Math.sin(stepI*1.0/1000)*0.2;//sine and cosine based on timestep cause circular advection
// double advectionY=Math.cos(stepI*1.0/1000)*0.2;

if (stepI !=0 && stepI % 1000 == 0) {
    rotator += 1.0 / 4.0;

}

double angle = Math.PI * rotator;

if (stepI == 0 || stepI % 1000 == 0) {
    System.out.println("rotator: " + rotator + " angle: " + angle%Math.PI + "pi");
}

double advectionX=Math.cos(angle)*0.2;//sine and cosine based on timestep cause circular advection
double advectionY=Math.sin(angle)*0.2;

(I initialized rotator as a field=0)

After further investigation I found that the Advection2 that is called from that example uses what seems to be the wrong calculation for DeltaY2D:

public static void Advection2(double[]field,double[]deltas,double xVel,double yVel,int xDim,int yDim,boolean wrapX,boolean wrapY,Coords2DDouble BC) {
        int disp = 0;
        for (int x = 0; x < xDim; x++) {
            for (int y = 0; y < yDim; y++) {
                int i = x * yDim + y;
                double centerVal=field[i];
                if (xVel > 0) {
                    deltas[i] += DeltaX2D(field,centerVal, x - 1,y, xDim,yDim, wrapX, BC)*xVel;
                } else if(xVel < 0){
                    deltas[i] -= DeltaX2D(field,centerVal, x + 1,y, xDim,yDim, wrapX, BC)*xVel;
                }
                if (yVel > 0) {
                    deltas[i] += DeltaY2D(field,centerVal, x ,y-1, xDim,yDim, wrapY, BC)*yVel;
                } else if(yVel < 0){
                    deltas[i] -= DeltaY2D(field,centerVal, x ,y-1, xDim,yDim, wrapY, BC)*yVel;
                }
            }
        }
    }

For negative x velocities a positive pixel displacement is used, DeltaX2D(field,centerVal, x + 1,y, xDim,yDim, wrapX, BC), while negative y velocities use a negative displacement, DeltaY2D(field,centerVal, x ,y-1, xDim,yDim, wrapY, BC). Flipping the displacement there fixed the simulation

The other Advection solvers should also be checked