CODEWITHSUNDEEP
fortransolverode
Adriana Huerta
Adriana Huerta

Reputation: 23

ODE solver in F90

I decided to write my own solver (RK4th order) in F90 for a NPZD box model. It seems to run fine for the first 120 time steps (timestep = 1 day) and then it grows unrealistically. I reduced the time step and have the same problem. At this point I'm wondering how can I make stable the solver. Or if a better solution is to use one of the ODEPACK, which I'm not sure how to properly use. I would like to keep the timestep = 1 day, and also would like to learn how to make it stable without reducing the time step.
I am actually rewriting in Fortran a code done already in python (for technical reasons we need it in Fortran now) my college runs it without problems, but uses the available ODEs from python. In this same forum someone suggested (to someone else) to write her own solver in fortran, therefore I thought I can do that as well.

Below are the module of the solver and the main program: 

subroutine solver(stateSim1E)
 type(ModelState) :: stateSim1E,stateSim1E1
type(ModelState) :: stateSim1E2,stateSim1E3,stateSim1E4

stateSim1E1 = dModelState(stateSim1E) !k1
!!$
stateSim1E2 = stateSim1E + 0.5*DTS*stateSim1E1 !cc1 K2
stateSim1E2 = dModelState(stateSim1E2)     

stateSim1E3 = stateSim1E + 0.5*DTS*(stateSim1E2)!k3
stateSim1E3 = dModelState(stateSim1E3) 

stateSim1E4 = stateSim1E + DTS*(stateSim1E3)!k4
stateSim1E4 = dModelState(stateSim1E4) 

stateSim1E = stateSim1E + (DTS/6)*(stateSim1E1 + 2._RK*(stateSim1E2 + stateSim1E3) +stateSim1E4)

end subroutine solver

main program

program main_npzdsv
use iso_fortran_env
use NumericKinds
use ModelState_class_npzdsv  
  integer :: i, outFileID,ti
  type(ModelState) :: stateSim1E, statePerturbation !,stateSim1E1

   !initial conditions
   stateSim1E = ModelState(N=N0,P=P0,Z=Z0,D=D0,S = s_mean0, V = s_var0)     


ti = 0 
do i = 1, NO_T_STEP
   ti = i     
   call getindex(ti)
   call solver(stateSim1E)
   print*,'new', index, stateSim1E

   if( mod(i-1, OUTPUT_FREQUENCY) == 0 ) then
      write(outFileID,'(i4)')index
      write(outFileID, '(18('// RK_FMT // ', 1x))') stateSim1E%getCurrModelState()
 end if
 end do
 close(outFileID) ! Clean-up for output
 close(outFileID2) ! Clean-up for output


end program main_npzdsv

!

Upvotes: 1

Views: 877

Answers (1)

innoSPG
innoSPG

Reputation: 4656

The first think to check is the CFL. Explicit methods sometime require very small time step as opposed to implicit methods that are sometime unconditionally stable. Checking the CFL can actually be tricky for problems where the velocity evolves in time. It is very challenging for spinup problem. By spinup, I mean when you want to run a long time simulation from very small perturbation to get to a realistic state. For example simulating the water circulation in a gulf from rest. The initial transient can generate very high velocity that makes it very difficult to satisfy the CFL and break the system.

If you absolutely want your time step of at least one day, you will possibly need to increase your step in space. Another simple approach is to read the doc of the toolbox used in python to get information about the underlying method.

(This is not the solution, it is an indication that is too long to be put as comment)

Upvotes: 1

Related Questions