!  Fortran-90 single precision example of solving 3D Helmholtz problem in a
!  parallelepiped domain using MKL Poisson Library
!
!*******************************************************************************

! ED: the following files must be added in order to include mkl modules
include 'mkl_dfti.f90'
include 'mkl_poisson.f90'

program Helmholtz_3D_single_precision

! Include modules defined by mkl_poisson.f90 and mkl_dfti.f90 header files
use mkl_poisson

implicit none

! Number of mesh intervals along x, y and z axis
integer nx, ny, nz
parameter(nx=201, ny=201, nz=201)

! Parameters of the poisson equation
real pi
parameter(pi=3.14159265358979324E0)
real eps_0 ! permittivity (C^2/(N m^2))
parameter(eps_0=8.85418782E-12)
real charge ! electron charge (C)
parameter(charge=1.60217662E-19)
integer q_ix, q_iy, q_iz ! charge position (in voxel units)
parameter(q_ix=int((nx+1)/2),q_iy=int((ny+1)/2),q_iz=int((nz+1)/2))
real eps
parameter(eps=1.0E-20)
real a_factor ! allows to transform length units to angstrom
parameter(a_factor=1.0E-10)

integer stat ! MKL routines completion status

! Boundary coordinates of the domain along x, y and z axis
real ax, bx, ay, by, az, bz

! Integer and single-precision data to be used by Fast Helmholtz Solver
integer ipar(128)
real spar(13*(nx+ny)/2+9)

! Note that proper packing of data in right-hand side array f and boundary 
! arrays bd_* is automatically provided by the following
! declarations of arrays
real f(nx+1,ny+1,nz+1), u(nx+1,ny+1,nz+1)
real bd_ax(ny+1,nz+1), bd_bx(ny+1,nz+1), bd_ay(nx+1,nz+1), bd_by(nx+1,nz+1)
real bd_az(nx+1,ny+1), bd_bz(nx+1,ny+1)

real q ! constant Helmholtz coefficient

! Data structures used by the Intel MKL FFT interface
type(DFTI_DESCRIPTOR), pointer :: xhandle, yhandle

character(6) BCtype ! type of boundary conditions on each boundary

! File unit for output results
integer iout
parameter(iout=8)

integer ix, iy, iz, i
real lx, ly, lz, hx, hy, hz, xi, yi, zi, c1

! Open output file and print header
open(unit=iout,file="output.txt",status='unknown',form='formatted')
write(iout,'(A,I3,A,I3)') 'iy= ',q_iy,', iz= ',q_iz
write(iout,'(A,A,A)') 'ix ', 'Calculated ', 'True '

! Printing the header for the example
print *, ''
print *, ' MKL Poisson Library - potential due to a point charge'
print *, ' *****************************************************'
print *, ''
print *, ' This program gives the solution of 3D Poisson problem'
print *, ' with the equation -u_xx-u_yy-u_zz=f(x,y,z), '
print *, ' f(x,y,z)=rho(x,y,z)/eps_0,'
print *, ' with Dirichlet boundary conditions '
print *, ' --------------------------------------------------------------------'
print *, ' In general, the error should be of order ', &
            'O(1.0/nx^2+1.0/ny^2+1.0/nz^2)'
print '(1x,a,I2)', ' For this example, the value of nx=ny=nz is ', nx
print *, ' The approximation error should be of order 0.5E+0,',&
         ' if everything is OK'
print *, ' --------------------------------------------------------------------'
print *, '                      SINGLE PRECISION COMPUTATIONS                  '
print *, ' ===================================================================='
print *, ''

! Defining the parallelepiped domain for 3D Helmholtz Solver
ax=-10.0E0
bx= 10.0E0
ay=-10.0E0
by= 10.0E0
az=-10.0E0
bz= 10.0E0

! Setting the coefficient q to 0, needed to use Helmholtz solver to solve 
! Poisson problem
q=0.0E0

! Computing the mesh size hx in x-direction
lx=bx-ax
hx=lx/nx
! Computing the mesh size hy in y-direction
ly=by-ay
hy=ly/ny
! Computing the mesh size hx in z-direction
lz=bz-az
hz=lz/nz

write (*,'(A,F5.3)') ' The value of the mesh size hx=hy=hz is ', hx
write (*,'(A,F7.3,A,F7.3)') ' Domain boundaries in x: ', ax, ', ', bx
write (*,'(A,F7.3,A,F7.3)') ' Domain boundaries in y: ', ay, ', ', by
write (*,'(A,F7.3,A,F7.3)') ' Domain boundaries in z: ', az, ', ', bz

! Filling in the values of the TRUE solution
! u(x,y,z) = q / (4*pi*eps_0 r)
! in the mesh points into the array u
! Filling in the right-hand side
! f(x,y,z) = rho(x,y,z)/eps_0
! in the mesh points into the array f.
! We choose the right-hand side to correspond to the TRUE solution of Helmholtz 
! equation.

! Here we are using the mesh sizes hx, hy, and hz computed before to compute 
! the coordinates (xi,yi,zi) of the mesh points. The
! coordinates are placed in the middle point of each voxel
do iz=1,nz+1
   do iy=1,ny+1
      do ix=1,nx+1
         xi=ax + hx*(ix-0.5)
         yi=ay + hy*(iy-0.5)
         zi=az + hz*(iz-0.5)
         u(ix,iy,iz)=charge/(4*pi*eps_0)/sqrt(xi**2+yi**2+zi**2+eps)/a_factor
         if((iz.eq.q_iz).and.(iy.eq.q_iy).and.(ix.eq.q_ix)) then
            f(ix,iy,iz)=charge/(hx*hy*hz)/eps_0/(a_factor)
         else
            f(ix,iy,iz)=0.0E0
         end if
      enddo
   enddo
