Shallow water flow - equation system unstable

[MartinV]

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]]]);
];
];

 

[MartinV]

I'm using this article as a reference:

http://www.engr.colostate.edu/~ramirez/ce_old/projects/Fiedler_Ramirez.pdf