Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rank mapping confusion #132

Open
YL-codehub opened this issue Aug 26, 2024 · 2 comments
Open

Rank mapping confusion #132

YL-codehub opened this issue Aug 26, 2024 · 2 comments

Comments

@YL-codehub
Copy link
Contributor

YL-codehub commented Aug 26, 2024

Hello,
I believe there might be a mistake in the SWFFT/Poisson example.
There is a block of code dedicated to remapping the ranks from Fortran to C ones, which is what is used by the DFFT. However, I’m not sure that the original array is Fortran-ordered. In fact, if I change the sinusoidal/periodic input built in the SWFFT_Test_F.F90 file to a random number drawn on each pixel (see attached modified .F90 and input picture), the final outcome looks like something went wrong in the ranks mapping (see output picture, clearly showing that a F <-> C reordering was done but should not have been).
When keeping the same mapping of the ranks for the DFFT quantities (i.e. uncommenting that line instead of the one below), the box is what we expect; the boxes of the result end up in the right places.
My interpretation is:

  • There might have been some confusion in the use of a Fortran code vs a Fortran order for the input.
  • The example with the periodic input looks like it is working because it is periodic and so the wrong ordering doesn’t show up if some boxes are not in the right places.

I would really like to hear what you (@cgilet , @etpalmer63 ?) think as I’m not as advanced in my AMReX learning as I would like and so I wonder if this isn’t just me...
Thank you so much! Looking forward to FFTing in AMReX
Yoann.

CODE MODIF: fortran fort_init_rhs routine in SWFFT_Test_F.F90

 subroutine fort_init_rhs (lo, hi, rhs, rlo, rhi, domlo, domhi, problo, probhi, dx) &
       bind(c,name="fort_init_rhs")

    integer, intent(in) :: lo(3), hi(3), rlo(3), rhi(3), domlo(3), domhi(3)
    real(amrex_real), intent(inout) :: rhs(rlo(1):rhi(1),rlo(2):rhi(2),rlo(3):rhi(3))
    real(amrex_real), intent(in) :: problo(3), probhi(3)
    real(amrex_real), intent(in) :: dx(3)

    integer :: i,j,k
    real(amrex_real) :: x,y,z
    real(amrex_real) :: Lx,Ly,Lz
    real(amrex_real) :: pi, fpi, tpi, fac
    real(amrex_real) :: noise                       ! Added
    real(amrex_real) :: noise_amplitude     ! Added

    ! Noise amplitude can be adjusted as needed
    noise_amplitude = 0.1                         ! Added

    Lx = probhi(1)-problo(1)
    Ly = probhi(2)-problo(2)
    Lz = probhi(3)-problo(3)

    pi = 4.d0 * atan(1.d0)
    tpi = 2.0d0 * pi
    fpi = 4.0d0 * pi
    fac = 3.0d0 * tpi**2 / (Lx**2 * Ly**2 * Lz**2)

    ! Seed the random number generator (if not already seeded elsewhere)
    call random_seed()                              ! Added

    do k = lo(3), hi(3)
       z = (dble(k)+0.5d0)*dx(3)/Lz

       do j = lo(2), hi(2)
          y = (dble(j)+0.5d0)*dx(2)/Ly

          do i = lo(1), hi(1)
             x = (dble(i)+0.5d0)*dx(1)/Lx

             ! Compute the deterministic part
             ! rhs(i,j,k) = -fac * (sin(tpi*x) * sin(tpi*y) * sin(tpi*z))  &     ! Commented
             !     &       -fac * (sin(fpi*x) * sin(fpi*y) * sin(fpi*z))             ! Commented

             ! Add random noise
             call random_number(noise)       ! Added
             rhs(i,j,k) = noise_amplitude * (2.0d0 * noise - 1.0d0)        ! Added
          end do
       end do
    end do

  end subroutine fort_init_rhs

INPUT:
inputSWFFT
OUTPUT:
outputSWFFT

Compilation details:

  • COMP = intel
  • USE_MPI = True
@WeiqunZhang
Copy link
Member

@YL-codehub I believe you are correct that the tutorial is wrong.

This is how rank_mapping is used.

    int Ndims[3] = { nbz, nby, nbx };
    int     n[3] = {domain.length(2), domain.length(1), domain.length(0)};
    hacc::Distribution d(MPI_COMM_WORLD,n,Ndims,&rank_mapping[0]);

The Fortran ordering has already been handled by the ordering in Ndims[3] and n[3]. So we should provide rank_mapping such that rank_mapping[0] is for box (0,0,0), and rank_mapping[1] for box (0,0,1). Because of the ordering passed to hacc functions, (0,0,1) is actually i=1, j=0, k=0 in that loop for (int ib = 0; ib < nboxes; ++ib).

@WeiqunZhang
Copy link
Member

@YL-codehub Please submit a PR to fix it. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants
@WeiqunZhang @YL-codehub and others