ode_monodomain_2d.f90

Go to the documentation of this file.

!!\f$ ~~~~~~~~~~~~~~~~\displaystyle{

program ode_monodomain_2d

  use choral_constants
  use choral

  implicit none

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

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

  !! MODEL EESCIPTION
  !!
  !! IONIC MODEL DEF.
  !!   im_type = type of ionic model
  !!   im      = definition of the ionic model 
  !!
  integer, parameter :: im_type = ionic_br  
  type(ionicmodel)   :: im
  !!
  !! CARDIAC TISSUE MODEL DEF.
  !!   cd_type = type of conductivities
  !!   am      = cell surface to volume ratio
  !!   cm      = definition of the cardiac tissue model
  !!
  integer , parameter :: cd_type = le_guyader
  real(RP), parameter :: am      = 500.0_rp
  type(cardiomodel)   :: cm

  !!
  !! STIMULATION DEF.
  !!
  !!   stim_time         = stimulation mid-time
  !!   stim_time_radius  = stimulation mid-duration
  !!   stim_space_radius = radius of the stimulation area
  !!
  type(stimulation)   :: stim
  real(RP), parameter :: stim_time         =  3._rp
  real(RP), parameter :: stim_time_radius  =  1.2_rp
  real(RP), parameter :: stim_space_radius =  0.15_rp


  !! TIME DISCRETISATION
  !!
  !!   t0      = initial time
  !!   T       = final time
  !!   dt      = time step
  !!   SL_meth = method for the semilinear eq.
  !!   NL_meth = method for the non-ilinear system
  !!
  real(RP), parameter :: t0      = 0.00_rp
  real(RP), parameter :: T       = 30.0_rp
  real(RP), parameter :: dt      = 0.05_rp
  integer , parameter :: SL_meth = ode_bdfsbdf3
  integer , parameter :: NL_meth = ode_rl3
  !!
  !!   pb      = definition of the ode problem
  !!   slv     = definition of the ode solver
  !!   sol     = data structure for the ode solution
  !!   out     = definition of the output
  !!
  type(ode_solution):: sol
  type(ode_problem) :: pb 
  type(ode_solver)  :: slv
  type(ode_output)  :: out


  !! SPACE DISCRETISATION
  !!
  !!   fe_type = finite element method
  !!   qd_type = quadrature method    
  !!
  integer, parameter :: fe_type = fe_p3_2d
  integer, parameter :: qd_type = quad_gauss_trg_12
  !! msh      = mesh
  !! X_h      = finite element space
  !! M, S     = mass and stiffness matrices
  !!
  type(mesh)    :: msh
  type(fespace) :: X_h
  type(csr)     :: M, S

  
  !! LINEAR SYSTEM AND PRECONDITIONING
  !!
  !!   pc_type = preconditionner type
  !!   kry     = krylov method def.
  !!   K       = linear system matrix: K = M + Cs*S
  !!   Cs      = stiffness matrix prefactor
  !!   pc      = preconditionner for 'Kx = y'
  !!
  integer, parameter :: pc_type = pc_jacobi
  real(RP)     :: Cs
  type(krylov) :: kry
  type(csr)    :: K  
  type(precond):: pc


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


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


  !!!!!!!!!!!!!!!!!!!!!!!  MODEL DEFINITION
  !!
  write(*,*) ""
  write(*,*) "==================================  MODEL DEFINITION"
  !!
  !!  Cardiac tissue model
  !!
  cm = cardiomodel(vector_field_e_x, am, le_guyader)
  if (verb>1) call print(cm)
  !!
  !!  ionic model
  !!
  im = ionicmodel(im_type)
  if (verb>1) call print(im)
  !!
  !! stimulation
  !!
  stim = stimulation(f_c3, l=im%Ist*1.2_rp, t0=stim_time, &
       & rt = stim_time_radius, rx = stim_space_radius)


  !!!!!!!!!!!!!!!!!!!!!!!  SPACE DISCRETISATION
  !!
  write(*,*) ""
  write(*,*) "==================================  SPACE DISCRETISATION"
  !!
  !! MESH AND FINITE ELEMENT MESH
  !!
  msh = mesh(trim(gmsh_dir)//"square/square_1.msh", 'gmsh')
  !!  
  x_h = fespace(msh)
  call set(x_h, fe_type)
  call assemble(x_h)
  !!
  !! MONODOMAIN MASS AND STIFFNESS MATRICES
  !!
  call monodomain_assemble(m, s, cm, x_h, qd_type, qd_type)
  !!
  !! SYSTEM SYSTEM AND PRECONDITIONING
  !!
  cs = s_prefactor(sl_meth, dt)
  call add(k, m, 1._rp, s, cs)
  pc  = precond(k, pc_type)
  kry = krylov(kry_cg, tol=1e-8_rp, itmax=100)



  !! !!!!!!!!!!!!!!!!!!!!! TIME DISCRETISATION
  !!
  write(*,*) ""
  write(*,*) "==================================  TIME DISCRETISATION"
  !!
  !! ODE PROBLEM DEF.
  !!
  pb = ode_problem(ode_pb_sl_nl, &
       &   dof = x_h%nbDof, x=x_h%dofCoord, &
       &   m=massmat, s=stiffmat, ab=pde_reaction, n=im%N, na=im%Na) 
  if (verb>1) call print(pb)
  !!
  !! ODE SOLVER DEF.
  !!
  slv = ode_solver(pb, ode_slv_ms, &
       & sl_meth=sl_meth, nl_meth=nl_meth)
  if (verb>1) call print(slv)
  !!
  !! ODE SOLUTION DEF.
  !!
  sol = ode_solution(slv, pb)
  if (verb>1) call print(sol)
  !!
  !! ODE OUTPUT DEF.
  !!
  out = ode_output(t0, t, im%N)
  call set(out, verb=2)
  call set(out, vtn_rec_prd = 1._rp, vtn_plot=.true.)
  call set(out, act_type=act_4)
  call set(out, act_rec=.true., act_plot=.true.)
  call set(out, pos=pos_gmsh)
  if (verb>1) call print(out)


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

  write(*,*) "monodomain_2d   : end"

contains


  !! !!!!!!!!!!!!!!!!!!!!!  REACTION TERM
  !!
  subroutine pde_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

    !! source term (stimulation current)
    !!
    i_app = stim%I_app(x, t)

    !! total reaction term
    !!
    call im%AB(a, b, i_app, y, n, na)

  end subroutine pde_reaction



  !! !!!!!!!!!!!!!!!!!!!!!  DIFFUSION TERMS
  !!
  !!  Mass matrix product
  !!
  subroutine massmat(yy, xx)
    real(RP), dimension(:), intent(out) :: yy
    real(RP), dimension(:), intent(in)  :: xx

    call matvecprod(yy, m, xx)

  end subroutine massmat
  !!
  !!  Stiffness matrix product
  !!
  subroutine stiffmat(yy, xx)
    real(RP), dimension(:), intent(out) :: yy
    real(RP), dimension(:), intent(in)  :: xx
  
    call matvecprod(yy, s, xx)

  end subroutine stiffmat


  !! !!!!!!!!!!!!!!!!!!!!!  LINEAR SYSTEM INVERSION
  !!
  !!  x = K**{-1}b
  !!
  subroutine kinv(xx, ierr, bb)
    real(RP), dimension(:), intent(inout) :: xx
    logical               , intent(out)   :: ierr
    real(RP), dimension(:), intent(in)    :: bb

    call solve(xx, kry, bb, k, pc)

    ierr =  kry%ierr
    
  end subroutine kinv

end program ode_monodomain_2d