Adding variables causes odd issue 2

[kbrownk]

Hi,

I've been stuck on a project for about 3 years and it has really stalled my hopes of becoming an adequate programmer. Any help would be so appreciated!

In short, when I have 33 variables y(1) through y(33) the code runs great. However, once I add in a 34th variable y(34) the code does something odd. I am using a runge-kutta method as 1 of 2 subroutines, which determines the size of time steps (= dt)for each step iteration. This time step becomes increasingly smaller for each time step causing the code to never finish.

It doesn't matter what the 34th variable is as I can set it to 0 and never use it.

I would be happy to condense and post the code if anyone feels it would help, but wanted to see if there was a general explanation first. I'm such a novice that I'm not sure what other info would be of use, so please let me know and I will gratefully post.

I posted this issue as well as the code on google groups a couple of years back but no one could offer an explanation!

I am using cygwin X on a Windows XP PC w/ a single 3GHz Intel pentium processor.

Thank you!
kbrownk

 

[mikrom]

y(1) through y(33) are not the variables, but the elements of an array y. Have you declared the size of the array properly?

 

[kbrownk]

I believe so, the beginning of the main routine is as follows:

implicit double precision (a-h,o-z)
parameter(neq=34,nflux=25,nstepmax=10000)
dimension y(neq),dy(neq)
dimension flux(nflux)


neq has been changed from less than 30 up to 44 in the past and it never works once it gets beyond 33. neq < OR = 33 gives reasonable/expected values.

 

[xwb]

There is possibly something hardcoded or a bit pattern somewhere. Another possibility is a variable name going over column 72. Since you have implicit DP, anything not declared will be DP, even misspelt variables.

 

[kbrownk]

I am confident I am not going over 72 columns and have not misspelled anything

Does hardcoding refer to the machine language used in cygwin or the fortran compiler I use in my cygwin version? Maybe something w/ my computer or OS? I have tried running the code in Salford Plato software which is an editor/compiler that runs directly on Windows, and I get the same problem. I'll try setting up Cygwin on a different PC and see what happens. I had tried on a UNIX machine a couple years ago, but I forget if I had the exact same issue (I feel like I did, though).

Thanks,
K

 

[xwb]

By hardcoding, I mean there is a 32, 33 or 34 somewhere in the code. If you are getting the same result on two compilers, the problem must be in the code. If you could post the code, maybe someone will have a look (you may get lots of answers on Christmas day when the MVPs get bored and look to TT for entertainment).

 

[mikrom]

Explain us please: Try you to solve one ODE or a system of ODEs and what Runge-Kutta Method are you using - explicit or implicit?
It's possible that the ODE or the method is causing the problem.
xwb is right - it would be best when you post the code. Maybe I will look at it at Christmas day or later

 

[kbrownk]

Thanks! I'll clean it up and post ASAP.
-K

 

[kbrownk]

Hi,

I've copied all the main and subroutines all to a single text file that I uploaded to

http://www.sendspace.com/file/0bnrye

I hope that is an OK way to provide it. There is nothing proprietary in the code, but I still didn't want it just hanging around on the internet in its current state.

I copied the makefile I use below.

Anyway, these 3 routines were initially 3 .f files and I did not provide the output files they write to since they change and are quite large.

I placed the non-working (i.e 34 array variables) version, where if I change it to 33 and character out the 34th variable in the main routine and fnc subroutine, the code works fine. Specifically, it will end after writing about 20,000 lines of output on the fort. files.

Please let me know if anything is unclear. I am obviously a novice programmer, so I apologize if I could have made the code easier to read, and would be happy to clear up anything.

Thanks again, I'm pretty excited that there's a chance this issue could finally get cleared up! -kbk

Here's the Makefile:

FFLAGS=-O3 +extend_source
F77=f77
fortkbk: main.o fnc.o rk4m.o
$(F77) $(FFLAGS) *.o -o fortkbk

.f.o:
$(F77) $(FFLAGS) -c $<

 

[mikrom]

I downloaded the file an compiled it with

g77 Fortran_kbk.f -o Fortran_kbk

It runs and produces these files:

fort.8
...
fort.14
fort.17
...
fort.21
rfinal.dat

All the files are filled wit values, e.g. rfinal.dat contain this

 -5.14101E+01  9.90000E-05  1.89623E+00  5.26480E-04  1.50000E-03  9.90000E-04  9.90000E-04  9.90000E-04  9.99900E-01  9.99900E-01  3.17131E-01  2.48843E-01  1.96440E-04  8.22870E-01  9.02719E-01  2.20932E-02  9.45389E-01  9.99900E-01  1.53463E-04  9.99992E-01  1.37206E-01  3.25527E-01  2.44761E-01  3.54351E-01  9.64725E-01  2.40000E-02  2.60939E+01  1.31057E-03  1.17801E+00  2.61958E-02
  2.00169E-30  4.42039E-02  1.84683E-01  2.05250E+00

