gravity_spectral.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: gravity_spectral.f90                                              #
!#                                                                           #
!# Copyright (C) 2011-2024                                                   #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!# Jannes Klee <jklee@astrophysik.uni-kiel.de>                               #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE gravity_spectral_mod
  USE gravity_base_mod
  USE fluxes_base_mod
  USE physics_base_mod
  USE mesh_base_mod
  USE logging_base_mod
  USE marray_base_mod
  USE marray_compound_mod
  USE functions
  USE common_dict
  USE fftw
#ifdef PARALLEL
  USE mpi
#endif
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  REAL, PARAMETER              :: sqrttwopi &
    = 2.50662827463100050241576528481104525300698674
 
  TYPE, EXTENDS(gravity_base) :: gravity_spectral
#if defined(HAVE_FFTW)
 
    TYPE(c_ptr)                      :: pfdensity, pfphi
    TYPE(c_ptr)                      :: plan_r2c
    TYPE(c_ptr)                      :: plan_c2r
    REAL(c_double), DIMENSION(:,:), POINTER &
                                     :: fdensity,fphi,block
    COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:,:), POINTER &
                                     :: cfdensity, cfphi, cblock
    REAL(c_double), DIMENSION(:,:,:), POINTER &
                                     :: fi
    COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(:,:,:), POINTER &
                                     :: cfi
    TYPE(c_ptr)                      :: p_fi
    REAL, DIMENSION(:,:,:), POINTER  :: phi
    REAL, DIMENSION(:,:), POINTER    :: tmp2d
    REAL, DIMENSION(:,:), POINTER    :: phi2d
    REAL, DIMENSION(:), POINTER      :: height1d
    INTEGER                          :: green
    REAL                             :: sigma, ecut
    INTEGER, DIMENSION(:), POINTER   :: sizes
    INTEGER, POINTER                 :: mcut
    INTEGER                          :: inum,knum
    INTEGER                          :: mnum
    INTEGER                          :: imax, kmax
#if defined(PARALLEL)
 
    INTEGER,DIMENSION(:), POINTER    :: displ, num
    INTEGER                          :: mpi_error
    REAL,DIMENSION(:,:), POINTER     :: sbuf1,sbuf2,& !< send buffers
                                        rbuf1,rbuf2
#endif
#endif
  CONTAINS
    PROCEDURE :: initgravity_spectral
    PROCEDURE :: setoutput
    PROCEDURE :: updategravity_single
    PROCEDURE :: calcdiskheight_single
    final :: finalize
#ifdef HAVE_FFTW
    PROCEDURE :: calcpotential
    PROCEDURE :: precomputei
    PROCEDURE :: calcmcut
#endif
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       gravity_spectral
  !--------------------------------------------------------------------------!
  CONTAINS
 
  SUBROUTINE initgravity_spectral(this,Mesh,Physics,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_spectral), INTENT(INOUT) :: this
    CLASS(mesh_base),            INTENT(IN) :: Mesh
    CLASS(physics_base),         INTENT(IN) :: Physics
    TYPE(dict_typ),              POINTER    :: config,IO
    !------------------------------------------------------------------------!
    INTEGER           :: err, valwrite
    CHARACTER(LEN=32) :: info_str
    !------------------------------------------------------------------------!
    CALL this%InitGravity(mesh,physics,"spectral",config,io)
 
#if !defined(HAVE_FFTW)
    CALL this%Error("InitGravity_spectral", &
         "Mandatory requirement fft has not been enabled. "//&
         "Please add --with-fftw=$FFTWDIR or similar to configure call.")
#endif
#if defined(HAVE_FFTW)
#ifdef PARALLEL
    ! Check domains a annular rings
    IF(.NOT.(mesh%dims(2).EQ.1)) &
      CALL this%Error("InitGravity_spectral", &
                 "Only domains shaped as annular rings are allowed (N x 1 decompositions).")
#endif
 
