! ***************************************************** !
! ****************** ftcs-ref-ref.f90 ***************** !
! ***************************************************** !
!
!    Leonardo Dagdug, Ivan Pompa-García, Jason Peña
!
! From the book:
! **************************************************** !
!      Diffusion Under Confinement:
!               A Journey Through Counterintuition    
! **************************************************** !

! **************************************************** !
! This file contains the source code to numerically 
! solve the 1D diffusion equation for a channel
! of length L.
! It has two reflecting points at x=0 and x=L.
! 
! The initial condition is such that
!                 p(x,t=0) = \delta(x-x0),
! where x0 is the initial position (default = 1/2).
!
! The solver uses the Forward Time-Centered Space
! (FTCS) Finite Difference Method (FDM).
!
! As the FTCS is stable only with:
!         D \Delta t / \Delta x^2 <= 1/2,
! those parameters must be chosen carefully.
!
! **************************************************** !

! **************************************************** !
! To compile with GFortran, you must include the
! helpers.f90 module.
!         gfortran helpers.f90 ftcs_ref_ref.f90
! Then you can run the program:
!         ./a.out
!
! After running the program, a file is generated:
!   - «ftcs-aa.dat», a file with time, position,
!                    and concentration
! It can be plotted with any software
! of your choice. To do this using gnuplot, execute:
!  gnuplot -p -e '
!     set title "FTCS FDM for the Diffusion Equation";
!     set xlabel "x (position)" ;
!     set ylabel "c(x,t) (concentration)" ;
!times = system("awk '"'"'{print $1}'"'"' ftcs-aa.dat | sort -u");
!     plot for [t in times] "<awk '"'"'$1==".t."'"'"' ftcs-aa.dat"
!       u 2:3 with linespoints title "t=".sprintf("%.3f", t*1.0)'
!
! -p means that the plot must persist until the
!    user closes it explicitly.
! -e is used to specify the instructions
!    for using  gnuplot directly in the commandline.
! **************************************************** !

program ftcsrefref

! Load the helpers module which contains functions,
! constants, and more...
use helpers

! Mandatory declaration of data type variables
! and constants.
implicit none

! **************************************************** !
!               Declaration of constants
! **************************************************** !

!!! FDM parameters
! Number of pieces in the partition of «x»
integer, parameter :: NX = 10
! Number of pieces in the partition of «t»
integer, parameter :: NT = 30
! Total system evolution time
! (assuming that it starts at t=0)
real(kind=dp), parameter :: TF = 0.1_dp
! Start of channel
real(kind=dp), parameter :: XI = 0.0_dp
! End of channel
real(kind=dp), parameter :: XF = 1.0_dp
! Initial Condition continous position
real(kind=dp), parameter :: X0 = 0.5_dp
! Temporal step size
real(kind=dp), parameter :: DT = TF / NT
! Spatial step size
real(kind=dp), parameter :: DX = (XF - XI) / NX
! Diffusion constant (bulk)
real(kind=dp), parameter :: D0 = 1.0_dp
! «r» factor
real(kind=dp), parameter :: R = D0 * DT / DX**2
! Discrete initial condition position
! Substract 1 from the index because our array
! will be defined to start at 0.
integer, parameter :: J0 = int(X0 / DX + 1)


! **************************************************** !
!                Declaration of variables
! **************************************************** !

! General counter
integer :: i
! Time counter
integer :: n
! Spatial counter
integer :: j
! Current continous position
real(kind=dp) :: x
! Current continous time
real(kind=dp) :: t
! Array to hold the current time concentrations «c»
real(kind=dp) :: c(NX+1) = 0.0_dp
! Array to hold the previous time concentrations «c»
real(kind=dp) :: c0(NX+1) = 0.0_dp

! A unit number for the file where 
! time, position, and concentration data
! will be saved.
! Fortran >= 2003 will
! assign this automatically in the open call.
integer(kind=i32) :: FUNIT



! **************************************************** !
!              Main calculation program
! **************************************************** !

! Open a file to write the data
open(newunit=FUNIT, file='ftcs-aa.dat')

! Show a warning if R > 1/2
if(R > 0.5) then
  print *, 'r > 0.5. Under this condition the method ' &
            // 'is not stable! Try to choose ' &
            // ' different values of DX and DT.'
end if

! Set the discretized Initial Condition
c(J0) = 1.0_dp / DX

! Print the values for the time zero
do j=1, NX+1
  write(FUNIT, *) 0.0_dp, (j-1) * DX, c(j)
end do

! The loop iterates over each time step until NT is
! reached.
! In other words, it walks along every point on the mesh
! for the respective amount of time.
! Keeping the notation of the book, the current time 
! step is held in «n».
do n=1, NT
  ! Obtain the current continous time
  t = n * DT

  ! Save the previous time concentration before
  ! calculating the next one.
  c0 = c

  ! The first element of «c» is calculated separately
  ! because it is a boundary and uses another equation.
  c(1) = c0(1) + 2.0_dp * R * ( c0(2) - c0(1) )

  ! Write out the time, space position, and
  ! concentration
  write(FUNIT, *) t, x, c(j)

  ! Now we go through every point on the mesh for the
  ! spatial axis, holding the time «n» constant.
  ! The final point in j=Nx is a boundary and it is
  ! excluded.
  do j=2, NX
    ! The continous position is calculated
    x = (j-1) * DX

    ! Obtain the «c»
    c(j) = c0(j) &
            + R * ( c0(j+1) - 2.0_dp * c0(j) + c0(j-1) )

    ! Write out the time, space position, and «c»
    write(FUNIT, *) t, x, c(j)
  end do ! End of the spatial loop for time «n»

  x = (j-1) * DX

  ! The last element (a boundary) is calculated outside
  ! the loop for the inner points of the mesh, because
  ! it obeys a different equation.
  c(j) = c0(j) + 2.0_dp * R * ( c0(j-1) - c0(j) )

  ! Write out the time, space position, and «c»
  write(FUNIT, *) t, x, c(j)
end do ! Ends the loop for all timepoints.

! Close the file.
close(FUNIT)

!!! At this point, the calculation is done.
! As a reminder, print out some numbers:
print *
print *, TAB // '=== FTCS FDM finished with params ==='
! ! The number of slides for space
print *, TAB // 'Nx = ' // nstr(NX)
! ! The spatial step
print *, TAB // 'dx = ' // nstr(DX, 3)
! ! The number of slides for time
print *, TAB // 'Nt = ' // nstr(NT)
! ! The time step, too
print *, TAB // 'dt = ' // nstr(DT, 3)
! ! r factor
print *, TAB // 'r = ' // nstr(R, 2)

end program ftcsrefref