1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
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...