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