How to Integrate a Single ODE in F90 for Water Flow Simulation?

[amberf90]

Hi all,

I've written a simulation of water flow in two dimensions in F90 but I'm having some trouble with it. Water flows from one cell to another using an equation for rate of change of depth and an algorithm for assigning flow direction.

The flow direction bit is fine but the dD/dt equation is giving me some trouble. I've tried an Euler forward method of integrating this equation through time but more water wants to move than is in the cell and I have to use a very fine time-step to allow the water to get where it wants to go, which is computationally expensive. A colleague suggested I adopt an implicit integrator as this will allow the hydraulic head (which depth is linked with) to reduce gradually, and presumably slow down the flow.

I've tried using two subroutines I found online; bsstep.f and ode.f but I can't manage to get either to work successfully. One of the problems is that the r.h.s includes values that need to be passed through from the main program. Also they talk about a system of ODEs but as far as I'm aware I only have 1. I'm not the best at maths sadly...

Does anybody have any ideas about simple ODE integrators, for one equation only, that aren't to difficult/complicated to use?

Thanks in advance...

 

[rcgldr]

Try Runge Kutta 4, wiki article:

wiki_rk4.htm

 

[amberf90]

Thanks for the suggestion; I was really looking for an implicit solution though...

 

[rcgldr]

Thanks for the suggestion; I was really looking for an implicit solution though...

You could consider RK4 to be "semi-implicit", since it calculates 3 intermediate increments using feedback from the initial euler increment and the first 2 intermediate increments. RK4 is more accurate than a simple predictor - corrector method such as the euler - trapezoidal example in this wiki article:

wiki_predictor_corrector_method.htm

Note that repeating correction step will quickly converge to a specific value, but since it's based on trapezoidal rule, it's a linear aproximation, versus RK4's 4th order approximation. If you did 1 euler and 3 trapezoidal correction steps, the overhead would be similar to RK4, but the result would not be as accurate (you'd need a much smaller step size).

 

[amberf90]

Thanks,

I'll code up the RK4 and see what happens!

A