SolverSetup.F90

 
Subroutine solversetup

Use basic_variables, Only: ncomp

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

Implicit None

#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscvec.h90>
#include <finclude/petscdmda.h>
#include <finclude/petscdm.h>
#include <finclude/petscdmda.h90>
#include <finclude/petscis.h>
#include <finclude/petscmat.h>
#include <finclude/petscsnes.h>
#include <finclude/petscsnes.h90>


dm             ::  da
petscerrorcode ::  ierr
type(userctx) ::  user
mat            ::  j
petscreal      ::  erel
petscbool      ::  flg
INTEGER        ::  jac_id
petscint       ::  jaclocal

character(80) :: filename=''

!external subroutines associated with Solver
external forminitialguess
external formfunction
external jac_shell_ad
external jac_matrix_empty
external monitortimer


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


!  Create solver context
     call snescreate(petsc_comm_world,snes,ierr) 

!  Create distributed array (DMDA) to manage parallel grid and vectors
     call dmdacreate1d(petsc_comm_world,                                &  !MPI communicator
     &     dmda_boundary_ghosted,                                          &  !Boundary type at boundary of physical domain
     &     ngrid,                                                       &  !global dimension of array (if negative number, then value can be changed by user via command line!)
     &     ncomp,                                                           &  !number of degrees of freedom per grid point (number of unknowns at each grid point)
     &     ngp,                                                           &  !number of ghost points accessible for local vectors
     &     petsc_null_integer,                                          &  !could be an array to specify the number of grid points per processor
     &     da,                                                          &  !the resulting distributed array object 
     &     ierr)   

!  Extract global arrays from DMDA: x: unknowns, r: residual
      call dmcreateglobalvector(da,x,ierr)
      call vecduplicate(x,r,ierr)

!  Get local grid boundaries (for 1-dimensional DMDA)
 call dmdagetcorners(da,       & !the distributed array 
     &   user%xs,                   & !corner index in x direction
     &   petsc_null_integer,        & !corner index in y direction
     &   petsc_null_integer,        & !corner index in z direction
     &   user%xm,                   & !width of locally owned part in x direction
     &   petsc_null_integer,        & !width of locally owned part in y direction
     &   petsc_null_integer,        & !width of locally owned part in z direction 
     &   ierr)                        !error check 
       
      call dmdagetghostcorners(da,  & !the distributed array
     &   user%gxs,                  & !corner index in x direction (but now counting includes ghost points) 
     &   petsc_null_integer,        & !corner index in y direction (but now counting includes ghost points)
     &   petsc_null_integer,        & !corner index in z direction (but now counting includes ghost points)
     &   user%gxm,                  & !width of locally owned part in x direction (but now including ghost points)
     &   petsc_null_integer,        & !width of locally owned part in y direction (but now including ghost points)
     &   petsc_null_integer,        & !width of locally owned part in z direction (but now including ghost points)
     &   ierr)                        !error check

!  Here we shift the starting indices up by one so that we can easily
!  use the Fortran convention of 1-based indices (rather 0-based indices).
      user%xs  = user%xs+1
      user%gxs = user%gxs+1
      user%xe  = user%xs+user%xm-1
      user%gxe = user%gxs+user%gxm-1

      call snessetfunction(snes,r,formfunction,user,ierr)
      call snessetapplicationcontext(snes,user,ierr)
      call snessetdm(snes,da,ierr)

! !Set up matrix free jacobian
!       erel = 1e-08
!       call MatCreateSNESMF(snes,J,ierr)  !matrix free jacobi matrix
!       call MatMFFDSetFunctionError(J,erel,ierr)
!       call SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,PETSC_NULL_OBJECT,ierr)


! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Set up nonlinear solver depending on option set in makefile (with -jac)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      call petscoptionsgetint(petsc_null_character,'-jac',jac_id,flg,ierr)
      If(flg) Then
                
      Else         
             jac_id = 0
             write(*,*)'Option -jac not set in makefile, using matrix-free with numerical approximations.'
      End If
      
      Select Case(jac_id)
       
            Case(-1) !for options that dont need a jacobian: Anderson mixing, Picard

            Case (0) !matrix-free with numerical approximation of J(x)v        
               !default value of erel
               erel = 1e-08
               !check if value for erel is specified in makefile
               call petscoptionsgetreal(petsc_null_character,'-erel',erel,flg,ierr)      
               If(user%rank == 0) write(*,*)'Using matrix-free Jacobian with numerical approximations. erel:',erel
               call matcreatesnesmf(snes,j,ierr)  !matrix free jacobi matrix
               call matmffdsetfunctionerror(j,erel,ierr)
               call snessetjacobian(snes,j,j,matmffdcomputejacobian,petsc_null_object,ierr)
 
            Case (1) !matrix-free with AD-calculated J(x)v
               If(user%rank == 0) write(*,*)'Using matrix-free AD Jacobian.'
               !determine local size of shell matrix
               jaclocal = int(ngrid / user%num_procs) 
               If(mod(ngrid,user%num_procs) /= 0 ) Then
                   write(*,*)'Surface Tension Code: Dimensions and number of cores dont match! mod(ngrid,n_cores) must equal 0'
                   stop 5                        
               End if
               
               !Make sure that local parts of JacShell have the same size as the local DMDA-Arrays (local size JacShell != user%xm!!)
               If(jaclocal /= user%xm) Then
                   write(*,*)'Surface Tension Code: Shell-Jacobi-Matrix and DMDA have to be parallelized accordingly.'
!               write(*,*) 'Surface Tension Code: Shell-Jacobi-Matrix and DMDA have to be parallelized accordingly!!(JacLocal = user%xm)'
                   stop 6 
               End If
               
               call matcreateshell(petsc_comm_world,ncomp*jaclocal,ncomp*jaclocal,ncomp*ngrid,ncomp*ngrid,petsc_null_object,j,ierr)
               call matshellsetoperation( j, matop_mult, jac_shell_ad, ierr )
               call snessetjacobian( snes, j, j, jac_matrix_empty, petsc_null_object, ierr)
              
           Case (2) !build complete Jacobi matrix via finite-differences
              If(user%rank == 0) write(*,*)'Using finite-difference Jacobian.'
!Now in makefile: -snes_fd
!               call MatCreate(PETSC_COMM_WORLD,J,ierr)
!               call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,ngrid,ngrid,ierr)  
!               call MatSetFromOptions(J,ierr)                            
!               call MatSetUp(J,ierr)      !!sets up the internal matrix data structure 
!               call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,PETSC_NULL_OBJECT,ierr)

               
           Case (3) !build the complete Jacobi matrix via AD
               If(user%rank == 0) write(*,*)'Surface Tension Code: Using AD-generated Jacobian.'
               write(*,*)'Not working yet!' 
               stop 6              
               !create Matrix that has the appropriate sparsity pattern for the da
               !But that also means that when values are inserted in positions that
               !dont fit this structure, an error occurs!
               !i.e. when calculating the jacobi matrix it is not possible anymore
               !to ignore the sparsity pattern and just create a dense jacobian!!                 
               
               !call DMCreateMatrix(da,MATMPIAIJ,J,ierr)   !CAUTION: in newer PETSc Versions DMCreateMatrix does not contain the MatType argument anymore!
               !call SNESSetJacobian(snes,J,J,Jac_AD,PETSC_NULL_OBJECT,ierr)
            
          Case default
              write(*,*) 'No valid option set for -jac in makefile, using full Jacobi via finite-differences.'  
              call snessetjacobian(snes,j,j,snescomputejacobiandefault,petsc_null_object,ierr)
 
      End Select

