32 #include "../base/base_uses.f90"
36 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rs_methods'
44 REAL(dp),
PARAMETER,
PRIVATE :: small_value = 1.0e-8_dp
61 TYPE(pw_r3d_rs_type),
INTENT(IN) :: f
62 TYPE(pw_r3d_rs_type),
DIMENSION(3),
INTENT(INOUT) :: df
63 TYPE(realspace_grid_type),
INTENT(IN) :: rs_grid
65 CHARACTER(LEN=*),
PARAMETER :: routinen =
'derive_fdm_cd3'
67 INTEGER :: handle, i, j, k
68 INTEGER,
DIMENSION(3) :: lb, ub
69 REAL(kind=dp),
DIMENSION(3) :: h
70 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: drdx, drdy, drdz, r
71 TYPE(realspace_grid_desc_type),
POINTER :: rs_desc
72 TYPE(realspace_grid_type),
DIMENSION(3) :: drs_grid
74 CALL timeset(routinen, handle)
77 rs_desc => rs_grid%desc
84 lb(1:3) = rs_grid%lb_real(1:3)
85 ub(1:3) = rs_grid%ub_real(1:3)
92 h(1:3) = 2.0_dp*f%pw_grid%dr(1:3)
99 drdx(i, j, k) = (r(i + 1, j, k) - r(i - 1, j, k))/h(1)
100 drdy(i, j, k) = (r(i, j + 1, k) - r(i, j - 1, k))/h(2)
101 drdz(i, j, k) = (r(i, j, k + 1) - r(i, j, k - 1))/h(3)
113 CALL timestop(handle)
130 TYPE(pw_r3d_rs_type),
INTENT(IN) :: f
131 TYPE(pw_r3d_rs_type),
DIMENSION(3),
INTENT(INOUT) :: df
132 TYPE(realspace_grid_type),
INTENT(IN) :: rs_grid
134 CHARACTER(LEN=*),
PARAMETER :: routinen =
'derive_fdm_cd5'
136 INTEGER :: handle, i, j, k
137 INTEGER,
DIMENSION(3) :: lb, ub
138 REAL(kind=dp),
DIMENSION(3) :: h
139 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: drdx, drdy, drdz, r
140 TYPE(realspace_grid_desc_type),
POINTER :: rs_desc
141 TYPE(realspace_grid_type),
DIMENSION(3) :: drs_grid
143 CALL timeset(routinen, handle)
146 rs_desc => rs_grid%desc
153 lb(1:3) = rs_grid%lb_real(1:3)
154 ub(1:3) = rs_grid%ub_real(1:3)
156 drdx => drs_grid(1)%r
157 drdy => drs_grid(2)%r
158 drdz => drs_grid(3)%r
161 h(1:3) = 12.0_dp*f%pw_grid%dr(1:3)
168 drdx(i, j, k) = (r(i - 2, j, k) - r(i + 2, j, k) + 8.0_dp*(r(i + 1, j, k) - r(i - 1, j, k)))/h(1)
169 drdy(i, j, k) = (r(i, j - 2, k) - r(i, j + 2, k) + 8.0_dp*(r(i, j + 1, k) - r(i, j - 1, k)))/h(2)
170 drdz(i, j, k) = (r(i, j, k - 2) - r(i, j, k + 2) + 8.0_dp*(r(i, j, k + 1) - r(i, j, k - 1)))/h(3)
182 CALL timestop(handle)
199 TYPE(pw_r3d_rs_type),
INTENT(IN) :: f
200 TYPE(pw_r3d_rs_type),
DIMENSION(3),
INTENT(INOUT) :: df
201 TYPE(realspace_grid_type),
INTENT(IN) :: rs_grid
203 CHARACTER(LEN=*),
PARAMETER :: routinen =
'derive_fdm_cd7'
205 INTEGER :: handle, i, j, k
206 INTEGER,
DIMENSION(3) :: lb, ub
207 REAL(kind=dp),
DIMENSION(3) :: h
208 REAL(kind=dp),
DIMENSION(:, :, :),
POINTER :: drdx, drdy, drdz, r
209 TYPE(realspace_grid_desc_type),
POINTER :: rs_desc
210 TYPE(realspace_grid_type),
DIMENSION(3) :: drs_grid
212 CALL timeset(routinen, handle)
215 rs_desc => rs_grid%desc
222 lb(1:3) = rs_grid%lb_real(1:3)
223 ub(1:3) = rs_grid%ub_real(1:3)
225 drdx => drs_grid(1)%r
226 drdy => drs_grid(2)%r
227 drdz => drs_grid(3)%r
230 h(1:3) = 60.0_dp*f%pw_grid%dr(1:3)
237 drdx(i, j, k) = (r(i + 3, j, k) - r(i - 3, j, k) + 9.0_dp*(r(i - 2, j, k) - r(i + 2, j, k)) + &
238 45.0_dp*(r(i + 1, j, k) - r(i - 1, j, k)))/h(1)
239 drdy(i, j, k) = (r(i, j + 3, k) - r(i, j - 3, k) + 9.0_dp*(r(i, j - 2, k) - r(i, j + 2, k)) + &
240 45.0_dp*(r(i, j + 1, k) - r(i, j - 1, k)))/h(2)
241 drdz(i, j, k) = (r(i, j, k + 3) - r(i, j, k - 3) + 9.0_dp*(r(i, j, k - 2) - r(i, j, k + 2)) + &
242 45.0_dp*(r(i, j, k + 1) - r(i, j, k - 1)))/h(3)
254 CALL timestop(handle)
275 TYPE(pw_grid_type),
INTENT(IN) :: pw_grid
276 REAL(dp),
ALLOCATABLE,
DIMENSION(:),
INTENT(OUT) :: x_glbl, y_glbl, z_glbl, x_locl, y_locl, &
279 CHARACTER(LEN=*),
PARAMETER :: routinen =
'setup_grid_axes'
281 INTEGER :: glb1, glb2, glb3, gub1, gub2, gub3, &
282 handle, i, lb1, lb2, lb3, ub1, ub2, ub3
283 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: gindx, gindy, gindz, lindx, lindy, lindz
284 INTEGER,
DIMENSION(2, 3) :: bounds, bounds_local
285 INTEGER,
DIMENSION(3) :: npts, npts_local
286 REAL(dp),
DIMENSION(3) :: dr
288 CALL timeset(routinen, handle)
291 bounds = pw_grid%bounds
292 bounds_local = pw_grid%bounds_local
294 npts_local = pw_grid%npts_local
297 glb1 = bounds(1, 1); gub1 = bounds(2, 1)
298 glb2 = bounds(1, 2); gub2 = bounds(2, 2)
299 glb3 = bounds(1, 3); gub3 = bounds(2, 3)
300 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
301 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
302 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
304 ALLOCATE (lindx(lb1:ub1), lindy(lb2:ub2), lindz(lb3:ub3))
305 ALLOCATE (gindx(glb1:gub1), gindy(glb2:gub2), gindz(glb3:gub3))
306 ALLOCATE (x_locl(lb1:ub1), y_locl(lb2:ub2), z_locl(lb3:ub3))
307 ALLOCATE (x_glbl(glb1:gub1), y_glbl(glb2:gub2), z_glbl(glb3:gub3))
309 gindx(:) = (/(i, i=0, npts(1) - 1)/)
310 gindy(:) = (/(i, i=0, npts(2) - 1)/)
311 gindz(:) = (/(i, i=0, npts(3) - 1)/)
312 lindx(:) = (/(i, i=0, npts_local(1) - 1)/)
313 lindy(:) = (/(i, i=0, npts_local(2) - 1)/)
314 lindz(:) = (/(i, i=0, npts_local(3) - 1)/)
316 x_glbl(:) = gindx*dr(1)
317 y_glbl(:) = gindy*dr(2)
318 z_glbl(:) = gindz*dr(3)
321 IF (lb1 .EQ. ub1)
THEN
324 lindx(:) = lindx(:)*((ub1 - lb1)/(npts_local(1) - 1)) + lb1
326 IF (lb2 .EQ. ub2)
THEN
329 lindy(:) = lindy(:)*((ub2 - lb2)/(npts_local(2) - 1)) + lb2
331 IF (lb3 .EQ. ub3)
THEN
334 lindz(:) = lindz(:)*((ub3 - lb3)/(npts_local(3) - 1)) + lb3
337 x_locl(:) = x_glbl(lindx)
338 y_locl(:) = y_glbl(lindy)
339 z_locl(:) = z_glbl(lindz)
341 CALL timestop(handle)
365 SUBROUTINE pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, pw_in, pw_out)
367 TYPE(pw_pool_type),
INTENT(IN),
POINTER :: pw_pool
368 REAL(dp),
INTENT(IN) :: zeta
369 REAL(dp),
ALLOCATABLE,
DIMENSION(:),
INTENT(IN) :: x_glbl, y_glbl, z_glbl
370 TYPE(pw_r3d_rs_type),
INTENT(IN) :: pw_in
371 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: pw_out
373 CHARACTER(LEN=*),
PARAMETER :: routinen =
'pw_mollifier'
375 INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
377 INTEGER,
DIMENSION(2, 3) :: bounds, bounds_local
378 REAL(dp) :: normfact, xi, xmax, xmin, yj, ymax, &
380 REAL(dp),
DIMENSION(3, 3) :: dh
381 TYPE(pw_c1d_gs_type) :: g_gs, pw_in_gs, pw_out_gs
382 TYPE(pw_grid_type),
POINTER :: pw_grid
383 TYPE(pw_r3d_rs_type) :: g
385 CALL timeset(routinen, handle)
387 pw_grid => pw_in%pw_grid
389 bounds_local = pw_grid%bounds_local
390 bounds = pw_grid%bounds
392 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
393 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
394 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
396 CALL pw_pool%create_pw(g)
397 CALL pw_pool%create_pw(g_gs)
398 CALL pw_pool%create_pw(pw_in_gs)
399 CALL pw_pool%create_pw(pw_out_gs)
402 xmin = x_glbl(bounds(1, 1)); xmax = x_glbl(bounds(2, 1))
403 ymin = y_glbl(bounds(1, 2)); ymax = y_glbl(bounds(2, 2))
404 zmin = z_glbl(bounds(1, 3)); zmax = z_glbl(bounds(2, 3))
409 xi = x_glbl(i); yj = y_glbl(j); zk = z_glbl(k)
410 IF (norm2((/(xi - xmin), (yj - ymin), (zk - zmin)/)) .LT. zeta - small_value)
THEN
411 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmin), (yj - ymin), (zk - zmin)/)/zeta)**2 - 1))
412 ELSE IF (norm2((/(xi - xmax), (yj - ymax), (zk - zmax)/)) .LT. zeta - small_value)
THEN
413 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmax), (yj - ymax), (zk - zmax)/)/zeta)**2 - 1))
414 ELSE IF (norm2((/(xi - xmin), (yj - ymax), (zk - zmax)/)) .LT. zeta - small_value)
THEN
415 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmin), (yj - ymax), (zk - zmax)/)/zeta)**2 - 1))
416 ELSE IF (norm2((/(xi - xmax), (yj - ymin), (zk - zmax)/)) .LT. zeta - small_value)
THEN
417 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmax), (yj - ymin), (zk - zmax)/)/zeta)**2 - 1))
418 ELSE IF (norm2((/(xi - xmax), (yj - ymax), (zk - zmin)/)) .LT. zeta - small_value)
THEN
419 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmax), (yj - ymax), (zk - zmin)/)/zeta)**2 - 1))
420 ELSE IF (norm2((/(xi - xmin), (yj - ymin), (zk - zmax)/)) .LT. zeta - small_value)
THEN
421 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmin), (yj - ymin), (zk - zmax)/)/zeta)**2 - 1))
422 ELSE IF (norm2((/(xi - xmin), (yj - ymax), (zk - zmin)/)) .LT. zeta - small_value)
THEN
423 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmin), (yj - ymax), (zk - zmin)/)/zeta)**2 - 1))
424 ELSE IF (norm2((/(xi - xmax), (yj - ymin), (zk - zmin)/)) .LT. zeta - small_value)
THEN
425 g%array(i, j, k) = exp(1.0_dp/(norm2((/(xi - xmax), (yj - ymin), (zk - zmin)/)/zeta)**2 - 1))
430 CALL pw_scale(g, (1.0_dp/zeta)**3)
431 normfact = pw_integrate_function(g)
432 CALL pw_scale(g, 1.0_dp/normfact)
434 CALL pw_transfer(g, g_gs)
435 CALL pw_transfer(pw_in, pw_in_gs)
436 pw_out_gs%array = g_gs%array*pw_in_gs%array
437 CALL pw_transfer(pw_out_gs, pw_out)
440 CALL pw_scale(pw_out, real(pw_grid%ngpts, kind=dp))
442 CALL pw_scale(pw_out, pw_grid%dvol)
447 IF (abs(pw_out%array(i, j, k)) .LE. 1.0e-10_dp) pw_out%array(i, j, k) = 0.0_dp
452 CALL pw_pool%give_back_pw(g)
453 CALL pw_pool%give_back_pw(g_gs)
454 CALL pw_pool%give_back_pw(pw_in_gs)
455 CALL pw_pool%give_back_pw(pw_out_gs)
456 CALL timestop(handle)
Defines the basic variable types.
integer, parameter, public dp
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public rs_grid_create(rs, desc)
...
subroutine, public transfer_pw2rs(rs, pw)
...
subroutine, public transfer_rs2pw(rs, pw)
...
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
numerical operations on real-space grid
subroutine, public setup_grid_axes(pw_grid, x_glbl, y_glbl, z_glbl, x_locl, y_locl, z_locl)
returns the global axes and the portion of the axes that are local to the current mpi rank
subroutine, public derive_fdm_cd7(f, df, rs_grid)
6th order finite difference derivative of a function on realspace grid
subroutine, public pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, pw_in, pw_out)
convolutes a function with a smoothing kernel K_\zeta v * K_\zeta K_\zeta is the standard mollifier d...
subroutine, public derive_fdm_cd5(f, df, rs_grid)
4th order finite difference derivative of a function on realspace grid
subroutine, public derive_fdm_cd3(f, df, rs_grid)
2nd order finite difference derivative of a function on realspace grid