MetricOps.f90

Go to the documentation of this file.

MODULE metricops

  USE operators
  USE grid
  
  IMPLICIT NONE
  
CONTAINS
  
  
  SUBROUTINE alphaweight(                          &
    numDim,numPointsBuffer,bufferSizes,opInterval, &
         gridType,gridMetric,alphaDir,alphaW)
    
    IMPLICIT NONE 

    INTEGER(KIND=4), INTENT(IN)         :: numDim, gridType, alphaDir
    INTEGER(KIND=8), INTENT(IN)         :: bufferSizes(numdim),numPointsBuffer
    INTEGER(KIND=8), INTENT(IN)         :: opInterval(2*numdim)
    REAL(KIND=8),    INTENT(IN), TARGET :: gridMetric(numdim*numdim*numpointsbuffer)
    REAL(KIND=8),    INTENT(OUT)        :: alphaW(numpointsbuffer)

    REAL(KIND=8)    :: gridScale
    INTEGER(KIND=8) :: metricOffset

    REAL(KIND=8),    DIMENSION(:), POINTER :: metricPtr

    IF(gridtype < rectilinear) THEN

       gridscale = abs(gridmetric(alphadir))

       CALL assignmentxa(numdim,numpointsbuffer,buffersizes,opinterval, &
            gridscale,alphaw)

    ELSE IF(gridtype == rectilinear) THEN


       metricoffset =  (alphadir-1)*numpointsbuffer
       metricptr    => gridmetric(metricoffset+1:metricoffset+numpointsbuffer)

       CALL assignmentyabsx(numdim,numpointsbuffer,buffersizes,opinterval, &
            metricptr,alphaw);
       
    ELSE ! Curvilinear

       metricoffset =  (alphadir-1)*numdim*numpointsbuffer
       metricptr    => gridmetric(metricoffset+1:metricoffset+(numdim*numpointsbuffer))

       CALL veclen(numdim,numpointsbuffer,buffersizes,opinterval, &
            numdim,metricptr,alphaw)
       
    END IF
    
  END SUBROUTINE alphaweight

  SUBROUTINE vhatcomponent(                          &
       numDim,numPointsBuffer,bufferSizes,opInterval, &
       gridType,gridMetric,velDir,velocity,velHatComponent)
    
    IMPLICIT NONE 

    INTEGER(KIND=4), INTENT(IN)         :: numDim, gridType,velDir
    INTEGER(KIND=8), INTENT(IN)         :: bufferSizes(numdim),numPointsBuffer
    INTEGER(KIND=8), INTENT(IN)         :: opInterval(2*numdim)
    REAL(KIND=8),    INTENT(IN), TARGET :: gridMetric(numdim*numdim*numpointsbuffer)
    REAL(KIND=8),    INTENT(IN), TARGET :: velocity(numdim*numpointsbuffer)
    REAL(KIND=8),    INTENT(OUT)        :: velHatComponent(numpointsbuffer)


    INTEGER         :: iDim, jDim
    INTEGER(KIND=8) :: metricOffset, velOffset
    REAL(KIND=8)    :: gridScale

    REAL(KIND=8),    DIMENSION(:), POINTER :: velPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: metricPtr

    ! Get vHat component for given grid type and direction
    !
    ! If grid < rectilinear,   vhat   = constant*V(idim)
    !                          metric = (constant_x, constant_y, constant_z)
    ! If grid == rectilinear,  vhat   = metric(idim)*V(idim)
    !                          metric = <M_x,M_y,M_z> [i.e. diagonal per point]
    ! If grid == curvilinear,  vhat   = Metric(idim)*V(idim)
    !                          metric = <Mx_i Mx_j Mx_k My_i ... Mz_k> [tensor per point]
    !
