AD_Routines.F90

 
 
 
 
 
 
 
!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
!
!  Differentiation of formfunctionlocal in forward (tangent) mode:
!   variations   of useful results: f rhop
!   with respect to varying inputs: rhop
!   RW status of diff variables: f:out rhop:in-out
SUBROUTINE formfunctionlocal_d(rhop, rhopd_nogp, fd, user, ierr)


!PETSc modules
 Use petscmanagement


!DFT modules
  USE mod_dft_fmt
  USE mod_dft_chain
  USE mod_dft_fmt_d
  USE mod_dft_chain_d
  USE mod_dft_disp_wda_d


  USE basic_variables, ONLY : ncomp 
  USE eos_variables, Only:dhs, rho
  USE vle_var, Only: rhob
  USE mod_dft, ONLY : fa, zp
  USE dft_fcn_module, ONLY : chempot_res
  USE global_x, ONLY : ngrid, ngp
  IMPLICIT NONE

#include <finclude/petscsys.h>

!  Input/output variables:
      type(userctx) user
      petscscalar  rhop(ncomp,user%gxs:user%gxe)
      petscscalar  rhopd_nogp(ncomp,user%xs:user%xe)
      !PetscScalar  f(ncomp,user%xs:user%xe)
      petscscalar  fd(ncomp,user%xs:user%xe)
      petscerrorcode ierr

!  Local variables:
      petscint  i
      petscint  k
      REAL,dimension(user%gxs:user%gxe)          :: n0,n1,n2,n3,nv1,nv2  !ngp muss groesser als fa+fa/2 sein!!
      REAL,dimension(user%gxs:user%gxe)          :: phi_dn0,phi_dn1,phi_dn2,phi_dn3,phi_dnv1,phi_dnv2
      REAL,dimension(user%gxs:user%gxe)          :: phi_dn0d,phi_dn1d,phi_dn2d,phi_dn3d,phi_dnv1d,phi_dnv2d
      REAL,dimension(user%xs:user%xe,ncomp)      :: df_drho_fmt,df_drho_fmtd
      REAL,dimension(user%gxs:user%gxe,ncomp)    :: rhobar,lambda,rhobard,lambdad  
      REAL,dimension(user%xs:user%xe,ncomp)      :: df_drho_chain,df_drho_chaind    
      REAL :: vext(ncomp)

      !DISP VAR
      REAL,dimension(user%gxs:user%gxe,ncomp)    :: rhop_wd,my_disp,df_disp_drk
      REAL,dimension(user%gxs:user%gxe,ncomp)    :: rhop_wdd,my_dispd
      REAL,dimension(user%gxs:user%gxe)          :: f_disp
      REAL, dimension(ncomp)                     :: df_drho_disp
      REAL, dimension(ncomp)                     :: df_drho_dispd



      INTRINSIC exp
      REAL :: arg1
      REAL :: arg1d
      petscscalar :: rhopd(ncomp,user%gxs:user%gxe)

      rhopd = 0.0
      rhopd(1:ncomp,user%xs:user%xe) = rhopd_nogp(1:ncomp,user%xs:user%xe)



!calculate weighted densities
  CALL fmt_weighted_densities_d(rhop, rhopd, n0, n1, n2, n3, nv1, nv2, &
&                         phi_dn0, phi_dn0d, phi_dn1, phi_dn1d, phi_dn2&
&                         , phi_dn2d, phi_dn3, phi_dn3d, phi_dnv1, &
&                         phi_dnv1d, phi_dnv2, phi_dnv2d, user)
!calculate averaged density rhobar and lambda (both needed for chain term)
  CALL chain_aux_d(rhop, rhopd, rhobar, rhobard, lambda, lambdad, user)
!HIER DISP!!!
  CALL disp_weighted_densities_d(rhop, rhopd, rhop_wd, rhop_wdd, user)
  CALL disp_mu_d(rhop_wd, rhop_wdd, f_disp, my_disp, my_dispd, &
&          df_disp_drk, user)

  fd = 0.0
  df_drho_chaind = 0.0
  df_drho_fmtd = 0.0
  df_drho_dispd = 0.0

  DO i = user%xs,user%xe
    CALL fmt_dfdrho_d(i, fa, user, phi_dn0, phi_dn0d, phi_dn1, phi_dn1d&
&               , phi_dn2, phi_dn2d, phi_dn3, phi_dn3d, phi_dnv1, &
&               phi_dnv1d, phi_dnv2, phi_dnv2d, df_drho_fmt, &
&               df_drho_fmtd)
    CALL chain_dfdrho_d(i, rhop, rhopd, lambda, lambdad, rhobar, rhobard&
&                 , df_drho_chain, df_drho_chaind, user)

