Choral
ode_monodomain_2d.f90
Go to the documentation of this file.
1 !!\f$ ~~~~~~~~~~~~~~~~\displaystyle{
42 
44 
46  use choral
47 
48  implicit none
49 
50  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
51  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
52  !!
53  !! VARAIBLE DEFINITION : BEGIN
54  !!
55  !!
56 
57  !! verbosity level
58  integer, parameter :: verb = 2
59 
60  !! MODEL EESCIPTION
61  !!
62  !! IONIC MODEL DEF.
63  !! im_type = type of ionic model
64  !! im = definition of the ionic model
65  !!
66  integer, parameter :: im_type = ionic_br
67  type(ionicmodel) :: im
68  !!
69  !! CARDIAC TISSUE MODEL DEF.
70  !! cd_type = type of conductivities
71  !! am = cell surface to volume ratio
72  !! cm = definition of the cardiac tissue model
73  !!
74  integer , parameter :: cd_type = le_guyader
75  real(RP), parameter :: am = 500.0_rp
76  type(cardiomodel) :: cm
77 
78  !!
79  !! STIMULATION DEF.
80  !!
81  !! stim_time = stimulation mid-time
82  !! stim_time_radius = stimulation mid-duration
83  !! stim_space_radius = radius of the stimulation area
84  !!
85  type(stimulation) :: stim
86  real(RP), parameter :: stim_time = 3._rp
87  real(RP), parameter :: stim_time_radius = 1.2_rp
88  real(RP), parameter :: stim_space_radius = 0.15_rp
89 
90 
91  !! TIME DISCRETISATION
92  !!
93  !! t0 = initial time
94  !! T = final time
95  !! dt = time step
96  !! SL_meth = method for the semilinear eq.
97  !! NL_meth = method for the non-ilinear system
98  !!
99  real(RP), parameter :: t0 = 0.00_rp
100  real(RP), parameter :: T = 30.0_rp
101  real(RP), parameter :: dt = 0.05_rp
102  integer , parameter :: SL_meth = ode_bdfsbdf3
103  integer , parameter :: NL_meth = ode_rl3
104  !!
105  !! pb = definition of the ode problem
106  !! slv = definition of the ode solver
107  !! sol = data structure for the ode solution
108  !! out = definition of the output
109  !!
110  type(ode_solution):: sol
111  type(ode_problem) :: pb
112  type(ode_solver) :: slv
113  type(ode_output) :: out
114 
115 
116  !! SPACE DISCRETISATION
117  !!
118  !! fe_type = finite element method
119  !! qd_type = quadrature method
120  !!
121  integer, parameter :: fe_type = fe_p3_2d
122  integer, parameter :: qd_type = quad_gauss_trg_12
123  !! msh = mesh
124  !! X_h = finite element space
125  !! M, S = mass and stiffness matrices
126  !!
127  type(mesh) :: msh
128  type(fespace) :: X_h
129  type(csr) :: M, S
130 
131 
132  !! LINEAR SYSTEM AND PRECONDITIONING
133  !!
134  !! pc_type = preconditionner type
135  !! kry = krylov method def.
136  !! K = linear system matrix: K = M + Cs*S
137  !! Cs = stiffness matrix prefactor
138  !! pc = preconditionner for 'Kx = y'
139  !!
140  integer, parameter :: pc_type = pc_jacobi
141  real(RP) :: Cs
142  type(krylov) :: kry
143  type(csr) :: K
144  type(precond):: pc
145 
146 
147  !!
148  !! VARAIBLE DEFINITION : END
149  !!
150  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
152 
153 
154  call choral_init(verb=verb)
155  write(*,*) "ode_monodomain_2d : start"
156 
157 
158  !!!!!!!!!!!!!!!!!!!!!!! MODEL DEFINITION
159  !!
160  write(*,*) ""
161  write(*,*) "================================== MODEL DEFINITION"
162  !!
163  !! Cardiac tissue model
164  !!
165  cm = cardiomodel(vector_field_e_x, am, le_guyader)
166  if (verb>1) call print(cm)
167  !!
168  !! ionic model
169  !!
170  im = ionicmodel(im_type)
171  if (verb>1) call print(im)
172  !!
173  !! stimulation
174  !!
175  stim = stimulation(f_c3, l=im%Ist*1.2_rp, t0=stim_time, &
176  & rt = stim_time_radius, rx = stim_space_radius)
177 
178 
179  !!!!!!!!!!!!!!!!!!!!!!! SPACE DISCRETISATION
180  !!
181  write(*,*) ""
182  write(*,*) "================================== SPACE DISCRETISATION"
183  !!
184  !! MESH AND FINITE ELEMENT MESH
185  !!
186  msh = mesh(trim(gmsh_dir)//"square/square_1.msh", 'gmsh')
187  !!
188  x_h = fespace(msh)
189  call set(x_h, fe_type)
190  call assemble(x_h)
191  !!
192  !! MONODOMAIN MASS AND STIFFNESS MATRICES
193  !!
194  call monodomain_assemble(m, s, cm, x_h, qd_type, qd_type)
195  !!
196  !! SYSTEM SYSTEM AND PRECONDITIONING
197  !!
198  cs = s_prefactor(sl_meth, dt)
199  call add(k, m, 1._rp, s, cs)
200  pc = precond(k, pc_type)
201  kry = krylov(kry_cg, tol=1e-8_rp, itmax=100)
202 
203 
204 
205  !! !!!!!!!!!!!!!!!!!!!!! TIME DISCRETISATION
206  !!
207  write(*,*) ""
208  write(*,*) "================================== TIME DISCRETISATION"
209  !!
210  !! ODE PROBLEM DEF.
211  !!
212  pb = ode_problem(ode_pb_sl_nl, &
213  & dof = x_h%nbDof, x=x_h%dofCoord, &
214  & m=massmat, s=stiffmat, ab=pde_reaction, n=im%N, na=im%Na)
215  if (verb>1) call print(pb)
216  !!
217  !! ODE SOLVER DEF.
218  !!
219  slv = ode_solver(pb, ode_slv_ms, &
220  & sl_meth=sl_meth, nl_meth=nl_meth)
221  if (verb>1) call print(slv)
222  !!
223  !! ODE SOLUTION DEF.
224  !!
225  sol = ode_solution(slv, pb)
226  if (verb>1) call print(sol)
227  !!
228  !! ODE OUTPUT DEF.
229  !!
230  out = ode_output(t0, t, im%N)
231  call set(out, verb=2)
232  call set(out, vtn_rec_prd = 1._rp, vtn_plot=.true.)
233  call set(out, act_type=act_4)
234  call set(out, act_rec=.true., act_plot=.true.)
235  call set(out, pos=pos_gmsh)
236  if (verb>1) call print(out)
237 
238 
239  !!!!!!!!!!!!!!!!!!!!!!! NUMERICAL RESOLUTION
240  !!
241  write(*,*) ""
242  write(*,*) "================================== NUMERICAL RESOLUTION"
243  !!
244  !! finalise the output definition
245  call assemble(out, dt, x_h)
246  !!
247  !! initial condition
248  call initialcond(sol, pb, slv, t0, im%y0)
249  !!
250  !! numerical resolution
251  call solve(sol, slv, pb, t0, t, dt, kinv, output=out)
252 
253  write(*,*) "monodomain_2d : end"
254 
255 contains
256 
257 
258  !! !!!!!!!!!!!!!!!!!!!!! REACTION TERM
259  !!
260  subroutine pde_reaction(a, b, x, t, y, N, Na)
261  real(RP), dimension(Na), intent(out) :: a
262  real(RP), dimension(N) , intent(out) :: b
263  real(RP), dimension(3) , intent(in) :: x
264  real(RP) , intent(in) :: t
265  real(RP), dimension(N) , intent(in) :: y
266  integer , intent(in) :: N, Na
267 
268  real(RP) :: I_app
269 
270  !! source term (stimulation current)
271  !!
272  i_app = stim%I_app(x, t)
273 
274  !! total reaction term
275  !!
276  call im%AB(a, b, i_app, y, n, na)
277 
278  end subroutine pde_reaction
279 
280 
281 
282  !! !!!!!!!!!!!!!!!!!!!!! DIFFUSION TERMS
283  !!
284  !! Mass matrix product
285  !!
286  subroutine massmat(yy, xx)
287  real(RP), dimension(:), intent(out) :: yy
288  real(RP), dimension(:), intent(in) :: xx
289 
290  call matvecprod(yy, m, xx)
291 
292  end subroutine massmat
293  !!
294  !! Stiffness matrix product
295  !!
296  subroutine stiffmat(yy, xx)
297  real(RP), dimension(:), intent(out) :: yy
298  real(RP), dimension(:), intent(in) :: xx
299 
300  call matvecprod(yy, s, xx)
301 
302  end subroutine stiffmat
303 
304 
305  !! !!!!!!!!!!!!!!!!!!!!! LINEAR SYSTEM INVERSION
306  !!
307  !! x = K**{-1}b
308  !!
309  subroutine kinv(xx, ierr, bb)
310  real(RP), dimension(:), intent(inout) :: xx
311  logical , intent(out) :: ierr
312  real(RP), dimension(:), intent(in) :: bb
313 
314  call solve(xx, kry, bb, k, pc)
315 
316  ierr = kry%ierr
317 
318  end subroutine kinv
319 
320 end program ode_monodomain_2d
subroutine stiffmat(yy, xx)
integer, parameter ode_bdfsbdf3
BDF / SBDF 3.
subroutine, public choral_init(verb)
Initialisation for the CHORAL library
Definition: choral.F90:166
TOP-LEVEL MODULE FOR THE LIBRARY CHORAL
Definition: choral.F90:27
subroutine pde_reaction(a, b, x, t, y, N, Na)
program ode_monodomain_2d
RESOLUTION OF: the monodomain model in cardiac-electrophysiology.
integer, parameter ionic_br
integer, parameter ode_slv_ms
multistep
subroutine kinv(x, bool, b)
Solver for K*x = b.
real(rp) function e(x, v1, v2)
CHORAL CONSTANTS
integer, parameter pos_gmsh
integer, parameter pc_jacobi
Jacobi diagonal preconditioning.
subroutine massmat(yy, xx)
integer, parameter ode_pb_sl_nl
SemiLinear ODE coupled with a non-linear ODE system for with .
integer, parameter le_guyader
integer, parameter ode_rl3
Rush Larsen 3.
integer, parameter na
integer, parameter fe_p3_2d
integer, parameter act_4
bi-quadratic interpolation
integer, parameter quad_gauss_trg_12
integer, parameter kry_cg
CG linear solver.