Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Shallow water flow - equation system unstable

  1. Apr 12, 2012 #1
    Hello, everyone.

    I'm working on simulating shallow water flow. I'm using the MacCormack scheme and fractional predictor-corrector steps. The Dxy factor should make the equation system stable but no matter what numbers I input, the results skyrocket, up to 10^80 and more.

    Has anyone worked with this before? If so, could you tell me what I'm doing wrong?


    Thank you,
    Martin





    Below is the code I've written in Mathematica (since I doubt this is an error by Mathematica, I think this doesn't belong in Computing part of the forum):

    g = 9.81; q1 = 0.0001; dzdx = -0.005; dzdy = -0.01; kov = 1.5*10^-3; dt = 0.0000005; dx =
    dy = 0.0005; dim = 20; qx = 0.00009; h = (q1 - qx/dx)*dt

    U[1] = Table[{h, qx, 0}, {y, 1, dim}, {z, 1, dim}];
    G[x_] := {x[[2]], x[[2]]*x[[2]]/x[[1]] + g/2*x[[1]]*x[[1]], x[[2]]*x[[3]]/x[[1]]};
    H[x_] := {x[[3]], x[[2]]*x[[3]]/x[[1]], x[[3]]*x[[3]]/x[[1]] + g/2*x[[1]]*x[[1]]};
    S[x_] := {q1, -g*dzdx*x[[1]] - kov/8*x[[2]]/x[[1]]/x[[1]], g*dzdy*x[[1]] - kov/8*x[[3]]/x[[1]]/x[[1]]};
    Dxy[x_] := {1, 1/(1 - dt*kov/8/x[[1]]/x[[1]]/x[[1]]), 1/(1 - dt*kov/8/x[[1]]/x[[1]]/x[[1]])};

    UU2 = ConstantArray[0, {dim, dim}];
    UU3 = ConstantArray[0, {dim, dim}]; UU4 = ConstantArray[0, {dim, dim}];
    UU5 = ConstantArray[0, {dim, dim}]; UU6 = ConstantArray[0, {dim, dim}];
    UU7 = ConstantArray[0, {dim, dim}]; UU8 = ConstantArray[0, {dim, dim}];
    UU9 = ConstantArray[0, {dim, dim}];
    UU = U[1];

    (*Predictor Lx1*)
    For[j = 1, j <= dim, j++,
    (*Calc. first line*)
    UU2[[1, j]] =
    UU[[1, j]] -
    Dxy[UU[[1, j]]]*dt/dx/2*(G[UU[[1, j]]] - G[UU[[2, j]]]) +
    Dxy[UU[[1, j]]]*dt/2*S[UU[[1, j]]];
    For[i = 2, i <= dim, i++,
    (*Calc. other lines*)
    UU2[[i, j]] =
    UU[[i, j]] -
    Dxy[UU[[i, j]]]*dt/dx/2*(G[UU[[i, j]]] - G[UU[[i - 1, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU[[i, j]]];
    ];
    ];

    (*Corrector Lx1*)
    For[j = 1, j <= dim, j++,
    For[i = 1, i <= dim - 1, i++,
    (*Calc. other lines*)
    UU3[[i, j]] =
    0.5*(UU[[i, j]] - UU2[[i, j]] -
    Dxy[UU[[i, j]]]*
    dt/2/dx*(G[UU2[[i + 1, j]]] - G[UU2[[i, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU2[[i, j]]]);
    ];
    (*Calc. last line*)
    UU3[[dim, j]] =
    0.5*(UU[[dim, j]] - UU2[[dim, j]] -
    Dxy[UU[[dim, j]]]*
    dt/2/dx*(G[UU2[[dim - 1, j]]] - G[UU2[[dim, j]]]) +
    Dxy[UU[[dim, j]]]*dt/2*S[UU2[[dim, j]]]);
    ];

    (*Predictor Ly1*)
    For[i = 1, i <= dim, i++,
    (*Calc. first collumn*)
    UU4[[i, 1]] =
    UU3[[i, 1]] -
    Dxy[UU[[i, 1]]]*dt/2/dx*(H[UU3[[i, 1]]] - H[UU3[[i, 2]]]) +
    Dxy[UU[[i, 1]]]*dt/2*S[UU3[[i, 1]]];
    (*Calc. other collumns*)
    For[j = 2, j <= dim, j++,
    UU4[[i, j]] =
    UU3[[i, j]] -
    Dxy[UU[[i, j]]]*dt/2/dx*(H[UU3[[i, j]]] - H[UU3[[i, j - 1]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU3[[i, j]]];
    ];
    ];

    (*Corrector Ly1*)
    For[i = 1, i <= dim, i++,
    (*Calc. other collumns*)
    For[j = 1, j <= dim - 1, j++,
    UU5[[i, j]] =
    0.5*(UU3[[i, j]] - UU4[[i, j]] -
    Dxy[UU[[i, j]]]*
    dt/2/dx*(H[UU4[[i, j + 1]]] - H[UU4[[i, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU4[[i, j]]]);
    ];
    (*Calc. last collumn*)
    UU5[[i, dim]] =
    0.5*(UU3[[i, dim]] - UU4[[i, dim]] -
    Dxy[UU[[i, dim]]]*
    dt/2/dx*(H[UU4[[i, dim - 1]]] - H[UU4[[i, dim]]]) +
    Dxy[UU[[i, dim]]]*dt/2*S[UU4[[i, dim]]]);
    ];

    (*Predictor Ly2*)
    For[i = 1, i <= dim, i++,
    (*Calc. other collumns*)
    For[j = 1, j <= dim - 1, j++,
    UU6[[i, j]] =
    UU5[[i, j]] -
    Dxy[UU[[i, j]]]*dt/2/dx*(H[UU5[[i, j + 1]]] - H[UU5[[i, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU5[[i, j]]];
    ];
    (*Calc. last collumn*)
    UU6[[i, dim]] =
    UU5[[i, dim]] -
    Dxy[UU[[i, dim]]]*
    dt/2/dx*(H[UU5[[i, dim - 1]]] - H[UU5[[i, dim]]]) +
    Dxy[UU[[i, dim]]]*dt/2*S[UU5[[i, dim]]];
    ];

    (*Corrector Ly2*)
    For[i = 1, i <= dim, i++,
    (*Calc. first collumn*)
    UU7[[i, 1]] =
    0.5*(UU5[[i, 1]] + UU6[[i, 1]] -
    Dxy[UU[[i, 1]]]*dt/2/dx*(H[UU6[[i, 1]]] - H[UU6[[i, 2]]]) +
    Dxy[UU[[i, 1]]]*dt/2*S[UU6[[i, 1]]]);
    (*Calc. other collumns*)
    For[j = 2, j <= dim, j++,
    UU7[[i, j]] =
    0.5*(UU5[[i, j]] + UU6[[i, j]] -
    Dxy[UU[[i, j]]]*
    dt/2/dx*(H[UU6[[i, j]]] - H[UU6[[i, j - 1]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU6[[i, j]]]);
    ];
    ];

    (*Predictor Lx2*)
    For[j = 1, j <= dim, j++,
    (*Calc. other lines*)
    For[i = 1, i <= dim - 1, i++,
    UU8[[i, j]] =
    UU7[[i, j]] -
    Dxy[UU[[i, j]]]*dt/2/dx*(G[UU7[[i + 1, j]]] - G[UU7[[i, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU7[[i, j]]];
    ];
    (*Calc. last line*)
    UU8[[dim, j]] =
    UU7[[dim, j]] -
    Dxy[UU[[dim, j]]]*
    dt/2/dx*(G[UU7[[dim - 1, j]]] - G[UU7[[dim, j]]]) +
    Dxy[UU[[dim, j]]]*dt/2*S[UU7[[dim, j]]];
    ];

    (*Corrector Lx2*)
    For[j = 1, j <= dim, j++,
    (*Calc. first line*)
    UU9[[1, j]] =
    0.5*(UU7[[1, j]] + UU8[[1, j]] -
    Dxy[UU[[1, j]]]*dt/2/dx*(G[UU8[[1, j]]] - G[UU8[[2, j]]]) +
    Dxy[UU[[1, j]]]*dt/2*S[UU8[[1, j]]]);
    (*Calc. other lines*)
    For[i = 2, i <= dim, i++,
    UU9[[i, j]] =
    0.5*(UU7[[i, j]] + UU8[[i, j]] -
    Dxy[UU[[i, j]]]*
    dt/2/dx*(G[UU8[[i, j]]] - G[UU8[[i - 1, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU8[[i, j]]]);
    ];
    ];
     
  2. jcsd
  3. Apr 12, 2012 #2
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Shallow water flow - equation system unstable
Loading...