LLNL / sundials

Official development repository for SUNDIALS - a SUite of Nonlinear and DIfferential/ALgebraic equation Solvers. Pull requests are welcome for bug fixes and minor changes.
https://computing.llnl.gov/projects/sundials
BSD 3-Clause "New" or "Revised" License
515 stars 125 forks source link

Run error: Process killed. The memory footprint is getting higher and higher #78

Closed yushimin13 closed 1 year ago

yushimin13 commented 3 years ago

Hello everyone, I called the latest version of SUNDIALS CVode in the program, during the running of my program, the process took up more and more memory, until it was killed. As you can see below. Cache_-483ecab8900799e7

Because in the system of differential equations in some parameters may change over time, so every time the node needs to be updated on the equation, so that each run a time step will need to update equations, and to call the interface, don't know is it because the mentor of memory, and if so, would you please tell me what the good under the solution? Thank you.

The code can see below. `! 三阶的可以修改输入参数的

module ode_mod

!======= Inclusions ===========
 use, intrinsic :: iso_c_binding   ! Fortran 调用 C 函数

 use fcvode_mod                    ! Fortran interface to CVODE
 use fsundials_nvector_mod         ! Fortran interface to generic N_Vector
 use fnvector_serial_mod           ! Fortran interface to serial N_Vector
 use fsundials_nonlinearsolver_mod ! Fortran interface to generic SUNNonlinearSolver
 use fsunnonlinsol_fixedpoint_mod  ! Fortran interface to fixed point SUNNonlinearSolver
 Use ModuleControlFlow, only : Inputdt
!======= Declarations =========
implicit none

! number of equations
integer(c_long), parameter :: neq = 4

! ODE parameters
double precision, parameter :: lamda = -100.0d0

! ODE parameters
double precision, parameter :: PI=3.1415926
double precision, parameter :: Rs=50.d0
double precision, parameter :: C1=150d-12,Cb=200.d-12,Cc=1.390814556207367d-11
double precision, parameter :: L=4.3d-6
double precision, parameter :: Rm=0.5d0
double precision, parameter :: V0=50.d0
double precision, parameter :: f=13.56d6

double precision, parameter :: Radius=0.1d0
double precision, parameter :: A=PI*Radius*Radius
double precision, parameter :: Phi_nc=0.0d0

real(c_double) :: Vt=0.d0
real(c_double) :: Phi_n0=0.d0
real(c_double) :: J_conv=0.d0  

contains

! ----------------------------------------------------------------
! RhsFn provides the right hand side function for the
! ODE: dy/dt = f(t,y)
!
! Return values:
!    0 = success,
!    1 = recoverable error,
!   -1 = non-recoverable error
! ----------------------------------------------------------------
integer(c_int) function RhsFn(tn, sunvec_y, sunvec_f, user_data) &
     result(ierr2) bind(C,name='RhsFn')

  !======= Inclusions ===========
  !use, intrinsic :: iso_c_binding
  !use fsundials_nvector_mod

  !======= Declarations =========
  implicit none

  ! calling variables
  real(c_double), value :: tn        ! current time
  type(N_Vector)        :: sunvec_y  ! solution N_Vector
  type(N_Vector)        :: sunvec_f  ! rhs N_Vector
  type(c_ptr), value    :: user_data ! user-defined data

  ! pointers to data in SUNDIALS vectors
  real(c_double), pointer :: yvec(:)
  real(c_double), pointer :: fvec(:)

  ! ODE system matrix
  real(c_double) :: Amat(neq,neq),Bmat(1,neq)

  !======= Internals ============

  ! get data arrays from SUNDIALS vectors
  yvec => FN_VGetArrayPointer(sunvec_y)
  fvec => FN_VGetArrayPointer(sunvec_f)

  fvec(1) = -yvec(1)/(Rs*C1)+yvec(2)/(Rs*C1)+(Vt/Rs)
  fvec(2) = yvec(3)
  fvec(3) = (yvec(1)-yvec(2))/(L*C1)-yvec(2)/(L*Cb)-(Rm*yvec(3))/L-Phi_n0/(L)
  fvec(4) = yvec(3)/A + J_conv

  ! return success
  ierr2 = 0
  return

end function RhsFn

subroutine odeout(yvec,Gloab)

    implicit none

    integer(c_int) x;
    ! solution vector, neq is set in the ode_mod module
    !real(c_double) :: yvec(neq)

    ! local variables
    real(c_double) :: tstart     ! initial time
    real(c_double) :: tend       ! final time
    real(c_double) :: rtol, atol ! relative and absolute tolerance
    real(c_double) :: dtout      ! output time interval
    real(c_double) :: tout       ! output time
    real(c_double) :: tcur(1)    ! current time
    integer(c_int) :: ierr       ! error flag from C functions
    integer(c_int) :: nout       ! number of outputs
    integer(c_int) :: outstep    ! output loop counter
    type(c_ptr)    :: cvode_mem  ! CVODE memory

    type(N_Vector), pointer           :: sunvec_y ! sundials vector
    type(SUNNonlinearSolver), pointer :: sunnls   ! sundials fixed-point nonlinear solver

    real(c_double),intent(inout) :: yvec(neq)
    real(c_double),intent(inout) :: Gloab(neq)

    !======= Internals ============

    ! initialize ODE
    tstart = 0.0d0
    tend   = Inputdt
    tcur   = tstart
    tout   = tstart
    dtout  = Inputdt
    nout   = ceiling(tend/dtout)

  ! initialize solution vector
    !yvec(1:neq) = 0

    ! create SUNDIALS N_Vector
    sunvec_y => FN_VMake_Serial(neq, yvec)
    if (.not. associated(sunvec_y)) then
       print *, 'ERROR: sunvec = NULL'
       stop 1
    end if

    ! create CVode memory
    cvode_mem = FCVodeCreate(CV_ADAMS)
    if (.not. c_associated(cvode_mem)) then
       print *, 'ERROR: cvode_mem = NULL'
       stop 1
    end if

    ! initialize CVode
    ierr = FCVodeInit(cvode_mem, c_funloc(RhsFn), tstart, sunvec_y)
    if (ierr /= 0) then
       print *, 'Error in FCVodeInit, ierr = ', ierr, '; halting'
       stop 1
    end if

    ! set relative and absolute tolerances
    rtol = 1.0d-6
    atol = 1.0d-10

    ierr = FCVodeSStolerances(cvode_mem, rtol, atol)
    if (ierr /= 0) then
       print *, 'Error in FCVodeSStolerances, ierr = ', ierr, '; halting'
       stop 1
    end if

    ! create fixed point nonlinear solver object
    sunnls => FSUNNonlinSol_FixedPoint(sunvec_y, 0)
    if (.not. associated(sunnls)) then
       print *,'ERROR: sunnls = NULL'
       stop 1
    end if

    ! attache nonlinear solver object to CVode
    ierr = FCVodeSetNonlinearSolver(cvode_mem, sunnls)
    if (ierr /= 0) then
      print *, 'Error in FCVodeSetNonlinearSolver, ierr = ', ierr, '; halting'
      stop 1
    end if

! start time stepping
    ! print *, '   '
    ! print *, 'Finished initialization, starting time steps'
    ! print *, '   '
    ! print *, '       t           y1           y2           y3       '
    ! print *, '------------------------------------------------------'
    ! open(10,position='append',file='solution6.dat')
    ! write(10,'(2x,7(es12.5,1x))') tcur, yvec(1:neq)
    ! print '(2x,4(es12.5,1x))',  tcur,yvec(1), yvec(2), yvec(3)
    ! close(10)
    do outstep = 1,nout

       ! call CVode
       tout = min(tout + dtout, tend)
       ierr = FCVode(cvode_mem, tout, sunvec_y, tcur, CV_NORMAL)
       if (ierr /= 0) then
          print *, 'Error in FCVODE, ierr = ', ierr, '; halting'
          stop 1
       endif

       ! output current solution
    !    print '(2x,4(es12.5,1x))', tcur, yvec(1), yvec(2), yvec(3)
       Gloab(1:neq) = yvec(1:neq)

    enddo

    ! clean up
    !call FCVodeFree(cvode_mem)
    !ierr = FSUNNonLinSolFree(sunnls)
    !call FN_VDestroy(sunvec_y)

    ! return success
    ierr = 0

end subroutine odeout

end module ode_mod ` The above interface is called in another module.