!    DO iDim = 1,numDim 
!       velPtr => velocity((iDim-1)*numPointsBuffer+1:(iDim-1)*numPointsBuffer+numPointsBuffer)
!       WRITE(*,*) 'Velocity(',iDim,') = ',velPtr
!    END DO
    
    IF(gridtype < rectilinear) THEN

       gridscale = gridmetric(veldir)
       veloffset = (veldir-1)*numpointsbuffer
       velptr    => velocity(veloffset+1:veloffset+numpointsbuffer)

       CALL yax(numdim,numpointsbuffer,buffersizes,opinterval, &
            gridscale,velptr,velhatcomponent)

    ELSE IF(gridtype == rectilinear) THEN


       metricoffset =  (veldir-1)*numpointsbuffer
       metricptr    => gridmetric(metricoffset+1:metricoffset+numpointsbuffer)
       velptr       => velocity(metricoffset+1:metricoffset+numpointsbuffer)

       CALL zxy(numdim,numpointsbuffer,buffersizes,opinterval, &
            metricptr,velptr,velhatcomponent);

    ELSE ! Curvilinear

       metricoffset =  (veldir-1)*numdim*numpointsbuffer
       metricptr    => gridmetric(metricoffset+1:metricoffset+(numdim*numpointsbuffer))

       CALL zxdoty(numdim,numpointsbuffer,buffersizes,opinterval, &
            numdim,metricptr,velocity,velhatcomponent);

    END IF

