problem with runge-kutta algorithm solving system of ODEs


Ocean
11-14-2009, 06:27 AM
Dear Forum,

I tried to use a Runge-Kutta method to solve a system of three ODEs.
The code is splitted in the main program, the ODE subroutine and a solver subroutine.
I am very new to Fortran, therefore my problem could be very basic and resulting from misunderstanding of basic fortran concepts, even though I read quite a few Introductions explaining the exchange of values between the main program and subroutines.
My problem is that the calculated values of the state variables (N,P,Z) in the ODE and the solver subroutines are not returned to the main program and hence the state variables stay constant over the time loop.
Attached is the basic code outline:

Thank you very much in advance for your help!








PROGRAM NPZ

USE solver

(Declarations...)

(Initial Cond.)

! combine all state variables:
y(1) = N
y(2) = P
y(3) = Z


DO k = 1,tend,dt


CALL rk4(y,dydt,yout) ! Call the solver subroutine

N = y(1)
P = y(2)
Z = y(3)

END DO

END PROGRAM




MODULE solver
CONTAINS



SUBROUTINE ode(y,f)

(Declarations)


! De-combine state variables
N = y(1)
P = y(2)
Z = y(3)

dNdt = ...
dPdt = ...
dZdt = ...

f(1) = dNdt
f(2) = dPdt
f(3) = dZdt

dydt = f


RETURN
END SUBROUTINE ode







SUBROUTINE rk4(y,dydt,yout)

(Declarations)

! De-combine state variables
N = y1(1)
P = y1(2)
Z = y1(3)



CALL ode(y,dydt)
DO i=1,imax
s1(i) = dydt(i)
END DO


DO i=1,imax
ytemp1(i) = y(i) + (dt/2.)*s1(i)
END DO


CALL ode(ytemp1,dydt)
DO i=1,imax
s2(i) = dydt(i)
END DO


DO i=1,imax
ytemp2(i) = y(i)+(dt/2.)*s2(i)
END DO

CALL ode(ytemp2,dydt)
DO i=1,imax
s3(i) = dydt(i)
END DO



DO i=1,imax
ytemp3(i) = y(i)+dt*s3(i)
END DO


CALL ode(ytemp3,dydt)
DO i=1,imax
s4(i) = dydt(i)
END DO


DO i=1,imax
yout(i) = y(i) + (dt/6.)*(s1(i) + 2*s2(i) + 2*s3(i) + s4(i))
y(i) = yout(i)
dydt(i) = s1(i)
END DO



RETURN
END SUBROUTINE rk4




END MODULE solver

Ocean
11-14-2009, 11:08 PM
Another problem is that a namelist that contains all the parameters is not global available to the subroutines.

I define and read the namelist in the main program as follows:


REAL :: alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero
INTEGER :: dt, tend

NAMELIST /NPZNAME/ alpha, vmaxNP, vmaxPZ, Ks, mP, mZ, Phalf, Pzero, dt, tend

OPEN(UNIT=41,FILE='npz.nml',STATUS='OLD')
READ(41,NML=NPZNAME)


What statement do I have to use in the subroutines that they can access those parameters?
Thanks!

davekw7x
11-15-2009, 12:46 PM
...
My problem is ...
If you have studied textbooks, course material, class examples, and/or other reference material, and if you have compared your code with examples that you found, and if you still can't see why things don't happen the way that you expect, then a brute-force debugging technique is:

In the main program put a print statement for the values of the variables just before the call to the subprogram.


Put a print statement at the beginning of the subprogram that shows values of the variables that it is working with.


Put a print statement just before returning from the subprogram that shows the newly computed values of the variables.


Back in the main program put a print statement just after the statement that calls the subprogram so that you can see the results.


If you have already done this, how about showing us exactly what happened? (Show the code and show the printed output values, and tell us exactly what expected and what you don't understand.)



Attached is the basic code outline...
How the heck can anyone give meaningful suggestions without seeing exactly how the variables are declared and what calculations are executed?

If you have performed the steps that I enumerated above, and if you still can't see where your problem is then I suggest that you do the following:

Instead of posting your entire program, create a minimal subset that makes your problem show up. Maybe a single subprogram that alters a single variable from the main program. If it's still broken, then post the entire (simple) program.

If the simpler program seems to work, then carefully add stuff a bit at a time and see where it breaks.


What statement do I have to use in the subroutines..


The method preferred by many (most?) modern programmers (and, I think, just about all programmers that have to debug and/or maintain Other People's Code):
Pass them as arguments to the subprogram.


A Golden Oldie, used by Fortran programmers for something like 50 years (maybe more):
Put the variables in a common block. Declare the common block in the main program and in all called subprograms that need access to those variables.



Regards,

Dave