ode_output_mod.F90

Go to the documentation of this file.

module ode_output_mod

  use choral_constants
  use choral_variables
  use choral_env
  use real_type
  use basic_tools
  use abstract_interfaces, only: r3tor
  use io
  use r1d_mod, only: copy
  use r2d_mod, only: copy
  use fespace_mod
  use integral
  use quadmesh_mod
  use ode_solution_mod
 
  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  public :: ode_output
  public :: clear, set, assemble, print
  public :: output_proc

  !! Wave-front analysis 
  public :: celerity_l2_dist, celerity_l2_norm


  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type :: ode_output

     logical :: created   = .false.

     logical :: assembled = .false.

     !! Constructor parameter: provided by the user with
     !! 'out = ode-output(...)'
     !!
     integer  :: dim       
     integer  :: n    
     real(RP) :: t0
     real(RP) :: t

     !! Parameters provided by the user with 
     !! 'call set(out, ...)'
     !!
     integer :: verb
     real(RP) :: ytn_period
     real(RP) :: vtn_period
     real(RP) :: vtn_rec_prd
     logical  :: vtn_plot
     integer, dimension(:), allocatable :: xn_list
     real(RP) :: vxn_period
     character(len=75) :: vxn_file 
     logical  :: vxn_plot
     integer  :: act_type 
     real(RP) :: act_vth 
     logical  :: act_plot
     logical  :: act_rec
     logical  :: stop_if_dep
     integer  :: pos 
     character(len=150) :: outdir
     character(len=50)  :: vtn_file
     character(len=50)  :: act_file


     !! Parameters provided by the user with 
     !! 'call assemble(out, ...)'
     !!
     real(RP) :: dt
     type(fespace), pointer :: x_h=>null()


     !! The following data are initialised by the code 
     !! when assembling
     !!
     real(RP), dimension(:,:,:), allocatable :: ytn 
     real(RP), dimension(:,:)  , allocatable :: vtn 
     real(RP), dimension(:,:)  , allocatable :: vxn 
     real(RP), dimension(:,:)  , allocatable :: v   
     real(RP), dimension(:)    , allocatable :: act 
     integer :: nbdof 
     integer :: nv 
     integer :: ytn_cpt 
     integer :: vtn_cpt
     integer :: vtn_rec_cpt
     integer :: vxn_cpt

   contains

     final :: ode_output_clear

  end type ode_output


  !       %----------------------------------------%
  !       |                                        |
  !       |          GENERIC INTERFACES            |
  !       |                                        |
  !       %----------------------------------------%
  interface clear
     module procedure ode_output_clear
  end interface clear

  interface ode_output
     module procedure ode_output_create
  end interface ode_output

  interface set
     module procedure ode_output_set
  end interface set

  interface assemble
     module procedure ode_output_assemble
  end interface assemble

  interface print
     module procedure ode_output_print
  end interface print