!    IF(.NOT.(GetType(Boundary(NORTH)).EQ.PERIODIC .AND. &
!             GetType(Boundary(SOUTH)).EQ.PERIODIC)) THEN
!      CALL this%Error(,"InitGravity_spectral", &
!                 "The boundary conditions in north and south direction have " // &
!                 "to be periodic! Don't forget: This kind of " // &
!                 "self-gravitation only works for polar-like coordinate " // &
!                 "systems.")
!    END IF
 
 
    IF(.NOT.(mod(mesh%JNUM,2)==0)) THEN
      CALL this%Error("InitGravity_spectral", &
                 "The spectral poisson solver needs an even number of cells " // &
                 "in the phi direction due to the discrete cosinus transform.")
    END IF
 
 
! If this is the innermost domain include 1 ghost cell,
! if this is the outermost domain include also 1 ghost cell.
#ifdef PARALLEL
    IF(mesh%IMAX.EQ.mesh%INUM) THEN
      this%IMAX = mesh%IMAX+1
    ELSE
      this%IMAX = mesh%IMAX
    END IF
#else
    !always innermost and outermost domain
    this%IMAX = mesh%IMAX+1
#endif
 
    this%INUM = this%IMAX - mesh%IMIN + 1
 
    ! number of fourier modes
    this%MNUM = mesh%JNUM/2+1
 
    ALLOCATE(this%phi2D(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX), &
             this%tmp2D(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX), &
             this%pot, &
             this%height1D(1:mesh%INUM), &
             this%sizes(1),&
             this%mcut, &
#ifdef PARALLEL
             this%sbuf1(mesh%GNUM,mesh%JNUM), &
             this%rbuf1(mesh%GNUM,mesh%JNUM), &
             this%sbuf2(mesh%GNUM,mesh%JNUM), &
             this%rbuf2(mesh%GNUM,mesh%JNUM), &
             this%displ(0:this%GetNumProcs()-1),&
             this%num(0:this%GetNumProcs()-1),&
#endif
             stat=err)
 
    IF (err.NE.0) &
      CALL this%Error("gravity_spectral::InitGravity_spectral","memory allocation failed")
 
    ! initialize potential
    this%pot = marray_base(4)
    this%pot%data1d(:) = 0.0
 
    CALL getattr(config, "green", this%green, 1)
    SELECT CASE(this%green)
    CASE(1,2,3)
      ! ok -> do nothing
    CASE DEFAULT
      CALL this%Error("gravity_spectral::InitGravity_spectral","type of Green's function should be one of 1,2,3")
    END SELECT
 
    CALL getattr(config, "sigma", this%sigma, 0.05)
    this%height1D(:) = this%sigma
 
    ! mode cut off.
    !  0: disabled
    ! >0: constant cutoff number
    !     (can be overwritten by ecut>0.)
    CALL getattr(config, "mcut", this%mcut, this%MNUM)
    this%mcut = min(this%mcut,this%MNUM)
 
    ! mode cutoff relative error boundary
    ! <=0: no automatic cutoff calculation
    ! >0 : automatic setting of mcut
    CALL getattr(config, "ecut", this%ecut, 0.)
 
    CALL this%Info("            Initializing")
 
    WRITE (info_str, '(I8)') this%green
    CALL this%Info("            green-fn type:     " // trim(info_str))
    WRITE (info_str, '(ES8.2)') this%sigma
    CALL this%Info("            sigma:             " // trim(info_str))
 
 
    !\todo only 2D variables
    this%p_FI = fftw_alloc_real(int(2*this%MNUM * (this%INUM)*(mesh%INUM), c_size_t))
    CALL c_f_pointer(this%p_FI, this%FI, &
                     [2*this%MNUM,this%INUM,mesh%INUM])
    CALL c_f_pointer(this%p_FI, this%cFI, &
                     [this%MNUM,this%INUM,mesh%INUM])
 
 
    this%pFdensity = fftw_alloc_real(int(2*this%MNUM*mesh%INUM, c_size_t))
    CALL c_f_pointer(this%pFdensity, this%Fdensity, &
                     [2*this%MNUM,mesh%INUM])
    CALL c_f_pointer(this%pFdensity, this%cFdensity, &
                     [mesh%JNUM/2+1,mesh%INUM])
 
    this%pFphi = fftw_alloc_real(int(2*this%MNUM*this%INUM,c_size_t))
    CALL c_f_pointer(this%pFphi, this%Fphi, &
                     [2*this%MNUM, this%INUM])
    CALL c_f_pointer(this%pFphi, this%cFphi, &
                     [mesh%JNUM/2+1, this%INUM])
 
    IF (err.NE.0) &
        CALL this%Error("InitGravity_spectral","Memory allocation failed.")
 
    this%block => this%Fdensity(:,mesh%IMIN:mesh%IMAX)
    this%cblock => this%cFdensity(:,mesh%IMIN:mesh%IMAX)
 
    this%phi2D(:,:) = 0.
    this%tmp2D(:,:) = 0.
 
    ! Create plans for fftw
 
    ! Use FFTW_MEASURE for calculating the fastest plan, but this
    ! costs some extra seconds
    this%sizes(1) = mesh%JNUM
    this%plan_r2c = fftw_plan_many_dft_r2c(1, this%sizes, (mesh%IMAX-mesh%IMIN+1),&
                                           this%block, this%sizes, &
                                           1, 2*this%MNUM, &
                                           this%cblock, this%sizes, &
                                           1, this%MNUM, &
                                           fftw_measure)
 
    this%sizes(1) = mesh%JNUM
    this%plan_c2r = fftw_plan_many_dft_c2r(1, this%sizes, this%INUM, &
                                           this%cFphi, this%sizes, &
                                           1, this%MNUM, &
                                           this%Fphi, this%sizes, &
                                           1, 2*this%MNUM, &
                                           fftw_measure)
 
    IF(this%ecut.GT.0.) &
      CALL setattr(io, "mcut",ref(this%mcut))
 
#ifdef PARALLEL
    this%displ((this%GetRank())) = (mesh%IMIN-1)*2*this%mcut
    this%num((this%GetRank())) = (mesh%IMAX-mesh%IMIN+1)*2*this%mcut
    CALL mpi_allgather(mpi_in_place, 0, mpi_datatype_null, &
                       this%displ, 1, mpi_integer, mpi_comm_world, this%mpi_error)
    CALL mpi_allgather(mpi_in_place, 0, mpi_datatype_null, &
                       this%num, 1, mpi_integer, mpi_comm_world, this%mpi_error)
 
#endif
 
    CALL this%PrecomputeI(mesh, physics)
 
    CALL this%Info("            .. done initializing")
#endif
 
  END SUBROUTINE initgravity_spectral
 
  SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_spectral), INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN)    :: Mesh
    CLASS(physics_base), INTENT(IN)    :: Physics
    TYPE(dict_typ),      POINTER       :: config,IO
    !------------------------------------------------------------------------!
    INTEGER          :: valwrite
    !------------------------------------------------------------------------!
    valwrite = 0
    CALL getattr(config, "output/potential", valwrite, 0)
    IF (valwrite .EQ. 1.AND.ALLOCATED(this%pot)) THEN
      IF (ASSOCIATED(this%pot%data4d)) &
        CALL setattr(io, "potential", &
              this%pot%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
    END IF
  END SUBROUTINE setoutput
 
 
#ifdef HAVE_FFTW
  FUNCTION calcmcut(this,Mesh) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_spectral)      :: this
    CLASS(mesh_base), INTENT(IN) :: mesh
    INTEGER                      :: res
    !------------------------------------------------------------------------!
    INTEGER           :: m,i
