BEGIN_PROVIDER [ integer, grid_x_num ] &BEGIN_PROVIDER [ integer, grid_y_num ] &BEGIN_PROVIDER [ integer, grid_z_num ] &BEGIN_PROVIDER [ real , grid_step , (3) ] &BEGIN_PROVIDER [ real , grid_origin, (3) ] real, parameter :: UNDEFINED=123456789.123456789 BEGIN_DOC ! Number of grid points in x, y, z directions END_DOC integer :: Npoints(3) real :: step_size(3) real :: origin(3) real :: opposite(3) Npoints (:) = 80 origin (:) = UNDEFINED opposite (:) = UNDEFINED step_size(:) = UNDEFINED call get_grid_point_num(Npoints) call get_grid_origin(origin) call get_grid_opposite(opposite) call get_grid_step_size(step_size) if (origin(1) == UNDEFINED) then integer :: i,l do l=1,3 origin(l) = nucl_coord(1,l) do i=2,nucl_num origin(l) = min(origin(l),nucl_coord(i,l)) enddo origin(l) = origin(l) - 4. enddo endif if (opposite(1) == UNDEFINED) then do l=1,3 opposite(l) = nucl_coord(1,l) do i=2,nucl_num opposite(l) = max(opposite(l),nucl_coord(i,l)) enddo opposite(l) = opposite(l) + 4. enddo endif if (step_size(1) == UNDEFINED) then do l=1,3 step_size(l) = (opposite(l) - origin(l))/float(Npoints(l)) enddo endif do l=1,3 grid_origin(l) = origin(l) grid_step(l) = step_size(l) enddo grid_x_num = Npoints(1) grid_y_num = Npoints(2) grid_z_num = Npoints(3) call set_grid_point_num(Npoints) call set_grid_origin(origin) call set_grid_opposite(opposite) call set_grid_step_size(step_size) END_PROVIDER BEGIN_TEMPLATE BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ] implicit none BEGIN_DOC ! $X on a grid END_DOC IRP_IF MPI include 'mpif.h' IRP_ENDIF integer :: ix, iy, iz integer :: ibegin, iend real, parameter :: UNDEFINED=123456789.123456789 grid_$X(1,1,1) = UNDEFINED call get_grid_data_$X(grid_$X) if (grid_$X(1,1,1) /= UNDEFINED) then return endif do iz=1,grid_z_num do iy=1,grid_y_num do ix=1,grid_x_num grid_$X(ix,iy,iz) = 0. enddo enddo enddo if (mpi_master) then print *, 'Computing $X' endif integer :: icount icount = mpi_size do iz=1,grid_z_num if (mpi_master) then print *, int(100*dble(iz)/dble(grid_z_num)), '%' endif point(3) = grid_origin(3)+(iz-1)*grid_step(3) do iy=1,grid_y_num point(2) = grid_origin(2)+(iy-1)*grid_step(2) do ix=1,grid_x_num icount -= 1 if (icount == mpi_rank) then point(1) = grid_origin(1)+(ix-1)*grid_step(1) TOUCH point grid_$X(ix,iy,iz) = $X_value_p endif if (icount == 0) then icount = mpi_size endif enddo enddo enddo IRP_IF MPI integer :: dim, ierr do iz=1,grid_z_num real :: buffer(grid_x_num*grid_y_num+1) icount = 0 do iy=1,grid_y_num do ix=1,grid_x_num buffer(icount+ix) = grid_$X(ix,iy,iz) enddo icount = icount + grid_x_num enddo dim = grid_x_num * grid_y_num call MPI_REDUCE(buffer,grid_$X(1,1,iz),dim,mpi_real, & mpi_sum,0,MPI_COMM_WORLD,ierr) if (ierr /= MPI_SUCCESS) then call abrt(irp_here,'Unable to fetch buffer') endif call barrier enddo IRP_ENDIF call set_grid_data_$X(grid_$X) END_PROVIDER BEGIN_PROVIDER [ real, grid_$X_grad, (grid_x_num,grid_y_num,grid_z_num,4) ] &BEGIN_PROVIDER [ real, grid_$X_lapl, (grid_x_num,grid_y_num,grid_z_num) ] BEGIN_DOC ! Laplacian of $X on a grid END_DOC implicit none BEGIN_DOC ! Gradient and lapacian of $X on a grid. 4th dimension of the grad is its norm. END_DOC IRP_IF MPI include 'mpif.h' IRP_ENDIF integer :: ix, iy, iz, it integer :: ibegin, iend real, parameter :: UNDEFINED=123456789.123456789 grid_$X_lapl(1,1,1) = UNDEFINED call get_grid_data_$X_lapl(grid_$X) if (grid_$X_lapl(1,1,1) /= UNDEFINED) then return endif do iz=1,grid_z_num do iy=1,grid_y_num do ix=1,grid_x_num do it=1,4 grid_$X_grad(ix,iy,iz,it) = 0. enddo grid_$X_lapl(ix,iy,iz) = 0. enddo enddo enddo if (mpi_master) then print *, 'Computing $X' endif integer :: icount icount = mpi_size do iz=1,grid_z_num if (mpi_master) then print *, int(100*dble(iz)/dble(grid_z_num)), '%' endif point(3) = grid_origin(3)+(iz-1)*grid_step(3) do iy=1,grid_y_num point(2) = grid_origin(2)+(iy-1)*grid_step(2) do ix=1,grid_x_num icount = icount-1 if (icount == mpi_rank) then point(1) = grid_origin(1)+(ix-1)*grid_step(1) TOUCH point do it=1,3 grid_$X_grad(ix,iy,iz,it) = $X_grad_p(it) enddo grid_$X_grad(ix,iy,iz,4) = $X_grad_p(1)**2 + $X_grad_p(2)**2 + $X_grad_p(3)**2 grid_$X_lapl(ix,iy,iz) = $X_lapl_p grid_$X_grad(ix,iy,iz,4) = sqrt(grid_$X_grad(ix,iy,iz,4)) endif if (icount == 0) then icount = mpi_size endif enddo enddo enddo IRP_IF MPI integer :: dim, ierr do it=1,4 do iz=1,grid_z_num real :: buffer(grid_x_num*grid_y_num) icount = 0 do iy=1,grid_y_num do ix=1,grid_x_num buffer(icount+ix) = grid_$X_grad(ix,iy,iz,it) enddo icount = icount + grid_x_num enddo dim = grid_x_num * grid_y_num call MPI_REDUCE(buffer,grid_$X_grad(1,1,iz,it),dim,mpi_real, & mpi_sum,0,MPI_COMM_WORLD,ierr) icount = 0 do iy=1,grid_y_num do ix=1,grid_x_num buffer(icount+ix) = grid_$X_lapl(ix,iy,iz) enddo icount = icount + grid_x_num enddo dim = grid_x_num * grid_y_num call MPI_REDUCE(buffer,grid_$X_lapl(1,1,iz),dim,mpi_real, & mpi_sum,0,MPI_COMM_WORLD,ierr) enddo enddo IRP_ENDIF call set_grid_data_$X_grad(grid_$X_grad) call set_grid_data_$X_lapl(grid_$X_lapl) END_PROVIDER SUBST [ X ] eplf;; elf;; density;; END_TEMPLATE BEGIN_TEMPLATE BEGIN_PROVIDER [ real, grid_$X_partition, (grid_x_num,grid_y_num,grid_z_num) ] implicit none BEGIN_DOC ! Create the topological partition of $X END_DOC integer :: ix, iy, iz real, parameter :: UNDEFINED=123456789.123456789 if (mpi_master) then grid_$X_partition(1,1,1) = UNDEFINED call get_grid_data_$X_partition(grid_$X_partition) if (grid_$X_partition(1,1,1) /= UNDEFINED) then return endif do iz=1,grid_z_num do iy=1,grid_y_num do ix=1,grid_x_num grid_$X_partition(ix,iy,iz) = 0. enddo enddo enddo integer :: icount, imax imax = 100 do iz=1,grid_z_num do iy=1,grid_y_num do ix=1,grid_x_num call find_attribution(ix,iy,iz,grid_$X,grid_$X_partition,imax) enddo enddo enddo call set_grid_data_$X_partition(grid_$X_partition) endif print *, int(current_partition_index), ' basins found' END_PROVIDER SUBST [ X ] elf;; eplf;; density;; END_TEMPLATE BEGIN_PROVIDER [ real, current_partition_index ] implicit none BEGIN_DOC ! Current number for the space partitioning END_DOC current_partition_index = 0. END_PROVIDER recursive subroutine find_attribution(ix,iy,iz,g,a,depth) implicit none integer, intent(in) :: ix, iy, iz real, intent(in) :: g(grid_x_num,grid_y_num,grid_z_num) real, intent(inout) :: a(grid_x_num,grid_y_num,grid_z_num) integer, intent(in) :: depth real :: buffer(-1:1,-1:1,-1:1) integer :: i,j,k if (depth == 0) then return endif if (a(ix,iy,iz) /= 0.) then return else if (ix == 1) then call find_attribution(2,iy,iz,g,a,depth-1) a(ix,iy,iz) = a(2,iy,iz) else if (ix == grid_x_num) then call find_attribution(ix-1,iy,iz,g,a,depth-1) a(ix,iy,iz) = a(ix-1,iy,iz) else if (iy == 1) then call find_attribution(ix,2,iz,g,a,depth-1) a(ix,iy,iz) = a(ix,2,iz) else if (iy == grid_y_num) then call find_attribution(ix,iy-1,iz,g,a,depth-1) a(ix,iy,iz) = a(ix,iy-1,iz) else if (iz == 1) then call find_attribution(ix,iy,2,g,a,depth-1) a(ix,iy,iz) = a(ix,iy,2) else if (iz == grid_z_num) then call find_attribution(ix,iy,iz-1,g,a,depth-1) a(ix,iy,iz) = a(ix,iy,iz-1) else ! Copy the environment of the current point do i=-1,1 do j=-1,1 buffer(-1,j,i) = g(ix-1, iy+j, iz+i) buffer( 0,j,i) = g(ix , iy+j, iz+i) buffer( 1,j,i) = g(ix+1, iy+j, iz+i) enddo enddo ! Find the index of the maximum value real :: m integer :: im(3) m = buffer(0,0,0) real :: average, variance im(1) = 0 im(2) = 0 im(3) = 0 average = 0. variance = 0. do i=-1,1 do j=-1,1 do k=-1,1 if (buffer(k,j,i) > m) then m = buffer(k,j,i) im(1) = k im(2) = j im(3) = i endif average += buffer(k,j,i) variance += buffer(k,j,i)**2 enddo enddo enddo variance = variance / 27. average = average / 27. variance = (variance - average**2) ! Avoid strict equality between points buffer(0,0,0) += sqrt(variance/26.) if ( (im(1) == 0).and.(im(2) == 0).and.(im(3) == 0) ) then ! If the maximum is at the center, a new maximum is found if (a(ix,iy,iz) == 0.) then current_partition_index += 1. a(ix,iy,iz) = current_partition_index real :: x,y,z call int_to_cart(ix,iy,iz,x,y,z) print '(3(2X,F10.7))', x,y,z endif else ! Otherwise, get the partition index of the max value call find_attribution(ix+im(1), iy+im(2), iz+im(3),g,a,depth-1) a(ix,iy,iz) = a(ix+im(1), iy+im(2), iz+im(3)) endif endif end subroutine int_to_cart(ix,iy,iz,x,y,z) implicit none integer, intent(in) :: ix, iy, iz real, intent(out) :: x, y, z x = grid_origin(1)+(ix-1)*grid_step(1) y = grid_origin(2)+(iy-1)*grid_step(2) z = grid_origin(3)+(iz-1)*grid_step(3) end