Choral
ode_output_mod.F90
Go to the documentation of this file.
1 
44 
45 
47 
50  use choral_env
51  use real_type
52  use basic_tools
53  use abstract_interfaces, only: r3tor
54  use io
55  use r1d_mod, only: copy
56  use r2d_mod, only: copy
57  use fespace_mod
58  use integral
59  use quadmesh_mod
61 
62  implicit none
63  private
64 
65 
66  ! %----------------------------------------%
67  ! | |
68  ! | PUBLIC DATA |
69  ! | |
70  ! %----------------------------------------%
71  public :: ode_output
72  public :: clear, set, assemble, print
73  public :: output_proc
74 
75  !! Wave-front analysis
77 
78 
79  ! %----------------------------------------%
80  ! | |
81  ! | DERIVED TYPE |
82  ! | |
83  ! %----------------------------------------%
87  type :: ode_output
88 
90  logical :: created = .false.
91 
93  logical :: assembled = .false.
94 
95  !! Constructor parameter: provided by the user with
96  !! 'out = ode-output(...)'
97  !!
99  integer :: dim
101  integer :: n
103  real(RP) :: t0
105  real(RP) :: t
106 
107  !! Parameters provided by the user with
108  !! 'call set(out, ...)'
109  !!
111  integer :: verb
113  real(RP) :: ytn_period
115  real(RP) :: vtn_period
117  real(RP) :: vtn_rec_prd
119  logical :: vtn_plot
121  integer, dimension(:), allocatable :: xn_list
123  real(RP) :: vxn_period
125  character(len=75) :: vxn_file
127  logical :: vxn_plot
129  integer :: act_type
131  real(RP) :: act_vth
133  logical :: act_plot
135  logical :: act_rec
137  logical :: stop_if_dep
139  integer :: pos
141  character(len=150) :: outdir
143  character(len=50) :: vtn_file
145  character(len=50) :: act_file
146 
147 
148  !! Parameters provided by the user with
149  !! 'call assemble(out, ...)'
150  !!
152  real(RP) :: dt
154  type(fespace), pointer :: x_h=>null()
155 
156 
157  !! The following data are initialised by the code
158  !! when assembling
159  !!
161  real(RP), dimension(:,:,:), allocatable :: ytn
163  real(RP), dimension(:,:) , allocatable :: vtn
165  real(RP), dimension(:,:) , allocatable :: vxn
169  real(RP), dimension(:,:) , allocatable :: v
171  real(RP), dimension(:) , allocatable :: act
173  integer :: nbdof
177  integer :: nv
179  integer :: ytn_cpt
181  integer :: vtn_cpt
183  integer :: vtn_rec_cpt
185  integer :: vxn_cpt
186 
187  contains
188 
191 
192  end type ode_output
193 
194 
195  ! %----------------------------------------%
196  ! | |
197  ! | GENERIC INTERFACES |
198  ! | |
199  ! %----------------------------------------%
200  interface clear
201  module procedure ode_output_clear
202  end interface clear
203 
204  interface ode_output
205  module procedure ode_output_create
206  end interface ode_output
207 
209  interface set
210  module procedure ode_output_set
211  end interface set
212 
214  interface assemble
215  module procedure ode_output_assemble
216  end interface assemble
217 
219  interface print
220  module procedure ode_output_print
221  end interface print
222 
223 contains
224 
226  subroutine ode_output_clear(out)
227  type(ode_output), intent(inout) :: out
228 
229  out%created = .false.
230  out%assembled = .false.
231 
232  call freemem(out%xn_list)
233 
234  call freemem(out%Ytn )
235  call freemem(out%Vtn )
236  call freemem(out%Vxn )
237  call freemem(out%V )
238  call freemem(out%act )
239 
240  out%X_h=>null()
241 
242  end subroutine ode_output_clear
243 
244 
255  function ode_output_create(t0, T, N) result(out)
256  type(ode_output) :: out
257  real(RP), intent(in) :: t0, T
258  integer , intent(in) :: N
259 
260  call clear(out)
261  out%t0 = t0
262  out%T = t
263  out%N = n
264 
265  !! default parameters
266  call allocmem(out%xn_list, 1)
267  out%xn_list = (/1/)
268  out%verb = 1
269  out%Ytn_period = -1._rp
270  out%Vtn_period = -1._rp
271  out%Vtn_rec_prd = -1._rp
272  out%Vtn_plot = .false.
273  out%Vxn_period = -1._rp
274  out%Vxn_file = 'V_xn.dat'
275  out%Vxn_plot = .false.
276  out%act_type = act_no
277  out%act_Vth = -35._rp
278  out%act_plot = .false.
279  out%act_rec = .false.
280  out%stop_if_dep = .false.
281  out%pos = pos_gnuplot
282  out%outDir = "./"
283  out%Vtn_FILE = 'v'
284  out%ACT_FILE = 'act'
285 
286  out%created = .true.
287 
288  end function ode_output_create
289 
290 
308  subroutine ode_output_set(out, &
309  & verb, &
310  & Ytn_period, &
311  & Vtn_period, Vtn_rec_prd, Vtn_plot, &
312  & act_type , act_rec, act_plot, stop_if_dep, act_Vth, &
313  & Vtn_file, act_file, outDir, pos, &
314  & xn_list, Vxn_period, Vxn_plot, Vxn_file)
316  type(ode_output), intent(inOut):: out
317  real(RP), intent(in), optional :: Ytn_period, act_Vth, Vxn_period
318  real(RP), intent(in), optional :: Vtn_period, Vtn_rec_prd
319  integer , intent(in), optional :: act_type, pos, verb
320  logical , intent(in), optional :: Vtn_plot, Vxn_plot
321  logical , intent(in), optional :: act_rec, act_plot, stop_if_dep
322  integer, dimension(:), intent(in), optional :: xn_list
323  character(len=*), intent(in), optional :: Vtn_file, act_file
324  character(len=*), intent(in), optional :: outDir, Vxn_file
325 
326  if (.NOT. out%created) call quit(&
327  & "ode_output_mod: ode_output_set: 'out' not created yet")
328 
329  if (present(verb )) out%verb = verb
330  if (present(ytn_period )) out%Ytn_period = ytn_period
331  if (present(vtn_period )) out%Vtn_period = vtn_period
332  if (present(vtn_rec_prd)) out%Vtn_rec_prd = vtn_rec_prd
333  if (present(act_type )) out%act_type = act_type
334  if (present(outdir )) out%outDir = trim(outdir)
335  if (present(vtn_file )) out%Vtn_file = trim(vtn_file)
336  if (present(act_file )) out%act_file = trim(act_file)
337  if (present(vtn_plot)) out%Vtn_plot = vtn_plot
338  if (present(vxn_period )) out%Vxn_period = vxn_period
339  if (present(vxn_plot)) out%Vxn_plot = vxn_plot
340  if (present(vxn_file)) out%Vxn_file = vxn_file
341 
342  !! Vtn_plot possible only if Vtn_rec_prd>0
343  if (out%Vtn_plot) then
344  if (out%Vtn_rec_prd<0) out%Vtn_plot = .false.
345  end if
346 
347  !! Vxn_plot possible only if Vxn_period>0
348  if (out%Vxn_plot) then
349  if (out%Vxn_period<0) out%Vxn_plot = .false.
350  end if
351 
352  if (present(xn_list)) then
353  call allocmem(out%xn_list, size(xn_list, 1))
354  out%xn_list = xn_list
355  end if
356 
357  if (out%act_type /= act_no) then
358  if (present(stop_if_dep)) out%stop_if_dep = stop_if_dep
359  if (present(act_vth)) out%act_Vth = act_vth
360  if (present(act_rec)) out%act_rec = act_rec
361  if (out%act_rec) then
362  if (present(act_plot)) out%act_plot = act_plot
363  end if
364  end if
365 
366  if (present(pos)) then
367 
368  select case(pos)
369  case(pos_gnuplot)
370  out%pos = pos
371  case(pos_gmsh)
372  out%pos = pos
373  case default
374  call quit("ode_output_mod: ode_output_set:&
375  & unknown post-treatment type")
376  end select
377 
378  end if
379 
380  end subroutine ode_output_set
381 
382 
391  subroutine ode_output_assemble(out, dt, X_h)
392  type(ode_output), intent(inout) :: out
393  real(RP) , intent(in) :: dt
394  type(fespace) , target, optional :: X_h
395 
396  integer :: ii
397 
398  if (choral_verb>1) write(*,*) &
399  & "ode_output_mod : assemble"
400 
401  if (.NOT. out%created) call quit(&
402  & "ode_output_mod: ode_output_assemble: 'out' not created yet")
403 
404  out%dt = dt
405  if (present(x_h)) then
406  out%X_h => x_h
407  out%nbDof = x_h%nbDof
408  out%dim = x_h%mesh%dim
409  else
410  out%nbDof = 1
411  out%dim = 0
412  end if
413 
414  if (choral_verb>1) then
415  write(*,*)" Dimension in space =", out%dim
416  write(*,*)" nbDof =", out%nbDof
417  write(*,*)" dt =", out%dt
418  end if
419 
420  ! counters initialization
421  out%Ytn_cpt = 0
422  out%Vtn_cpt = 0
423  out%Vxn_cpt = 0
424  out%Vtn_rec_cpt = 0
425 
426  ! Y(.,tn)
427  if (out%Ytn_period>0._rp) then
428  ii = int( (out%T - out%t0) / out%Ytn_period) + 1
429  call allocmem(out%Ytn, out%N, out%nbDof, ii)
430  if (choral_verb>1) write(*,*) &
431  & " Array Ytn shape =", &
432  & shape(out%Ytn)
433  end if
434 
435  ! V(.,tn)
436  if (out%Vtn_period>0._rp) then
437  ii = int( (out%T - out%t0) / out%Vtn_period) + 1
438  call allocmem(out%Vtn, out%nbDof, ii)
439  if (choral_verb>1) write(*,*) &
440  & " Array Vtn shape =", &
441  & shape(out%Vtn)
442  end if
443 
444  ! V(xn,.)
445  if (out%Vxn_period>0._rp) then
446 
447  ii = int( (out%T - out%t0) / out%Vxn_period) + 1
448  call allocmem( out%Vxn, size(out%xn_list, 1)+1, ii )
449  if (choral_verb>1) write(*,*) &
450  & " Array Vxn shape =", &
451  & shape(out%Vxn)
452 
453  end if
454 
455  ! activation times
456  out%nv = def_nv(out%act_type)
457  if (out%nv > 0) then
458 
459  call allocmem( out%act, out%nbDof)
460  call allocmem( out%V, out%nV, out%nbDof)
461  if (choral_verb>1) write(*,*) &
462  &" Array V shape =", &
463  &shape(out%V)
464  if (choral_verb>1) write(*,*) &
465  &" Array act shape =", &
466  & shape(out%act)
467 
468  out%act = -1e150_rp
469  ! by default, the initial state is assumed to be depolarised
470  out%V = out%act_Vth - 1._rp
471 
472  end if
473 
474  out%assembled = .true.
475 
476  end subroutine ode_output_assemble
477 
480  function def_nv(method) result(nv)
481  integer :: nv, method
482 
483  select case(method)
484  case(act_no)
485  nv = 0
486  case(act_0)
487  nv = 1
488  case(act_1)
489  nv = 2
490  case(act_2)
491  nv = 3
492  case(act_4)
493  nv = 5
494  case default
495  nv = 0
496  call quit( "ode_output_mod: def_nv:&
497  & unknown computation method for isochrones" )
498  end select
499 
500  end function def_nv
501 
502 
505  subroutine ode_output_print(out)
506  type(ode_output), intent(in) :: out
507 
508  write(*,*)"ode_output_mod : print"
509  write(*,*)" Starting time =", out%t0
510  write(*,*)" Final time =", out%T
511  write(*,*)" Y(tn,xn) \in R**N, N =", out%N
512  if (out%Ytn_period>0.0_rp) then
513  write(*,*)" Storing period for Y(.,tn) =", out%Ytn_period
514  end if
515  if (out%Vtn_period>0.0_rp) then
516  write(*,*)" Storing period for V(.,tn) =", out%Vtn_period
517  end if
518  if (out%Vtn_rec_prd>0.0_rp) then
519  write(*,*)" Saving period for V(.,tn) =", out%Vtn_rec_prd
520  write(*,*)" File name for saving V(tn,.) = ",trim(out%Vtn_file)
521  end if
522  if (out%Vtn_plot) then
523  write(*,*)" Plot of V(tn,:) =", out%Vtn_plot
524  end if
525  if (out%Vxn_period>0.0_rp) then
526  write(*,*)" Storing period for V(xn,.) =", out%Vxn_period
527  write(*,*)" Number of saving nodes xn =", size(out%xn_list, 1)
528  end if
529  if (out%Vxn_plot) then
530  write(*,*)" Plot of V(xn,.) =", out%Vxn_plot
531  end if
532  if(out%act_type >= 0) then
533  write(*,*)" Activation time comp. method =", out%act_type
534  write(*,*)" Threshold for act. time comp. =", out%act_Vth
535  write(*,*)" Stop comp if fully depolarised =", out%stop_if_dep
536  write(*,*)" Saving activation times =", out%act_rec
537  write(*,*)" Plot activation times =", out%act_plot
538  write(*,*)" File name for saving act times = ",trim(out%act_file)
539  end if
540  select case(out%pos)
541  case(pos_gnuplot)
542  write(*,*)" Graphical post process tool = GNUPLOT"
543  case(pos_gmsh)
544  write(*,*)" Graphical post process tool = GMSH"
545  case default
546  write(*,*)" Graphical post process tool = unknown"
547  end select
548  write(*,*)" Directory for output savings = ",trim(out%outDir)
549 
550  end subroutine ode_output_print
551 
552 
569  subroutine output_proc(out, time, s, final)
570  class(ode_output), intent(inout) :: out
571 
572  real(RP) , intent(in) :: time
573  type(ode_solution), intent(in) :: s
574  logical , intent(out) :: final
575 
576  logical :: exit_comp
577  integer :: ii
578  real(RP) :: v1, v2, tn
579  character(LEN=100) :: fic
580 
581  exit_comp=.false.
582  select case(out%verb)
583 
584  ! prints time, potential min max at each ms
585  case(2)
586 
587  if ( abs(anint(time)-time)<1e-4_rp ) then
588  if (allocated(s%V)) then
589  v1 = minval( s%V(:, s%v_i(1)) )
590  v2 = maxval( s%V(:, s%v_i(1)) )
591  else
592  v1 = minval( s%Y(s%N, s%y_i(1), :) )
593  v2 = maxval( s%Y(s%N, s%y_i(1), :) )
594  end if
595  write(*,*)"ode_output : t, V_min, V_Max=", &
596  & real(time,SP), real(v1,SP), real(v2,SP)
597  end if
598 
599  ! prints time, potential min max at each time step
600  case (3)
601  if (allocated(s%V)) then
602  v1 = minval( s%V(:, s%v_i(1)) )
603  v2 = maxval( s%V(:, s%v_i(1)) )
604  else
605  v1 = minval( s%Y(s%N, s%y_i(1), :) )
606  v2 = maxval( s%Y(s%N, s%y_i(1), :) )
607  end if
608  write(*,*)"ode_output : t, V_m, V_M =", &
609  & real(time,SP), real(v1,SP), real(v2,SP)
610 
611  ! prints full component Y min and max
612  case (4)
613  if (allocated(s%Y)) then
614  do ii=1, s%N
615  write(*,*) " Y(ii,:) ii, min, max =", &
616  & int(ii,4), &
617  & real( minval(s%Y(ii, s%y_i(1), :)), SP ), &
618  & real( maxval(s%Y(ii, s%y_i(1), :)), SP )
619  end do
620  else
621  v1 = minval( s%V(:, s%v_i(1)) )
622  v2 = maxval( s%V(:, s%v_i(1)) )
623  write(*,*)"ode_output : t, V_m, V_M =", &
624  & real(time,SP), real(v1,SP), real(v2,SP)
625  end if
626  end select
627 
628  ! storing of Y(., tn)
629  if (out%Ytn_period>0._rp) then
630 
631  tn = out%t0 + re(out%Ytn_cpt) * out%Ytn_period
632  if ( (time + real_tol) >= tn ) then
633 
634  out%Ytn_cpt = out%Ytn_cpt + 1
635  call copy(out%Ytn(:,:,out%Ytn_cpt), s%Y(:,s%y_i(1), :) )
636 
637  if (out%verb>0) then
638  write(*,*)"ode_output : Y(.,tn) stored = tn =", &
639  &real(tn,sp)
640  end if
641 
642  end if
643  end if
644 
645  ! storing of V(., tn)
646  if (out%Vtn_period>0._rp) then
647 
648  tn = out%t0 + re(out%Vtn_cpt) * out%Vtn_period
649  if ( (time + real_tol) >= tn ) then
650 
651  out%Vtn_cpt = out%Vtn_cpt + 1
652  if (allocated(s%V)) then
653  call copy(out%Vtn(:,out%Vtn_cpt), s%v(:,s%v_i(1)))
654  else
655  call copy(out%Vtn(:,out%Vtn_cpt), s%Y(s%N,s%y_i(1),:))
656  end if
657  if (out%verb>0) then
658  write(*,*)"ode_output : V(.,tn) stored = tn =", &
659  & real(tn,sp)
660  end if
661 
662  end if
663  end if
664 
665 
666  ! storing of V(xn, .)
667  if (out%Vxn_period>0._rp) then
668 
669  tn = out%t0 + re(out%Vxn_cpt) * out%Vxn_period
670  if ( (time + real_tol) >= tn ) then
671 
672  out%Vxn_cpt = out%Vxn_cpt + 1
673  if ( out%Vxn_cpt > size(out%Vxn,2) ) return
674  ii = size(out%Vxn,1)
675 
676  out%Vxn(1, out%Vxn_cpt) = time
677  if (allocated(s%V)) then
678  out%Vxn(2:ii,out%Vxn_cpt) = s%v(out%xn_list, s%v_i(1))
679  else
680  out%Vxn(2:ii,out%Vxn_cpt) = s%Y(s%N,s%y_i(1),out%xn_list)
681  end if
682  if (out%verb>2) then
683  write(*,*)"ode_output : V(xn,.) stored = tn =", &
684  & real(tn,sp)
685  end if
686 
687  end if
688  end if
689 
690  ! Rcording of V(., tn) to file vn.dat
691  if (out%Vtn_rec_prd>0._rp) then
692 
693  tn = out%t0 + re(out%Vtn_rec_cpt) * out%Vtn_rec_prd
694  if ( (time + real_tol) >= tn ) then
695 
696  if (out%verb>0) then
697  write(*,"(A35,F10.3)")&
698  & " ode_output : V(., tn) saved = tn = ", tn
699  end if
700 
701  call inttostring(fic, out%Vtn_rec_cpt)
702  fic = trim(out%OUTDIR)//trim(out%Vtn_FILE)//trim(fic)//".dat"
703 
704  if (allocated(s%V)) then
705  call write(s%V(:, s%v_i(1)), trim(fic))
706  else
707  call write(s%Y(s%N, s%y_i(1), :), trim(fic))
708  end if
709 
710  out%Vtn_rec_cpt = out%Vtn_rec_cpt + 1
711  end if
712 
713  end if
714 
715 
716  ! Activation times
717  if (out%nV>0) then
718 
719  ! update potential
720  do ii=out%nV, 2, -1
721  out%V(ii, :) = out%V(ii-1, :)
722  end do
723  if (allocated(s%V)) then
724  call copy(out%V(1,:), s%V(:, s%v_i(1)) )
725  else
726  call copy(out%V(1,:), s%Y(s%N, s%y_i(1),:) )
727  end if
728  select case(out%act_type)
729  case(act_0)
730  call comp_act_0(out%act, out%V, time, out%act_Vth)
731 
732  case(act_1)
733  call comp_act_1(out%act, out%V, time, out%dt, out%act_Vth)
734 
735  case(act_2)
736  call comp_act_2(out%act, out%V, time, out%dt, out%act_Vth)
737 
738  case(act_4)
739  call comp_act_4(out%act, out%V, time, out%dt, out%act_Vth)
740 
741  end select
742 
743  if (out%stop_if_dep) then
744 
745  if (minval(out%act)>-1e50_rp) then
746  exit_comp = .true.
747  if (out%verb>0) write(*,*)"ode_output : &
748  &time =", real(time,sp)
749  if (out%verb>0) write(*,*)" fully depolarised domain"
750  end if
751  end if
752 
753  end if
754 
755  !! Computation is over
756  if (final) then
757  call ode_output_last_output(out)
758  return
759  end if
760 
761  final = exit_comp
762 
763  end subroutine output_proc
764 
765 
768  subroutine ode_output_last_output(out)
769  type(ode_output), intent(in) :: out
770 
771  if (choral_verb>1) write(*,*) &
772  & "ode_output_mod : last_output"
773 
774  if (out%act_rec ) then
775  call write(out%act, trim(out%OUTDIR)//trim(out%ACT_FILE)//".dat")
776  end if
777 
778  if (out%Vtn_plot) call plot_vtn(out)
779  if (out%act_plot) call plot_act(out)
780  if (out%Vxn_plot) call plot_vxn_gnuplot(out)
781 
782  end subroutine ode_output_last_output
783 
784 
785 
788  subroutine plot_vtn(out)
789  type(ode_output), intent(in) :: out
790 
791  if (choral_verb>0) write(*,*) &
792  & "ode_output_mod : plot_vtn"
793 
794  if (out%pos==pos_gnuplot) call plot_vtn_gnuplot(out)
795  if (out%pos==pos_gmsh ) call plot_vtn_gmsh(out)
796 
797  end subroutine plot_vtn
801  subroutine plot_vtn_gnuplot(out)
802  type(ode_output), intent(in) :: out
803 
804  character(LEN=300) :: com
805 
806  if (out%dim>0) then
807  call write(out%X_h%dofCoord(1:out%dim,:), &
808  & trim(out%OUTDIR)//"X.dat", transpose=.true.)
809  end if
810 
811  select case(out%dim)
812 
813  case(1)
814  call inttostring(com, out%Vtn_rec_cpt-1)
815  com = 'gnuplot -e "N=' // trim(com) // '" ' &
816  & // ' -e "root='''// trim(out%outDir) // '''" '&
817  & // trim(gnuplot_dir) // "plot_v_1d.gp"
818 
819  if (choral_verb>2) write(*,*) &
820  & " Shell command > "// trim(com)
821  call system(com)
822 
823  case(2)
824  call inttostring(com, out%Vtn_rec_cpt-1)
825  com = 'gnuplot -e "N=' // trim(com) // '" ' &
826  & // ' -e "root='''// trim(out%outDir) // '''" '&
827  & // trim(gnuplot_dir) // "2d_v.gp"
828 
829  if (choral_verb>2) write(*,*) &
830  & " Shell command > "// trim(com)
831  call system(com)
832 
833  case(3)
834  call warning("ode_output_mod : plot_v_gnuplot=&
835  & dim 3 not implemented", 1)
836 
837  end select
838 
839  end subroutine plot_vtn_gnuplot
840 
841  subroutine plot_vtn_gmsh(out)
842  type(ode_output), intent(in) :: out
843 
844  real(RP), dimension(out%nbDof) :: v
845  character(LEN=300) :: com
846  integer :: ii
847  real(RP) :: time, min, max, val
848 
849  if (out%dim<=1) then
850  call warning('ode_output_mod : plot_v_gmsh = &
851  & dim 2, 3 only', 1)
852  return
853  end if
854 
855  call write(out%X_h, trim(out%OUTDIR)//'v_plot.msh', 'gmsh')
856 
857  min = 0.0_rp
858  max = 0.0_rp
859  do ii=0, out%Vtn_rec_cpt-1
860  call inttostring(com, ii)
861  com = trim(out%OUTDIR)//trim(out%Vtn_FILE)//trim(com)//".dat"
862  call read(v, trim(com))
863 
864  val = minval(v)
865  if (val < min ) min = val
866  val = maxval(v)
867  if (val > max ) max = val
868 
869  time = out%t0 + out%Vtn_rec_prd*re(ii)
870 
871  call gmsh_addview(out%X_h, v, &
872  & trim(out%OUTDIR)//'v_plot.msh', &
873  & 'V', time, ii)
874  end do
875 
876  !! add an option file
877  !! set the range
878  open(unit=60,file=trim(out%OUTDIR)//'v_plot.msh.opt')
879  write(60,*) "View[0].RangeType = 2;"
880  write(60,*) "View[0].CustomMax = ", REAL(max, SP),";"
881  write(60,*) "View[0].CustomMin = ", REAL(min, SP),";"
882  close(60)
883 
884  com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'
885  com = trim(com)//' '//trim(out%OUTDIR)//'v_plot.msh'
886  if (choral_verb>2) write(*,*) &
887  & " Shell command > "// trim(com)
888  call system(com)
889 
890  end subroutine plot_vtn_gmsh
891 
892 
893 
894  subroutine plot_vxn_gnuplot(out)
895  type(ode_output), intent(in) :: out
896 
897  integer :: ii, nn
898  character(LEN=300) :: com, fic
899 
900  !! Recording file
901  !!
902  fic = trim(out%outDir)//trim(out%vxn_file)
903 
904  !! print the transmembrane potential V(xi, tn) to a file
905  open(unit=60,file=trim(fic))
906  do ii=1, out%Vxn_cpt
907  write(60,*) out%Vxn(:,ii)
908  end do
909  close(60)
910 
911  do nn=2, size(out%Vxn,1)
912 
913  !! plot
914  call inttostring(com, nn)
915  com = 'gnuplot -e "fic=''' // trim(fic)//'''" ' &
916  & // ' -e col=' // trim(com)//' '&
917  & // trim(gnuplot_dir)//"plot_v_0d.gp"
918 
919  if (choral_verb>2) write(*,*) &
920  & " Shell command $> "//trim(com)
921  call system(com)
922 
923  end do
924 
925  end subroutine plot_vxn_gnuplot
926 
927 
930  subroutine plot_act(out)
931  type(ode_output), intent(in) :: out
932 
933  if (choral_verb>0) write(*,*) &
934  & "ode_output_mod : plot_act"
935  if (out%pos==pos_gnuplot) call plot_act_gnuplot(out)
936  if (out%pos==pos_gmsh ) call plot_act_gmsh(out)
937 
938  end subroutine plot_act
939 
940  subroutine plot_act_gnuplot(out)
941  type(ode_output), intent(in) :: out
942 
943  character(LEN=300) :: com
944 
945  call write(out%X_h%dofCoord(1:out%dim,:), &
946  & trim(out%OUTDIR)//"X.dat", transpose=.true.)
947 
948  select case(out%dim)
949  case(1)
950  com = 'gnuplot -e "root=''' // trim(out%outDir) // '''" '&
951  & // trim(gnuplot_dir) // "plot_act_1d.gp"
952 
953  if (choral_verb>2) write(*,*) &
954  & " Shell command > "// trim(com)
955  call system(com)
956 
957  case(2)
958  com = 'gnuplot -e "root=''' // trim(out%outDir) // '''" '&
959  & // trim(gnuplot_dir) // "2d_act.gp"
960 
961  if (choral_verb>2) write(*,*) &
962  & " Shell command > "// trim(com)
963  call system(com)
964 
965  case(3)
966  write(*,*)"ode_output_mod : plot_act: &
967  &dim 3 not implemented"
968 
969  end select
970 
971  end subroutine plot_act_gnuplot
972 
973  subroutine plot_act_gmsh(out)
974  type(ode_output), intent(in) :: out
975 
976  real(RP), dimension(out%nbDof) :: v
977  character(LEN=300) :: com
978 
979  if (out%dim==1) then
980  call warning('ode_output_mod : plot_act_gmsh = &
981  & dim 1 not implemented', 1)
982  return
983  end if
984 
985  call write(out%X_h, trim(out%OUTDIR)//'act_plot.msh', 'gmsh')
986 
987  com = trim(out%OUTDIR)//trim(out%ACT_FILE)//".dat"
988  call read(v, trim(com))
989 
990  call gmsh_addview(out%X_h, v, &
991  & trim(out%OUTDIR)//'act_plot.msh', &
992  & 'Activation times in ms', 0._rp, 0)
993 
994  com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'
995  com = trim(com)//' '//trim(out%OUTDIR)//'act_plot.msh'
996  if (choral_verb>2) write(*,*) &
997  & " Shell command > "// trim(com)
998  call system(com)
999 
1000  end subroutine plot_act_gmsh
1001 
1002 
1003 
1004  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1005  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1006  ! !!
1007  ! !!
1008  ! !! ACTIVATION TIMES COMPUTATION
1009  ! !!
1010  ! !!
1011  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1012  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1013 
1016  subroutine comp_act_0(act, V, time, Vth)
1017  real(RP), dimension(:) , intent(out) :: act
1018  real(RP), dimension(:,:), intent(in) :: V
1019  real(RP) , intent(in) :: time, Vth
1020 
1021  integer :: ii
1022 
1023  do ii=1, size(act,1)
1024 
1025  if (act(ii)<-1e50_rp) then
1026  if (v(1,ii)>vth) act(ii) = time
1027  end if
1028 
1029  end do
1030 
1031  end subroutine comp_act_0
1032 
1033 
1034 
1035 
1036 
1039  subroutine comp_act_1(act, V, time, dt, Vth)
1040  real(RP), dimension(:) , intent(out) :: act
1041  real(RP), dimension(:,:), intent(in) :: V
1042  real(RP) , intent(in) :: time, dt, Vth
1043 
1044  integer :: ii
1045  real(RP) :: v2, v1
1046 
1047  do ii=1, size(act, 1)
1048 
1049  if (act(ii)<-1e50_rp) then
1050  v1 = v(1,ii)
1051  if (v1>vth) then
1052  v2 = v(2,ii)
1053  act(ii) = time + ( vth - v1 )/( v1 - v2 ) * dt
1054  end if
1055  end if
1056 
1057  end do
1058 
1059  end subroutine comp_act_1
1060 
1061 
1064  subroutine comp_act_2(act, V, time, dt, Vth)
1065  real(RP), dimension(:) , intent(out) :: act
1066  real(RP), dimension(:,:), intent(in) :: V
1067  real(RP) , intent(in) :: time, dt, Vth
1068 
1069  integer :: ii
1070  real(RP) :: v3, v2, v1, a, b, delta
1071 
1072  do ii=1, size(act, 1)
1073 
1074  if (act(ii)<-1e50_rp) then
1075  v1 = v(1,ii)
1076  if (v1>vth) then
1077  v2 = v(2,ii)
1078  v3 = v(3,ii)
1079 
1080  a = (v1 - 2._rp*v2 + v3)*0.5_rp
1081  b = (v1 - v3)*0.5_rp
1082  delta = b**2 - 4._rp*a*(v2-vth)
1083 
1084  ! highest root
1085  a = ( sqrt(delta) - b ) / (2._rp * a)
1086 
1087  act(ii) = time + dt * (a - 1._rp)
1088 
1089  end if
1090  end if
1091 
1092  end do
1093 
1094  end subroutine comp_act_2
1095 
1096 
1097 
1100  subroutine comp_act_4(act, V, time, dt, Vth)
1101  real(RP), dimension(:) , intent(out) :: act
1102  real(RP), dimension(:,:), intent(in) :: V
1103  real(RP) , intent(in) :: time, dt, Vth
1104 
1105  integer :: ii
1106  real(RP) :: v5, v4, v3, v2, v1, a, b, c, d, e, tol
1107 
1108  tol = real_tol * 1e4_rp
1109 
1110  do ii=1, size(act, 1)
1111 
1112  if (act(ii)<-1e50_rp) then
1113  v2 = v(2,ii)
1114  if (v2>vth) then
1115  v1 = v(1,ii)
1116  v3 = v(3,ii)
1117  v4 = v(4,ii)
1118  v5 = v(5,ii)
1119 
1120  a = v1 - 4._rp*v2 + 6._rp*v3 - 4._rp*v4 + v5
1121  a = a/24._rp
1122 
1123  c =-v1 + 16._rp*v2 - 30._rp*v3 + 16._rp*v4 - v5
1124  c = c/24._rp
1125 
1126  b = v1 - 2._rp*v2 + 2._rp*v4 - v5
1127  b = b/12._rp
1128 
1129  d =-v1 + 8._rp*v2 - 8._rp*v4 + v5
1130  d = d/12._rp
1131 
1132  e = v3 - vth
1133 
1134  a = local_newton()
1135 
1136  act(ii) = time + dt * (a - 2._rp)
1137 
1138  end if
1139  end if
1140 
1141  end do
1142 
1143  contains
1144 
1145  function p(tau)
1146  real(RP) :: p
1147  real(RP), intent(in) :: tau
1148 
1149  p = a*tau**4 + b*tau**3 + c*tau**2 + d*tau + e
1150 
1151  end function p
1152 
1153  function d_p(tau)
1154  real(RP) :: d_p
1155  real(RP), intent(in) :: tau
1156 
1157  d_p = 4._rp*a*tau**3 + 3._rp*b*tau**2 + 2._rp*c*tau + d
1158 
1159  end function d_p
1160 
1161  function local_newton() result(s)
1162  real(RP) :: s
1163 
1164  real(RP):: ps
1165  integer :: cpt
1166 
1167  cpt = 0
1168  s = 0._rp !(vth-v3) / (v2-v3)
1169  ps = e !p(s)
1170 
1171  do while( abs(ps) > tol )
1172 
1173  s = s - ps / d_p(s)
1174  ps = p(s)
1175 
1176  cpt = cpt + 1
1177  if (cpt==1000) then
1178  call warning('ode_output_mod: comp_isoch_4 =&
1179  & Newton error', 1)
1180  return
1181  end if
1182  end do
1183 
1184  end function local_newton
1185 
1186  end subroutine comp_act_4
1187 
1188 
1189  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1190  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1191  ! !!
1192  ! !!
1193  ! !! WAVE FRONT ANALYSIS
1194  ! !!
1195  ! !!
1196  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1197  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1198 
1210  function celerity_l2_dist(u1, u2, X_h, qdm, weight) result(dist)
1211  real(RP) :: dist
1212  real(RP), dimension(:), intent(in) :: u1, u2
1213  type(fespace) , intent(in) :: X_h
1214  type(quadmesh) , intent(in) :: qdm
1215  procedure(r3tor) :: weight
1216 
1217  dist = integ(e, u1, u2, x_h, qdm)
1218  dist = sqrt(dist)
1219 
1220  contains
1221 
1222  function e(x, v1, v2)
1223  real(RP) :: E
1224  real(RP), dimension(3), intent(in) :: x, v1, v2
1225 
1226  real(RP), dimension(3) :: yy
1227 
1228  yy = v1 / sum(v1*v1) - v2 / sum(v2*v2)
1229  e = sum(yy*yy)
1230  e = weight(x) * e
1231 
1232  end function e
1233 
1234  end function celerity_l2_dist
1235 
1236 
1248  function celerity_l2_norm(u1, X_h, qdm, weight) result(dist)
1249  real(RP) :: dist
1250  real(RP) , dimension(:), intent(in) :: u1
1251  type(fespace) , intent(in) :: X_h
1252  type(quadmesh) , intent(in) :: qdm
1253  procedure(r3tor) :: weight
1254 
1255 
1256  dist = integ(e, u1, u1, x_h, qdm)
1257  dist = sqrt(dist)
1258 
1259  contains
1260 
1261  function e(x, v1, v2)
1262  real(RP) :: E
1263  real(RP), dimension(3), intent(in) :: x, v1, v2
1264 
1265  real(RP), dimension(3) :: yy
1266 
1267  yy = v1 / sum(v1*v1)
1268  e = sum(yy*yy)
1269  e = weight(x) * e
1270 
1271  end function e
1272 
1273  end function celerity_l2_norm
1274 
1275 
1276 end module ode_output_mod
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
Type ode_solution: data structure to solve ODE/PDE problems.
real(rp) function local_newton()
BASIC TOOLS
Definition: basic_tools.F90:8
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
subroutine comp_act_2(act, V, time, dt, Vth)
Computation of activation times : quadradic interp.
assemble the output (finalise construction)
integer, parameter act_1
linear interpolation
subroutine ode_output_last_output(out)
ends computation output
subroutine plot_vtn_gmsh(out)
integer, parameter pos_gnuplot
subroutine comp_act_0(act, V, time, Vth)
Computation of activation times : cut off.
subroutine, public quit(message)
Stop program execution, display an error messahe.
real(rp) function d_p(tau)
write to file
Definition: io.f90:33
real(rp) function, public celerity_l2_dist(u1, u2, X_h, qdm, weight)
returns || (c(u1) - c(u2)) * weight(x) ||_L2()
character(len=150), parameter, public gmsh_dir
path to &#39;gmsh&#39; directory = CHORAL_DIR/ress/gmsh
Definition: choral_env.f90:33
type(ode_output) function ode_output_create(t0, T, N)
Constructor for the type ode_output
print a short description
conversion integers or rational to real
Definition: real_type.F90:153
character(len=150), parameter, public gnuplot_dir
path to &#39;gnuplot&#39; directory = CHORAL_DIR/ress/gmsh
Definition: choral_env.f90:39
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
DERIVED TYPE quadMesh: integration methods on meshes
IO: module for input/output
Definition: io.f90:8
subroutine ode_output_clear(out)
destructor for ode_output
DERIVED TYPE feSpace: finite element spaces
Definition: feSpace_mod.F90:58
DERIVED TYPE ode_solution: data straucture to solve ODEs
y = x
Definition: R1d_mod.F90:44
subroutine ode_output_print(out)
print ode_output
real(rp) function e(x, v1, v2)
CHORAL CONSTANTS
integer, parameter pos_gmsh
real(rp) function p(tau)
subroutine, public output_proc(out, time, s, final)
Pre-defined output for ODE resolution
set the output parameters
real(rp) function, public celerity_l2_norm(u1, X_h, qdm, weight)
returns || c(u1) * weight(x) ||_L2
Integration of scalar functions.
Definition: integral.F90:54
subroutine plot_act_gnuplot(out)
DERIVED TYPE ode_output: handles output for PDE/ODE simulations
subroutine comp_act_4(act, V, time, dt, Vth)
Computation of activation times : bi-quadradic interp.
allocate memory for real(RP) arrays
Definition: real_type.F90:158
Derived type feSpace: finite element space on a mesh.
COMPUTATION OF INTEGRALS using a quadrature methods on a mesh.
Definition: integral.F90:21
subroutine comp_act_1(act, V, time, dt, Vth)
Computation of activation times : linear interp.
integer choral_verb
Verbosity level.
subroutine plot_vtn(out)
Plot V(x, tn) "co$outDir/vn.dat", n=1..Ytn_cpt.
subroutine plot_act_gmsh(out)
DEFINITION OF DIRECTORY PATHS
Definition: choral_env.f90:16
The type ode_output handles output for ODE simulations.
read from file
Definition: io.f90:39
real(rp), parameter, public real_tol
Definition: real_type.F90:91
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
The type quadMesh defines integration methods on meshes.
subroutine ode_output_assemble(out, dt, X_h)
assemble the output, follows set(out, ...)
integer, parameter act_no
no computation
integer, parameter act_4
bi-quadratic interpolation
OPENMP OPERATIONS ON 2-DIMENSIONAL REAL ARRAYS
Definition: R2d_mod.F90:24
subroutine, public inttostring(str, ii)
convert an integer to a string
Definition: basic_tools.F90:42
subroutine plot_vxn_gnuplot(out)
integer, parameter act_0
Heaviside.
integer function def_nv(method)
size of V required to compute act
integer, parameter act_2
quadratic interpolation
subroutine, public warning(message, verb)
Warning message.
subroutine plot_act(out)
plot activation times
subroutine, public gmsh_addview(X_h, v, fileName, vname, time, idx)
add a view to X_h output file gmesh file format
subroutine ode_output_set(out, verb, Ytn_period, Vtn_period, Vtn_rec_prd, Vtn_plot, act_type, act_rec, act_plot, stop_if_dep, act_Vth, Vtn_file, act_file, outDir, pos, xn_list, Vxn_period, Vxn_plot, Vxn_file)
set the following parameters = optional see ode_output_mod::ode_output for the default values ...
subroutine plot_vtn_gnuplot(out)
plot ../temp/vn.dat