` real(c_double),save :: y(1,neq) ! real(c_double),save :: Gloab(neq) !

 Iconv = EC%Qconv(1)/CF%Dt

    Vt  = FB%V1
        J_conv=Iconv/EC%A0
        Phi_n0=FG%Phi(1)

        if(CF%Timer==737000) then
           y = reshape([&
           8.42082d-09,  6.30835d-09,  4.20469d-01, -9.33230d-08], &
           [1,4])
        end if

        open(10,position='append',file='input.dat')
            write(10,FMt="(*(es21.14,1x))") CF%Timer, Vt, Iconv, J_conv, Phi_n0
         close(10)

     call odeout(y,Gloab)

         y = reshape([&
           Gloab(1), Gloab(2), Gloab(3), Gloab(4)], &
           [1,4])

        open(10,position='append',file='solution.dat')
           write(10,'(2x,7(es12.5,1x))') CF%Timer, Gloab(1:neq)
        close(10)

            EC%Qc(1) = Gloab(2)  !带返回电荷量Q
            EC%Sigma(1) = Gloab(4)  !带返回的西格玛

`

My running environment is WSL+VSCODE+oneAPI. Compile with Cmake and run the program under LINUX.

The complete code is shown below, where an Killed occurs when the run time NRUN is changed to large(For example, NRUN=5000) code.zip

balos1 commented 1 year ago

Closing as stale. Open a new issue if this is still a problem for you.