LCPQ
/
eplf
mirror of https://gitlab.com/scemama/eplf
10
0
Fork 0
Code Issues Releases Wiki Activity
eplf/src/grid.irp.f

426 lines
9.4 KiB
Fortran
Raw Permalink Blame History

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
print *, 'Attribution : $X'
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)
print *, int(current_partition_index), ' basins found'
endif
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
Reference in New Issue View Git Blame Copy Permalink