#ifdef PARALLEL
    INTEGER           :: ier
#endif
    REAL              :: cumsum(this%mnum)
    INTEGER           :: mcut(mesh%imin:mesh%imax)
    !------------------------------------------------------------------------!
    DO i=mesh%IMIN,mesh%IMAX
      cumsum(1) = this%Fdensity(2*1-1,i)**2 + this%Fdensity(2*1,i)**2
      DO m=2,this%MNUM
        cumsum(m) = cumsum(m-1) + this%Fdensity(2*m-1,i)**2 + this%Fdensity(2*m,i)**2
      END DO
      DO m=1,this%MNUM
        IF(abs(cumsum(m)/cumsum(this%MNUM)-1.).LT.this%ecut) THEN
          mcut(i) = m
          EXIT
        END IF
      END DO
    END DO
    res = maxval(mcut)
 
#ifdef PARALLEL
    CALL mpi_allreduce(mpi_in_place,res,1,mpi_integer,mpi_max,mpi_comm_world,ier)
#endif
  END FUNCTION calcmcut
#endif
 
#if defined(HAVE_FFTW)
  SUBROUTINE calcpotential(this,Mesh,Physics,time,pvar)
    USE physics_eulerisotherm_mod, ONLY : statevector_eulerisotherm
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_spectral), INTENT(INOUT):: this
    CLASS(mesh_base),    INTENT(IN) :: Mesh
    CLASS(physics_base), INTENT(IN) :: Physics
    CLASS(marray_compound),   INTENT(INOUT) :: pvar
    REAL, INTENT(IN)                :: time
    !------------------------------------------------------------------------!
    INTEGER           :: i,k,m,i0,oldmcut
