void setupInitialConditions(Array<float,3>& P1, Array<float,3>& P2,
    Array<float,3>& P3, Array<float,3>& c, int N)
{
    // Set the velocity field
    c(Range(0,N/2-1), Range::all(), Range::all()) = 0.05;
    c(Range(N/2,N-1), Range::all(), Range::all()) = 0.3;

    int cavityLeft = 3*N/7.0-1;
    int cavityRight = 4*N/7.0-1;
    int cavityFront = 3*N/7.0-1;
    int cavityBack = 4*N/7.0-1;
    int cavityTop = 5*N/7.0-1;
    int cavityBottom = 6*N/7.0-1;

    c(Range(cavityTop,cavityBottom),Range(cavityLeft,cavityRight), 
        Range(cavityFront,cavityBack)) = 0.02;

    int cavityTop2 = 1*N/7.0-1;
    int cavityBottom2 = 2*N/7.0-1;
    c(Range(cavityTop2,cavityBottom2),Range(cavityLeft,cavityRight),
        Range(cavityFront,cavityBack)) = 0.001;

    // Initial pressure distribution
    using namespace blitz::tensor;
    float ci = N/2-1;
    float cj = N/2-1;
    float ck = N/2-1;
    float s2 = 64.0 * 9.0 / pow2(N/2.0);
    P1 = 0.0;
    P2 = exp(-(pow2(i-ci)+pow2(j-cj)+pow2(k-ck)) * s2);
    P3 = 0.0;
}