! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! Maintainers : support@fluidnumerics.com ! Official Repository : https://github.com/FluidNumerics/self/ ! ! Copyright © 2024 Fluid Numerics LLC ! ! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: ! ! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in ! the documentation and/or other materials provided with the distribution. ! ! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF ! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! module SELF_Vector_3D_t use SELF_Constants use SELF_Lagrange use SELF_Metadata use FEQParse use SELF_HDF5 use SELF_Data use SELF_DataPool use HDF5 use iso_c_binding implicit none type,extends(SELF_DataObj),public :: Vector3D_t real(prec),pointer,contiguous,dimension(:,:,:,:,:,:) :: interior real(prec),pointer,contiguous,dimension(:,:,:,:,:,:) :: boundary real(prec),pointer,contiguous,dimension(:,:,:,:,:,:) :: extBoundary real(prec),pointer,contiguous,dimension(:,:,:,:,:,:) :: avgBoundary real(prec),pointer,contiguous,dimension(:,:,:,:,:) :: boundaryNormal !! High-water-mark backing store for the arrays above (AMR Stage 6b; see SELF_DataPool). !! The public members are pointers remapped onto the leading part of these pools at the !! exact logical shape, so every extent, stride and shape() is unchanged while Resize can !! reuse the storage across an adaptation epoch. Nothing should index the pools directly. real(prec),pointer,contiguous :: pool_interior(:) => null() real(prec),pointer,contiguous :: pool_boundary(:) => null() real(prec),pointer,contiguous :: pool_extBoundary(:) => null() real(prec),pointer,contiguous :: pool_avgBoundary(:) => null() real(prec),pointer,contiguous :: pool_boundaryNormal(:) => null() contains procedure,public :: Init => Init_Vector3D_t procedure,public :: Resize => Resize_Vector3D_t procedure,public :: MapArrays => MapArrays_Vector3D_t procedure,public :: Free => Free_Vector3D_t procedure,public :: UpdateHost => UpdateHost_Vector3D_t procedure,public :: UpdateDevice => UpdateDevice_Vector3D_t procedure,public :: BoundaryInterp => BoundaryInterp_Vector3D_t procedure,public :: AverageSides => AverageSides_Vector3D_t generic,public :: GridInterp => GridInterp_Vector3D_t procedure,private :: GridInterp_Vector3D_t generic,public :: Gradient => Gradient_Vector3D_t procedure,private :: Gradient_Vector3D_t generic,public :: Curl => Curl_Vector3D_t procedure,private :: Curl_Vector3D_t generic,public :: Divergence => Divergence_Vector3D_t procedure,private :: Divergence_Vector3D_t generic,public :: SetEquation => SetEquation_Vector3D_t procedure,private :: SetEquation_Vector3D_t generic,public :: WriteHDF5 => WriteHDF5_MPI_Vector3D_t,WriteHDF5_Vector3D_t procedure,private :: WriteHDF5_MPI_Vector3D_t procedure,private :: WriteHDF5_Vector3D_t endtype Vector3D_t contains subroutine Init_Vector3D_t(this,interp,nVar,nElem) implicit none class(Vector3D_t),intent(out) :: this type(Lagrange),target,intent(in) :: interp integer,intent(in) :: nVar integer,intent(in) :: nElem ! local integer :: i this%interp => interp this%nVar = nVar this%nElem = nElem this%N = interp%N this%M = interp%M call this%MapArrays(interp%N+1,nVar,nElem) allocate(this%meta(1:nVar)) allocate(this%eqn(1:3*nVar)) ! Initialize equation parser ! This is done to prevent segmentation faults that arise ! when building with amdflang that are traced back to ! feqparse_functions.f90 : finalize routine ! When the equation parser is not initialized, the ! functions are not allocated, which I think are the ! source of the segfault - joe@fluidnumerics.com do i = 1,3*nvar this%eqn(i) = EquationParser('f=0',(/'x','y','z','t'/)) enddo this%interior = 0.0_prec this%boundary = 0.0_prec this%boundarynormal = 0.0_prec this%extBoundary = 0.0_prec this%avgBoundary = 0.0_prec endsubroutine Init_Vector3D_t subroutine Free_Vector3D_t(this) implicit none class(Vector3D_t),intent(inout) :: this this%interp => null() this%nVar = 0 this%nElem = 0 ! The public arrays point into the pools rather than owning storage (AMR Stage 6b), so they ! are nullified and the pools released. this%interior => null() this%boundary => null() this%boundaryNormal => null() this%extBoundary => null() this%avgBoundary => null() if(associated(this%pool_interior)) deallocate(this%pool_interior) if(associated(this%pool_boundary)) deallocate(this%pool_boundary) if(associated(this%pool_boundaryNormal)) deallocate(this%pool_boundaryNormal) if(associated(this%pool_extBoundary)) deallocate(this%pool_extBoundary) if(associated(this%pool_avgBoundary)) deallocate(this%pool_avgBoundary) deallocate(this%meta) deallocate(this%eqn) endsubroutine Free_Vector3D_t subroutine MapArrays_Vector3D_t(this,Np,nVar,nElem) !! Size the backing pools for (Np,nVar,nElem) and remap the public arrays onto them at the !! exact logical shape. See SELF_DataPool for why this is done with pointer remapping rather !! than by over-allocating the element dimension. implicit none class(Vector3D_t),intent(inout) :: this integer,intent(in) :: Np integer,intent(in) :: nVar integer,intent(in) :: nElem ! Local integer :: nInt,nBnd,nNrm nInt = Np*Np*Np*nElem*nVar*3 nBnd = Np*Np*6*nElem*nVar*3 nNrm = Np*Np*6*nElem*nVar call EnsurePool(this%pool_interior,nInt) call EnsurePool(this%pool_boundary,nBnd) call EnsurePool(this%pool_extBoundary,nBnd) call EnsurePool(this%pool_avgBoundary,nBnd) call EnsurePool(this%pool_boundaryNormal,nNrm) this%interior(1:Np,1:Np,1:Np,1:nElem,1:nVar,1:3) => this%pool_interior(1:nInt) this%boundary(1:Np,1:Np,1:6,1:nElem,1:nVar,1:3) => this%pool_boundary(1:nBnd) this%extBoundary(1:Np,1:Np,1:6,1:nElem,1:nVar,1:3) => this%pool_extBoundary(1:nBnd) this%avgBoundary(1:Np,1:Np,1:6,1:nElem,1:nVar,1:3) => this%pool_avgBoundary(1:nBnd) this%boundaryNormal(1:Np,1:Np,1:6,1:nElem,1:nVar) => this%pool_boundaryNormal(1:nNrm) endsubroutine MapArrays_Vector3D_t subroutine Resize_Vector3D_t(this,interp,nVar,nElem) !! Rebind a live object to a new element count, reusing existing storage when it fits (AMR !! Stage 6b). Unlike Init - which is intent(out), so it resets the object, reallocates, !! zeroes and reconstructs the equation parsers - this preserves metadata and parsers !! (neither depends on nElem) and touches storage only when it must grow. Preserving the !! parsers also avoids repeating the EquationParser construction that Init performs 3*nVar !! times as an amdflang workaround. !! !! All arrays are zeroed, exactly as Init leaves them - see Resize_Scalar3D_t for why that is !! required rather than merely tidy. implicit none class(Vector3D_t),intent(inout) :: this type(Lagrange),target,intent(in) :: interp integer,intent(in) :: nVar integer,intent(in) :: nElem if(nVar /= this%nVar) then print*,__FILE__,':',__LINE__, & ' : Error : Resize cannot change nVar. Use Free followed by Init.' stop 1 endif this%interp => interp this%nElem = nElem this%N = interp%N this%M = interp%M call this%MapArrays(interp%N+1,nVar,nElem) this%interior = 0.0_prec this%boundary = 0.0_prec this%boundarynormal = 0.0_prec this%extBoundary = 0.0_prec this%avgBoundary = 0.0_prec endsubroutine Resize_Vector3D_t subroutine UpdateHost_Vector3D_t(this) implicit none class(Vector3D_t),intent(inout) :: this if(.false.) this%N = this%N ! CPU stub; suppress unused-dummy-argument warning endsubroutine UpdateHost_Vector3D_t subroutine UpdateDevice_Vector3D_t(this) implicit none class(Vector3D_t),intent(inout) :: this if(.false.) this%N = this%N ! CPU stub; suppress unused-dummy-argument warning endsubroutine UpdateDevice_Vector3D_t subroutine SetEquation_Vector3D_t(this,idir,ivar,eqnChar) !! Sets the equation parser for the `idir` direction and `ivar-th` variable implicit none class(Vector3D_t),intent(inout) :: this integer,intent(in) :: idir,ivar character(*),intent(in) :: eqnChar this%eqn(idir+3*(ivar-1)) = EquationParser(trim(eqnChar), & (/'x','y','z','t'/)) endsubroutine SetEquation_Vector3D_t subroutine GridInterp_Vector3D_t(this,f) implicit none class(Vector3D_t),intent(in) :: this real(prec),intent(out) :: f(1:this%M+1,1:this%M+1,1:this%M+1,1:this%nelem,1:this%nvar,1:3) !! (Output) Array of function values, defined on the target grid ! Local integer :: i,j,k,ii,jj,kk,iel,ivar,idir real(prec) :: fi,fij,fijk do concurrent(i=1:this%M+1,j=1:this%M+1, & k=1:this%M+1,iel=1:this%nelem,ivar=1:this%nvar,idir=1:3) fijk = 0.0_prec do kk = 1,this%N+1 fij = 0.0_prec do jj = 1,this%N+1 fi = 0.0_prec do ii = 1,this%N+1 fi = fi+this%interior(ii,jj,kk,iel,ivar,idir)*this%interp%iMatrix(ii,i) enddo fij = fij+fi*this%interp%iMatrix(jj,j) enddo fijk = fijk+fij*this%interp%iMatrix(kk,k) enddo f(i,j,k,iel,ivar,idir) = fijk enddo endsubroutine GridInterp_Vector3D_t subroutine AverageSides_Vector3D_t(this) implicit none class(Vector3D_t),intent(inout) :: this ! Local integer :: iel integer :: iside integer :: ivar integer :: i,j integer :: idir do concurrent(i=1:this%N+1,j=1:this%N+1, & iside=1:6,iel=1:this%nelem,ivar=1:this%nvar, & idir=1:3) this%boundary(i,j,iside,iel,ivar,idir) = 0.5_prec*( & this%boundary(i,j,iside,iel,ivar,idir)+ & this%extBoundary(i,j,iside,iel,ivar,idir)) enddo endsubroutine AverageSides_Vector3D_t subroutine BoundaryInterp_Vector3D_t(this) implicit none class(Vector3D_t),intent(inout) :: this ! Local integer :: i,j,ii,idir,iel,ivar real(prec) :: fbb,fbs,fbe,fbn,fbw,fbt do concurrent(i=1:this%N+1,j=1:this%N+1, & iel=1:this%nelem,ivar=1:this%nvar,idir=1:3) fbb = 0.0_prec fbs = 0.0_prec fbe = 0.0_prec fbn = 0.0_prec fbw = 0.0_prec fbt = 0.0_prec do ii = 1,this%N+1 fbb = fbb+this%interp%bMatrix(ii,1)*this%interior(i,j,ii,iel,ivar,idir) ! Bottom fbs = fbs+this%interp%bMatrix(ii,1)*this%interior(i,ii,j,iel,ivar,idir) ! South fbe = fbe+this%interp%bMatrix(ii,2)*this%interior(ii,i,j,iel,ivar,idir) ! East fbn = fbn+this%interp%bMatrix(ii,2)*this%interior(i,ii,j,iel,ivar,idir) ! North fbw = fbw+this%interp%bMatrix(ii,1)*this%interior(ii,i,j,iel,ivar,idir) ! West fbt = fbt+this%interp%bMatrix(ii,2)*this%interior(i,j,ii,iel,ivar,idir) ! Top enddo this%boundary(i,j,1,iel,ivar,idir) = fbb this%boundary(i,j,2,iel,ivar,idir) = fbs this%boundary(i,j,3,iel,ivar,idir) = fbe this%boundary(i,j,4,iel,ivar,idir) = fbn this%boundary(i,j,5,iel,ivar,idir) = fbw this%boundary(i,j,6,iel,ivar,idir) = fbt enddo endsubroutine BoundaryInterp_Vector3D_t subroutine Gradient_Vector3D_t(this,df) implicit none class(Vector3D_t),intent(in) :: this real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:3,1:3) ! Local integer :: i,j,k,ii,idir,iel,ivar real(prec) :: dfds1,dfds2,dfds3 do concurrent(i=1:this%N+1,j=1:this%N+1, & k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar, & idir=1:3) dfds1 = 0.0_prec dfds2 = 0.0_prec dfds3 = 0.0_prec do ii = 1,this%N+1 dfds1 = dfds1+this%interp%dMatrix(ii,i)*this%interior(ii,j,k,iel,ivar,idir) dfds2 = dfds2+this%interp%dMatrix(ii,j)*this%interior(i,ii,k,iel,ivar,idir) dfds3 = dfds3+this%interp%dMatrix(ii,k)*this%interior(i,j,ii,iel,ivar,idir) enddo df(i,j,k,iel,ivar,idir,1) = dfds1 df(i,j,k,iel,ivar,idir,2) = dfds2 df(i,j,k,iel,ivar,idir,3) = dfds3 enddo endsubroutine Gradient_Vector3D_t subroutine Curl_Vector3D_t(this,curlf) implicit none class(Vector3D_t),intent(in) :: this real(prec),intent(out) :: curlf(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:3) ! Local integer :: i,j,k,ii,idir,iel,ivar real(prec) :: dfds1,dfds2,dfds3 real(prec) :: df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar,1:3,1:3) do concurrent(i=1:this%N+1,j=1:this%N+1, & k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar, & idir=1:3) dfds1 = 0.0_prec dfds2 = 0.0_prec dfds3 = 0.0_prec do ii = 1,this%N+1 dfds1 = dfds1+this%interp%dMatrix(ii,i)*this%interior(ii,j,k,iel,ivar,idir) dfds2 = dfds2+this%interp%dMatrix(ii,j)*this%interior(i,ii,k,iel,ivar,idir) dfds3 = dfds3+this%interp%dMatrix(ii,k)*this%interior(i,j,ii,iel,ivar,idir) enddo df(i,j,k,iel,ivar,idir,1) = dfds1 df(i,j,k,iel,ivar,idir,2) = dfds2 df(i,j,k,iel,ivar,idir,3) = dfds3 enddo do concurrent(i=1:this%N+1,j=1:this%N+1, & k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar) curlf(i,j,k,iel,ivar,1) = (df(i,j,k,iel,ivar,3,2)-df(i,j,k,iel,ivar,2,3)) curlf(i,j,k,iel,ivar,2) = (df(i,j,k,iel,ivar,1,3)-df(i,j,k,iel,ivar,3,1)) curlf(i,j,k,iel,ivar,3) = (df(i,j,k,iel,ivar,2,1)-df(i,j,k,iel,ivar,1,2)) enddo endsubroutine Curl_Vector3D_t subroutine Divergence_Vector3D_t(this,df) implicit none class(Vector3D_t),intent(in) :: this real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nelem,1:this%nvar) ! Local integer :: i,j,k,ii,iel,ivar real(prec) :: dfLoc do concurrent(i=1:this%N+1,j=1:this%N+1, & k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar) dfLoc = 0.0_prec do ii = 1,this%N+1 dfLoc = dfLoc+this%interp%dMatrix(ii,i)*this%interior(ii,j,k,iel,ivar,1) enddo dF(i,j,k,iel,ivar) = dfLoc enddo do concurrent(i=1:this%N+1,j=1:this%N+1, & k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar) dfLoc = 0.0_prec do ii = 1,this%N+1 dfLoc = dfLoc+this%interp%dMatrix(ii,j)*this%interior(i,ii,k,iel,ivar,2) enddo dF(i,j,k,iel,ivar) = dF(i,j,k,iel,ivar)+dfLoc enddo do concurrent(i=1:this%N+1,j=1:this%N+1, & k=1:this%N+1,iel=1:this%nelem,ivar=1:this%nvar) dfLoc = 0.0_prec do ii = 1,this%N+1 dfLoc = dfLoc+this%interp%dMatrix(ii,k)*this%interior(i,j,ii,iel,ivar,3) enddo dF(i,j,k,iel,ivar) = dF(i,j,k,iel,ivar)+dfLoc enddo endsubroutine Divergence_Vector3D_t subroutine WriteHDF5_MPI_Vector3D_t(this,fileId,group,elemoffset,nglobalelem) implicit none class(Vector3D_t),intent(in) :: this character(*),intent(in) :: group integer(HID_T),intent(in) :: fileId integer,intent(in) :: elemoffset integer,intent(in) :: nglobalelem ! Local integer(HID_T) :: offset(1:4) integer(HID_T) :: globalDims(1:4) integer :: ivar,idir character(4) :: dimvar offset(1:4) = (/0,0,0,elemoffset/) globalDims(1:4) = (/this%interp%N+1, & this%interp%N+1, & this%interp%N+1, & nglobalelem/) call CreateGroup_HDF5(fileId,trim(group)) do idir = 1,3 write(dimvar,'(I1)') idir dimvar = "dim"//trim(dimvar) do ivar = 1,this%nVar call WriteArray_HDF5(fileId, & trim(group)//"/"//trim(this%meta(ivar)%name)//"_"//dimvar, & this%interior(:,:,:,:,ivar,idir),offset,globalDims) enddo enddo endsubroutine WriteHDF5_MPI_Vector3D_t subroutine WriteHDF5_Vector3D_t(this,fileId,group) implicit none class(Vector3D_t),intent(in) :: this integer(HID_T),intent(in) :: fileId character(*),intent(in) :: group ! Local integer :: ivar,idir character(4) :: dimvar call CreateGroup_HDF5(fileId,trim(group)) do ivar = 1,this%nVar call this%meta(ivar)%WriteHDF5(group,ivar,fileId) enddo do idir = 1,3 write(dimvar,'(I1)') idir dimvar = "dim"//trim(dimvar) do ivar = 1,this%nVar call WriteArray_HDF5(fileId, & trim(group)//"/"//trim(this%meta(ivar)%name)//"_"//dimvar, & this%interior(:,:,:,:,ivar,idir)) enddo enddo endsubroutine WriteHDF5_Vector3D_t endmodule SELF_Vector_3D_t