#ifdef PARALLEL
    INTEGER           :: status(MPI_STATUS_SIZE)
#endif
    !------------------------------------------------------------------------!
    ! Fourier transform the density with respect to Phi
    SELECT TYPE(p => pvar)
    CLASS IS(statevector_eulerisotherm)
      DO k=mesh%KMIN, mesh%KMAX
        DO i=mesh%IMIN, mesh%IMAX
          this%Fdensity(1:mesh%JNUM,i) = p%density%data3d(i,mesh%JMIN:mesh%JMAX,k)
        END DO
      END DO
    CLASS DEFAULT
      CALL this%Error("gravity_spectral::CalcPotential","unsupported state vector")
    END SELECT
 
    CALL fftw_execute_dft_r2c(this%plan_r2c, this%block, this%cblock)
 
    IF(this%ecut.GT.0.) THEN
      oldmcut = this%mcut
 
      this%mcut = this%CalcMcut(mesh)
 
#ifdef PARALLEL
      this%displ = this%displ / oldmcut * this%mcut
      this%num = this%num / oldmcut * this%mcut
#endif
    END IF
 
#ifdef PARALLEL
    ! Distribute Fdensity to all processes
      CALL mpi_allgatherv(mpi_in_place, 0, mpi_datatype_null,&
                          this%Fdensity(1:2*this%mcut,1:mesh%INUM), &
                          this%num, this%displ, default_mpi_real, &
                          mpi_comm_world, this%mpi_error)
#endif
 
    ! Integrate in radial direction by numerical quadrature
    ! Fphi_m(t,r) = int_rmin^rmax Fdensity_m(t,r') * I_m(r,r') dr'
    ! dr' is already multiplicated into FI.
 
    ! Direct summation and multiplication with complex numbers:
    !DO m=1, this%MNUM
    !  ! integrate
    !  DO i=1, this%INUM
    !    this%cFPhi(m,i) = SUM(this%cFdensity(m,:) * this%cFI(m,i,:))
    !  END DO
    !END DO
 
    ! or with real numbers, real and imaginary part seperate:
    this%FPhi(:,:) = 0.
!NEC$ IVDEP
    DO i0=1,mesh%INUM
!NEC$ IVDEP
      DO i=1,this%INUM
!NEC$ IVDEP
        DO m=1,this%mcut
           ! real part
           this%FPhi(2*m-1,i) = this%FPhi(2*m-1,i) &
             + this%Fdensity(2*m-1,i0) * this%FI(2*m-1,i,i0)
           ! imag part
           this%FPhi(2*m,i) = this%FPhi(2*m,i) &
             + this%Fdensity(2*m,i0) * this%FI(2*m-1,i,i0)
        END DO
      END DO
    END DO
 
    ! Inverse fourier transform Fphi with respect to Phi
    CALL fftw_execute_dft_c2r(this%plan_c2r, this%cFphi, this%Fphi)
 
    DO i=1,this%INUM
      this%phi2D(i+mesh%IMIN-1,mesh%JMIN:mesh%JMAX) = this%FPhi(1:mesh%JNUM,i)
    END DO
 
    ! calculate boundary data (only one ghost cell is enough)
