(git:495eafe)
Loading...
Searching...
No Matches
ps_wavelet_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Definition and initialisation of the ps_wavelet data type.
10!> \history 01.2014 Renamed from ps_wavelet_types to disentangle dependencies (Ole Schuett)
11!> \author Florian Schiffmann (09.2007,fschiff)
12! **************************************************************************************************
13MODULE ps_wavelet_methods
14
15 USE bibliography, ONLY: genovese2006,&
16 genovese2007,&
17 cite_reference
18 USE kinds, ONLY: dp
19 USE ps_wavelet_kernel, ONLY: createkernel
20 USE ps_wavelet_types, ONLY: wavelet0d,&
21 wavelet2d,&
22 ps_wavelet_release,&
23 ps_wavelet_type
24 USE ps_wavelet_util, ONLY: f_fft_dimensions,&
25 psolver,&
26 p_fft_dimensions,&
27 s_fft_dimensions
28 USE pw_grid_types, ONLY: pw_grid_type
29 USE pw_poisson_types, ONLY: pw_poisson_parameter_type
30 USE pw_types, ONLY: pw_r3d_rs_type
31 USE util, ONLY: get_limit
32#include "../base/base_uses.f90"
33
34 IMPLICIT NONE
35
36 PRIVATE
37
38 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_methods'
39
40! *** Public data types ***
41
42 PUBLIC :: ps_wavelet_create, &
43 cp2k_distribution_to_z_slices, &
44 z_slices_to_cp2k_distribution, &
45 ps_wavelet_solve
46
47CONTAINS
48
49! **************************************************************************************************
50!> \brief creates the ps_wavelet_type which is needed for the link to
51!> the Poisson Solver of Luigi Genovese
52!> \param poisson_params ...
53!> \param wavelet wavelet to create
54!> \param pw_grid the grid that is used to create the wavelet kernel
55!> \author Flroian Schiffmann
56! **************************************************************************************************
57 SUBROUTINE ps_wavelet_create(poisson_params, wavelet, pw_grid)
58 TYPE(pw_poisson_parameter_type), INTENT(IN) :: poisson_params
59 TYPE(ps_wavelet_type), POINTER :: wavelet
60 TYPE(pw_grid_type), POINTER :: pw_grid
61
62 CHARACTER(len=*), PARAMETER :: routinen = 'ps_wavelet_create'
63
64 INTEGER :: handle
65 REAL(kind=dp) :: hx, hy, hz
66
67 CALL timeset(routinen, handle)
68
69 CALL cite_reference(genovese2006)
70 CALL cite_reference(genovese2007)
71
72 IF (ASSOCIATED(wavelet)) THEN
73 CALL ps_wavelet_release(wavelet)
74 NULLIFY (wavelet)
75 END IF
76
77 ALLOCATE (wavelet)
78
79 NULLIFY (wavelet%karray, wavelet%rho_z_sliced)
80
81 wavelet%geocode = poisson_params%wavelet_geocode
82 wavelet%method = poisson_params%wavelet_method
83 wavelet%special_dimension = poisson_params%wavelet_special_dimension
84 wavelet%itype_scf = poisson_params%wavelet_scf_type
85 wavelet%datacode = "D"
86
87 CALL set_wavelet_axis(wavelet)
88 hx = pw_grid%dr(wavelet%axis(1))
89 hy = pw_grid%dr(wavelet%axis(2))
90 hz = pw_grid%dr(wavelet%axis(3))
91
92 IF (poisson_params%wavelet_method == wavelet0d) THEN
93 IF (hx /= hy) &
94 cpabort("Poisson solver for non cubic cells not yet implemented")
95 IF (hz /= hy) &
96 cpabort("Poisson solver for non cubic cells not yet implemented")
97 END IF
98
99 CALL rs_z_slice_distribution(wavelet, pw_grid)
100
101 CALL timestop(handle)
102 END SUBROUTINE ps_wavelet_create
103
104! **************************************************************************************************
105!> \brief ...
106!> \param wavelet ...
107!> \param pw_grid ...
108! **************************************************************************************************
109 SUBROUTINE rs_z_slice_distribution(wavelet, pw_grid)
110
111 TYPE(ps_wavelet_type), POINTER :: wavelet
112 TYPE(pw_grid_type), POINTER :: pw_grid
113
114 CHARACTER(len=*), PARAMETER :: routinen = 'RS_z_slice_distribution'
115
116 CHARACTER(LEN=1) :: geocode
117 INTEGER :: handle, iproc, m1, m2, m3, md1, md2, &
118 md3, n1, n2, n3, nd1, nd2, nd3, nproc, &
119 nx, ny, nz, z_dim
120 REAL(kind=dp) :: hx, hy, hz
121
122 CALL timeset(routinen, handle)
123 nproc = product(pw_grid%para%group%num_pe_cart)
124 iproc = pw_grid%para%group%mepos
125 geocode = wavelet%geocode
126 CALL get_wavelet_grid(wavelet, pw_grid, nx, ny, nz, hx, hy, hz)
127
128 !calculate Dimensions for the z-distributed density and for the kernel
129
130 IF (geocode == 'P') THEN
131 CALL p_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
132 ELSE IF (geocode == 'S') THEN
133 CALL s_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
134 ELSE IF (geocode == 'F') THEN
135 CALL f_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
136 END IF
137
138 wavelet%PS_grid(1) = md1
139 wavelet%PS_grid(2) = md3
140 wavelet%PS_grid(3) = md2
141 z_dim = md2/nproc
142 !!!!!!!!! indices y and z are interchanged !!!!!!!
143 ALLOCATE (wavelet%rho_z_sliced(md1, md3, z_dim))
144
145 CALL createkernel(geocode, nx, ny, nz, hx, hy, hz, wavelet%itype_scf, iproc, nproc, wavelet%karray, &
146 pw_grid%para%group)
147
148 CALL timestop(handle)
149 END SUBROUTINE rs_z_slice_distribution
150
151! **************************************************************************************************
152!> \brief ...
153!> \param density ...
154!> \param wavelet ...
155!> \param pw_grid ...
156! **************************************************************************************************
157 SUBROUTINE cp2k_distribution_to_z_slices(density, wavelet, pw_grid)
158
159 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
160 TYPE(ps_wavelet_type), POINTER :: wavelet
161 TYPE(pw_grid_type), POINTER :: pw_grid
162
163 CHARACTER(len=*), PARAMETER :: routinen = 'cp2k_distribution_to_z_slices'
164
165 INTEGER :: dest, handle, i, ii, iproc, j, k, l, &
166 local_z_dim, loz, m, m2, md2, nproc, &
167 should_warn
168 INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl, scount, sdispl, tmp
169 INTEGER, DIMENSION(2) :: cart_pos, lox, loy
170 INTEGER, DIMENSION(3) :: lb, ub
171 REAL(kind=dp) :: max_val_low, max_val_up
172 REAL(kind=dp), DIMENSION(:), POINTER :: rbuf, sbuf
173
174 CALL timeset(routinen, handle)
175
176 cpassert(ASSOCIATED(wavelet))
177
178 IF (.NOT. wavelet_axis_is_identity(wavelet)) THEN
179 CALL warn_density_edges(density, wavelet, pw_grid)
180 CALL cp2k_distribution_to_z_slices_permuted(density, wavelet, pw_grid)
181 CALL timestop(handle)
182 RETURN
183 END IF
184
185 nproc = product(pw_grid%para%group%num_pe_cart)
186 iproc = pw_grid%para%group%mepos
187 md2 = wavelet%PS_grid(3)
188 m2 = pw_grid%npts(3)
189 lb(:) = pw_grid%bounds_local(1, :)
190 ub(:) = pw_grid%bounds_local(2, :)
191 local_z_dim = max((md2/nproc), 1)
192
193 ALLOCATE (sbuf(product(pw_grid%npts_local)))
194 ALLOCATE (rbuf(product(wavelet%PS_grid)/nproc))
195 ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), tmp(nproc))
196
197 rbuf = 0.0_dp
198 ii = 1
199 DO k = lb(3), ub(3)
200 DO j = lb(2), ub(2)
201 DO i = lb(1), ub(1)
202 sbuf(ii) = density%array(i, j, k)
203 ii = ii + 1
204 END DO
205 END DO
206 END DO
207
208 should_warn = 0
209 IF (wavelet%geocode == 'S' .OR. wavelet%geocode == 'F') THEN
210 max_val_low = 0._dp
211 max_val_up = 0._dp
212 IF (lb(2) == pw_grid%bounds(1, 2)) max_val_low = maxval(abs(density%array(:, lb(2), :)))
213 IF (ub(2) == pw_grid%bounds(2, 2)) max_val_up = maxval(abs(density%array(:, ub(2), :)))
214 IF (max_val_low >= 0.0001_dp) should_warn = 1
215 IF (max_val_up >= 0.0001_dp) should_warn = 1
216 IF (wavelet%geocode == 'F') THEN
217 max_val_low = 0._dp
218 max_val_up = 0._dp
219 IF (lb(1) == pw_grid%bounds(1, 1)) max_val_low = maxval(abs(density%array(lb(1), :, :)))
220 IF (ub(1) == pw_grid%bounds(2, 1)) max_val_up = maxval(abs(density%array(ub(1), :, :)))
221 IF (max_val_low >= 0.0001_dp) should_warn = 1
222 IF (max_val_up >= 0.0001_dp) should_warn = 1
223 max_val_low = 0._dp
224 max_val_up = 0._dp
225 IF (lb(3) == pw_grid%bounds(1, 3)) max_val_low = maxval(abs(density%array(:, :, lb(3))))
226 IF (ub(3) == pw_grid%bounds(2, 3)) max_val_up = maxval(abs(density%array(:, :, ub(3))))
227 IF (max_val_low >= 0.0001_dp) should_warn = 1
228 IF (max_val_up >= 0.0001_dp) should_warn = 1
229 END IF
230 END IF
231
232 CALL pw_grid%para%group%max(should_warn)
233 IF (should_warn > 0 .AND. iproc == 0) THEN
234 cpwarn("Density non-zero on the edges of the unit cell: wrong results in WAVELET solver")
235 END IF
236 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
237 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
238 cart_pos = [i, j]
239 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
240 IF ((ub(1) >= lb(1)) .AND. (ub(2) >= lb(2))) THEN
241 IF (dest*local_z_dim <= m2) THEN
242 IF ((dest + 1)*local_z_dim <= m2) THEN
243 scount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*local_z_dim)
244 ELSE
245 scount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*mod(m2, local_z_dim))
246 END IF
247 ELSE
248 scount(dest + 1) = 0
249 END IF
250 ELSE
251 scount(dest + 1) = 0
252 END IF
253 lox = get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), i)
254 loy = get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), j)
255 IF ((lox(2) >= lox(1)) .AND. (loy(2) >= loy(1))) THEN
256 IF (iproc*local_z_dim <= m2) THEN
257 IF ((iproc + 1)*local_z_dim <= m2) THEN
258 rcount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*local_z_dim)
259 ELSE
260 rcount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*mod(m2, local_z_dim))
261 END IF
262 ELSE
263 rcount(dest + 1) = 0
264 END IF
265 ELSE
266 rcount(dest + 1) = 0
267 END IF
268
269 END DO
270 END DO
271 sdispl(1) = 0
272 rdispl(1) = 0
273 DO i = 2, nproc
274 sdispl(i) = sdispl(i - 1) + scount(i - 1)
275 rdispl(i) = rdispl(i - 1) + rcount(i - 1)
276 END DO
277 CALL pw_grid%para%group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
278 !!!! and now, how to put the right cubes to the right position!!!!!!
279
280 wavelet%rho_z_sliced = 0.0_dp
281
282 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
283 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
284 cart_pos = [i, j]
285 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
286
287 lox = get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), i)
288 loy = get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), j)
289 IF (iproc*local_z_dim <= m2) THEN
290 IF ((iproc + 1)*local_z_dim <= m2) THEN
291 loz = local_z_dim
292 ELSE
293 loz = mod(m2, local_z_dim)
294 END IF
295 ii = 1
296 DO k = 1, loz
297 DO l = loy(1), loy(2)
298 DO m = lox(1), lox(2)
299 wavelet%rho_z_sliced(m, l, k) = rbuf(ii + rdispl(dest + 1))
300 ii = ii + 1
301 END DO
302 END DO
303 END DO
304 END IF
305 END DO
306 END DO
307
308 DEALLOCATE (sbuf, rbuf, scount, sdispl, rcount, rdispl, tmp)
309
310 CALL timestop(handle)
311
312 END SUBROUTINE cp2k_distribution_to_z_slices
313
314! **************************************************************************************************
315!> \brief ...
316!> \param density ...
317!> \param wavelet ...
318!> \param pw_grid ...
319! **************************************************************************************************
320 SUBROUTINE z_slices_to_cp2k_distribution(density, wavelet, pw_grid)
321
322 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
323 TYPE(ps_wavelet_type), POINTER :: wavelet
324 TYPE(pw_grid_type), POINTER :: pw_grid
325
326 INTEGER :: dest, i, ii, iproc, j, k, l, &
327 local_z_dim, loz, m, m2, md2, nproc
328 INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl, scount, sdispl, tmp
329 INTEGER, DIMENSION(2) :: cart_pos, lox, loy, min_x, min_y
330 INTEGER, DIMENSION(3) :: lb, ub
331 REAL(kind=dp), DIMENSION(:), POINTER :: rbuf, sbuf
332
333 cpassert(ASSOCIATED(wavelet))
334
335 IF (.NOT. wavelet_axis_is_identity(wavelet)) THEN
336 CALL z_slices_to_cp2k_distribution_permuted(density, wavelet, pw_grid)
337 RETURN
338 END IF
339
340 nproc = product(pw_grid%para%group%num_pe_cart)
341 iproc = pw_grid%para%group%mepos
342 md2 = wavelet%PS_grid(3)
343 m2 = pw_grid%npts(3)
344
345 lb(:) = pw_grid%bounds_local(1, :)
346 ub(:) = pw_grid%bounds_local(2, :)
347
348 local_z_dim = max((md2/nproc), 1)
349
350 ALLOCATE (rbuf(product(pw_grid%npts_local)))
351 ALLOCATE (sbuf(product(wavelet%PS_grid)/nproc))
352 ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), tmp(nproc))
353 scount = 0
354 rcount = 0
355 rbuf = 0.0_dp
356 ii = 1
357 IF (iproc*local_z_dim <= m2) THEN
358 IF ((iproc + 1)*local_z_dim <= m2) THEN
359 loz = local_z_dim
360 ELSE
361 loz = mod(m2, local_z_dim)
362 END IF
363 ELSE
364 loz = 0
365 END IF
366
367 min_x = get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), 0)
368 min_y = get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), 0)
369 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
370 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
371 cart_pos = [i, j]
372 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
373 IF ((ub(1) >= lb(1)) .AND. (ub(2) >= lb(2))) THEN
374 IF (dest*local_z_dim <= m2) THEN
375 IF ((dest + 1)*local_z_dim <= m2) THEN
376 rcount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*local_z_dim)
377 ELSE
378 rcount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*mod(m2, local_z_dim))
379 END IF
380 ELSE
381 rcount(dest + 1) = 0
382 END IF
383 ELSE
384 rcount(dest + 1) = 0
385 END IF
386 lox = get_limit(pw_grid%npts(1), pw_grid%para%group%num_pe_cart(1), i)
387 loy = get_limit(pw_grid%npts(2), pw_grid%para%group%num_pe_cart(2), j)
388 IF ((lox(2) >= lox(1)) .AND. (loy(2) >= loy(1))) THEN
389 scount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*loz)
390 DO k = lox(1) - min_x(1) + 1, lox(2) - min_x(1) + 1
391 DO l = loy(1) - min_y(1) + 1, loy(2) - min_y(1) + 1
392 DO m = 1, loz
393 sbuf(ii) = wavelet%rho_z_sliced(k, l, m)
394 ii = ii + 1
395 END DO
396 END DO
397 END DO
398 ELSE
399 scount(dest + 1) = 0
400 END IF
401 END DO
402 END DO
403 sdispl(1) = 0
404 rdispl(1) = 0
405 DO i = 2, nproc
406 sdispl(i) = sdispl(i - 1) + scount(i - 1)
407 rdispl(i) = rdispl(i - 1) + rcount(i - 1)
408 END DO
409 CALL pw_grid%para%group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
410
411 !!!! and now, how to put the right cubes to the right position!!!!!!
412
413 DO i = 0, pw_grid%para%group%num_pe_cart(1) - 1
414 DO j = 0, pw_grid%para%group%num_pe_cart(2) - 1
415 cart_pos = [i, j]
416 CALL pw_grid%para%group%rank_cart(cart_pos, dest)
417 IF (dest*local_z_dim <= m2) THEN
418 IF ((dest + 1)*local_z_dim <= m2) THEN
419 loz = local_z_dim
420 ELSE
421 loz = mod(m2, local_z_dim)
422 END IF
423 ii = 1
424 IF (lb(3) + (dest*local_z_dim) <= ub(3)) THEN
425 DO m = lb(1), ub(1)
426 DO l = lb(2), ub(2)
427 DO k = lb(3) + (dest*local_z_dim), lb(3) + (dest*local_z_dim) + loz - 1
428 density%array(m, l, k) = rbuf(ii + rdispl(dest + 1))
429 ii = ii + 1
430 END DO
431 END DO
432 END DO
433 END IF
434 END IF
435 END DO
436 END DO
437 DEALLOCATE (sbuf, rbuf, scount, sdispl, rcount, rdispl, tmp)
438
439 END SUBROUTINE z_slices_to_cp2k_distribution
440
441! **************************************************************************************************
442!> \brief Set the internal Wavelet solver axes for the requested boundary conditions.
443!> \param wavelet ...
444! **************************************************************************************************
445 SUBROUTINE set_wavelet_axis(wavelet)
446
447 TYPE(ps_wavelet_type), POINTER :: wavelet
448
449 wavelet%axis = [1, 2, 3]
450
451 IF (wavelet%method == wavelet2d) THEN
452 SELECT CASE (wavelet%special_dimension)
453 CASE (1)
454 wavelet%axis = [2, 1, 3]
455 CASE (2)
456 wavelet%axis = [1, 2, 3]
457 CASE (3)
458 wavelet%axis = [1, 3, 2]
459 CASE DEFAULT
460 cpabort("Invalid isolated dimension for WAVELET 2D")
461 END SELECT
462 END IF
463
464 END SUBROUTINE set_wavelet_axis
465
466! **************************************************************************************************
467!> \brief Return grid sizes and spacings in the internal Wavelet solver axes.
468!> \param wavelet ...
469!> \param pw_grid ...
470!> \param nx ...
471!> \param ny ...
472!> \param nz ...
473!> \param hx ...
474!> \param hy ...
475!> \param hz ...
476! **************************************************************************************************
477 SUBROUTINE get_wavelet_grid(wavelet, pw_grid, nx, ny, nz, hx, hy, hz)
478
479 TYPE(ps_wavelet_type), POINTER :: wavelet
480 TYPE(pw_grid_type), POINTER :: pw_grid
481 INTEGER, INTENT(OUT) :: nx, ny, nz
482 REAL(kind=dp), INTENT(OUT) :: hx, hy, hz
483
484 nx = pw_grid%npts(wavelet%axis(1))
485 ny = pw_grid%npts(wavelet%axis(2))
486 nz = pw_grid%npts(wavelet%axis(3))
487 hx = pw_grid%dr(wavelet%axis(1))
488 hy = pw_grid%dr(wavelet%axis(2))
489 hz = pw_grid%dr(wavelet%axis(3))
490
491 END SUBROUTINE get_wavelet_grid
492
493! **************************************************************************************************
494!> \brief Return whether CP2K and internal Wavelet solver axes are identical.
495!> \param wavelet ...
496!> \return ...
497! **************************************************************************************************
498 FUNCTION wavelet_axis_is_identity(wavelet) RESULT(is_identity)
499
500 TYPE(ps_wavelet_type), POINTER :: wavelet
501 LOGICAL :: is_identity
502
503 is_identity = all(wavelet%axis == [1, 2, 3])
504
505 END FUNCTION wavelet_axis_is_identity
506
507! **************************************************************************************************
508!> \brief Warn if the density is non-zero at isolated Wavelet cell edges.
509!> \param density ...
510!> \param wavelet ...
511!> \param pw_grid ...
512! **************************************************************************************************
513 SUBROUTINE warn_density_edges(density, wavelet, pw_grid)
514
515 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
516 TYPE(ps_wavelet_type), POINTER :: wavelet
517 TYPE(pw_grid_type), POINTER :: pw_grid
518
519 INTEGER :: idir, iproc, should_warn
520
521 should_warn = 0
522 iproc = pw_grid%para%group%mepos
523
524 IF (wavelet%geocode == 'S') THEN
525 CALL update_edge_warning(density, pw_grid, wavelet%special_dimension, should_warn)
526 ELSE IF (wavelet%geocode == 'F') THEN
527 DO idir = 1, 3
528 CALL update_edge_warning(density, pw_grid, idir, should_warn)
529 END DO
530 END IF
531
532 CALL pw_grid%para%group%max(should_warn)
533 IF (should_warn > 0 .AND. iproc == 0) THEN
534 cpwarn("Density non-zero on the edges of the unit cell: wrong results in WAVELET solver")
535 END IF
536
537 END SUBROUTINE warn_density_edges
538
539! **************************************************************************************************
540!> \brief Update the density edge warning for one real-space direction.
541!> \param density ...
542!> \param pw_grid ...
543!> \param direction ...
544!> \param should_warn ...
545! **************************************************************************************************
546 SUBROUTINE update_edge_warning(density, pw_grid, direction, should_warn)
547
548 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
549 TYPE(pw_grid_type), POINTER :: pw_grid
550 INTEGER, INTENT(IN) :: direction
551 INTEGER, INTENT(INOUT) :: should_warn
552
553 INTEGER, DIMENSION(3) :: lb, ub
554 REAL(kind=dp) :: max_val_low, max_val_up
555
556 lb(:) = pw_grid%bounds_local(1, :)
557 ub(:) = pw_grid%bounds_local(2, :)
558 IF (.NOT. all(ub >= lb)) RETURN
559
560 max_val_low = 0._dp
561 max_val_up = 0._dp
562 SELECT CASE (direction)
563 CASE (1)
564 IF (lb(1) == pw_grid%bounds(1, 1)) max_val_low = maxval(abs(density%array(lb(1), :, :)))
565 IF (ub(1) == pw_grid%bounds(2, 1)) max_val_up = maxval(abs(density%array(ub(1), :, :)))
566 CASE (2)
567 IF (lb(2) == pw_grid%bounds(1, 2)) max_val_low = maxval(abs(density%array(:, lb(2), :)))
568 IF (ub(2) == pw_grid%bounds(2, 2)) max_val_up = maxval(abs(density%array(:, ub(2), :)))
569 CASE (3)
570 IF (lb(3) == pw_grid%bounds(1, 3)) max_val_low = maxval(abs(density%array(:, :, lb(3))))
571 IF (ub(3) == pw_grid%bounds(2, 3)) max_val_up = maxval(abs(density%array(:, :, ub(3))))
572 CASE DEFAULT
573 cpabort("Invalid WAVELET isolated dimension")
574 END SELECT
575
576 IF (max_val_low >= 0.0001_dp) should_warn = 1
577 IF (max_val_up >= 0.0001_dp) should_warn = 1
578
579 END SUBROUTINE update_edge_warning
580
581! **************************************************************************************************
582!> \brief Convert counts into zero-based displacements for all-to-all communication.
583!> \param counts ...
584!> \param displacements ...
585! **************************************************************************************************
586 SUBROUTINE set_displacements(counts, displacements)
587
588 INTEGER, DIMENSION(:), INTENT(IN) :: counts
589 INTEGER, DIMENSION(:), INTENT(OUT) :: displacements
590
591 INTEGER :: i
592
593 displacements(1) = 0
594 DO i = 2, SIZE(counts)
595 displacements(i) = displacements(i - 1) + counts(i - 1)
596 END DO
597
598 END SUBROUTINE set_displacements
599
600! **************************************************************************************************
601!> \brief Return the distributed grid owner of a compact one-based grid index.
602!> \param index ...
603!> \param npts ...
604!> \param nparts ...
605!> \return ...
606! **************************************************************************************************
607 FUNCTION grid_owner(index, npts, nparts) RESULT(owner)
608
609 INTEGER, INTENT(IN) :: index, npts, nparts
610 INTEGER :: owner
611
612 INTEGER :: ipart
613 INTEGER, DIMENSION(2) :: limits
614
615 DO ipart = 0, nparts - 1
616 limits = get_limit(npts, nparts, ipart)
617 IF (index >= limits(1) .AND. index <= limits(2)) THEN
618 owner = ipart
619 RETURN
620 END IF
621 END DO
622 cpabort("Grid index outside distributed bounds")
623
624 END FUNCTION grid_owner
625
626! **************************************************************************************************
627!> \brief Return the Wavelet z-slice owner of a compact one-based grid index.
628!> \param index ...
629!> \param local_z_dim ...
630!> \param nproc ...
631!> \return ...
632! **************************************************************************************************
633 FUNCTION z_slice_owner(index, local_z_dim, nproc) RESULT(owner)
634
635 INTEGER, INTENT(IN) :: index, local_z_dim, nproc
636 INTEGER :: owner
637
638 owner = min((index - 1)/local_z_dim, nproc - 1)
639
640 END FUNCTION z_slice_owner
641
642! **************************************************************************************************
643!> \brief Return the CP2K real-space rank owning a real-space x/y grid point.
644!> \param ix ...
645!> \param iy ...
646!> \param pw_grid ...
647!> \return ...
648! **************************************************************************************************
649 FUNCTION cp2k_rank_owner(ix, iy, pw_grid) RESULT(rank)
650
651 INTEGER, INTENT(IN) :: ix, iy
652 TYPE(pw_grid_type), POINTER :: pw_grid
653 INTEGER :: rank
654
655 INTEGER, DIMENSION(2) :: cart_pos
656
657 cart_pos(1) = grid_owner(ix - pw_grid%bounds(1, 1) + 1, pw_grid%npts(1), &
658 pw_grid%para%group%num_pe_cart(1))
659 cart_pos(2) = grid_owner(iy - pw_grid%bounds(1, 2) + 1, pw_grid%npts(2), &
660 pw_grid%para%group%num_pe_cart(2))
661 CALL pw_grid%para%group%rank_cart(cart_pos, rank)
662
663 END FUNCTION cp2k_rank_owner
664
665! **************************************************************************************************
666!> \brief Transfer a CP2K real-space grid into a permuted Wavelet z-sliced layout.
667!> \param density ...
668!> \param wavelet ...
669!> \param pw_grid ...
670! **************************************************************************************************
671 SUBROUTINE cp2k_distribution_to_z_slices_permuted(density, wavelet, pw_grid)
672
673 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
674 TYPE(ps_wavelet_type), POINTER :: wavelet
675 TYPE(pw_grid_type), POINTER :: pw_grid
676
677 INTEGER :: dest, i, idir, ii, j, k, local_z_dim, &
678 md2, nproc, nrecv, nsend
679 INTEGER, ALLOCATABLE, DIMENSION(:) :: coord_rbuf, coord_rcount, coord_rdispl, coord_sbuf, &
680 coord_scount, coord_sdispl, rcount, rdispl, scount, sdispl, send_pos
681 INTEGER, DIMENSION(3) :: lb, q, u, ub
682 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rbuf, sbuf
683
684 nproc = product(pw_grid%para%group%num_pe_cart)
685 md2 = wavelet%PS_grid(3)
686 local_z_dim = max(md2/nproc, 1)
687 lb(:) = pw_grid%bounds_local(1, :)
688 ub(:) = pw_grid%bounds_local(2, :)
689
690 ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), send_pos(nproc))
691 scount = 0
692 DO k = lb(3), ub(3)
693 DO j = lb(2), ub(2)
694 DO i = lb(1), ub(1)
695 u = [i, j, k]
696 DO idir = 1, 3
697 q(idir) = u(wavelet%axis(idir)) - pw_grid%bounds(1, wavelet%axis(idir)) + 1
698 END DO
699 dest = z_slice_owner(q(3), local_z_dim, nproc)
700 scount(dest + 1) = scount(dest + 1) + 1
701 END DO
702 END DO
703 END DO
704
705 CALL pw_grid%para%group%alltoall(scount, rcount, 1)
706 CALL set_displacements(scount, sdispl)
707 CALL set_displacements(rcount, rdispl)
708 nsend = sum(scount)
709 nrecv = sum(rcount)
710
711 ALLOCATE (sbuf(max(nsend, 1)), rbuf(max(nrecv, 1)))
712 ALLOCATE (coord_sbuf(max(3*nsend, 1)), coord_rbuf(max(3*nrecv, 1)))
713 ALLOCATE (coord_scount(nproc), coord_sdispl(nproc), coord_rcount(nproc), coord_rdispl(nproc))
714 coord_scount(:) = 3*scount(:)
715 coord_rcount(:) = 3*rcount(:)
716 coord_sdispl(:) = 3*sdispl(:)
717 coord_rdispl(:) = 3*rdispl(:)
718 send_pos(:) = sdispl(:) + 1
719
720 DO k = lb(3), ub(3)
721 DO j = lb(2), ub(2)
722 DO i = lb(1), ub(1)
723 u = [i, j, k]
724 DO idir = 1, 3
725 q(idir) = u(wavelet%axis(idir)) - pw_grid%bounds(1, wavelet%axis(idir)) + 1
726 END DO
727 dest = z_slice_owner(q(3), local_z_dim, nproc)
728 ii = send_pos(dest + 1)
729 sbuf(ii) = density%array(i, j, k)
730 coord_sbuf(3*ii - 2) = q(1)
731 coord_sbuf(3*ii - 1) = q(2)
732 coord_sbuf(3*ii) = q(3)
733 send_pos(dest + 1) = ii + 1
734 END DO
735 END DO
736 END DO
737
738 CALL pw_grid%para%group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
739 CALL pw_grid%para%group%alltoall(coord_sbuf, coord_scount, coord_sdispl, &
740 coord_rbuf, coord_rcount, coord_rdispl)
741
742 wavelet%rho_z_sliced = 0.0_dp
743 DO ii = 1, nrecv
744 q(1) = coord_rbuf(3*ii - 2)
745 q(2) = coord_rbuf(3*ii - 1)
746 q(3) = coord_rbuf(3*ii)
747 wavelet%rho_z_sliced(q(1), q(2), q(3) - pw_grid%para%group%mepos*local_z_dim) = rbuf(ii)
748 END DO
749
750 DEALLOCATE (sbuf, rbuf, coord_sbuf, coord_rbuf, coord_scount, coord_sdispl, coord_rcount, &
751 coord_rdispl, scount, sdispl, rcount, rdispl, send_pos)
752
753 END SUBROUTINE cp2k_distribution_to_z_slices_permuted
754
755! **************************************************************************************************
756!> \brief Transfer a permuted Wavelet z-sliced layout back to a CP2K real-space grid.
757!> \param density ...
758!> \param wavelet ...
759!> \param pw_grid ...
760! **************************************************************************************************
761 SUBROUTINE z_slices_to_cp2k_distribution_permuted(density, wavelet, pw_grid)
762
763 TYPE(pw_r3d_rs_type), INTENT(IN) :: density
764 TYPE(ps_wavelet_type), POINTER :: wavelet
765 TYPE(pw_grid_type), POINTER :: pw_grid
766
767 INTEGER :: dest, i, idir, ii, j, k, local_z, &
768 local_z_dim, md2, n1, n2, n3, nproc, &
769 nrecv, nsend, z_end, z_start
770 INTEGER, ALLOCATABLE, DIMENSION(:) :: coord_rbuf, coord_rcount, coord_rdispl, coord_sbuf, &
771 coord_scount, coord_sdispl, rcount, rdispl, scount, sdispl, send_pos
772 INTEGER, DIMENSION(3) :: lb, q, u, ub
773 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rbuf, sbuf
774
775 nproc = product(pw_grid%para%group%num_pe_cart)
776 md2 = wavelet%PS_grid(3)
777 local_z_dim = max(md2/nproc, 1)
778 n1 = pw_grid%npts(wavelet%axis(1))
779 n2 = pw_grid%npts(wavelet%axis(2))
780 n3 = pw_grid%npts(wavelet%axis(3))
781 z_start = pw_grid%para%group%mepos*local_z_dim + 1
782 z_end = min((pw_grid%para%group%mepos + 1)*local_z_dim, n3)
783 lb(:) = pw_grid%bounds_local(1, :)
784 ub(:) = pw_grid%bounds_local(2, :)
785
786 ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), send_pos(nproc))
787 scount = 0
788 DO k = z_start, z_end
789 DO j = 1, n2
790 DO i = 1, n1
791 q = [i, j, k]
792 DO idir = 1, 3
793 u(wavelet%axis(idir)) = q(idir) + pw_grid%bounds(1, wavelet%axis(idir)) - 1
794 END DO
795 dest = cp2k_rank_owner(u(1), u(2), pw_grid)
796 scount(dest + 1) = scount(dest + 1) + 1
797 END DO
798 END DO
799 END DO
800
801 CALL pw_grid%para%group%alltoall(scount, rcount, 1)
802 CALL set_displacements(scount, sdispl)
803 CALL set_displacements(rcount, rdispl)
804 nsend = sum(scount)
805 nrecv = sum(rcount)
806
807 ALLOCATE (sbuf(max(nsend, 1)), rbuf(max(nrecv, 1)))
808 ALLOCATE (coord_sbuf(max(3*nsend, 1)), coord_rbuf(max(3*nrecv, 1)))
809 ALLOCATE (coord_scount(nproc), coord_sdispl(nproc), coord_rcount(nproc), coord_rdispl(nproc))
810 coord_scount(:) = 3*scount(:)
811 coord_rcount(:) = 3*rcount(:)
812 coord_sdispl(:) = 3*sdispl(:)
813 coord_rdispl(:) = 3*rdispl(:)
814 send_pos(:) = sdispl(:) + 1
815
816 DO k = z_start, z_end
817 DO j = 1, n2
818 DO i = 1, n1
819 q = [i, j, k]
820 DO idir = 1, 3
821 u(wavelet%axis(idir)) = q(idir) + pw_grid%bounds(1, wavelet%axis(idir)) - 1
822 END DO
823 dest = cp2k_rank_owner(u(1), u(2), pw_grid)
824 ii = send_pos(dest + 1)
825 local_z = k - pw_grid%para%group%mepos*local_z_dim
826 sbuf(ii) = wavelet%rho_z_sliced(i, j, local_z)
827 coord_sbuf(3*ii - 2) = u(1)
828 coord_sbuf(3*ii - 1) = u(2)
829 coord_sbuf(3*ii) = u(3)
830 send_pos(dest + 1) = ii + 1
831 END DO
832 END DO
833 END DO
834
835 CALL pw_grid%para%group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
836 CALL pw_grid%para%group%alltoall(coord_sbuf, coord_scount, coord_sdispl, &
837 coord_rbuf, coord_rcount, coord_rdispl)
838
839 DO ii = 1, nrecv
840 u(1) = coord_rbuf(3*ii - 2)
841 u(2) = coord_rbuf(3*ii - 1)
842 u(3) = coord_rbuf(3*ii)
843 IF (u(1) >= lb(1) .AND. u(1) <= ub(1) .AND. u(2) >= lb(2) .AND. u(2) <= ub(2) .AND. &
844 u(3) >= lb(3) .AND. u(3) <= ub(3)) THEN
845 density%array(u(1), u(2), u(3)) = rbuf(ii)
846 END IF
847 END DO
848
849 DEALLOCATE (sbuf, rbuf, coord_sbuf, coord_rbuf, coord_scount, coord_sdispl, coord_rcount, &
850 coord_rdispl, scount, sdispl, rcount, rdispl, send_pos)
851
852 END SUBROUTINE z_slices_to_cp2k_distribution_permuted
853
854! **************************************************************************************************
855!> \brief ...
856!> \param wavelet ...
857!> \param pw_grid ...
858! **************************************************************************************************
859 SUBROUTINE ps_wavelet_solve(wavelet, pw_grid)
860
861 TYPE(ps_wavelet_type), POINTER :: wavelet
862 TYPE(pw_grid_type), POINTER :: pw_grid
863
864 CHARACTER(len=*), PARAMETER :: routinen = 'ps_wavelet_solve'
865
866 CHARACTER(LEN=1) :: geocode
867 INTEGER :: handle, iproc, nproc, nx, ny, nz
868 REAL(kind=dp) :: hx, hy, hz
869
870 CALL timeset(routinen, handle)
871 nproc = product(pw_grid%para%group%num_pe_cart)
872 iproc = pw_grid%para%group%mepos
873 geocode = wavelet%geocode
874 CALL get_wavelet_grid(wavelet, pw_grid, nx, ny, nz, hx, hy, hz)
875
876 CALL psolver(geocode, iproc, nproc, nx, ny, nz, hx, hy, hz, &
877 wavelet%rho_z_sliced, wavelet%karray, pw_grid)
878 CALL timestop(handle)
879 END SUBROUTINE ps_wavelet_solve
880
881END MODULE ps_wavelet_methods
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::genovese2006
integer, save, public genovese2006
Definition bibliography.F:36
bibliography::genovese2007
integer, save, public genovese2007
Definition bibliography.F:36
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
ps_wavelet_kernel
Creates the wavelet kernel for the wavelet based poisson solver.
Definition ps_wavelet_kernel.F:12
ps_wavelet_kernel::createkernel
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 ps_wavelet_kernel.F:80
ps_wavelet_methods
Definition and initialisation of the ps_wavelet data type. \history 01.2014 Renamed from ps_wavelet_t...
Definition ps_wavelet_methods.F:13
ps_wavelet_methods::ps_wavelet_create
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
Definition ps_wavelet_methods.F:58
ps_wavelet_methods::ps_wavelet_solve
subroutine, public ps_wavelet_solve(wavelet, pw_grid)
...
Definition ps_wavelet_methods.F:860
ps_wavelet_methods::z_slices_to_cp2k_distribution
subroutine, public z_slices_to_cp2k_distribution(density, wavelet, pw_grid)
...
Definition ps_wavelet_methods.F:321
ps_wavelet_methods::cp2k_distribution_to_z_slices
subroutine, public cp2k_distribution_to_z_slices(density, wavelet, pw_grid)
...
Definition ps_wavelet_methods.F:158
ps_wavelet_types
Definition and initialisation of the ps_wavelet data type.
Definition ps_wavelet_types.F:12
ps_wavelet_types::wavelet0d
integer, parameter, public wavelet0d
Definition ps_wavelet_types.F:23
ps_wavelet_types::ps_wavelet_release
subroutine, public ps_wavelet_release(wavelet)
...
Definition ps_wavelet_types.F:54
ps_wavelet_types::wavelet2d
integer, parameter, public wavelet2d
Definition ps_wavelet_types.F:23
ps_wavelet_util
Performs a wavelet based solution of the Poisson equation.
Definition ps_wavelet_util.F:12
ps_wavelet_util::p_fft_dimensions
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...
Definition ps_wavelet_util.F:238
ps_wavelet_util::psolver
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$,...
Definition ps_wavelet_util.F:83
ps_wavelet_util::s_fft_dimensions
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...
Definition ps_wavelet_util.F:334
ps_wavelet_util::f_fft_dimensions
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.
Definition ps_wavelet_util.F:441
pw_grid_types
Definition pw_grid_types.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_types
Definition pw_types.F:24
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
ps_wavelet_types::ps_wavelet_type
Definition ps_wavelet_types.F:36
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:52
pw_poisson_types::pw_poisson_parameter_type
parameters for the poisson solver independet of input_section
Definition pw_poisson_types.F:97
pw_types::pw_r3d_rs_type
Definition pw_types.F:62