31#include "../base/base_uses.f90"
37 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ps_wavelet_methods'
61 CHARACTER(len=*),
PARAMETER :: routinen =
'ps_wavelet_create'
63 INTEGER :: handle, iproc, nproc, nx, ny, nz
64 REAL(kind=
dp) :: hx, hy, hz
66 CALL timeset(routinen, handle)
71 IF (
ASSOCIATED(wavelet))
THEN
86 nproc = product(pw_grid%para%group%num_pe_cart)
88 iproc = pw_grid%para%group%mepos
90 NULLIFY (wavelet%karray, wavelet%rho_z_sliced)
92 wavelet%geocode = poisson_params%wavelet_geocode
93 wavelet%method = poisson_params%wavelet_method
94 wavelet%special_dimension = poisson_params%wavelet_special_dimension
95 wavelet%itype_scf = poisson_params%wavelet_scf_type
96 wavelet%datacode =
"D"
98 IF (poisson_params%wavelet_method ==
wavelet0d)
THEN
100 cpabort(
"Poisson solver for non cubic cells not yet implemented")
102 cpabort(
"Poisson solver for non cubic cells not yet implemented")
105 CALL rs_z_slice_distribution(wavelet, pw_grid)
107 CALL timestop(handle)
115 SUBROUTINE rs_z_slice_distribution(wavelet, pw_grid)
120 CHARACTER(len=*),
PARAMETER :: routinen =
'RS_z_slice_distribution'
122 CHARACTER(LEN=1) :: geocode
123 INTEGER :: handle, iproc, m1, m2, m3, md1, md2, &
124 md3, n1, n2, n3, nd1, nd2, nd3, nproc, &
126 REAL(kind=
dp) :: hx, hy, hz
128 CALL timeset(routinen, handle)
129 nproc = product(pw_grid%para%group%num_pe_cart)
130 iproc = pw_grid%para%group%mepos
131 geocode = wavelet%geocode
141 IF (geocode ==
'P')
THEN
142 CALL p_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
143 ELSE IF (geocode ==
'S')
THEN
144 CALL s_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
145 ELSE IF (geocode ==
'F')
THEN
146 CALL f_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
149 wavelet%PS_grid(1) = md1
150 wavelet%PS_grid(2) = md3
151 wavelet%PS_grid(3) = md2
154 ALLOCATE (wavelet%rho_z_sliced(md1, md3, z_dim))
156 CALL createkernel(geocode, nx, ny, nz, hx, hy, hz, wavelet%itype_scf, iproc, nproc, wavelet%karray, &
159 CALL timestop(handle)
160 END SUBROUTINE rs_z_slice_distribution
174 CHARACTER(len=*),
PARAMETER :: routinen =
'cp2k_distribution_to_z_slices'
176 INTEGER :: dest, handle, i, ii, iproc, j, k, l, &
177 local_z_dim, loz, m, m2, md2, nproc, &
179 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: rcount, rdispl, scount, sdispl, tmp
180 INTEGER,
DIMENSION(2) :: cart_pos, lox, loy
181 INTEGER,
DIMENSION(3) :: lb, ub
182 REAL(kind=
dp) :: max_val_low, max_val_up
183 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rbuf, sbuf
185 CALL timeset(routinen, handle)
187 cpassert(
ASSOCIATED(wavelet))
189 nproc = product(pw_grid%para%group%num_pe_cart)
190 iproc = pw_grid%para%group%mepos
191 md2 = wavelet%PS_grid(3)
193 lb(:) = pw_grid%bounds_local(1, :)
194 ub(:) = pw_grid%bounds_local(2, :)
195 local_z_dim = max((md2/nproc), 1)
197 ALLOCATE (sbuf(product(pw_grid%npts_local)))
198 ALLOCATE (rbuf(product(wavelet%PS_grid)/nproc))
199 ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), tmp(nproc))
206 sbuf(ii) = density%array(i, j, k)
213 IF (wavelet%geocode ==
'S' .OR. wavelet%geocode ==
'F')
THEN
216 IF (lb(2) == pw_grid%bounds(1, 2)) max_val_low = maxval(abs(density%array(:, lb(2), :)))
217 IF (ub(2) == pw_grid%bounds(2, 2)) max_val_up = maxval(abs(density%array(:, ub(2), :)))
218 IF (max_val_low .GE. 0.0001_dp) should_warn = 1
219 IF (max_val_up .GE. 0.0001_dp) should_warn = 1
220 IF (wavelet%geocode ==
'F')
THEN
223 IF (lb(1) == pw_grid%bounds(1, 1)) max_val_low = maxval(abs(density%array(lb(1), :, :)))
224 IF (ub(1) == pw_grid%bounds(2, 1)) max_val_up = maxval(abs(density%array(ub(1), :, :)))
225 IF (max_val_low .GE. 0.0001_dp) should_warn = 1
226 IF (max_val_up .GE. 0.0001_dp) should_warn = 1
229 IF (lb(3) == pw_grid%bounds(1, 3)) max_val_low = maxval(abs(density%array(:, :, lb(3))))
230 IF (ub(3) == pw_grid%bounds(2, 3)) max_val_up = maxval(abs(density%array(:, :, ub(3))))
231 IF (max_val_low .GE. 0.0001_dp) should_warn = 1
232 IF (max_val_up .GE. 0.0001_dp) should_warn = 1
236 CALL pw_grid%para%group%max(should_warn)
237 IF (should_warn > 0 .AND. iproc == 0)
THEN
238 cpwarn(
"Density non-zero on the edges of the unit cell: wrong results in WAVELET solver")
240 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
241 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
243 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
244 IF ((ub(1) .GE. lb(1)) .AND. (ub(2) .GE. lb(2)))
THEN
245 IF (dest*local_z_dim .LE. m2)
THEN
246 IF ((dest + 1)*local_z_dim .LE. m2)
THEN
247 scount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*local_z_dim)
249 scount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*mod(m2, local_z_dim))
257 lox =
get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), i)
258 loy =
get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), j)
259 IF ((lox(2) .GE. lox(1)) .AND. (loy(2) .GE. loy(1)))
THEN
260 IF (iproc*local_z_dim .LE. m2)
THEN
261 IF ((iproc + 1)*local_z_dim .LE. m2)
THEN
262 rcount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*local_z_dim)
264 rcount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*mod(m2, local_z_dim))
278 sdispl(i) = sdispl(i - 1) + scount(i - 1)
279 rdispl(i) = rdispl(i - 1) + rcount(i - 1)
281 CALL pw_grid%para%group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
284 wavelet%rho_z_sliced = 0.0_dp
286 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
287 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
289 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
291 lox =
get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), i)
292 loy =
get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), j)
293 IF (iproc*local_z_dim .LE. m2)
THEN
294 IF ((iproc + 1)*local_z_dim .LE. m2)
THEN
297 loz = mod(m2, local_z_dim)
301 DO l = loy(1), loy(2)
302 DO m = lox(1), lox(2)
303 wavelet%rho_z_sliced(m, l, k) = rbuf(ii + rdispl(dest + 1))
312 DEALLOCATE (sbuf, rbuf, scount, sdispl, rcount, rdispl, tmp)
314 CALL timestop(handle)
330 INTEGER :: dest, i, ii, iproc, j, k, l, &
331 local_z_dim, loz, m, m2, md2, nproc
332 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: rcount, rdispl, scount, sdispl, tmp
333 INTEGER,
DIMENSION(2) :: cart_pos, lox, loy, min_x, min_y
334 INTEGER,
DIMENSION(3) :: lb, ub
335 REAL(kind=
dp),
DIMENSION(:),
POINTER :: rbuf, sbuf
337 cpassert(
ASSOCIATED(wavelet))
339 nproc = product(pw_grid%para%group%num_pe_cart)
340 iproc = pw_grid%para%group%mepos
341 md2 = wavelet%PS_grid(3)
344 lb(:) = pw_grid%bounds_local(1, :)
345 ub(:) = pw_grid%bounds_local(2, :)
347 local_z_dim = max((md2/nproc), 1)
349 ALLOCATE (rbuf(product(pw_grid%npts_local)))
350 ALLOCATE (sbuf(product(wavelet%PS_grid)/nproc))
351 ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), tmp(nproc))
356 IF (iproc*local_z_dim .LE. m2)
THEN
357 IF ((iproc + 1)*local_z_dim .LE. m2)
THEN
360 loz = mod(m2, local_z_dim)
366 min_x =
get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), 0)
367 min_y =
get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), 0)
368 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
369 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
371 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
372 IF ((ub(1) .GE. lb(1)) .AND. (ub(2) .GE. lb(2)))
THEN
373 IF (dest*local_z_dim .LE. m2)
THEN
374 IF ((dest + 1)*local_z_dim .LE. m2)
THEN
375 rcount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*local_z_dim)
377 rcount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*mod(m2, local_z_dim))
385 lox =
get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), i)
386 loy =
get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), j)
387 IF ((lox(2) .GE. lox(1)) .AND. (loy(2) .GE. loy(1)))
THEN
388 scount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*loz)
389 DO k = lox(1) - min_x(1) + 1, lox(2) - min_x(1) + 1
390 DO l = loy(1) - min_y(1) + 1, loy(2) - min_y(1) + 1
392 sbuf(ii) = wavelet%rho_z_sliced(k, l, m)
405 sdispl(i) = sdispl(i - 1) + scount(i - 1)
406 rdispl(i) = rdispl(i - 1) + rcount(i - 1)
408 CALL pw_grid%para%group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
412 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
413 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
415 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
416 IF (dest*local_z_dim .LE. m2)
THEN
417 IF ((dest + 1)*local_z_dim .LE. m2)
THEN
420 loz = mod(m2, local_z_dim)
423 IF (lb(3) + (dest*local_z_dim) .LE. ub(3))
THEN
426 DO k = lb(3) + (dest*local_z_dim), lb(3) + (dest*local_z_dim) + loz - 1
427 density%array(m, l, k) = rbuf(ii + rdispl(dest + 1))
436 DEALLOCATE (sbuf, rbuf, scount, sdispl, rcount, rdispl, tmp)
450 CHARACTER(len=*),
PARAMETER :: routinen =
'ps_wavelet_solve'
452 CHARACTER(LEN=1) :: geocode
453 INTEGER :: handle, iproc, nproc, nx, ny, nz
454 REAL(kind=
dp) :: hx, hy, hz
456 CALL timeset(routinen, handle)
457 nproc = product(pw_grid%para%group%num_pe_cart)
458 iproc = pw_grid%para%group%mepos
459 geocode = wavelet%geocode
467 CALL psolver(geocode, iproc, nproc, nx, ny, nz, hx, hy, hz, &
468 wavelet%rho_z_sliced, wavelet%karray, pw_grid)
469 CALL timestop(handle)
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public genovese2006
integer, save, public genovese2007
Defines the basic variable types.
integer, parameter, public dp
Creates the wavelet kernel for the wavelet based poisson solver.
subroutine, public createkernel(geocode, n01, n02, n03, hx, hy, hz, itype_scf, iproc, nproc, kernel, mpi_group)
Allocate a pointer which corresponds to the zero-padded FFT slice needed for calculating the convolut...
Definition and initialisation of the ps_wavelet data type. \history 01.2014 Renamed from ps_wavelet_t...
subroutine, public ps_wavelet_create(poisson_params, wavelet, pw_grid)
creates the ps_wavelet_type which is needed for the link to the Poisson Solver of Luigi Genovese
subroutine, public ps_wavelet_solve(wavelet, pw_grid)
...
subroutine, public z_slices_to_cp2k_distribution(density, wavelet, pw_grid)
...
subroutine, public cp2k_distribution_to_z_slices(density, wavelet, pw_grid)
...
Definition and initialisation of the ps_wavelet data type.
integer, parameter, public wavelet0d
subroutine, public ps_wavelet_release(wavelet)
...
Performs a wavelet based solution of the Poisson equation.
subroutine, public p_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the convolution for the periodic syste...
subroutine, public psolver(geocode, iproc, nproc, n01, n02, n03, hx, hy, hz, rhopot, karray, pw_grid)
Calculate the Poisson equation $\nabla^2 V(x,y,z)=-4 \pi \rho(x,y,z)$ from a given $\rho$,...
subroutine, public s_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the convolution for the surface system...
subroutine, public f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the zero-padded convolution.
functions related to the poisson solver on regular grids
All kind of helpful little routines.
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
parameters for the poisson solver independet of input_section