!    WRITE(*,*) 'VelHat(',velDir,') = ',velHatComponent

  END SUBROUTINE vhatcomponent

 
  SUBROUTINE ijkgradtoxyzdiv(                         &
       numDim, numPoints, gridSizes, opInterval,      & 
       gridType, gridJacobian,gridMetric,gradVBuffer, & 
       divBuffer)
    
    IMPLICIT NONE

    INTEGER(KIND=4), INTENT(IN)         :: numDim, gridType
    INTEGER(KIND=8), INTENT(IN)         :: gridSizes(numdim),  numPoints
    INTEGER(KIND=8), INTENT(IN)         :: opInterval(2*numdim)
    REAL(KIND=8),    INTENT(IN)         :: gridJacobian(numpoints)
    REAL(KIND=8),    INTENT(IN), TARGET :: gridMetric(numdim*numdim*numpoints)
    REAL(KIND=8),    INTENT(IN), TARGET :: gradVBuffer(numdim*numdim*numpoints)
    REAL(KIND=8),    INTENT(OUT)        :: divBuffer(numpoints)

    INTEGER         :: iDim, jDim
    INTEGER(KIND=8) :: metricOffset, gradOffset
    REAL(KIND=8)    :: gridScale

    REAL(KIND=8),    DIMENSION(:), POINTER :: gradPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: metricPtr

    IF(gridtype > rectilinear) THEN ! full curvilinear
       DO idim = 1,numdim
          DO jdim = 1,numdim
             gradoffset   =  ((idim-1)*numdim+(jdim-1))*numpoints
             metricoffset =  ((idim-1)+(jdim-1)*numdim)*numpoints
             metricptr    => gridmetric(metricoffset+1:metricoffset+numpoints)
             gradptr      => gradvbuffer(gradoffset+1:gradoffset+numpoints)
             CALL ywxpy(numdim,numpoints,gridsizes,opinterval, &
                  metricptr,gradptr,divbuffer)
          END DO
       END DO
       CALL yxy(numdim,numpoints,gridsizes,opinterval,gridjacobian,divbuffer)
    ELSE IF(gridtype < rectilinear) THEN ! uniform rectangular
       DO idim = 1,numdim
          gridscale  = gridmetric(idim)
          gradoffset = ((idim-1)*(numdim+1))*numpoints
          gradptr    => gradvbuffer(gradoffset+1:gradoffset+numpoints)
          CALL yaxpy(numdim,numpoints,gridsizes,opinterval, &
               gridscale,gradptr, divbuffer)
       END DO
       CALL xax(numdim,numpoints,gridsizes,opinterval, &
            gridjacobian(1),divbuffer)
    ELSE ! rectilinear (weirdo)
       DO idim = 1,numdim
          metricoffset =  (idim-1)*numpoints
          metricptr    => gridmetric(metricoffset+1:metricoffset+numpoints)
          gradoffset   = ((idim-1)*(numdim+1))*numpoints
          gradptr      => gradvbuffer(gradoffset+1:gradoffset+numpoints)
          CALL ywxpy(numdim,numpoints,gridsizes,opinterval, &
               metricptr,gradptr,divbuffer)          
       END DO
       CALL yxy(numdim,numpoints,gridsizes,opinterval,gridjacobian,divbuffer)
    END IF

  END SUBROUTINE ijkgradtoxyzdiv

  SUBROUTINE gradijktogradxyz(                             &
       numDim, numPoints, gridSizes, opInterval, gridType, &
       gridJacobian,gridMetric,gradIJK, gradXYZ)
    
    IMPLICIT NONE

    INTEGER(KIND=4), INTENT(IN)         :: numDim, gridType
    INTEGER(KIND=8), INTENT(IN)         :: gridSizes(numdim),  numPoints
    INTEGER(KIND=8), INTENT(IN)         :: opInterval(2*numdim)
    REAL(KIND=8),    INTENT(IN)         :: gridJacobian(numpoints)
    REAL(KIND=8),    INTENT(IN), TARGET :: gridMetric(numdim*numdim*numpoints)
    REAL(KIND=8),    INTENT(IN), TARGET :: gradIJK(numdim*numpoints)
    REAL(KIND=8),    INTENT(OUT),TARGET :: gradXYZ(numdim*numpoints)

    INTEGER         :: iDim, jDim
    INTEGER(KIND=8) :: metricOffset, gradIJKOffset, gradXYZOffset
    REAL(KIND=8)    :: gridScale
    REAL(KIND=8)    :: zero

    REAL(KIND=8),    DIMENSION(:), POINTER :: gradIJKPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: gradXYZPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: metricPtr

    zero = 0.0_8

    IF(gridtype > rectilinear) THEN ! full curvilinear
       
       ! Note: only need to zero for curvilinear, other metric types
       ! utilize assignment-type kernels
       DO idim = 1,numdim
          gradxyzoffset =  (idim-1)*numpoints
          gradxyzptr    => gradxyz(gradxyzoffset+1:gradxyzoffset+numpoints)
          CALL assignmentxa(numdim,numpoints,gridsizes,opinterval,zero,gradxyzptr)
       END DO

       DO idim = 1,numdim
          gradxyzoffset   =  (idim-1)*numpoints
          gradxyzptr      => gradxyz(gradxyzoffset+1:gradxyzoffset+numpoints)
          DO jdim = 1,numdim
             gradijkoffset   =  (jdim-1)*numpoints
             gradijkptr      => gradijk(gradijkoffset+1:gradijkoffset+numpoints)
             metricoffset    =  ((jdim-1)*numdim+(idim-1))*numpoints
             metricptr       => gridmetric(metricoffset+1:metricoffset+numpoints)
             CALL ywxpy(numdim,numpoints,gridsizes,opinterval, &
                  metricptr,gradijkptr,gradxyzptr)
          END DO
          CALL yxy(numdim,numpoints,gridsizes,opinterval,gridjacobian,gradxyzptr)
       END DO

    ELSE IF(gridtype < rectilinear) THEN ! uniform rectangular
       DO idim = 1,numdim
          gridscale     = gridmetric(idim)*gridjacobian(1)
          gradijkoffset = (idim-1)*numpoints
          gradijkptr    => gradijk(gradijkoffset+1:gradijkoffset+numpoints)
          gradxyzptr    => gradxyz(gradijkoffset+1:gradijkoffset+numpoints)
          CALL yax(numdim,numpoints,gridsizes,opinterval, &
               gridscale,gradijkptr, gradxyzptr)
       END DO
    ELSE ! rectilinear (weirdo)
       DO idim = 1,numdim
          metricoffset =  (idim-1)*numpoints
          metricptr    => gridmetric(metricoffset+1:metricoffset+numpoints)
          gradijkptr   => gradijk(metricoffset+1:metricoffset+numpoints)
          gradxyzptr   => gradxyz(metricoffset+1:metricoffset+numpoints)
          CALL zxy(numdim,numpoints,gridsizes,opinterval, &
               metricptr,gradijkptr,gradxyzptr)          
          CALL yxy(numdim,numpoints,gridsizes,opinterval,gridjacobian,gradxyzptr)
       END DO
    END IF

  END SUBROUTINE gradijktogradxyz


END MODULE metricops