contains

  subroutine ode_output_clear(out)
    type(ode_output), intent(inout) :: out

    out%created   = .false.
    out%assembled = .false.
    
    call freemem(out%xn_list)
 
    call freemem(out%Ytn )
    call freemem(out%Vtn )
    call freemem(out%Vxn )
    call freemem(out%V )
    call freemem(out%act )

    out%X_h=>null()

  end subroutine ode_output_clear


  function ode_output_create(t0, T, N) result(out)
    type(ode_output)     :: out
    real(RP), intent(in) :: t0, T
    integer , intent(in) :: N

    call clear(out)
    out%t0  = t0
    out%T   = t
    out%N   = n

    !! default parameters
    call allocmem(out%xn_list, 1)
    out%xn_list     = (/1/)
    out%verb        = 1
    out%Ytn_period  = -1._rp
    out%Vtn_period  = -1._rp
    out%Vtn_rec_prd = -1._rp
    out%Vtn_plot    = .false.
    out%Vxn_period  = -1._rp
    out%Vxn_file    = 'V_xn.dat' 
    out%Vxn_plot    = .false.
    out%act_type    =  act_no     
    out%act_Vth     = -35._rp   
    out%act_plot    = .false.
    out%act_rec     = .false.
    out%stop_if_dep = .false.  
    out%pos         =  pos_gnuplot
    out%outDir      = "./"
    out%Vtn_FILE    = 'v'
    out%ACT_FILE    = 'act'

    out%created     = .true.

  end function ode_output_create


  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)    

    type(ode_output), intent(inOut):: out
    real(RP), intent(in), optional :: Ytn_period, act_Vth, Vxn_period
    real(RP), intent(in), optional :: Vtn_period, Vtn_rec_prd
    integer , intent(in), optional :: act_type, pos, verb
    logical , intent(in), optional :: Vtn_plot, Vxn_plot
    logical , intent(in), optional :: act_rec, act_plot, stop_if_dep
    integer, dimension(:), intent(in), optional :: xn_list
    character(len=*), intent(in), optional :: Vtn_file, act_file
    character(len=*), intent(in), optional :: outDir, Vxn_file 

    if (.NOT. out%created) call quit(&
         & "ode_output_mod: ode_output_set: 'out' not created yet")

    if (present(verb       )) out%verb        = verb
    if (present(ytn_period )) out%Ytn_period  = ytn_period
    if (present(vtn_period )) out%Vtn_period  = vtn_period
    if (present(vtn_rec_prd)) out%Vtn_rec_prd = vtn_rec_prd
    if (present(act_type   )) out%act_type    = act_type
    if (present(outdir     )) out%outDir      = trim(outdir)
    if (present(vtn_file   )) out%Vtn_file    = trim(vtn_file)
    if (present(act_file   )) out%act_file    = trim(act_file)
    if (present(vtn_plot)) out%Vtn_plot = vtn_plot
    if (present(vxn_period )) out%Vxn_period  = vxn_period
    if (present(vxn_plot)) out%Vxn_plot = vxn_plot
    if (present(vxn_file)) out%Vxn_file = vxn_file

    !! Vtn_plot possible only if Vtn_rec_prd>0
    if (out%Vtn_plot) then
       if (out%Vtn_rec_prd<0) out%Vtn_plot = .false.
    end if

    !! Vxn_plot possible only if Vxn_period>0
    if (out%Vxn_plot) then
       if (out%Vxn_period<0) out%Vxn_plot = .false.
    end if

    if (present(xn_list)) then
       call allocmem(out%xn_list, size(xn_list, 1)) 
       out%xn_list = xn_list
    end if

    if (out%act_type /= act_no) then
       if (present(stop_if_dep)) out%stop_if_dep = stop_if_dep
       if (present(act_vth))     out%act_Vth     = act_vth
       if (present(act_rec))     out%act_rec     = act_rec
       if (out%act_rec) then
          if (present(act_plot)) out%act_plot    = act_plot
       end if
    end if

    if (present(pos)) then

       select case(pos)
       case(pos_gnuplot)
          out%pos = pos
       case(pos_gmsh)
          out%pos = pos
       case default
          call quit("ode_output_mod: ode_output_set:&
            &  unknown post-treatment type")
       end select

    end if

  end subroutine ode_output_set


  subroutine ode_output_assemble(out, dt, X_h)
    type(ode_output), intent(inout)    :: out
    real(RP)        , intent(in)       :: dt
    type(fespace)   , target, optional :: X_h

    integer :: ii

    if (choral_verb>1) write(*,*) &
         & "ode_output_mod  : assemble"

    if (.NOT. out%created) call quit(&
         & "ode_output_mod: ode_output_assemble: 'out' not created yet")

    out%dt    =  dt
    if (present(x_h)) then
       out%X_h   => x_h
       out%nbDof =  x_h%nbDof
       out%dim   =  x_h%mesh%dim
    else
       out%nbDof =  1
       out%dim   =  0
    end if

    if (choral_verb>1) then
       write(*,*)"  Dimension in space             =", out%dim
       write(*,*)"  nbDof                          =", out%nbDof
       write(*,*)"  dt                             =", out%dt   
    end if

    ! counters initialization
    out%Ytn_cpt     = 0
    out%Vtn_cpt     = 0
    out%Vxn_cpt     = 0
    out%Vtn_rec_cpt = 0

    ! Y(.,tn) 
    if (out%Ytn_period>0._rp) then 
       ii = int( (out%T - out%t0) / out%Ytn_period) + 1
       call allocmem(out%Ytn, out%N, out%nbDof, ii)
       if (choral_verb>1) write(*,*) &
            & "  Array Ytn shape                =", &
            & shape(out%Ytn)
    end if

    ! V(.,tn) 
    if (out%Vtn_period>0._rp) then
       ii = int( (out%T - out%t0) / out%Vtn_period) + 1
       call allocmem(out%Vtn, out%nbDof, ii)
       if (choral_verb>1) write(*,*) &
            & "  Array Vtn shape                =", &
            & shape(out%Vtn)
    end if

    ! V(xn,.) 
    if (out%Vxn_period>0._rp) then      

       ii = int( (out%T - out%t0) / out%Vxn_period) + 1
       call allocmem( out%Vxn, size(out%xn_list, 1)+1, ii )
       if (choral_verb>1) write(*,*) &
            & "  Array Vxn shape                =", &
            & shape(out%Vxn)

    end if

    ! activation times
    out%nv = def_nv(out%act_type)
    if (out%nv > 0) then

       call allocmem( out%act, out%nbDof) 
       call allocmem( out%V, out%nV, out%nbDof) 
       if (choral_verb>1) write(*,*) &
            &"  Array V   shape                =", &
            &shape(out%V)
       if (choral_verb>1) write(*,*) &
            &"  Array act shape                =", &
            & shape(out%act)

       out%act = -1e150_rp
       ! by default, the initial state is assumed to be depolarised
       out%V   = out%act_Vth - 1._rp 

    end if

    out%assembled = .true.

  end subroutine ode_output_assemble

  function def_nv(method) result(nv)
    integer :: nv, method

    select case(method)
    case(act_no)
       nv       = 0
    case(act_0)
       nv       = 1
    case(act_1)
       nv       = 2
    case(act_2)
       nv       = 3
    case(act_4)
       nv       = 5
    case default
       nv = 0
       call quit( "ode_output_mod: def_nv:&
            & unknown computation method for isochrones" )
    end select

  end function def_nv


  subroutine ode_output_print(out)
    type(ode_output), intent(in) :: out

    write(*,*)"ode_output_mod  : print"
    write(*,*)"  Starting time                  =", out%t0
    write(*,*)"  Final time                     =", out%T 
    write(*,*)"  Y(tn,xn) \in R**N,  N          =", out%N
    if (out%Ytn_period>0.0_rp) then
       write(*,*)"  Storing period for Y(.,tn)     =", out%Ytn_period
    end if
    if (out%Vtn_period>0.0_rp) then
       write(*,*)"  Storing period for V(.,tn)     =", out%Vtn_period
    end if
    if (out%Vtn_rec_prd>0.0_rp) then
       write(*,*)"  Saving  period for V(.,tn)     =", out%Vtn_rec_prd
       write(*,*)"  File name for saving V(tn,.)   = ",trim(out%Vtn_file)
    end if
    if (out%Vtn_plot) then
       write(*,*)"  Plot of V(tn,:)                =", out%Vtn_plot
    end if
    if (out%Vxn_period>0.0_rp) then
       write(*,*)"  Storing period for V(xn,.)     =", out%Vxn_period
       write(*,*)"  Number of saving nodes xn      =", size(out%xn_list, 1)
    end if
    if (out%Vxn_plot) then
       write(*,*)"  Plot of V(xn,.)                =", out%Vxn_plot
    end if
    if(out%act_type >= 0) then
       write(*,*)"  Activation time comp. method   =", out%act_type
       write(*,*)"  Threshold for act. time comp.  =", out%act_Vth
       write(*,*)"  Stop comp if fully depolarised =", out%stop_if_dep
       write(*,*)"  Saving activation times        =", out%act_rec
       write(*,*)"  Plot activation times          =", out%act_plot
       write(*,*)"  File name for saving act times = ",trim(out%act_file)
    end if
    select case(out%pos)
    case(pos_gnuplot)
       write(*,*)"  Graphical post process tool    = GNUPLOT"
    case(pos_gmsh)
       write(*,*)"  Graphical post process tool    = GMSH"
    case default
       write(*,*)"  Graphical post process tool    = unknown"
    end select
    write(*,*)"  Directory for output savings   = ",trim(out%outDir)

  end subroutine ode_output_print


  subroutine output_proc(out, time, s, final)
    class(ode_output), intent(inout) :: out

    real(RP)          , intent(in)  :: time
    type(ode_solution), intent(in)  :: s
    logical           , intent(out) :: final
    
    logical            :: exit_comp
    integer            :: ii
    real(RP)           :: v1, v2, tn
    character(LEN=100) :: fic

    exit_comp=.false.
    select case(out%verb)

    ! prints time, potential min max at each ms
    case(2) 

       if ( abs(anint(time)-time)<1e-4_rp ) then
          if (allocated(s%V)) then
             v1 = minval( s%V(:, s%v_i(1)) )
             v2 = maxval( s%V(:, s%v_i(1)) )
          else
             v1 = minval( s%Y(s%N, s%y_i(1), :) )
             v2 = maxval( s%Y(s%N, s%y_i(1), :) )
          end if
          write(*,*)"ode_output      : t, V_min, V_Max=", &
               & real(time,SP), real(v1,SP), real(v2,SP) 
       end if

    ! prints time, potential min max at each time step
    case (3) 
       if (allocated(s%V)) then
          v1 = minval( s%V(:, s%v_i(1)) )
          v2 = maxval( s%V(:, s%v_i(1)) )
       else
          v1 = minval( s%Y(s%N, s%y_i(1), :) )
          v2 = maxval( s%Y(s%N, s%y_i(1), :) )
       end if
       write(*,*)"ode_output      : t, V_m, V_M    =", &
               & real(time,SP), real(v1,SP), real(v2,SP) 

    ! prints full component Y min and max
    case (4) 
       if (allocated(s%Y)) then
          do ii=1, s%N
             write(*,*) "  Y(ii,:) ii, min, max          =", &
                  & int(ii,4),  &
                  & real( minval(s%Y(ii, s%y_i(1), :)), SP ), &
                  & real( maxval(s%Y(ii, s%y_i(1), :)), SP )
          end do
       else 
          v1 = minval( s%V(:, s%v_i(1)) )
          v2 = maxval( s%V(:, s%v_i(1)) )
          write(*,*)"ode_output      : t, V_m, V_M    =", &
               & real(time,SP), real(v1,SP), real(v2,SP) 
       end if
    end select

    ! storing of Y(., tn)
    if (out%Ytn_period>0._rp) then

       tn = out%t0 + re(out%Ytn_cpt) * out%Ytn_period
       if ( (time + real_tol) >= tn ) then

          out%Ytn_cpt = out%Ytn_cpt + 1
          call copy(out%Ytn(:,:,out%Ytn_cpt),  s%Y(:,s%y_i(1), :) )

          if (out%verb>0) then
             write(*,*)"ode_output      : Y(.,tn) stored = tn =", &
                  &real(tn,sp)
          end if

       end if
    end if

    ! storing of V(., tn)
    if (out%Vtn_period>0._rp) then

       tn = out%t0 + re(out%Vtn_cpt) * out%Vtn_period
       if ( (time + real_tol) >= tn ) then

          out%Vtn_cpt = out%Vtn_cpt + 1
          if (allocated(s%V)) then
             call copy(out%Vtn(:,out%Vtn_cpt), s%v(:,s%v_i(1)))
          else
             call copy(out%Vtn(:,out%Vtn_cpt), s%Y(s%N,s%y_i(1),:))
          end if
          if (out%verb>0) then
             write(*,*)"ode_output      : V(.,tn) stored = tn =", &
                  & real(tn,sp)
          end if

       end if
    end if


    ! storing of V(xn, .)
    if (out%Vxn_period>0._rp) then

       tn = out%t0 + re(out%Vxn_cpt) * out%Vxn_period
       if ( (time + real_tol) >= tn ) then

          out%Vxn_cpt = out%Vxn_cpt + 1
          if ( out%Vxn_cpt > size(out%Vxn,2) ) return
          ii = size(out%Vxn,1)

          out%Vxn(1, out%Vxn_cpt) = time
          if (allocated(s%V)) then
             out%Vxn(2:ii,out%Vxn_cpt) = s%v(out%xn_list, s%v_i(1))
          else
             out%Vxn(2:ii,out%Vxn_cpt) = s%Y(s%N,s%y_i(1),out%xn_list)
          end if
          if (out%verb>2) then
             write(*,*)"ode_output      : V(xn,.) stored = tn =", &
                  & real(tn,sp)
          end if

       end if
    end if

    ! Rcording of V(., tn) to file vn.dat
    if (out%Vtn_rec_prd>0._rp) then

       tn = out%t0 + re(out%Vtn_rec_cpt) * out%Vtn_rec_prd
       if ( (time + real_tol) >= tn ) then

          if (out%verb>0) then
             write(*,"(A35,F10.3)")&
                  & " ode_output      : V(., tn) saved = tn = ", tn
          end if
          
          call inttostring(fic, out%Vtn_rec_cpt)
          fic = trim(out%OUTDIR)//trim(out%Vtn_FILE)//trim(fic)//".dat" 

          if (allocated(s%V)) then
             call write(s%V(:, s%v_i(1)), trim(fic))
          else
             call write(s%Y(s%N, s%y_i(1), :), trim(fic))
          end if

          out%Vtn_rec_cpt = out%Vtn_rec_cpt + 1          
       end if

    end if


    ! Activation times
    if (out%nV>0) then

       ! update potential
       do ii=out%nV, 2, -1
          out%V(ii, :) = out%V(ii-1, :)
       end do
       if (allocated(s%V)) then
          call copy(out%V(1,:), s%V(:, s%v_i(1)) ) 
       else
          call copy(out%V(1,:), s%Y(s%N, s%y_i(1),:) ) 
       end if
       select case(out%act_type)
       case(act_0)
          call comp_act_0(out%act, out%V, time, out%act_Vth)

       case(act_1)
          call comp_act_1(out%act, out%V, time, out%dt, out%act_Vth)

       case(act_2)
          call comp_act_2(out%act, out%V, time, out%dt, out%act_Vth)

       case(act_4)
          call comp_act_4(out%act, out%V, time, out%dt, out%act_Vth)

       end select

       if (out%stop_if_dep) then

          if (minval(out%act)>-1e50_rp) then
             exit_comp = .true.
             if (out%verb>0) write(*,*)"ode_output      : &
                  &time           =", real(time,sp)
             if (out%verb>0) write(*,*)"    fully depolarised domain"
          end if
       end if

    end if

    !! Computation is over 
    if (final) then
       call ode_output_last_output(out)
       return
    end if

    final = exit_comp

  end subroutine output_proc


  subroutine ode_output_last_output(out)
    type(ode_output), intent(in) :: out

    if (choral_verb>1) write(*,*) &
         & "ode_output_mod  : last_output"

    if (out%act_rec ) then
       call write(out%act, trim(out%OUTDIR)//trim(out%ACT_FILE)//".dat")
    end if

    if (out%Vtn_plot) call plot_vtn(out)
    if (out%act_plot) call plot_act(out)
    if (out%Vxn_plot) call plot_vxn_gnuplot(out)

  end subroutine ode_output_last_output



  subroutine plot_vtn(out)
    type(ode_output), intent(in) :: out

    if (choral_verb>0) write(*,*) &
         & "ode_output_mod  : plot_vtn"

    if (out%pos==pos_gnuplot) call plot_vtn_gnuplot(out)
    if (out%pos==pos_gmsh   ) call plot_vtn_gmsh(out)

  end subroutine plot_vtn
  subroutine plot_vtn_gnuplot(out)
    type(ode_output), intent(in) :: out

    character(LEN=300) :: com

    if (out%dim>0) then       
       call write(out%X_h%dofCoord(1:out%dim,:), &
            &     trim(out%OUTDIR)//"X.dat", transpose=.true.)
    end if

    select case(out%dim)

    case(1)
       call inttostring(com, out%Vtn_rec_cpt-1)
       com =   'gnuplot -e "N=' // trim(com)  // '" ' & 
            &  // ' -e "root='''// trim(out%outDir) // '''" '&
            &  // trim(gnuplot_dir) // "plot_v_1d.gp"
       
       if (choral_verb>2) write(*,*) &
            & "  Shell command > "// trim(com) 
       call system(com)

    case(2)
       call inttostring(com, out%Vtn_rec_cpt-1)
       com =   'gnuplot -e "N=' // trim(com)  // '" ' & 
            &  // ' -e "root='''// trim(out%outDir) // '''" '&
            &  // trim(gnuplot_dir) // "2d_v.gp"
       
       if (choral_verb>2) write(*,*) &
            & "  Shell command > "// trim(com) 
       call system(com)

    case(3)
       call warning("ode_output_mod  : plot_v_gnuplot=&
            & dim 3 not implemented", 1)

    end select

  end subroutine plot_vtn_gnuplot

  subroutine plot_vtn_gmsh(out)
    type(ode_output), intent(in) :: out

    real(RP), dimension(out%nbDof) :: v
    character(LEN=300) :: com
    integer  :: ii
    real(RP) :: time, min, max, val

    if (out%dim<=1) then
       call warning('ode_output_mod  : plot_v_gmsh    = &
            & dim 2, 3 only', 1)
       return
    end if
    
    call write(out%X_h, trim(out%OUTDIR)//'v_plot.msh', 'gmsh')

    min = 0.0_rp
    max = 0.0_rp
    do ii=0, out%Vtn_rec_cpt-1
       call inttostring(com, ii)
       com = trim(out%OUTDIR)//trim(out%Vtn_FILE)//trim(com)//".dat"
       call read(v, trim(com))
       
       val = minval(v)
       if (val < min ) min = val
       val = maxval(v)
       if (val > max ) max = val

       time = out%t0 + out%Vtn_rec_prd*re(ii)

       call gmsh_addview(out%X_h, v, &
            & trim(out%OUTDIR)//'v_plot.msh', &
            & 'V', time, ii)
    end do

    !! add an option file
    !! set the range
    open(unit=60,file=trim(out%OUTDIR)//'v_plot.msh.opt')   
    write(60,*) "View[0].RangeType = 2;"
    write(60,*) "View[0].CustomMax = ", REAL(max, SP),";"
    write(60,*) "View[0].CustomMin = ", REAL(min, SP),";"
    close(60)

    com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'
    com = trim(com)//' '//trim(out%OUTDIR)//'v_plot.msh'
    if (choral_verb>2) write(*,*) &
         & "  Shell command > "// trim(com) 
    call system(com)

  end subroutine plot_vtn_gmsh



  subroutine plot_vxn_gnuplot(out)
    type(ode_output), intent(in) :: out

    integer  :: ii, nn
    character(LEN=300) :: com, fic

    !! Recording file
    !!
    fic = trim(out%outDir)//trim(out%vxn_file)

    !! print the transmembrane potential V(xi, tn) to a file
    open(unit=60,file=trim(fic))   
    do ii=1, out%Vxn_cpt
       write(60,*) out%Vxn(:,ii)
    end do
    close(60)

    do nn=2, size(out%Vxn,1)

       !! plot
       call inttostring(com, nn)
       com = 'gnuplot -e "fic=''' // trim(fic)//'''"  '  &
            &  // ' -e col=' // trim(com)//'  '&
            &  // trim(gnuplot_dir)//"plot_v_0d.gp"
  
       if (choral_verb>2) write(*,*) &
            & "    Shell command $> "//trim(com)
       call system(com)

    end do

  end subroutine plot_vxn_gnuplot


  subroutine plot_act(out)
    type(ode_output), intent(in) :: out

    if (choral_verb>0) write(*,*) &
         & "ode_output_mod  : plot_act"
    if (out%pos==pos_gnuplot) call plot_act_gnuplot(out)
    if (out%pos==pos_gmsh   ) call plot_act_gmsh(out)

  end subroutine plot_act

  subroutine plot_act_gnuplot(out)
    type(ode_output), intent(in) :: out

    character(LEN=300) :: com

    call write(out%X_h%dofCoord(1:out%dim,:), &
         &     trim(out%OUTDIR)//"X.dat", transpose=.true.)

    select case(out%dim) 
    case(1)
       com =   'gnuplot -e "root=''' // trim(out%outDir) // '''" '&
            &  // trim(gnuplot_dir)   // "plot_act_1d.gp"

       if (choral_verb>2) write(*,*) &
            & "  Shell command > "// trim(com) 
       call system(com)

    case(2)
       com =   'gnuplot -e "root=''' // trim(out%outDir) // '''" '&
            &  // trim(gnuplot_dir)   // "2d_act.gp"

       if (choral_verb>2) write(*,*) &
            & "  Shell command > "// trim(com) 
       call system(com)

    case(3)
       write(*,*)"ode_output_mod  : plot_act: &
            &dim 3 not implemented" 

    end select

  end subroutine plot_act_gnuplot

  subroutine plot_act_gmsh(out)
    type(ode_output), intent(in) :: out

    real(RP), dimension(out%nbDof) :: v
    character(LEN=300) :: com
 
    if (out%dim==1) then
       call warning('ode_output_mod  : plot_act_gmsh  = &
            & dim 1 not implemented', 1)
       return
    end if

    call write(out%X_h, trim(out%OUTDIR)//'act_plot.msh', 'gmsh')

    com = trim(out%OUTDIR)//trim(out%ACT_FILE)//".dat"
    call read(v, trim(com))

    call gmsh_addview(out%X_h, v, &
            & trim(out%OUTDIR)//'act_plot.msh', &
            & 'Activation times in ms', 0._rp, 0)

    com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'  
    com = trim(com)//' '//trim(out%OUTDIR)//'act_plot.msh'
    if (choral_verb>2) write(*,*) &
         & "  Shell command > "// trim(com) 
    call system(com)

  end subroutine plot_act_gmsh



  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
  ! !!
  ! !!
  ! !!           ACTIVATION TIMES COMPUTATION
  ! !!
  ! !!
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine comp_act_0(act, V, time, Vth)
    real(RP), dimension(:)  , intent(out) :: act
    real(RP), dimension(:,:), intent(in)  :: V
    real(RP)                , intent(in)  :: time, Vth

    integer :: ii

    do ii=1, size(act,1)

       if (act(ii)<-1e50_rp) then
          if (v(1,ii)>vth) act(ii) = time
       end if

    end do

  end subroutine comp_act_0





  subroutine comp_act_1(act, V, time, dt, Vth)
    real(RP), dimension(:)  , intent(out) :: act
    real(RP), dimension(:,:), intent(in)  :: V
    real(RP)                , intent(in)  :: time, dt, Vth

    integer  :: ii
    real(RP) :: v2, v1

    do ii=1, size(act, 1)

       if (act(ii)<-1e50_rp) then
          v1 = v(1,ii)
          if (v1>vth) then
             v2 = v(2,ii)
             act(ii) = time + ( vth - v1 )/( v1 - v2 ) * dt
          end if
       end if

    end do

  end subroutine comp_act_1


  subroutine comp_act_2(act, V, time, dt, Vth)
    real(RP), dimension(:)  , intent(out) :: act
    real(RP), dimension(:,:), intent(in)  :: V
    real(RP)                , intent(in)  :: time, dt, Vth

    integer  :: ii
    real(RP) :: v3, v2, v1, a, b, delta

    do ii=1, size(act, 1)

       if (act(ii)<-1e50_rp) then
          v1 = v(1,ii)
          if (v1>vth) then
             v2 = v(2,ii)
             v3 = v(3,ii)

             a     = (v1 - 2._rp*v2 + v3)*0.5_rp
             b     = (v1 - v3)*0.5_rp          
             delta = b**2 - 4._rp*a*(v2-vth)

             ! highest root
             a = ( sqrt(delta) - b ) / (2._rp * a)

             act(ii) = time + dt * (a - 1._rp)

          end if
       end if

    end do

  end subroutine comp_act_2



  subroutine comp_act_4(act, V, time, dt, Vth)
    real(RP), dimension(:)  , intent(out) :: act
    real(RP), dimension(:,:), intent(in)  :: V
    real(RP)                , intent(in)  :: time, dt, Vth

    integer  :: ii
    real(RP) :: v5, v4, v3, v2, v1, a, b, c, d, e, tol

    tol = real_tol * 1e4_rp

    do ii=1, size(act, 1)

       if (act(ii)<-1e50_rp) then
          v2 = v(2,ii)
          if (v2>vth) then
             v1 = v(1,ii)
             v3 = v(3,ii)
             v4 = v(4,ii)
             v5 = v(5,ii)

             a = v1 - 4._rp*v2 + 6._rp*v3 - 4._rp*v4 + v5
             a = a/24._rp

             c =-v1 + 16._rp*v2 - 30._rp*v3 + 16._rp*v4 - v5
             c = c/24._rp

             b = v1 - 2._rp*v2 + 2._rp*v4 - v5
             b = b/12._rp

             d =-v1 + 8._rp*v2 - 8._rp*v4 + v5 
             d = d/12._rp

             e = v3 - vth

             a = local_newton()

             act(ii) = time + dt * (a - 2._rp)

          end if
       end if

    end do

    contains

      function p(tau)
        real(RP)             :: p
        real(RP), intent(in) :: tau 

        p = a*tau**4 + b*tau**3 + c*tau**2 + d*tau + e

      end function p

      function d_p(tau)
        real(RP)             :: d_p
        real(RP), intent(in) :: tau 

        d_p = 4._rp*a*tau**3 + 3._rp*b*tau**2 + 2._rp*c*tau + d

      end function d_p

      function local_newton() result(s)
        real(RP) :: s

        real(RP):: ps
        integer :: cpt

        cpt = 0
        s  = 0._rp   !(vth-v3) / (v2-v3)
        ps = e       !p(s)

        do while( abs(ps) > tol )

           s  = s - ps / d_p(s)
           ps = p(s)

           cpt = cpt + 1
           if (cpt==1000) then
              call warning('ode_output_mod: comp_isoch_4 =&
                   & Newton error', 1)
              return
           end if
        end do

      end function local_newton

  end subroutine comp_act_4


  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
  ! !!
  ! !!
  ! !!           WAVE FRONT ANALYSIS
  ! !!
  ! !!
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function celerity_l2_dist(u1, u2, X_h, qdm, weight) result(dist)
    real(RP)                           :: dist
    real(RP), dimension(:), intent(in) :: u1, u2
    type(fespace)         , intent(in) :: X_h
    type(quadmesh)        , intent(in) :: qdm
    procedure(r3tor)                   :: weight
        
    dist = integ(e, u1, u2, x_h, qdm)
    dist = sqrt(dist)

  contains

    function e(x, v1, v2)
      real(RP) :: E
      real(RP), dimension(3), intent(in) :: x, v1, v2

      real(RP), dimension(3) :: yy

      yy = v1 / sum(v1*v1) - v2 / sum(v2*v2)
      e  = sum(yy*yy)
      e  = weight(x) * e

    end function e

  end function celerity_l2_dist


  function celerity_l2_norm(u1, X_h, qdm, weight) result(dist)
    real(RP)                             :: dist
    real(RP)  , dimension(:), intent(in) :: u1
    type(fespace)           , intent(in) :: X_h
    type(quadmesh)          , intent(in) :: qdm
    procedure(r3tor)                     :: weight
    

    dist = integ(e, u1, u1, x_h, qdm)
    dist = sqrt(dist)

  contains

    function e(x, v1, v2)
      real(RP) :: E
      real(RP), dimension(3), intent(in) :: x, v1, v2

      real(RP), dimension(3) :: yy

      yy = v1 / sum(v1*v1) 
      e  = sum(yy*yy)
      e  = weight(x) * e

    end function e

  end function celerity_l2_norm
  

end module ode_output_mod