!> @author Google DeepMind Team & Ignis !> @brief DNS post-processing mathematical operations and derivatives, matching appropriate compact schemes in X, Y, Z directions for a defined flow domain and optimizing calculations from a cache perspective. !! !! This module calculates spatial derivatives (first and second spatial derivatives in X, Y, Z directions like ddx, ddy, ddz, etc.) !! using high-order compact finite difference schemes. It also calculates chemical reaction rates, !! positive/negative component extraction, and threshold operations, incorporating memory cache !! efficiency (cache blocking) and transpose optimization operations (`tp2`). module m_calculate use, intrinsic :: iso_fortran_env, only: real64 use Compact use m_parameters implicit none real(real64), allocatable, dimension(:,:) :: xsrc real(real64), allocatable, dimension(:,:) :: xdst real(real64), allocatable, dimension(:,:) :: rsrc real(real64), allocatable, dimension(:,:) :: rdst integer, parameter :: nb = BLOCKSIZE private :: nb private :: xsrc, xdst, rsrc, rdst contains !> Initializes the workspace arrays and matrices. !! !! Pre-allocates transposed workspace buffers (`xsrc`, `xdst`, `rsrc`, `rdst`) !! and performs the LU decomposition setup via `ludcmp` for the tridiagonal compact schemes. subroutine m_calculate_init integer :: ierr call ludcmp(nxp,nyp,nzp,1,0,0) ! 1,1,0 allocate(xsrc(nb, nxp), stat=ierr) if (ierr /= 0) then write(0,*) "Error: allocation of xsrc failed on process", myid call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err) end if allocate(xdst(nb, nxp), stat=ierr) if (ierr /= 0) then write(0,*) "Error: allocation of xdst failed on process", myid call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err) end if allocate(rsrc(nxp, nzp), stat=ierr) if (ierr /= 0) then write(0,*) "Error: allocation of rsrc failed on process", myid call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err) end if allocate(rdst(nxp, nzp), stat=ierr) if (ierr /= 0) then write(0,*) "Error: allocation of rdst failed on process", myid call MPI_ABORT(MPI_COMM_TASK, 1, mpi_err) end if end subroutine m_calculate_init subroutine m_calculate_finalize deallocate(xsrc) deallocate(xdst) deallocate(rsrc) deallocate(rdst) end subroutine m_calculate_finalize subroutine ddx1d(dst, src) real(real64), dimension(1,nxp), intent(in) :: src real(real64), dimension(1,nxp), intent(out) :: dst call dfnonp(nxp, hxp, src, dst, 1, 1) end subroutine ddx1d !> Computes the first-order derivative in the X-direction. !! !! @param src 3D input scalar field (nxp, nyp, nzp). !! @param dst 3D output derivative field (nxp, nyp, nzp). subroutine ddx(dst, src) real(real64), dimension(nxp,nyp,nzp), intent(in) :: src real(real64), dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k integer :: ju do k = 1,nzp do j = 1,nyp,nb ju = min(j+nb-1,nyp) call tp2(xsrc, src(:,j:ju,k), nb, nxp) call dfnonp(nxp, hxp, xsrc, xdst, nb, 1) call tp2(dst(:,j:ju,k), xdst, nxp, nb) end do end do end subroutine ddx !> Computes the first-order derivative in the Y-direction. !! !! @param src 3D input scalar field (nxp, nyp, nzp). !! @param dst 3D output derivative field (nxp, nyp, nzp). subroutine ddy(dst, src) real(real64), dimension(nxp,nyp,nzp), intent(in) :: src real(real64), dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k do k = 1,nzp call dfp(nyp, hyp, src(:,:,k), dst(:,:,k), nxp, 2) end do end subroutine ddy !> Computes the first-order derivative in the Z-direction. !! !! @param src 3D input scalar field (nxp, nyp, nzp). !! @param dst 3D output derivative field (nxp, nyp, nzp). subroutine ddz(dst, src) real(real64), dimension(nxp,nyp,nzp), intent(in) :: src real(real64), dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k do j = 1,nyp do k = 1,nzp rsrc(:,k) = src(:,j,k) end do call dfp(nzp, hzp, rsrc, rdst, nxp, 3) do k = 1,nzp dst(:,j,k) = rdst(:,k) end do end do end subroutine ddz subroutine d2dx1d(dst, src) real(real64), dimension(nxp), intent(in) :: src real(real64), dimension(nxp), intent(out) :: dst call d2fnonp(nxp, hxp, src, dst, 1, 1) end subroutine d2dx1d subroutine d2dx(dst, src) real(real64), dimension(nxp,nyp,nzp), intent(in) :: src real(real64), dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k integer :: ju do k = 1,nzp do j = 1,nyp,nb ju = min(j+nb-1,nyp) call tp2(xsrc, src(:,j:ju,k), nb, nxp) call d2fnonp(nxp, hxp, xsrc, xdst, nb, 1) call tp2(dst(:,j:ju,k), xdst, nxp, nb) end do end do end subroutine d2dx subroutine d2dy(dst, src) real(real64), dimension(nxp,nyp,nzp), intent(in) :: src real(real64), dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k do k = 1,nzp call d2fp(nyp, hyp, src(:,:,k), dst(:,:,k), nxp, 2) end do end subroutine d2dy subroutine d2dz(dst, src) real(real64), dimension(nxp,nyp,nzp), intent(in) :: src real(real64), dimension(nxp,nyp,nzp), intent(out) :: dst integer :: i, j ,k do j = 1,nyp do k = 1,nzp rsrc(:,k) = src(:,j,k) end do call d2fp(nzp, hzp, rsrc, rdst, nxp, 3) do k = 1,nzp dst(:,j,k) = rdst(:,k) end do end do end subroutine d2dz subroutine tp(a, b, nx) ! a(nb,nx) = transpose(b(nx,nb)) integer,intent(in) :: nx real(real64),intent(out) :: a(nb,nx) real(real64),intent(in) :: b(nx,nb) call tp2(a, b, nb, nx) end subroutine tp subroutine tp2 (a, b, n1, n2) ! a = transpose(b) implicit none integer,intent(in) :: n1, n2 real(real64),intent(out) :: a(n1,n2) real(real64),intent(in) :: b(n2,n1) integer :: i,j,ii,jj DO jj=1,n2,nb DO ii=1,n1,nb DO j=jj,min(jj+nb-1,n2) DO i=ii,min(ii+nb-1,n1) a(i,j) = b(j,i) ENDDO ENDDO ENDDO ENDDO end subroutine tp2 !> Computes the chemical reaction rate based on the progress variable c. !! !! This uses a piecewise exponential/Arrhenius model depending on whether the progress variable !! is below c_cut, above c_ref, or intermediate. !! !! @param c The progress variable (0.0 to 1.0). !! @return The computed reaction rate. real(real64) function rxn_rate (c) real(real64) :: c if(c.lt.0._real64) c=0._real64 if(c.gt.1._real64) c=1._real64 if (c.le.c_cut) then rxn_rate = min_wr else if (c.gt.c_ref) then rxn_rate = pre*(1.-c)*exp(-ac/(1.+bc*c)) else rxn_rate = & ((refwr-min_wr)*exp(prof_wr*(c-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) & / (1.-exp(prof_wr*(c_cut-c_ref))) endif end function rxn_rate real(real64) function threshold_min_max (c, minc, maxc) real(real64) :: c real(real64) :: minc, maxc if ((c.lt.minc) .or. (c.gt.maxc)) then threshold_min_max = 0._real64 else threshold_min_max = 1.0_real64 end if end function threshold_min_max real(real64) function positive (c) real(real64) :: c if (c > 0.0_real64) then positive = c else positive = 0._real64 end if end function positive real(real64) function negative (c) real(real64) :: c if (c < 0.0_real64) then negative = c else negative = 0._real64 end if end function negative end module m_calculate