Parallel Integral

[MarcoArtiano]

Hello, I writing a function who performs the easiest approximation possible to do 3D integral, but I'm getting some weird results: whenever I change the number of processor in the parallel computation I get quite different results and don't know why. I guess the problem is not in the 3D vector inputs but in the way I'm managing data in order to get the 3d scalar at each time step. The function I wrote looks like this:

subroutine Reeffcomp(epsni,epsrk,Reeff2)

    USE param
    USE variables
    USE var, ONLY: px1, py1, pz1
    USE decomp_2d
    USE MPI
    implicit none

    !! INPUTS
real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: epsni,epsrk
integer :: j,i,k,code
real(mytype) :: epsnism, epsrksm, epsrks, epsnis,dxv,dyv,dzv, Reeff2

epsnism = zero
epsnis = zero
epsrks = zero
epsrksm = zero
Reeff2 = zero

do k=xstart(3),xsize(3)

      if (k.ne.one.and.k.ne.nz) then
       dzv = dz
      else
       dzv = dz*half
      endif



        do j=xstart(2),xend(2)

         if (j.eq.one) then
               dyv = dyp(j)*half
         elseif (j.eq.ny) then
               dyv = dyp(j)*half
         else
               dyv = dyp(j)
         endif

         do i = xstart(1),xend(1)

      if (i.ne.one.and.i.ne.nx) then
       dxv = dx
      else
       dxv = dx*half
      endif

          epsnism = epsnism + epsni(i,j,k)*dxv*dyv*dzv
          epsrksm = epsrksm + epsrk(i,j,k)*dxv*dyv*dzv
          enddo
      enddo
  enddo

    call MPI_ALLREDUCE(epsnism,epsnis,1,real_type,MPI_SUM,MPI_COMM_WORLD,code)
    call MPI_ALLREDUCE(epsrksm,epsrks,1,real_type,MPI_SUM,MPI_COMM_WORLD,code)

epsnis = epsnis/dt
Reeff = epsnis/(-half*epsrks/dt+epsnis)

    if (mod(itime, ioutput).eq.0.and.nrank.eq.0) then
       print*, 'Reeff= ', Reeff
    endif


endsubroutine Reeffcomp

Maybe I'm wondering that do statement is wrong and I should write something like "do i = 1,xsize(1)" and same for j and k. Thanks for the help

ps: the function is called inside the main xcompact3d.f90, and I'm 100% sure that the inputs epsni and epsrk have no problem, then the scalar Reeff is saved as follows
if (mod(itime, ioutput).eq.0) then call Reeffcomp(epsni,epsrk,Reeff) Reeffvec(itime) = Reeff endif
and the write into a file to read at the end of the time cycle
if (nrank == 0) then open(10,file='Reeff.dat', form='formatted') do itime=ifirst,ilast write(10,*)Reeffvec(itime) enddo close(10) endif

[mathrack]

Hello,
Have you looked at the estimation of the bulk velocity in the channel flow subroutine ? See https://github.com/xcompact3d/Incompact3d/blob/master/src/BC-Channel-flow.f90#L320
In your case, some minor modification is probably needed at i=1, i=nx and when xstart(3)=1 or xend(3)=nz. But it should be a good starting point.

EDIT : avoid the comparison between i (k as well) and "one" because "one" is a floating point number

EDIT2 : the loops do k=xstart(3),xsize(3) and do j=xstart(2),xend(2) are not consistent with epsni(i,j,k). It should be epsni(i,j-xstart(2)+1,k-xstart(3)+1). The i loop is ok because xstart(1)=1 and xend(1)=xsize(1)=nx. Also, the ending index of the k loop is wrong.

[MarcoArtiano]

Thank you so much, now it works. I'm so grateful.