enddo

write (*,'(A,I3,A,I3,A,I3,A)') ' Coordinates of the point charge: (', q_ix, &
', ', q_iy, ', ', q_iz, ')'

! Setting the type of the boundary conditions on each surface of the 
! parallelepiped domain:
! On all boundaries Dirichlet boundary condition will be used
BCtype = 'DDDDDD'

! Setting the values of the boundary function G(x,y,z) that is equal to the 
! TRUE solution in the mesh points laying on Dirichlet
! boundaries
do iy = 1,ny+1
   do iz = 1,nz+1
      bd_ax(iy,iz) = charge/(4*pi*eps_0)/sqrt(ax**2+yi**2+zi**2+eps)/a_factor
      bd_bx(iy,iz) = charge/(4*pi*eps_0)/sqrt(bx**2+yi**2+zi**2+eps)/a_factor
   enddo
enddo
do ix = 1,nx+1
   do iz = 1,nz+1
      bd_ay(ix,iz) = charge/(4*pi*eps_0)/sqrt(xi**2+ay**2+zi**2+eps)/a_factor
      bd_by(ix,iz) = charge/(4*pi*eps_0)/sqrt(xi**2+by**2+zi**2+eps)/a_factor
   enddo
enddo
do ix = 1,nx+1
   do iy = 1,ny+1
      bd_az(ix,iy) = charge/(4*pi*eps_0)/sqrt(xi**2+yi**2+az**2+eps)/a_factor
      bd_bz(ix,iy) = charge/(4*pi*eps_0)/sqrt(xi**2+yi**2+bz**2+eps)/a_factor
   enddo
enddo

! Initializing ipar array to make it free from garbage
do i=1,128
   ipar(i)=0
enddo

! Initializing simple data structures of Poisson Library for 3D Helmholtz Solver
call s_init_Helmholtz_3D(ax, bx, ay, by, az, bz, nx, ny, nz, BCtype, q, ipar, &
                        spar, stat)
if (stat.ne.0) goto 999

! Initializing complex data structures of Poisson Library for 3D Helmholtz
! Solver
! NOTE: Right-hand side f may be altered after the Commit step. If you want to 
! keep it, you should save it in another memory
! location!
call s_commit_Helmholtz_3D(f, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, &
                            xhandle, yhandle, ipar, spar, stat)
if (stat.ne.0) goto 999

! Computing the approximate solution of 3D Helmholtz problem
! NOTE: Boundary data stored in the arrays bd_ax, bd_bx, bd_ay, bd_by, bd_az, 
! bd_bz should not be changed between the Commit step
! and the subsequent call to the Solver routine! Otherwise the results may be 
! wrong.
call s_Helmholtz_3D(f, bd_ax, bd_bx, bd_ay, bd_by, bd_az, bd_bz, xhandle, &
                    yhandle, ipar, spar, stat)
if (stat.ne.0) goto 999

! Cleaning the memory used by xhandle and yhandle
call free_Helmholtz_3D(xhandle, yhandle, ipar, stat)
if (stat.ne.0) goto 999
! Now we can use xhandle and yhandle to solve another 3D Helmholtz problem

! Printing the results
write(*,9)  nx
write(*,10) ny
write(*,11) nz

! Watching the error in the plane x=q_ix
ix=q_ix
c1 = 0.0
do iy=1,ny+1
   do iz=1,nz+1
!      write(*,12) ax + hx*(ix-0.5), ay + hy*(iy-0.5), az + hz*(iz-0.5), &
!                  f(ix,iy,iz)-u(ix,iy,iz)
      if (abs(f(ix,iy,iz)-u(ix,iy,iz)).ge.c1) c1 = abs(f(ix,iy,iz)-u(ix,iy,iz))
   enddo
enddo
print *, ''

! Print calculated and exact solution to Poisson equation
iy=q_iy
iz=q_iz
do ix=1,nx+1
    write(iout,*) ax + hx*(ix-0.5), f(ix,iy,iz), u(ix,iy,iz)
enddo

write(*,*) 'Maximum error found in plane x=hx : ', c1
if (c1.ge.0.5E+0) then
   print *, 'The computed solution seems to be inaccurate.'
   goto 999
endif

! Free MKL memory if any was allocated
call mkl_free_buffers
! Success message to print if everything is OK
print *, ' Single precision 3D Helmholtz example has successfully PASSED'
print *, ' through all steps of computation!'
stop 0
! Failure message to print if something went wrong
999 print *, &
'Single precision 3D Helmholtz example FAILED to compute the solution...'
    ! Free MKL memory if any was allocated
    call mkl_free_buffers
    stop 1

9     format(1x,'The number of mesh intervals in x-direction is nx=',I3)
10    format(1x,'The number of mesh intervals in y-direction is ny=',I3)
11    format(1x,'The number of mesh intervals in z-direction is nz=',I3)
12    format(1x,'In the mesh point (',F5.3,',',F5.3,',',F5.3,') the error', &
            'between the computed and the true solution is equal to ',sE10.3)

! End of the example code
close(iout)

end