But when I try to compile the program with the option -O3 it seems to run in an infinite loop.
I cannot debug it with the option -O3 which means code optimization, because it doesn't step right in the debugger.

So it seams, that your problem is caused by the optimization option -O3.

 

[xwb]

There is a problem around line 400

      dy(1) =-(aICaT+aICaN+aICaL+aINaf+aINas+aIK+aIA+aIDTX
     +aIBK+aISK+aIBCa+aIBNa+aINaCa+aICaP+aINaK+stmm)/Cm

Because you have implicit DP a-h, this one doesn't show up. Note that spaces are allowed in variable names in Fortran. The + is a continuation character so aIDTX is concatenated to aIBK to make aIDTXaIBK. I think the code needs to be

      dy(1) =-(aICaT+aICaN+aICaL+aINaf+aINas+aIK+aIA+aIDTX[b][COLOR=red]+[/color][/b]
     +aIBK+aISK+aIBCa+aIBNa+aINaCa+aICaP+aINaK+stmm)/Cm

Otherwise you have ... aIA + aIDTXaIBK + ...

 

[kbrownk]

Hi mikrom and xwb,

xwb, you're right about the concatenation issue. I fixed it and am currently running, but this a recent issue as I've had this right before and had the same issue, even after charactering out this equation and making = 0.0.