!HIER DISP!!!
    CALL disp_dfdrho_wda_d(i, rhop, rhop_wd, my_disp, my_dispd, f_disp, &
&                    df_disp_drk, df_drho_disp, df_drho_dispd, user)


    vext(1:ncomp) = 0.0
    DO k=1,ncomp
      IF (zp(i) .LT. dhs(k)/2.0) vext(k) = 100000.0
!  If( zp(ngrid) - zp(i) < dhs(k)/2.0  )   Vext(k) = 100000.0
       arg1d = -df_drho_fmtd(i,k)-df_drho_chaind(i,k)-df_drho_dispd(k)
       arg1 = chempot_res(k) - vext(k) - df_drho_fmt(i, k) - &
&       df_drho_chain(i, k) - df_drho_disp(k)

      fd(k, i) = rhob(1,k)*arg1d*exp(arg1) - rhopd(k, i)
      !f(k, i) = xx(k)*rho*EXP(arg1) - rhop(k, i)
    END DO
  END DO
END SUBROUTINE formfunctionlocal_d









!------------------------------------------------------------------------------------------------
!This Subroutine calculates the Jacobi-Vector product using derivatives obtained via AD
!------------------------------------------------------------------------------------------------
Subroutine jac_shell_ad(Jshell,v_in,v_out)

use global_x, only: x,snes
Use petscmanagement

#include "finclude/petsc.h90"  


!passed
        mat                  :: jshell
        vec                  :: v_in        !has global size discret (NOT discret +- ghost points!!)
        vec                  :: v_out
!local
        petscscalar, pointer :: xd(:), rhop_loc(:,:),fd(:) !for dof2, xd and fd are twice the size as for dof1
        petscerrorcode       :: ierr
        type(userctx)       :: user
        dm                   :: da
        vec                  :: rhop_local

       

!get the user context and DM which are associated with nonlinear solver 
        call snesgetapplicationcontext(snes,user,ierr)
        call snesgetdm(snes,da,ierr)

!get local vector for x (needed because we need the ghost point values of x)
        call dmgetlocalvector(da,rhop_local,ierr)

!copy global to local for x (then x_local also contains ghost points)
        call dmglobaltolocalbegin(da,x,insert_values,rhop_local,ierr)
        call dmglobaltolocalend(da,x,insert_values,rhop_local,ierr)

!get pointers to the vectors x_local,v_in and v_out
        call dmdavecgetarrayf90(da,rhop_local,rhop_loc,ierr)
        call vecgetarrayf90(v_in, xd, ierr )
        call vecgetarrayf90(v_out, fd, ierr )

!  Get local grid boundaries (dont know why this is neccessary again!)
 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 AD generated Routine 
        !x_loc: x
        !xd   : direction which the derivative is calculated for
        !fd   : the directional derivative in direction xd
       
         call formfunctionlocal_d(rhop_loc,xd,fd,user,ierr)
!---------------------------------------------------------        




!restore arrays
        call dmdavecrestorearrayf90(da,rhop_local,rhop_loc,ierr )
        call vecrestorearrayf90( v_in, xd, ierr )
        call vecrestorearrayf90( v_out, fd, ierr )   
        call dmrestorelocalvector(da,rhop_local,ierr)


End Subroutine jac_shell_ad




! empty subroutine for shell jacobian; probably should copy v_x to x, as they might not be same
Subroutine jac_matrix_empty(snes,v_x,jac,B,flag,dummy,ierr)
        implicit none
#include "finclude/petsc.h90"
        snes            :: snes
        vec             :: v_x
        mat             :: jac,b
        matstructure    :: flag
        petscerrorcode  :: ierr
        integer dummy(*)
End Subroutine jac_matrix_empty