! **************************************************** !
! ***************** btcs-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 Backward Time-Centered Space
! (BTCS) Finite Difference Method (FDM).
!
! The BTCS is inconditionally stable, so the increment
! values do not need to be adapted to a
! specific rule as in FTCS.
!
! **************************************************** !

! **************************************************** !
! To compile with GFortran, you must include the
! helpers.f90 module, and also indicate that the
! program needs to be linked with LAPACK:
!         gfortran helpers.f90 btcs_ref_ref.f90 -llapack
!
! LAPACK must be installed in your operating system's
! standard location.
!
! Then you can run the program:
!         ./a.out
!
! After running the program, a file is generated:
!   - «btcs-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 "BTCS FDM for the Diffusion Equation";
!     set xlabel "x (position)" ;
!     set ylabel "c(x,t) (concentration)" ;
!times = system("awk '"'"'{print $1}'"'"' btcs-aa.dat | sort -u");
!     plot for [t in times] "<awk '"'"'$1==".t."'"'"' btcs-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 btcsrefref

! Load the helpers module, which contains functions,
! constants, ...
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
! Cintinuous initial condition 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
! Column vector «b»
real(kind=dp) :: b(NX+1) = 0.0_dp
! Matrix A
real(kind=dp) :: A(NX+1,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

!!! Auxilliary variables for LAPACK
integer :: ipiv(NX+1) = 0
integer :: info = 0
integer :: iter = 0
real(kind=dp) :: work(NX+1) = 0.0_dp
real(kind=dp) :: swork((NX+1)*(NX+1+1)) = 0.0_dp


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

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

! Populate matrix «A» with the appropiate values.
do j=2, NX
  A(j,j-1) = -r
  A(j,j+1) = -r
  A(j,j) = 1 + 2*r
end do

! Set the boundary conditions.
A(1,2) = -2*r
A(NX+1,NX) = -2*r
A(1,1) = 1 + 2*r
A(NX+1,NX+1) = 1 + 2*r

! Set the discretized initial condition.
c(J0) = 1.0_dp / DX

! Print out the values for 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

  ! Vector «b» must be set to the values of
  ! «c» for the previous time step.
  b = c

  ! The equation «Ac=b» is solved using LAPACK.
  ! The arguments are (in order):
    ! Number of linear equations
    ! Number of columns in «b»
    ! «A»
    ! Leading dimension of A
    ! Pivot integer
    ! «b», also the output if the operation was 
    ! successful.
    ! Leading dimension of «b»: NX+1
    ! The result: «c»
    ! Leading dimension of the result: «NX+1»
    ! Residual hold vectors: «NX+1»
    ! To use a single precision matrix: «(NX+1)*(NX+1+1)»
    ! An indicator if iterative methods were used.
  call dsgesv(NX+1, 1, A, NX+1, ipiv, b, NX+1, &
                      c, NX+1, work, swork, iter, info)

  if(info /= 0) then
    print *, 'An error has ocurred solving the' &
              // ' problem Ac=b, verify!'

    call exit(info)
  end if

  ! Set the BCs
  c(1) = c(2)
  c(NX+1) = c(NX)

  ! Print out the whole «c» for this specific time step.
  do j=1, NX+1
    write(FUNIT, *) t, (j-1) * DX, c(j)
  end do
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 // '=== BTCS 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 btcsrefref