ode_cardio_0d.f90

Go to the documentation of this file.

program ode_cardio_0d

  use choral_constants
  use choral

  implicit none

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!      VARAIBLE DEFINITION : BEGIN
  !!
  !!  

  !! verbosity level
  integer, parameter :: verb = 3

  !!         MODEL EESCIPTION
  !!
  !! IONIC MODEL DEF.
  !!   im_type = type of ionic model
  !!   im      = definition of the ionic model 
  !!
  integer, parameter :: im_type = ionic_tnnp
  type(ionicmodel)   :: im
  !!
  !! STIMULATION DEF.
  !!   stim_base = base function for the stimulation
  !!   stim_time = mid-time for the stimulation
  !!   stim_time_radius  = mid stimulation duration
  !!
  procedure(rtor), pointer   :: stim_base         => f_c3
  real(RP)       , parameter :: stim_time         =  20._rp
  real(RP)       , parameter :: stim_time_radius  =  1._rp
  !!
  !! OUTPUT DEF.
  !!   co = definition of the output
  type(ode_output) :: co


  !!         TIME DISCRETISATION
  !!
  !! t0      = initial time
  !! T       = final time
  !! dt      = time step
  !! slv_type = solver type (multistep, onestep)
  !! NL_meth  = method for the non-ilinear system
  !!
  real(RP), parameter :: t0       = 0.0_rp
  real(RP), parameter :: T        = 400.0_rp
  !! first choice of method (stable)
  ! integer , parameter :: slv_type = ODE_SLV_MS
  ! integer , parameter :: NL_meth  = ODE_RL3
  ! real(RP), parameter :: dt       = .05_RP
  !! second choice of method (less stable)
  integer , parameter :: slv_type = ode_slv_1s
  integer , parameter :: NL_meth  = ode_modif_erk2_b
  real(RP), parameter :: dt       = 0.4_rp
  !!
  !! pb      = definition of the ode problem
  !! slv     = definition of the ode solver
  !! sol     = data structure for the ode solution
  !!
  type(ode_solution) :: sol
  type(ode_problem)  :: pb 
  type(ode_solver)   :: slv


  !!
  !!   VARAIBLE DEFINITION : END
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  call choral_init(verb=verb)
  write(*,*) "ode_cardio_0D       : start"


  !!!!!!!!!!!!!!!!!!!!!!!  MODEL DEFINITION
  !!
  write(*,*) ""
  write(*,*) "==================================  MODEL DEFINITION"
  !!
  !!  ionic model
  !!
  im = ionicmodel(im_type)
  if (verb>0) call print(im)


  !! !!!!!!!!!!!!!!!!!!!!! TIME DISCRETISATION
  !!
  write(*,*) ""
  write(*,*) "==================================  TIME DISCRETISATION"
  !!
  !! ODE PROBLEM DEF.
  !!
  pb = ode_problem(ode_pb_nl, &
       &   dim=0, ab=ode_reaction, n=im%N, na=im%Na) 
  if (verb>0) call print(pb)
  !!
  !! ODE SOLVER DEF.
  !!
  slv = ode_solver(pb, slv_type, nl_meth=nl_meth, verb=1)
  if (verb>0) call print(slv)
  !!
  !! ODE SOLUTION DEF.
  !!
  sol = ode_solution(slv, pb)
  if (verb>0) call print(sol)
  !!
  !! ODE OUTPUT DEF.
  !!
  co = ode_output(t0, t, im%N)
  call set(co, verb=1)
  call set(co, vxn_period = 1._rp, vxn_plot=.true.)
  if (verb>0) call print(co)


  !!!!!!!!!!!!!!!!!!!!!!!  NUMERICAL RESOLUTION
  !!
  write(*,*) ""
  write(*,*) "==================================  NUMERICAL RESOLUTION"
  !!
  !! finalise the output definition
  call assemble(co, dt)
  !!
  !! initial condition
  call initialcond(sol, pb, slv, t0, im%y0)
  !!
  !! numerical resolution
  call solve(sol, slv, pb, t0, t, dt, output=co)

  write(*,*) "ode_cardio_0D       : end"

contains


  !! !!!!!!!!!!!!!!!!!!!!!  SOURCE TERM
  !!
  function stim(x, t) result(res)
    real(RP)                           :: res
    real(RP), dimension(3), intent(in) :: x
    real(RP)              , intent(in) :: t

    real(RP) :: t2

    t2 = abs( t - stim_time)/stim_time_radius

    res = stim_base(t2) * im%Ist

  end function stim



  !! !!!!!!!!!!!!!!!!!!!!!  REACTION TERMS
  !!
  subroutine ode_reaction(a, b, x, t, y, N, Na)
    real(RP), dimension(Na), intent(out) :: a
    real(RP), dimension(N) , intent(out) :: b
    real(RP), dimension(3) , intent(in)  :: x
    real(RP)               , intent(in)  :: t
    real(RP), dimension(N) , intent(in)  :: y
    integer                , intent(in)  :: N, Na

    real(RP) :: I_app

    i_app = stim(x, t)
    call im%AB(a, b, i_app, y, n, na)

  end subroutine ode_reaction

end program ode_cardio_0d