Main.F90

!!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
!


Program dft

Use parameters, Only: np,nc,kbol
Use basic_variables, Only: ncomp,t,p,ensemble_flag,num,parame,eos,pol,xif,compna,mm, &
                           surfactant, wif_surfactant
Use eos_variables, Only: dhs,mseg,eta,dd_term,dq_term,qq_term
Use eos_constants, Only: ap,bp,dnm
Use mod_dft, Only: box,dzp,fa,zp,fa_disp,ab_disp,pbulk
Use vle_var, ONLY: tc,pc,rhob,density
USE dft_fcn_module, ONLY: chempot_total

!PETSc modules
Use petscmanagement
Use f90moduleinterfaces
Use global_x      !(snes,ngrid,ngp,x,r,timer)


Implicit None

#include <finclude/petscsys.h>


petscerrorcode     :: ierr
type(userctx)     :: user
petscbool          :: flg
INTEGER            :: i
character(80)      :: filename=''
REAL               :: zges(np)
REAL               :: dhs_save(nc)
REAL               :: cif(nc)
REAL               :: psi_factor
REAL               :: x_surfactant
INTEGER            :: position_surfactant

!external subroutines associated with Solver
external forminitialguess
external formfunction


      call petscinitialize(petsc_null_character,ierr)
      call mpi_comm_rank(petsc_comm_world,user%rank,ierr)
      call mpi_comm_size(petsc_comm_world,user%num_procs,ierr)


      filename='./in.txt'
      CALL file_open(filename,77)       ! open input file
      READ (77,*) t                     ! read temperature
      READ (77,*) ncomp                 ! read number of components in the system
      Do i = 1,ncomp                    ! read component names
        READ (77,*) compna(i)
      End Do
      Do i = 1,ncomp                    ! read component overall molar concentrations
       READ (77,*) xif(i)
     End Do

     !check whether surfactant present in system
     surfactant = .false.
     If(ncomp > 1) Then
        Do i = 1,ncomp
          If(trim(compna(i)) == 'surfactant') Then
             surfactant = .true.
             !reduce number of components
             ncomp = ncomp - 1
             x_surfactant = xif(i)
             position_surfactant = i
          End If
        End Do
     End If

     If(surfactant) Then
        !calculate new molar concentrations without surfactant
        If(position_surfactant == 1) Then
            Do i = 2,ncomp+1
               xif(i-1) = xif(i)
               compna(i-1) = compna(i)
            End Do
            xif(1:ncomp) = xif(1:ncomp) / sum(xif(1:ncomp))
        Else If(position_surfactant == ncomp + 1) Then
            xif(1:ncomp) = xif(1:ncomp) / sum(xif(1:ncomp))
        Else
            Do i = 1,ncomp+1
              If(i>position_surfactant) Then
                  xif(i-1) = xif(i)
                  compna(i-1) = compna(i)
              End If
            End Do
            xif(1:ncomp) = xif(1:ncomp) / sum(xif(1:ncomp))
        End If

     End If

     !In cases where only 0.0 are passed for xif, xif(1:ncomp) = xif(1:ncomp) / sum(xif(1:ncomp)) leads to NAN results
     If(xif(1) /= xif(1)) Then
        If(ncomp == 2) xif(1:ncomp) = 0.
        If(ncomp >  2) xif(1:ncomp) = 1./Real(ncomp)
     End If


!       !calculate molar fractions from molar concentrations
!       xif(1:ncomp) =cif(1:ncomp) / sum(cif(1:ncomp))


! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  GENERAL SIMULATION SET UP
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  num = 0                     ! (num=0: using analytical derivatives of EOS)
                              ! (num=1: using numerical derivatives of EOS)
                              ! (num=2: White's renormalization group theory)
  IF ( num /= 0 ) CALL set_default_eos_numerical

  eos = 1
  pol = 1

  p = 1.000e05

  CALL para_input             ! retriev pure comp. parameters

  ensemble_flag = 'tp'        ! this specifies, whether the eos-subroutines
                              ! are executed for constant 'tp' or for constant 'tv'



!calculate feed mass fraction of surfactant
If(surfactant) wif_surfactant = 0.01!x_surfactant*2655.24078 / (sum( xif(1:ncomp)*mm(1:ncomp) ) + x_surfactant*2655.24078 )


  CALL eos_const (ap, bp, dnm)

      dd_term = 'GV'              ! dipole-dipole term of Gross & Vrabec (2006)
      qq_term = 'JG'              ! quadrupole-quadrupole term of Gross (2005)
      dq_term = 'VG'              ! dipole-quadrupole term of Vrabec & Gross (2008)

  CALL vle_mix(rhob,density,chempot_total,user) !user is passed so user%rank is known and only node 0 prints out the reslts of the VLE calculation

  !determine critical point
    num = 0
    dhs_save = dhs               !needed because subroutine CRIT_POINT_MIX changes the value
    CALL crit_point_mix(tc,user) !of the global variable dhs! In cases where Tc doesnt converge
    dhs = dhs_save               !dhs can be NAN after CRIT_POINT_MIX!!
    !chance to overwite NAN results
!     IF(user%rank == 0) THEN
!         WRITE (*,*) 'Give final value of Tc:'
!         READ (*,*) tc
!     END IF
!     CALL MPI_Bcast(tc,1,MPI_DOUBLE_PRECISION,0,PETSC_COMM_WORLD,ierr)





   num = 1
   CALL set_default_eos_numerical

   !set default values (overwritten from makefile)
   box   = 45.0                !box length [A]
   ngrid = 614                 !grid points

   !overwrite box size and ngrid if options are set in makefile
   call petscoptionsgetreal(petsc_null_character,'-box',box,flg,ierr)
   call petscoptionsgetint(petsc_null_character,'-nx',ngrid,flg,ierr)


   dzp   = box / REAL(ngrid)   !grid spacing [A]
   Allocate(fa(ncomp))
   fa(1:ncomp) = nint( parame(1:ncomp,2) / dzp  ) !grid points per sigma [-]

   !FOR DISP
   ALLOCATE(fa_disp(ncomp),ab_disp(ncomp))
   psi_factor = 1.5
   fa_disp(1:ncomp) = nint( psi_factor * REAL(fa(1:ncomp)) )
   Do i=1,ncomp
     if( mod(fa_disp(i),2) /= 0 ) fa_disp(i) = fa_disp(i) + 1
   End Do


   ab_disp(1:ncomp) = psi_factor * dhs(1:ncomp)


   ngp = 2 * maxval(fa_disp(1:ncomp)) + 5   !number of ghost points (+5 kann eig weg)

   Allocate(zp(-ngp:ngrid+ngp))

   Do i=-ngp,ngrid+ngp
      zp(i) = REAL(i) * dzp          !z-distance from origin of grid points [A]
   End Do

   pbulk = ( p * 1.e-30 ) / ( kbol*t* rhob(1,0) ) * rhob(1,0) !p/kT


   !update z3t, the T-dependent quantity that relates eta and rho, as eta = z3t*rho
   CALL perturbation_parameter

     call solversetup




End Program dft