#ifdef PARALLEL
    ! send boundary data to western and receive from eastern neighbor
    IF (mesh%neighbor(west).NE.mpi_proc_null) &
         this%sbuf1 = this%phi2D(mesh%IMIN:mesh%IMIN+mesh%GNUM-1,mesh%JMIN:mesh%JMAX)
    CALL mpi_sendrecv(this%sbuf1,mesh%GNUM*mesh%JNUM, &
         default_mpi_real,mesh%neighbor(west),10+west,this%rbuf1, &
         mesh%GNUM*mesh%JNUM,default_mpi_real,mesh%neighbor(east), &
         mpi_any_tag,mesh%comm_cart,status,this%mpi_error)
    IF (mesh%neighbor(east).NE.mpi_proc_null) &
         this%phi2D(mesh%IMAX+1:mesh%IGMAX,mesh%JMIN:mesh%JMAX) = this%rbuf1
    ! send boundary data to western and receive from eastern neighbor
    IF (mesh%neighbor(east).NE.mpi_proc_null) &
         this%sbuf2 = this%phi2D(mesh%IMAX-mesh%GNUM+1:mesh%IMAX,mesh%JMIN:mesh%JMAX)
    CALL mpi_sendrecv(this%sbuf2,mesh%GNUM*mesh%JNUM, &
         default_mpi_real,mesh%neighbor(east),10+east,this%rbuf2, &
         mesh%GNUM*mesh%JNUM,default_mpi_real,mesh%neighbor(west), &
         mpi_any_tag,mesh%comm_cart,status,this%mpi_error)
    IF (mesh%neighbor(west).NE.mpi_proc_null) &
         this%phi2D(mesh%IGMIN:mesh%IMIN-1,mesh%JMIN:mesh%JMAX) = this%rbuf2
#endif
    this%phi2D(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JMIN-1) &
      = this%phi2D(mesh%IGMIN:mesh%IGMAX,mesh%JMAX-mesh%GNUM+1:mesh%JMAX)
    this%phi2D(mesh%IGMIN:mesh%IGMAX,mesh%JMAX+1:mesh%JGMAX) &
      = this%phi2D(mesh%IGMIN:mesh%IGMAX,mesh%JMIN:mesh%JMIN+mesh%GNUM-1)
 
  END SUBROUTINE calcpotential
#endif
 
#ifdef HAVE_FFTW
  ELEMENTAL FUNCTION greenfunction(dr2, green, sigma) RESULT(G)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL,    INTENT(IN) :: dr2
    INTEGER, INTENT(IN) :: green
    REAL,    INTENT(IN) :: sigma
    REAL                :: g
    !------------------------------------------------------------------------!
    SELECT CASE(green)
      CASE(1)  ! Razor sharp disc
        g = -1.0/sqrt(dr2)
      CASE(2)  ! gaussian spheres
        g = -1.0 *  bessel_k0e(dr2/(4.0*sigma*sigma)) &
            / (sqrttwopi*sigma)
      CASE(3)  ! Used for the orbiting cylinder example
        g = log(dr2)
      CASE DEFAULT
        ! should never happen
        g = 0.0
    END SELECT
  END FUNCTION greenfunction
#endif
 
#ifdef HAVE_FFTW
  SUBROUTINE precomputei(this, Mesh, Physics)
   IMPLICIT NONE
   !------------------------------------------------------------------------!
   CLASS(gravity_spectral)             :: this
   CLASS(mesh_base),        INTENT(IN) :: Mesh
   CLASS(physics_base)                 :: Physics
   !------------------------------------------------------------------------!
   TYPE(c_ptr)        :: plan_r2c
   INTEGER            :: i1, i, j, k
   REAL               :: r0, notzero
   CHARACTER(LEN=128) :: str
   REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX) &
                      :: phi
   REAL, DIMENSION(1:Mesh%INUM+1) :: r,hx,hy
   REAL, DIMENSION(Mesh%JMIN:Mesh%JMAX,Mesh%IMIN:this%IMAX) :: dr2
   INTEGER            :: rank, howmany, istride, idist, ostride, odist
   INTEGER, DIMENSION(1) &
                      :: n