!Set monitoring function that ouputs the surface tension at every iteration
      call snesmonitorset(snes,monitortimer,petsc_null_object,petsc_null_function,ierr)

!Enable modification from makefile
      call snessetfromoptions(snes,ierr)

!Evaluate Initial Guess
      call forminitialguess(snes,x,ierr)
      !write(filename,*)'0_initial_profile_global.xlo'
      !call PrintGlobalVec(x,filename)

      !start timer
      total_time = 0.0
      timer_old = mpi_wtime()

!Solve system
      call snessolve(snes,petsc_null_object,x,ierr)
      !write(filename,*)'1_final_profile_global.xlo'
      !call PrintGlobalVec(x,filename)
      !write(filename,'(a,I3.3,a)') '2_final_profile_local_proc_',user%rank,'.xlo' 
      !call PrintLocalVec(x,da,filename)  

      !plot residual vs time graph
      !call system('gnuplot gnuplot_script.srp')  
    

End Subroutine solversetup







subroutine monitortimer(snes,its,norm,dummy,ierr)

  Use global_x, Only: timer,timer_old,total_time,r
  Use mod_dft,  Only: free    
  Use parameters, Only: kbol,pi
  Use basic_variables, Only: t,parame, ncomp, surfactant, wif_surfactant
  Use vle_var, Only: rhob,tc 

! ! !PETSc modules
! !       use f90module   
! ! !DFT modules
! !       use DFT_MODULE, ONLY: free
! !       use VLE_VAR, ONLY: tc,rhob
! !       use BASIC_VARIABLES, ONLY:t,parame,ncomp
! !       use PARAMETERS, ONLY: PI,KBOL   

      implicit none

#include <finclude/petscsys.h>
#include <finclude/petscvec.h>

!  Input/output variables:
      snes           snes
      petscint   dummy
      Integer :: its
      REAL    :: norm
      petscerrorcode ierr
! local
      petscreal :: f2norm
      vec       :: current_solution
      DOUBLE PRECISION :: delta_time
      INTEGER   :: rank
      REAL      :: m_average,surftens,st_macro
      character(80) :: filename=''

      Real :: scale_factor !to scale surface tension in case surfactant is present


        !calculate interfacial tension 
          !sum up the variable free from all procs on proc 0
          call mpi_reduce(free,free,1,mpi_double_precision,mpi_sum,0,petsc_comm_world,ierr)  

          !Only proc 0 calculates surface tension
          call mpi_comm_rank(petsc_comm_world,rank,ierr)
          IF(rank == 0) THEN
    
             !calculate surface tension  
             surftens = kbol * t *1.e20*1000.0 *free 
             m_average = sum( rhob(1,1:ncomp)*parame(1:ncomp,1) ) / rhob(1,0)

             st_macro = surftens / ( 1.0 + 3.0/8.0/pi *t/tc  &
                                   * (1.0/2.55)**2  / (0.0674*m_average+0.0045) )

            write(*,*)'ST',st_macro
            
            !calculate scaling factor to implicitly account for surfactant
            If(surfactant) Then
               scale_factor = 0.8 - 2.4828*wif_surfactant  !adjusted to experimental results
            Else
               scale_factor = 1.
            End If
             
            If(surfactant) write(*,*)'ST including surfactant', st_macro*scale_factor

             
            !write result to outputfile
            filename='./out.txt'
            CALL file_open(filename,99)
            WRITE (99,*) st_macro*scale_factor
            close(99)
             
          End If
        

        !calculate elapsed time
        timer = mpi_wtime()  
        delta_time = timer - timer_old
        timer_old = timer
        total_time = total_time + delta_time
        
        !get norm of residual array
        call vecnorm(r,2,f2norm,ierr)

        IF(rank == 0) THEN
           !print results to file
           filename = 'ItsTimeNorm.dat'
           open(unit = 44, file = filename)
           write(44,*) its,total_time,f2norm
        End If 
           
end subroutine monitortimer