mikrom, I never considered tweaking the makefile, since I never really understood it (I didn't write it). I'm unsure how I can change the make file to something similar to what you used, but am looking in to it. In the meantime, after fixing what xwb pointed out, I am trying using -O2. Cygwin didn't like me avoiding the makefile altogether and just using g77 Fortran_kbk.f -o Fortran_kbk. So if -O2 doesn't work, I'll have some research to do.

Thanks to both for your generosity, I think a combination of these types of fixes is required and should get it up and running hopefully soon!

-kbk

 

[kbrownk]

Mikrom,

Do you know if you got more than 1 line of values in your fort. files when you ran it using:

g77 Fortran_kbk.f -o Fortran_kbk

I tried using this command in cygwin xterm (x11) instead of make but it completes immediately w/ only one iteration of values. When I use the makefile and only 33 variables I get over 20,000 lines of values with varying timestep sizes and value changes I would expect to see.

I tried changing the timestep-related variables in case the runge-kutta equations were creating too big of timesteps without the optimization being used, but it doesn't help.

Thanks again,
kbk

 

[xwb]

When I ran it on IVF with NEQ=39, it started going wrong at j=2, i=277. dy jumps from -5E-1 to -9E61. Thereafter, everything goes wrong.

I don't know a lot about RK so I can't help you there.

 

[mikrom]

Hi kbrownk,

I used g77 compiled with MinGW o Windows Vista.
Without -O3 the long loop

c     Start integration time loop
    2 do 20 j=1,nstepmax
    ....
   20 continue

was terminated by condition

if (y(33).gt.1.0) goto 255

The files produced fort.8,..,fort.21 are 31 lines long.

With the option -O3 the above loop was not terminated.

 

[kbrownk]

xwb, I've suspected the rk method for a long time but keep coming back to the fact that an extra variable that is never used should not have any effect on it. I have noticed using different programs that the calculations simply go haywire at some point once I have more than 33 variables, but I've never figured out why.

Thanks,
kbk

 

[mikrom]

Hi kbrownk,

Runge-Kutta is for integrating ODE of the form
y'(x) = g(x, y(x))
y(x0) = y0
For one ODE the program is something like this (without step control)

[COLOR=#0000ff]C-----Main Program for RK-method[/color]
      [COLOR=#2e8b57][b]EXTERNAL[/b][/color] G
      [COLOR=#804040][b]WRITE[/b][/color]([COLOR=#ff00ff]6[/color],[COLOR=#ff00ff]1[/color])
  [COLOR=#6a5acd]1[/color]   [COLOR=#804040][b]FORMAT[/b][/color]([COLOR=#008080]10X[/color],1HX,[COLOR=#008080]15X[/color],4HY(X))
      X[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]0[/color].
      Y[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color].
      [COLOR=#804040][b]WRITE[/b][/color]([COLOR=#ff00ff]6[/color],[COLOR=#ff00ff]3[/color])X,Y
  [COLOR=#6a5acd]2[/color]   [COLOR=#a020f0]CALL[/color] RK(G,X,Y,[COLOR=#ff00ff]0.1[/color])
      [COLOR=#804040][b]WRITE[/b][/color]([COLOR=#ff00ff]6[/color],[COLOR=#ff00ff]3[/color])X,Y
  [COLOR=#6a5acd]3[/color]   [COLOR=#804040][b]FORMAT[/b][/color]([COLOR=#008080]F13.2[/color],[COLOR=#008080]E18.6[/color])
      [COLOR=#804040][b]IF[/b][/color](X[COLOR=#804040][b].LT.[/b][/color][COLOR=#ff00ff]1[/color].) [COLOR=#804040][b]GO TO[/b][/color] [COLOR=#ff00ff]2[/color]
      [COLOR=#804040][b]STOP[/b][/color]
      [COLOR=#a020f0]END[/color]
[COLOR=#0000ff]C[/color]
[COLOR=#0000ff]C-----Method Runge-Kutta of 4.order [/color]
      [COLOR=#a020f0]SUBROUTINE[/color] RK(F,X,Y,H)
[COLOR=#2e8b57][b]      REAL[/b][/color] K([COLOR=#ff00ff]4[/color])
      K([COLOR=#ff00ff]1[/color])[COLOR=#804040][b]=[/b][/color]H[COLOR=#804040][b]*[/b][/color]F(X,Y)
      X[COLOR=#804040][b]=[/b][/color]X[COLOR=#804040][b]+[/b][/color]H[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]2[/color].
      K([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]=[/b][/color]H[COLOR=#804040][b]*[/b][/color]F(X,Y[COLOR=#804040][b]+[/b][/color]K([COLOR=#ff00ff]1[/color])[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]2[/color].)
      K([COLOR=#ff00ff]3[/color])[COLOR=#804040][b]=[/b][/color]H[COLOR=#804040][b]*[/b][/color]F(X,Y[COLOR=#804040][b]+[/b][/color]K([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]2[/color].)
      X[COLOR=#804040][b]=[/b][/color]X[COLOR=#804040][b]+[/b][/color]H[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]2[/color].
      K([COLOR=#ff00ff]4[/color])[COLOR=#804040][b]=[/b][/color]H[COLOR=#804040][b]*[/b][/color]F(X,Y[COLOR=#804040][b]+[/b][/color]K([COLOR=#ff00ff]3[/color]))
      Y[COLOR=#804040][b]=[/b][/color]Y[COLOR=#804040][b]+[/b][/color](K([COLOR=#ff00ff]1[/color])[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]2[/color].[COLOR=#804040][b]*[/b][/color]K([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]2[/color].[COLOR=#804040][b]*[/b][/color]K([COLOR=#ff00ff]3[/color])[COLOR=#804040][b]+[/b][/color]K([COLOR=#ff00ff]4[/color]))[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]6[/color].
      [COLOR=#804040][b]RETURN[/b][/color]
      [COLOR=#a020f0]END[/color]
[COLOR=#0000ff]C[/color]
[COLOR=#0000ff]C-----Function from the right side of ODE[/color]
      [COLOR=#a020f0]FUNCTION[/color] G(X,Y)
      G[COLOR=#804040][b]=[/b][/color][COLOR=#008080]SIN[/color](X[COLOR=#804040][b]+[/b][/color]Y)[COLOR=#804040][b]/[/b][/color][COLOR=#008080]SQRT[/color](X[COLOR=#804040][b]*[/b][/color]Y[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color].)
      [COLOR=#804040][b]RETURN[/b][/color]
      [COLOR=#a020f0]END[/color]

If I good understand from your example, you try to solve a system of 33 or 34 ODEs.
I don't understand in your program these things:
1. What is the right side of your ODE system? It must be the function which returns vector of dimension 33 or 34.
2. Are you using the correct RK-formula? I coudn't find such coefficient for method of order 4.

The method you are using is explicit. But not every ODE or system of ODEs could be solved by explicit methods. Some ODEs are bad preconditioned - so called stiff ODEs. These could not be solved with explicit RK-methods - nor the step control helps here. For solving stiff ODEs there are implicit RK-methods suitable.

So, even if your program should be error free, then the stiffness may cause a problem: Maybe the system of 33 ODEs is good preconditioned for solving with explicit RK, but when you add one more ODE, so the system of 34 ODEs may get bad preconditioned - stiff.

One case where stiff system of ODEs occur is by solving parabolic PDEs as a system of ODEs. By discretization of the space component of parabolic PDE we get a large system of time dependent ODEs which is usually stiff.

Programming implicit RK solver is not so simple as explicit. There is a good book by Hairer, Norset, Wanner on this topic available. You can download the fortran codes too:

http://www.unige.ch/~hairer/software.html

 

[kbrownk]

mikrom, I think you have hit upon what my target should be which is what is used for the rkm. A lot of the things you discuss here I was beginning to learn when I first started working on this program, but due to the long hiatus I am going to have to take some time to re-learn and beyond.

In the meantime, I wanted to respond and thank you for giving this advice and providing a way to try and move forward!

kbk