#ifdef PARALLEL
   INTEGER, DIMENSION(0:(this%GetNumProcs())-1) :: num,displ
#endif
   !------------------------------------------------------------------------!
    n(1)    = mesh%JNUM
    rank    = 1
    howmany = this%INUM * mesh%INUM
    istride = 1
    idist = 2*this%MNUM
    ostride = istride
    odist   = this%MNUM
    plan_r2c = fftw_plan_many_dft_r2c(rank, n, howmany, &
                                  this%FI, n, &
                                  istride, idist, &
                                  this%cFI, n, &
                                  ostride, odist, &
                                  fftw_measure &
                                  )
 
    phi = mesh%curv%center(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN,2)&
        - mesh%curv%center(mesh%IMIN,mesh%JMIN,mesh%KMIN,2)
 
    r(mesh%IMIN:this%IMAX) = mesh%radius%bcenter(mesh%IMIN:this%IMAX,mesh%JMIN,mesh%KMIN)
    hx(mesh%IMIN:this%IMAX) = mesh%hx%bcenter(mesh%IMIN:this%IMAX,mesh%JMIN,mesh%KMIN)
    hy(mesh%IMIN:this%IMAX) = mesh%hy%bcenter(mesh%IMIN:this%IMAX,mesh%JMIN,mesh%KMIN)
#ifdef PARALLEL
    displ(this%GetRank()) = (mesh%IMIN-1)
    num(this%GetRank()) = this%INUM
    CALL mpi_allgather(mpi_in_place, 0, mpi_datatype_null, &
                       displ, 1, mpi_integer, mpi_comm_world, this%mpi_error)
    CALL mpi_allgather(mpi_in_place, 0, mpi_datatype_null, &
                       num, 1, mpi_integer, mpi_comm_world, this%mpi_error)
    CALL mpi_allgatherv(mpi_in_place, 0, mpi_datatype_null, &
                        r, num, displ, &
                        default_mpi_real, mpi_comm_world, this%mpi_error)
    CALL mpi_allgatherv(mpi_in_place, 0, mpi_datatype_null, &
                        hx, num, displ, &
                        default_mpi_real, mpi_comm_world, this%mpi_error)
    CALL mpi_allgatherv(mpi_in_place, 0, mpi_datatype_null, &
                        hy, num, displ, &
                        default_mpi_real, mpi_comm_world, this%mpi_error)
#endif
 
    DO k=mesh%KMIN,mesh%KMAX
      DO i1=1,mesh%INUM
!NEC$ IVDEP
        DO i=mesh%IMIN,this%IMAX
          DO j=mesh%JMIN,mesh%JMAX
            r0 = mesh%radius%faces(i,j,k,1)
            dr2(j,i) = r(i1)**2 + r0**2 - 2.*r(i1)*r0*cos(phi(j))
          END DO
        END DO
 
        this%FI(1:mesh%JNUM,1:this%INUM,i1) &
          = 2.0 * pi * hx(i1) * hy(i1) * mesh%dx &
          * physics%Constants%GN &
          * greenfunction(dr2(:,:),this%green,this%height1D(i1)) &
          / (mesh%JNUM)**2
      END DO
    END DO
 
    CALL fftw_execute_dft_r2c(plan_r2c, this%FI, this%cFI)
 
    ! Free plan_r2c
    CALL fftw_destroy_plan(plan_r2c)
 
    ! I has an even symmetry around 0. Normally a cosinus tranform (r2r)
    ! would be sufficent. But this is not available on all platforms. There-
    ! fore do a standard r2c dft and ditch the imaginary part of the result,
    ! which is approx zero numerically (and exactly zero analytically).
    !this%cFI = REAL(this%cFI)
    notzero = sum(abs(this%FI(2:2*this%MNUM:2,:,:)))/(this%MNUM*mesh%JNUM*this%INUM)
