1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module poisson_corrections_oct_m 22 use derivatives_oct_m 23 use global_oct_m 24 use lalg_basic_oct_m 25 use loct_math_oct_m 26 use math_oct_m 27 use mesh_function_oct_m 28 use mesh_oct_m 29 use messages_oct_m 30 use namespace_oct_m 31 use nl_operator_oct_m 32 use parser_oct_m 33 use profiling_oct_m 34 use simul_box_oct_m 35 36 implicit none 37 38 private 39 public :: & 40 poisson_corrections_init, & 41 poisson_corrections_end, & 42 correct_rho, & 43 poisson_boundary_conditions, & 44 poisson_corr_t, & 45 internal_laplacian_op, & 46 internal_dotp, & 47 der_pointer, & 48 mesh_pointer 49 50 FLOAT, parameter :: alpha_ = M_FIVE 51 52 type poisson_corr_t 53 private 54 integer :: method 55 integer :: maxl 56 FLOAT, pointer :: phi(:, :) 57 FLOAT, pointer :: aux(:, :) 58 FLOAT, pointer :: gaussian(:) 59 end type poisson_corr_t 60 61 type(derivatives_t), pointer :: der_pointer 62 type(mesh_t), pointer :: mesh_pointer 63 64 integer, parameter :: & 65 CORR_MULTIPOLE = 1, & 66 CORR_EXACT = 3 67 68contains 69 70 ! --------------------------------------------------------- 71 subroutine poisson_corrections_init(this, namespace, ml, mesh) 72 type(poisson_corr_t), intent(out) :: this 73 type(namespace_t), intent(in) :: namespace 74 integer, intent(in) :: ml 75 type(mesh_t), intent(in) :: mesh 76 77 FLOAT :: alpha, gamma, ylm, rr, xx(MAX_DIM) 78 integer :: ip, ll, add_lm, lldfac, jj, mm 79 80 PUSH_SUB(poisson_corrections_init) 81 82 if(simul_box_is_periodic(mesh%sb)) & 83 call messages_not_implemented("Poisson boundary corrections for periodic systems") 84 85 !%Variable PoissonSolverBoundaries 86 !%Type integer 87 !%Section Hamiltonian::Poisson 88 !%Default multipole 89 !%Description 90 !% For finite systems, some Poisson solvers (<tt>multigrid</tt>, 91 !% <tt>cg_corrected</tt>, and <tt>fft</tt> with <tt>PoissonFFTKernel = multipole_correction</tt>) 92 !% require the calculation of the 93 !% boundary conditions with an auxiliary method. This variable selects that method. 94 !%Option multipole 1 95 !% A multipole expansion of the density is used to approximate the potential on the boundaries. 96 !%Option exact 3 97 !% An exact integration of the Poisson equation is done over the boundaries. This option is 98 !% experimental, and not implemented for domain parallelization. 99 !%End 100 call parse_variable(namespace, 'PoissonSolverBoundaries', CORR_MULTIPOLE, this%method) 101 102 select case(this%method) 103 case(CORR_MULTIPOLE) 104 this%maxl = ml 105 106 add_lm = 0 107 do ll = 0, this%maxl 108 do mm = -ll, ll 109 add_lm = add_lm + 1 110 end do 111 end do 112 113 SAFE_ALLOCATE(this%phi(1:mesh%np, 1:add_lm)) 114 SAFE_ALLOCATE(this%aux(1:mesh%np, 1:add_lm)) 115 SAFE_ALLOCATE(this%gaussian(1:mesh%np)) 116 117 alpha = alpha_ * mesh%spacing(1) 118 do ip = 1, mesh%np 119 call mesh_r(mesh, ip, rr, coords = xx) 120 this%gaussian(ip) = exp(-(rr/alpha)**2) 121 add_lm = 1 122 do ll = 0, this%maxl 123 lldfac = 1 124 do jj = 1, 2*ll+1, 2 125 lldfac = lldfac * jj 126 end do 127 gamma = ( sqrt(M_PI)*2**(ll+3) ) / lldfac 128 do mm = -ll, ll 129 call grylmr(xx(1), xx(2), xx(3), ll, mm, ylm) 130 if(rr > M_EPSILON) then 131 this%phi(ip, add_lm) = gamma*isubl(ll, rr/alpha)*ylm/rr**(ll+1) 132 this%aux(ip, add_lm) = rr**ll*ylm 133 else 134 this%phi(ip, add_lm) = gamma*ylm / alpha 135 if(ll == 0) then 136 this%aux(ip, add_lm) = ylm 137 else 138 this%aux(ip, add_lm) = M_ZERO 139 end if 140 end if 141 add_lm = add_lm + 1 142 end do 143 end do 144 end do 145 146 case(CORR_EXACT) 147 call messages_experimental('Exact Poisson solver boundaries') 148 if(mesh%parallel_in_domains) call messages_not_implemented('Exact Poisson solver boundaries with domain parallelization') 149 150 end select 151 152 POP_SUB(poisson_corrections_init) 153 154 contains 155 156 ! --------------------------------------------------------- 157 FLOAT function isubl( ll, xx) 158 integer, intent(in) :: ll 159 FLOAT, intent(in) :: xx 160 161 ! no push_sub, called too frequently 162 isubl = M_HALF*loct_gamma(ll + M_HALF)*(M_ONE - loct_incomplete_gamma(ll+M_HALF, xx**2) ) 163 164 end function isubl 165 166 end subroutine poisson_corrections_init 167 168 169 ! --------------------------------------------------------- 170 subroutine poisson_corrections_end(this) 171 type(poisson_corr_t), intent(inout) :: this 172 173 PUSH_SUB(poisson_corrections_end) 174 175 select case(this%method) 176 case(CORR_MULTIPOLE) 177 SAFE_DEALLOCATE_P(this%phi) 178 SAFE_DEALLOCATE_P(this%aux) 179 SAFE_DEALLOCATE_P(this%gaussian) 180 case(CORR_EXACT) 181 end select 182 183 POP_SUB(poisson_corrections_end) 184 end subroutine poisson_corrections_end 185 186 187 ! --------------------------------------------------------- 188 subroutine correct_rho(this, der, rho, rho_corrected, vh_correction) 189 type(poisson_corr_t), intent(in) :: this 190 type(derivatives_t), intent(in) :: der 191 FLOAT, intent(in) :: rho(:) 192 FLOAT, intent(out) :: rho_corrected(:) 193 FLOAT, intent(out) :: vh_correction(:) 194 195 integer :: ip, add_lm, ll, mm, lldfac, jj, ip2 196 FLOAT :: alpha, vv, rr 197 FLOAT, allocatable :: mult(:) 198 FLOAT, allocatable :: betal(:) 199 type(profile_t), save :: prof 200 201 PUSH_SUB(correct_rho) 202 call profiling_in(prof, "POISSON_CORRECT") 203 204 ASSERT(ubound(vh_correction, dim = 1) == der%mesh%np_part) 205 206 select case(this%method) 207 case(CORR_MULTIPOLE) 208 209 SAFE_ALLOCATE(mult(1:(this%maxl+1)**2)) 210 call get_multipoles(this, der%mesh, rho, this%maxl, mult) 211 212 alpha = alpha_ * der%mesh%spacing(1) 213 214 SAFE_ALLOCATE(betal(1:(this%maxl+1)**2)) 215 add_lm = 1 216 do ll = 0, this%maxl 217 do mm = -ll, ll 218 lldfac = 1 219 do jj = 1, 2*ll+1, 2 220 lldfac = lldfac*jj 221 end do 222 betal(add_lm) = (2**(ll+2))/( alpha**(2*ll+3) * sqrt(M_PI) * lldfac ) 223 add_lm = add_lm + 1 224 end do 225 end do 226 227 call lalg_copy(der%mesh%np, rho, rho_corrected) 228 vh_correction = M_ZERO 229 add_lm = 1 230 do ll = 0, this%maxl 231 do mm = -ll, ll 232 do ip = 1, der%mesh%np 233 rho_corrected(ip) = rho_corrected(ip) - mult(add_lm)*betal(add_lm)*this%aux(ip, add_lm)*this%gaussian(ip) 234 end do 235 call lalg_axpy(der%mesh%np, mult(add_lm), this%phi(:, add_lm), vh_correction) 236 add_lm = add_lm + 1 237 end do 238 end do 239 240 SAFE_DEALLOCATE_A(mult) 241 SAFE_DEALLOCATE_A(betal) 242 243 case(CORR_EXACT) 244 245 do ip = 1, der%mesh%np 246 vh_correction(ip) = M_ZERO 247 end do 248 249 do ip = der%mesh%np + 1, der%mesh%np_part 250 vv = M_ZERO 251 do ip2 = 1, der%mesh%np 252 rr = sqrt(sum((der%mesh%x(ip, 1:der%mesh%sb%dim) - der%mesh%x(ip2, 1:der%mesh%sb%dim))**2)) 253 vv = vv + rho(ip2)/rr 254 end do 255 vh_correction(ip) = der%mesh%volume_element*vv 256 end do 257 258 ASSERT(.not. nl_operator_compact_boundaries(der%lapl)) 259 260 call dderivatives_lapl(der, vh_correction, rho_corrected, set_bc = .false.) 261 262 do ip = 1, der%mesh%np 263 rho_corrected(ip) = rho(ip) + CNST(1.0)/(CNST(4.0)*M_PI)*rho_corrected(ip) 264 end do 265 266 end select 267 268 call profiling_out(prof) 269 POP_SUB(correct_rho) 270 end subroutine correct_rho 271 272 273 ! --------------------------------------------------------- 274 subroutine get_multipoles(this, mesh, rho, ml, mult) 275 type(poisson_corr_t), intent(in) :: this 276 type(mesh_t), intent(in) :: mesh 277 FLOAT, intent(in) :: rho(:) !< rho(mesh%np) 278 integer, intent(in) :: ml 279 FLOAT, intent(out) :: mult((ml+1)**2) 280 281 integer :: add_lm, ll, mm 282 283 PUSH_SUB(get_multipoles) 284 285 mult(:) = M_ZERO 286 add_lm = 1 287 do ll = 0, ml 288 do mm = -ll, ll 289 mult(add_lm) = dmf_dotp(mesh, rho, this%aux(:, add_lm)) 290 add_lm = add_lm + 1 291 end do 292 end do 293 294 POP_SUB(get_multipoles) 295 end subroutine get_multipoles 296 297 ! --------------------------------------------------------- 298 subroutine internal_laplacian_op(xx, yy) 299 FLOAT, intent(inout) :: xx(:) 300 FLOAT, intent(out) :: yy(:) 301 302 PUSH_SUB(internal_laplacian_op) 303 call dderivatives_lapl(der_pointer, xx, yy) 304 POP_SUB(internal_laplacian_op) 305 306 end subroutine internal_laplacian_op 307 308 309 ! --------------------------------------------------------- 310 FLOAT function internal_dotp(xx, yy) result(res) 311 FLOAT, intent(inout) :: xx(:) 312 FLOAT, intent(in) :: yy(:) 313 314 PUSH_SUB(internal_dotp) 315 316 res = dmf_dotp(mesh_pointer, xx, yy) 317 POP_SUB(internal_dotp) 318 end function internal_dotp 319 320 321 ! --------------------------------------------------------- 322 subroutine poisson_boundary_conditions(this, mesh, rho, pot) 323 type(poisson_corr_t), intent(in) :: this 324 type(mesh_t), intent(in) :: mesh 325 FLOAT, intent(in) :: rho(:) !< rho(mesh%np) 326 FLOAT, intent(inout) :: pot(:) !< pot(mesh%np_part) 327 328 integer :: ip, add_lm, ll, mm, bp_lower 329 FLOAT :: xx(MAX_DIM), rr, s1, sa 330 FLOAT, allocatable :: mult(:) 331 332 PUSH_SUB(poisson_boundary_conditions) 333 334 SAFE_ALLOCATE(mult(1:(this%maxl+1)**2)) 335 336 call get_multipoles(this, mesh, rho, this%maxl, mult) 337 338 bp_lower = mesh%np + 1 339 if(mesh%parallel_in_domains) then 340 bp_lower = bp_lower + mesh%vp%np_ghost 341 end if 342 343 pot(bp_lower:mesh%np_part) = M_ZERO 344 do ip = bp_lower, mesh%np_part ! boundary conditions 345 call mesh_r(mesh, ip, rr, coords = xx) 346 add_lm = 1 347 do ll = 0, this%maxl 348 s1 = M_FOUR*M_PI/((M_TWO*ll + M_ONE)*rr**(ll + 1)) 349 do mm = -ll, ll 350 call grylmr(xx(1), xx(2), xx(3), ll, mm, sa) 351 pot(ip) = pot(ip) + sa * mult(add_lm) * s1 352 add_lm = add_lm+1 353 end do 354 end do 355 end do 356 357 SAFE_DEALLOCATE_A(mult) 358 POP_SUB(poisson_boundary_conditions) 359 end subroutine poisson_boundary_conditions 360 361end module poisson_corrections_oct_m 362 363!! Local Variables: 364!! mode: f90 365!! coding: utf-8 366!! End: 367