(implicit ps) bugfix in setting up arbitrary planar dirichlet regions…

…, and improve tiling of such regions; consider forward fft scaling prefactor for smoothing via convolution with a mollifier function.

svn-origin-rev: 15268

src/pw/dirichlet_bc_methods.F

@@ -279,7 +279,7 @@ SUBROUTINE dirichlet_bc_partition(v_D, smooth, zeta, n_prtn, pw_pool, x_glbl, y_
     REAL(dp)                                 :: phi1, phi2
     REAL(dp), DIMENSION(3)                   :: rot_axis1, rot_axis2
     TYPE(cp_logger_type), POINTER            :: logger
-    TYPE(cs_rectangle), POINTER              :: rectangle_aa
+    TYPE(cs_rectangle), POINTER              :: rectangle_aa, rectangle_tmp
     TYPE(tile_p_type), DIMENSION(:), POINTER :: tiles
 
     CALL timeset(routineN,handle)
@@ -315,7 +315,10 @@ SUBROUTINE dirichlet_bc_partition(v_D, smooth, zeta, n_prtn, pw_pool, x_glbl, y_
       DO k = 1, n_tiles
          ALLOCATE(dirichlet_bc%tiles(k)%tile)
          ALLOCATE(dirichlet_bc%tiles(k)%tile%rectangle)
-         CALL rotate_rectangle(tiles(k)%tile%rectangle, phi1, rot_axis1, BWROT, &
+         ALLOCATE(rectangle_tmp)
+         CALL rotate_rectangle(tiles(k)%tile%rectangle, phi2, rot_axis2, BWROT, &
+                               rectangle_tmp, error)
+         CALL rotate_rectangle(rectangle_tmp, phi1, rot_axis1, BWROT, &
                                dirichlet_bc%tiles(k)%tile%rectangle, error)
 
          dirichlet_bc%tiles(k)%tile%tile_id = 8000 + k
@@ -341,6 +344,7 @@ SUBROUTINE dirichlet_bc_partition(v_D, smooth, zeta, n_prtn, pw_pool, x_glbl, y_
          tile_npts = NINT(SUM(dirichlet_bc%tiles(k)%tile%tile_pw%cr3d),KIND=KIND(tile_npts))
          CALL mp_sum(tile_npts, pw_pool%pw_grid%para%group)
          dirichlet_bc%tiles(k)%tile%npts = tile_npts
+         CALL cs_rectangle_release(rectangle_tmp, error)
       END DO
 
       IF ((unit_nr .GT. 0) .AND. debug) WRITE(unit_nr, '(T3,A)') REPEAT('=', 78) 

src/pw/dirichlet_bc_utils.F

@@ -878,9 +878,10 @@ SUBROUTINE partition_aa_rectangle_into_tiles(rectangle, x_glbl, y_glbl, z_glbl,
                                                 yprtn_ubind
     INTEGER, DIMENSION(n_prtn(3))            :: zprtn_lbind, zprtn_npts, &
                                                 zprtn_ubind
-    REAL(dp)                                 :: step
-    REAL(dp), ALLOCATABLE, DIMENSION(:)      :: xglbl_sym, yglbl_sym, &
-                                                zglbl_sym
+    REAL(dp)                                 :: dx, dy, dz, step
+    REAL(dp), ALLOCATABLE, DIMENSION(:)      :: xglbl_dbl, xglbl_sym, &
+                                                yglbl_dbl, yglbl_sym, &
+                                                zglbl_dbl, zglbl_sym
     REAL(dp), DIMENSION(2)                   :: x_xtnt, y_xtnt, z_xtnt
     REAL(dp), DIMENSION(3)                   :: A, B, C, D
     REAL(dp), DIMENSION(n_prtn(1))           :: xprtn_lb, xprtn_ub
@@ -889,12 +890,20 @@ SUBROUTINE partition_aa_rectangle_into_tiles(rectangle, x_glbl, y_glbl, z_glbl,
 
     CALL timeset(routineN,handle)
 
-    ALLOCATE(xglbl_sym(2*SIZE(x_glbl)-1), yglbl_sym(2*SIZE(y_glbl)-1), zglbl_sym(2*SIZE(z_glbl)-1))
+    ALLOCATE(xglbl_dbl(2*SIZE(x_glbl)), yglbl_dbl(2*SIZE(y_glbl)), zglbl_dbl(2*SIZE(z_glbl)))
+    ALLOCATE(xglbl_sym(2*SIZE(xglbl_dbl)-1), yglbl_sym(2*SIZE(yglbl_dbl)-1), zglbl_sym(2*SIZE(zglbl_dbl)-1))
 
 ! extend the global axes and make them symmetric
-    xglbl_sym(:) = (/ (- x_glbl(i), i = UBOUND(x_glbl,1), LBOUND(x_glbl,1)+1, -1), x_glbl /)
-    yglbl_sym(:) = (/ (- y_glbl(i), i = UBOUND(y_glbl,1), LBOUND(y_glbl,1)+1, -1), y_glbl /)
-    zglbl_sym(:) = (/ (- z_glbl(i), i = UBOUND(z_glbl,1), LBOUND(z_glbl,1)+1, -1), z_glbl /)
+    dx = x_glbl(LBOUND(x_glbl,1)+1) - x_glbl(LBOUND(x_glbl,1))
+    dy = y_glbl(LBOUND(y_glbl,1)+1) - y_glbl(LBOUND(y_glbl,1))
+    dz = z_glbl(LBOUND(z_glbl,1)+1) - z_glbl(LBOUND(z_glbl,1))
+
+    xglbl_dbl(:) = (/ x_glbl, x_glbl+SIZE(x_glbl)*dx /)
+    yglbl_dbl(:) = (/ y_glbl, y_glbl+SIZE(y_glbl)*dy /)
+    zglbl_dbl(:) = (/ z_glbl, z_glbl+SIZE(z_glbl)*dz /)
+    xglbl_sym(:) = (/ (- xglbl_dbl(i), i = UBOUND(xglbl_dbl,1), LBOUND(xglbl_dbl,1)+1, -1), xglbl_dbl /)
+    yglbl_sym(:) = (/ (- yglbl_dbl(i), i = UBOUND(yglbl_dbl,1), LBOUND(yglbl_dbl,1)+1, -1), yglbl_dbl /)
+    zglbl_sym(:) = (/ (- zglbl_dbl(i), i = UBOUND(zglbl_dbl,1), LBOUND(zglbl_dbl,1)+1, -1), zglbl_dbl /)
 
 ! find the extents of the rectangle in all three directions
     x_xtnt(1) = MINVAL(rectangle%vertices(1,:));    x_xtnt(2) = MAXVAL(rectangle%vertices(1,:))

src/pw/ps_implicit_methods.F

@@ -34,6 +34,7 @@ MODULE ps_implicit_methods
                                              pw_copy,&
                                              pw_derive,&
                                              pw_integral_ab,&
+                                             pw_integrate_function,&
                                              pw_scale,&
                                              pw_transfer,&
                                              pw_zero
@@ -2279,15 +2280,19 @@ SUBROUTINE pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, pw_in, pw_out, er
       END DO
     END DO
     CALL pw_scale(G, (1.0_dp/zeta)**3, error)
-    normfact = SUM ( G%cr3d ( :, :, : ) )
-    CALL mp_sum(normfact, pw_pool%pw_grid%para%group)
+    normfact = pw_integrate_function(G,error=error)
     CALL pw_scale(G, 1.0_dp/normfact, error)
 
     CALL pw_transfer(G, G_gs, error=error)
     CALL pw_transfer(pw_in, pw_in_gs, error=error)
     pw_out_gs%cc = G_gs%cc * pw_in_gs%cc
     CALL pw_transfer(pw_out_gs, pw_out, error=error)
 
+    ! multiply by the reciprocal of the forward Fourier transform normalization prefactor (here 1/N, by convention)
+    CALL pw_scale(pw_out, REAL(pw_grid%ngpts,KIND=dp), error=error)
+    ! from discrete convolution to continuous convolution
+    CALL pw_scale(pw_out, pw_grid%dvol, error=error)
+
     DO k = lb3, ub3
       DO j = lb2, ub2
         DO i = lb1, ub1

tests/QS/regtest-ps-implicit-1-2/TEST_FILES_RESET

@@ -9,3 +9,6 @@ Ar_mixed_periodic_planar.inp
 Ar_mixed_periodic_cuboidal.inp     
 Ar_mixed_periodic_cylindrical.inp  
 Ar_neumann.inp                     
+
+# bugfix in constructing arbitrary planar regions and improve tiling of such regions.
+Ar_mixed_periodic_planar.inp       

tests/QS/regtest-ps-implicit-1-3/TEST_FILES_RESET

@@ -7,3 +7,6 @@ Ar_mixed_aa_planar-ns_cell.inp
 Ar_mixed_aa_planar.inp           
 Ar_mixed_planar.inp              
 Ar_mixed_aa_planar-ns_cell.inp   
+
+# bugfix in constructing arbitrary planar regions and improve tiling of such regions.
+Ar_mixed_planar.inp              

tests/QS/regtest-ps-implicit-2-2/TEST_FILES_RESET

@@ -5,3 +5,6 @@ H2O_mixed_periodic_planar.inp
 # properly include the dielectric contribution to the Hamiltonian, v_eps.
 H2O_mixed_periodic_aa_planar.inp   
 H2O_mixed_periodic_planar.inp      
+
+# bugfix in constructing arbitrary planar regions and improve tiling of such regions.
+H2O_mixed_periodic_planar.inp