#ifdef PARALLEL
    CALL mpi_allreduce(mpi_in_place, notzero, 1, default_mpi_real, &
      mpi_max, &
      mpi_comm_world, this%mpi_error)
#endif
    IF(notzero.GT.1.) THEN
      WRITE(str,'(A,ES12.4,A)')"Imag(FI) should be zero, but ",notzero,&
        " > 1."
      CALL this%Warning("ProcomputeI_spectral",trim(str))
    END IF
 
    this%FI(2:2*this%MNUM:2,:,:) = 0.
  END SUBROUTINE precomputei
#endif
 
  SUBROUTINE updategravity_single(this,Mesh,Physics,Fluxes,pvar,time,dt)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_spectral),  INTENT(INOUT) :: this
    CLASS(mesh_base),            INTENT(IN) :: Mesh
    CLASS(physics_base),         INTENT(IN) :: Physics
    CLASS(fluxes_base),          INTENT(IN) :: Fluxes
    CLASS(marray_compound),   INTENT(INOUT) :: pvar
    REAL,                        INTENT(IN) :: time,dt
    !------------------------------------------------------------------------!
    INTEGER              :: i, j, k
    !------------------------------------------------------------------------!
#ifdef HAVE_FFTW
    ! calc potential first
    CALL this%CalcPotential(mesh,physics,time,pvar)
 
    ! compute gravitational acceleration using the gradient of the potential
    ! ATTENTION: ghost cell values and any other component of the gravitational
    !            acceleration are set to zero in InitGravity_spectral
    ! \todo This difference quotient is still 2D and only the third component is
    !       added. It will not work in 3D.
!NEC$ IVDEP
    DO k = mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
      DO j = mesh%JMIN-mesh%JP1,mesh%JMAX+mesh%JP1
!NEC$ IVDEP
        DO i = mesh%IMIN-mesh%IP1,mesh%IMAX
          this%accel%data4d(i,j,k,1) = -1.0*(this%phi2D(i+1,j)-this%phi2D(i,j))/mesh%dlx%data3d(i,j,k)
          this%accel%data4d(i,j,k,2) = -1.0*(this%phi2D(i+1,j+1)+this%phi2D(i,j+1) &
                           -this%phi2D(i+1,j-1)-this%phi2D(i,j-1))&
                            /(4.0*mesh%dly%data3d(i,j,k))
          this%pot%data4d(i,j,k,1) = 0.5*(this%phi2D(i,j)+this%phi2D(i+1,j))
          this%pot%data4d(i,j,k,2) = this%phi2D(i+1,j)
          this%pot%data4d(i,j,k,3) = 0.25*(this%phi2D(i,j)+this%phi2D(i+1,j)+this%phi2D(i,j+1)+this%phi2D(i+1,j+1))
        END DO
      END DO
    END DO
#endif
  END SUBROUTINE updategravity_single
 
 
  SUBROUTINE calcdiskheight_single(this,Mesh,Physics,pvar,bccsound,h_ext,height)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_spectral), INTENT(INOUT) :: this
    CLASS(mesh_base),        INTENT(IN)    :: Mesh
    CLASS(physics_base),     INTENT(IN)    :: Physics
    CLASS(marray_compound),  INTENT(INOUT) :: pvar
    TYPE(marray_base),       INTENT(INOUT) :: bccsound,h_ext,height
    !------------------------------------------------------------------------!
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,2) &
                      :: phi
    INTEGER           :: i,j,k
    REAL              :: cs2,p,q
    !------------------------------------------------------------------------!
