90 logical :: created = .false.
93 logical :: assembled = .false.
113 real(RP) :: ytn_period
115 real(RP) :: vtn_period
117 real(RP) :: vtn_rec_prd
121 integer,
dimension(:),
allocatable :: xn_list
123 real(RP) :: vxn_period
125 character(len=75) :: vxn_file
137 logical :: stop_if_dep
141 character(len=150) :: outdir
143 character(len=50) :: vtn_file
145 character(len=50) :: act_file
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
183 integer :: vtn_rec_cpt
229 out%created = .false.
230 out%assembled = .false.
257 real(RP),
intent(in) :: t0, T
258 integer ,
intent(in) :: N
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.
277 out%act_Vth = -35._rp
278 out%act_plot = .false.
279 out%act_rec = .false.
280 out%stop_if_dep = .false.
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)
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
326 if (.NOT. out%created)
call quit(&
327 &
"ode_output_mod: ode_output_set: 'out' not created yet")
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
343 if (out%Vtn_plot)
then 344 if (out%Vtn_rec_prd<0) out%Vtn_plot = .false.
348 if (out%Vxn_plot)
then 349 if (out%Vxn_period<0) out%Vxn_plot = .false.
352 if (
present(xn_list))
then 353 call allocmem(out%xn_list,
size(xn_list, 1))
354 out%xn_list = xn_list
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
366 if (
present(pos))
then 374 call quit(
"ode_output_mod: ode_output_set:& 375 & unknown post-treatment type")
393 real(RP) ,
intent(in) :: dt
394 type(
fespace) ,
target,
optional :: X_h
399 &
"ode_output_mod : assemble" 401 if (.NOT. out%created)
call quit(&
402 &
"ode_output_mod: ode_output_assemble: 'out' not created yet")
405 if (
present(x_h))
then 407 out%nbDof = x_h%nbDof
408 out%dim = x_h%mesh%dim
415 write(*,*)
" Dimension in space =", out%dim
416 write(*,*)
" nbDof =", out%nbDof
417 write(*,*)
" dt =", out%dt
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)
431 &
" Array Ytn shape =", &
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)
440 &
" Array Vtn shape =", &
445 if (out%Vxn_period>0._rp)
then 447 ii = int( (out%T - out%t0) / out%Vxn_period) + 1
448 call allocmem( out%Vxn,
size(out%xn_list, 1)+1, ii )
450 &
" Array Vxn shape =", &
456 out%nv =
def_nv(out%act_type)
460 call allocmem( out%V, out%nV, out%nbDof)
462 &
" Array V shape =", &
465 &
" Array act shape =", &
470 out%V = out%act_Vth - 1._rp
474 out%assembled = .true.
480 function def_nv(method)
result(nv)
481 integer :: nv, method
496 call quit(
"ode_output_mod: def_nv:& 497 & unknown computation method for isochrones" )
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
515 if (out%Vtn_period>0.0_rp)
then 516 write(*,*)
" Storing period for V(.,tn) =", out%Vtn_period
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)
522 if (out%Vtn_plot)
then 523 write(*,*)
" Plot of V(tn,:) =", out%Vtn_plot
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)
529 if (out%Vxn_plot)
then 530 write(*,*)
" Plot of V(xn,.) =", out%Vxn_plot
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)
542 write(*,*)
" Graphical post process tool = GNUPLOT" 544 write(*,*)
" Graphical post process tool = GMSH" 546 write(*,*)
" Graphical post process tool = unknown" 548 write(*,*)
" Directory for output savings = ",trim(out%outDir)
572 real(RP) ,
intent(in) :: time
574 logical ,
intent(out) :: final
578 real(RP) :: v1, v2, tn
579 character(LEN=100) :: fic
582 select case(out%verb)
587 if ( abs(anint(time)-time)<1
e-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)) )
592 v1 = minval( s%Y(s%N, s%y_i(1), :) )
593 v2 = maxval( s%Y(s%N, s%y_i(1), :) )
595 write(*,*)
"ode_output : t, V_min, V_Max=", &
596 &
real(time,SP),
real(v1,SP),
real(v2,SP) 601 if (
allocated(s%V))
then 602 v1 = minval( s%V(:, s%v_i(1)) )
603 v2 = maxval( s%V(:, s%v_i(1)) )
605 v1 = minval( s%Y(s%N, s%y_i(1), :) )
606 v2 = maxval( s%Y(s%N, s%y_i(1), :) )
608 write(*,*)
"ode_output : t, V_m, V_M =", &
609 &
real(time,SP),
real(v1,SP),
real(v2,SP) 613 if (
allocated(s%Y))
then 615 write(*,*)
" Y(ii,:) ii, min, max =", &
617 &
real( minval(s%Y(ii, s%y_i(1), :)), SP ), &
618 & real( maxval(s%Y(ii, s%y_i(1), :)), SP )
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) 629 if (out%Ytn_period>0._rp)
then 631 tn = out%t0 +
re(out%Ytn_cpt) * out%Ytn_period
634 out%Ytn_cpt = out%Ytn_cpt + 1
635 call copy(out%Ytn(:,:,out%Ytn_cpt), s%Y(:,s%y_i(1), :) )
638 write(*,*)
"ode_output : Y(.,tn) stored = tn =", &
646 if (out%Vtn_period>0._rp)
then 648 tn = out%t0 +
re(out%Vtn_cpt) * out%Vtn_period
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)))
655 call copy(out%Vtn(:,out%Vtn_cpt), s%Y(s%N,s%y_i(1),:))
658 write(*,*)
"ode_output : V(.,tn) stored = tn =", &
667 if (out%Vxn_period>0._rp)
then 669 tn = out%t0 +
re(out%Vxn_cpt) * out%Vxn_period
672 out%Vxn_cpt = out%Vxn_cpt + 1
673 if ( out%Vxn_cpt >
size(out%Vxn,2) )
return 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))
680 out%Vxn(2:ii,out%Vxn_cpt) = s%Y(s%N,s%y_i(1),out%xn_list)
683 write(*,*)
"ode_output : V(xn,.) stored = tn =", &
691 if (out%Vtn_rec_prd>0._rp)
then 693 tn = out%t0 +
re(out%Vtn_rec_cpt) * out%Vtn_rec_prd
697 write(*,
"(A35,F10.3)")&
698 &
" ode_output : V(., tn) saved = tn = ", tn
702 fic = trim(out%OUTDIR)//trim(out%Vtn_FILE)//trim(fic)//
".dat" 704 if (
allocated(s%V))
then 705 call write(s%V(:, s%v_i(1)), trim(fic))
707 call write(s%Y(s%N, s%y_i(1), :), trim(fic))
710 out%Vtn_rec_cpt = out%Vtn_rec_cpt + 1
721 out%V(ii, :) = out%V(ii-1, :)
723 if (
allocated(s%V))
then 724 call copy(out%V(1,:), s%V(:, s%v_i(1)) )
726 call copy(out%V(1,:), s%Y(s%N, s%y_i(1),:) )
728 select case(out%act_type)
730 call comp_act_0(out%act, out%V, time, out%act_Vth)
733 call comp_act_1(out%act, out%V, time, out%dt, out%act_Vth)
736 call comp_act_2(out%act, out%V, time, out%dt, out%act_Vth)
739 call comp_act_4(out%act, out%V, time, out%dt, out%act_Vth)
743 if (out%stop_if_dep)
then 745 if (minval(out%act)>-1e50_rp)
then 747 if (out%verb>0)
write(*,*)
"ode_output : & 748 &time =",
real(time,sp)
749 if (out%verb>0)
write(*,*)
" fully depolarised domain" 772 &
"ode_output_mod : last_output" 774 if (out%act_rec )
then 775 call write(out%act, trim(out%OUTDIR)//trim(out%ACT_FILE)//
".dat")
778 if (out%Vtn_plot)
call plot_vtn(out)
779 if (out%act_plot)
call plot_act(out)
792 &
"ode_output_mod : plot_vtn" 804 character(LEN=300) :: com
807 call write(out%X_h%dofCoord(1:out%dim,:), &
808 & trim(out%OUTDIR)//
"X.dat", transpose=.true.)
815 com =
'gnuplot -e "N=' // trim(com) //
'" ' &
816 & //
' -e "root='''// trim(out%outDir) //
'''" '&
820 &
" Shell command > "// trim(com)
825 com =
'gnuplot -e "N=' // trim(com) //
'" ' &
826 & //
' -e "root='''// trim(out%outDir) //
'''" '&
830 &
" Shell command > "// trim(com)
834 call warning(
"ode_output_mod : plot_v_gnuplot=& 835 & dim 3 not implemented", 1)
844 real(RP),
dimension(out%nbDof) :: v
845 character(LEN=300) :: com
847 real(RP) :: time, min, max, val
850 call warning(
'ode_output_mod : plot_v_gmsh = & 855 call write(out%X_h, trim(out%OUTDIR)//
'v_plot.msh',
'gmsh')
859 do ii=0, out%Vtn_rec_cpt-1
861 com = trim(out%OUTDIR)//trim(out%Vtn_FILE)//trim(com)//
".dat" 862 call read(v, trim(com))
865 if (val < min ) min = val
867 if (val > max ) max = val
869 time = out%t0 + out%Vtn_rec_prd*
re(ii)
872 & trim(out%OUTDIR)//
'v_plot.msh', &
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),
";" 884 com =
'gmsh -option '//trim(
gmsh_dir)//
'gmsh-options-view' 885 com = trim(com)//
' '//trim(out%OUTDIR)//
'v_plot.msh' 887 &
" Shell command > "// trim(com)
898 character(LEN=300) :: com, fic
902 fic = trim(out%outDir)//trim(out%vxn_file)
905 open(unit=60,file=trim(fic))
907 write(60,*) out%Vxn(:,ii)
911 do nn=2,
size(out%Vxn,1)
915 com =
'gnuplot -e "fic=''' // trim(fic)//
'''" ' &
916 & //
' -e col=' // trim(com)//
' '&
920 &
" Shell command $> "//trim(com)
934 &
"ode_output_mod : plot_act" 943 character(LEN=300) :: com
945 call write(out%X_h%dofCoord(1:out%dim,:), &
946 & trim(out%OUTDIR)//
"X.dat", transpose=.true.)
950 com =
'gnuplot -e "root=''' // trim(out%outDir) //
'''" '&
954 &
" Shell command > "// trim(com)
958 com =
'gnuplot -e "root=''' // trim(out%outDir) //
'''" '&
962 &
" Shell command > "// trim(com)
966 write(*,*)
"ode_output_mod : plot_act: & 967 &dim 3 not implemented" 976 real(RP),
dimension(out%nbDof) :: v
977 character(LEN=300) :: com
980 call warning(
'ode_output_mod : plot_act_gmsh = & 981 & dim 1 not implemented', 1)
985 call write(out%X_h, trim(out%OUTDIR)//
'act_plot.msh',
'gmsh')
987 com = trim(out%OUTDIR)//trim(out%ACT_FILE)//
".dat" 988 call read(v, trim(com))
991 & trim(out%OUTDIR)//
'act_plot.msh', &
992 &
'Activation times in ms', 0._rp, 0)
994 com =
'gmsh -option '//trim(
gmsh_dir)//
'gmsh-options-view' 995 com = trim(com)//
' '//trim(out%OUTDIR)//
'act_plot.msh' 997 &
" Shell command > "// trim(com)
1017 real(RP),
dimension(:) ,
intent(out) :: act
1018 real(RP),
dimension(:,:),
intent(in) :: V
1019 real(RP) ,
intent(in) :: time, Vth
1023 do ii=1,
size(act,1)
1025 if (act(ii)<-1e50_rp)
then 1026 if (v(1,ii)>vth) act(ii) = time
1040 real(RP),
dimension(:) ,
intent(out) :: act
1041 real(RP),
dimension(:,:),
intent(in) :: V
1042 real(RP) ,
intent(in) :: time, dt, Vth
1047 do ii=1,
size(act, 1)
1049 if (act(ii)<-1e50_rp)
then 1053 act(ii) = time + ( vth - v1 )/( v1 - v2 ) * dt
1065 real(RP),
dimension(:) ,
intent(out) :: act
1066 real(RP),
dimension(:,:),
intent(in) :: V
1067 real(RP) ,
intent(in) :: time, dt, Vth
1070 real(RP) :: v3, v2, v1, a, b, delta
1072 do ii=1,
size(act, 1)
1074 if (act(ii)<-1e50_rp)
then 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)
1085 a = ( sqrt(delta) - b ) / (2._rp * a)
1087 act(ii) = time + dt * (a - 1._rp)
1101 real(RP),
dimension(:) ,
intent(out) :: act
1102 real(RP),
dimension(:,:),
intent(in) :: V
1103 real(RP) ,
intent(in) :: time, dt, Vth
1106 real(RP) :: v5, v4, v3, v2, v1, a, b, c, d, e, tol
1110 do ii=1,
size(act, 1)
1112 if (act(ii)<-1e50_rp)
then 1120 a = v1 - 4._rp*v2 + 6._rp*v3 - 4._rp*v4 + v5
1123 c =-v1 + 16._rp*v2 - 30._rp*v3 + 16._rp*v4 - v5
1126 b = v1 - 2._rp*v2 + 2._rp*v4 - v5
1129 d =-v1 + 8._rp*v2 - 8._rp*v4 + v5
1136 act(ii) = time + dt * (a - 2._rp)
1147 real(RP),
intent(in) :: tau
1149 p = a*tau**4 + b*tau**3 + c*tau**2 + d*tau + e
1155 real(RP),
intent(in) :: tau
1157 d_p = 4._rp*a*tau**3 + 3._rp*b*tau**2 + 2._rp*c*tau + d
1171 do while( abs(ps) > tol )
1178 call warning(
'ode_output_mod: comp_isoch_4 =& 1212 real(RP),
dimension(:),
intent(in) :: u1, u2
1213 type(
fespace) ,
intent(in) :: X_h
1215 procedure(r3tor) :: weight
1217 dist =
integ(
e, u1, u2, x_h, qdm)
1222 function e(x, v1, v2)
1224 real(RP),
dimension(3),
intent(in) :: x, v1, v2
1226 real(RP),
dimension(3) :: yy
1228 yy = v1 / sum(v1*v1) - v2 / sum(v2*v2)
1250 real(RP) ,
dimension(:),
intent(in) :: u1
1251 type(
fespace) ,
intent(in) :: X_h
1253 procedure(r3tor) :: weight
1256 dist =
integ(
e, u1, u1, x_h, qdm)
1261 function e(x, v1, v2)
1263 real(RP),
dimension(3),
intent(in) :: x, v1, v2
1265 real(RP),
dimension(3) :: yy
1267 yy = v1 / sum(v1*v1)
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
Type ode_solution: data structure to solve ODE/PDE problems.
real(rp) function local_newton()
deallocate memory for real(RP) arrays
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.
real(rp) function d_p(tau)
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 'gmsh' directory = CHORAL_DIR/ress/gmsh
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
character(len=150), parameter, public gnuplot_dir
path to 'gnuplot' directory = CHORAL_DIR/ress/gmsh
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
DERIVED TYPE quadMesh: integration methods on meshes
IO: module for input/output
subroutine ode_output_clear(out)
destructor for ode_output
DERIVED TYPE feSpace: finite element spaces
DERIVED TYPE ode_solution: data straucture to solve ODEs
subroutine ode_output_print(out)
print ode_output
real(rp) function e(x, v1, v2)
integer, parameter pos_gmsh
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.
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
Derived type feSpace: finite element space on a mesh.
COMPUTATION OF INTEGRALS using a quadrature methods on a mesh.
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
The type ode_output handles output for ODE simulations.
real(rp), parameter, public real_tol
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
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 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