#ifdef HAVE_FFTW
    ! compute the laplacian of the gravitational potential in the disks
    ! equatorial plane; gravity_generic updates the gravitational acceleration
    DO k=mesh%KMIN-mesh%KP1,mesh%KMAX
      DO j=mesh%JMIN-mesh%JP1,mesh%JMAX
        DO i=mesh%IMIN-mesh%IP1,mesh%IMAX
          phi(i,j,1) = mesh%hy%faces(i,j,k,2)*mesh%hz%faces(i,j,k,2)/mesh%hx%faces(i,j,k,2)&
           * this%pot%data4d(i,j,k,2)
          phi(i,j,2) = mesh%hx%faces(i,j,k,4)*mesh%hz%faces(i,j,k,4)/mesh%hy%faces(i,j,k,4)&
           * this%pot%data4d(i,j,k,3)
        END DO
      END DO
    END DO
    ! div(-grad(phi)))
    DO k=mesh%KMIN,mesh%KMAX
      DO j=mesh%JMIN,mesh%JMAX
        DO i=mesh%IMIN,mesh%IMAX
          this%tmp2D(i,j) = - 1./mesh%sqrtg%bcenter(i,j,k) &
            * (  (phi(i,j,1) - phi(i-1,j,1)) / mesh%dlx%data3d(i,j,k) &
               + (phi(i,j,2) - phi(i,j-1,2)) / mesh%dly%data3d(i,j,k))
        END DO
      END DO
    END DO
 
    IF (ASSOCIATED(h_ext%data1d).AND.(h_ext.match.height).AND. &
       (minval(h_ext%data1d,mask=mesh%without_ghost_zones%mask1d(:)).GT.0.0)) THEN
      ! disk height with external potential
      DO k=mesh%KGMIN,mesh%KGMAX
        DO j=mesh%JGMIN,mesh%JGMAX
          DO i=mesh%IGMIN,mesh%IGMAX
            cs2 = bccsound%data3d(i,j,k)**2
            p = sqrttwopi*physics%Constants%GN*pvar%data4d(i,j,k,physics%DENSITY)*h_ext%data3d(i,j,k)/cs2
            q = 1. + this%tmp2D(i,j)*h_ext%data3d(i,j,k)**2/cs2
            ! return the new disk height
            height%data3d(i,j,k) = h_ext%data3d(i,j,k) / (p+sqrt(q+p*p))
          END DO
        END DO
      END DO
    ELSE
      ! pure self-gravitating disk without external potential
      DO k=mesh%KGMIN,mesh%KGMAX
        DO j=mesh%JGMIN,mesh%JGMAX
          DO i=mesh%IGMIN,mesh%IGMAX
            cs2 = bccsound%data3d(i,j,k)**2
            p = sqrttwopi*physics%Constants%GN*pvar%data4d(i,j,k,physics%DENSITY)*mesh%radius%bcenter(i,j,k)/cs2
            q = this%tmp2D(i,j)*h_ext%data3d(i,j,k)**2/cs2
            ! return the new disk height
            height%data3d(i,j,k) = mesh%radius%bcenter(i,j,k) / (p+sqrt(q+p*p))
          END DO
        END DO
      END DO
    END IF
#endif
  END SUBROUTINE calcdiskheight_single
 
 
  SUBROUTINE finalize(this)
   IMPLICIT NONE
   !------------------------------------------------------------------------!
   TYPE(gravity_spectral), INTENT(INOUT) :: this
   !------------------------------------------------------------------------!
#if defined(HAVE_FFTW)
    ! Destroy plans
    CALL fftw_destroy_plan(this%plan_r2c)
    CALL fftw_destroy_plan(this%plan_c2r)
 
    ! Free memomry
    DEALLOCATE(&
#ifdef PARALLEL
               this%sbuf1, &
               this%rbuf1, &
               this%sbuf2, &
               this%rbuf2, &
               this%displ,this%num,&
#endif
               this%phi2D, &
               this%pot, &
               this%tmp2D, &
               this%sizes, &
               this%mcut &
               )
 
    CALL fftw_free(this%p_FI)
    CALL fftw_free(this%pFdensity)
    CALL fftw_free(this%pFphi)
 
    ! Free all temporary memory of FFTW
    CALL fftw_cleanup()
#endif
    IF (ALLOCATED(this%pot)) DEALLOCATE(this%pot)
    CALL this%Finalize_base()
   END SUBROUTINE finalize
 
END MODULE gravity_spectral_mod