1!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
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!-----------------------------------------------------------------
20subroutine poisson_solve_direct(this, pot, rho)
21  type(poisson_t), intent(in)  :: this
22  FLOAT,           intent(out) :: pot(:)
23  FLOAT,           intent(in)  :: rho(:)
24
25  FLOAT                :: prefactor, aa1, aa2, aa3, aa4
26
27  integer              :: ip, jp
28  integer              :: dim_ele   !< physical dimensions  (= dimensions of simulation box      if electrons only)
29                                    !<                      (= dimensions of simulation box - 1  for dressed electrons)
30  integer              :: dim_eff   !< effective dimensions (= dim_ele,     if no soft coulomb)
31                                    !<                      (= dim_ele + 1  if soft coulomb)
32
33  FLOAT                :: xx1(1:MAX_DIM+1), xx2(1:MAX_DIM+1), xx3(1:MAX_DIM+1), xx4(1:MAX_DIM+1)
34  FLOAT                :: xx(1:MAX_DIM+1), yy(1:MAX_DIM+1)
35
36  logical              :: include_diag
37
38  FLOAT                :: xg(MAX_DIM)
39  integer, allocatable :: ip_v(:), part_v(:)
40  FLOAT, allocatable   :: pvec(:), tmp(:)
41
42  PUSH_SUB(poisson_solve_direct)
43
44  xx  = M_ZERO
45  xx1 = M_ZERO
46  xx2 = M_ZERO
47  xx3 = M_ZERO
48  xx4 = M_ZERO
49  yy  = M_ZERO
50
51  if (.not. this%is_dressed) then
52    dim_ele = this%der%mesh%sb%dim
53  else
54    dim_ele = this%der%mesh%sb%dim - 1
55  end if
56
57  if (this%poisson_soft_coulomb_param**2 > M_ZERO) then
58    dim_eff = dim_ele + 1
59    xx(dim_eff) = this%poisson_soft_coulomb_param
60    xx1(dim_eff) = this%poisson_soft_coulomb_param
61    xx2(dim_eff) = this%poisson_soft_coulomb_param
62    xx3(dim_eff) = this%poisson_soft_coulomb_param
63    xx4(dim_eff) = this%poisson_soft_coulomb_param
64    include_diag = .true.
65  else
66    dim_eff = dim_ele
67    include_diag = .false.
68  endif
69
70  ASSERT(this%method == POISSON_DIRECT_SUM)
71
72  select case(dim_ele)
73  case(3)
74    prefactor = M_TWO*M_PI*(M_THREE/(M_PI*M_FOUR))**(M_TWOTHIRD)
75  case(2)
76    prefactor = M_TWO*sqrt(M_PI)
77  case(1)
78    prefactor = M_ONE
79  case default
80    message(1) = "Internal error: poisson_solve_direct can only be called for 1D, 2D or 3D."
81    ! why not? all that is needed is the appropriate prefactors to be defined above, actually. then 1D, 4D etc. can be done
82    call messages_fatal(1)
83  end select
84
85  if(.not. this%der%mesh%use_curvilinear) then
86    prefactor = prefactor / (this%der%mesh%volume_element**(M_ONE/dim_ele))
87  end if
88
89  if(this%der%mesh%parallel_in_domains) then
90    SAFE_ALLOCATE(pvec(1:this%der%mesh%np))
91    SAFE_ALLOCATE(part_v(1:this%der%mesh%np_global))
92    SAFE_ALLOCATE(ip_v(1:this%der%mesh%np_global))
93    SAFE_ALLOCATE(tmp(1:this%der%mesh%np_global))
94    do ip = 1, this%der%mesh%np_global
95      ip_v(ip) = ip
96    end do
97    call partition_get_partition_number(this%der%mesh%inner_partition, this%der%mesh%np_global, ip_v, part_v)
98
99    pot = M_ZERO
100    do ip = 1, this%der%mesh%np_global
101      xg = mesh_x_global(this%der%mesh, ip)
102      xx(1:dim_ele) = xg(1:dim_ele)
103      if(this%der%mesh%use_curvilinear) then
104        do jp = 1, this%der%mesh%np
105          if(vec_global2local(this%der%mesh%vp, ip, this%der%mesh%vp%partno) == jp .and. .not. include_diag) then
106            pvec(jp) = rho(jp)*prefactor**(M_ONE - M_ONE/dim_ele)
107          else
108            yy(1:dim_ele) = this%der%mesh%x(jp, 1:dim_ele)
109            pvec(jp) = rho(jp)/sqrt(sum((xx(1:dim_eff) - yy(1:dim_eff))**2))
110          end if
111        end do
112      else
113        do jp = 1, this%der%mesh%np
114          if (vec_global2local(this%der%mesh%vp, ip, this%der%mesh%vp%partno) == jp .and. .not. include_diag) then
115            pvec(jp) = rho(jp)*prefactor
116          else
117            yy(1:dim_ele) = this%der%mesh%x(jp, 1:dim_ele)
118            pvec(jp) = rho(jp)/sqrt(sum((xx(1:dim_eff) - yy(1:dim_eff))**2))
119          end if
120        end do
121      end if
122      tmp(ip) = dmf_integrate(this%der%mesh, pvec, reduce = .false.)
123    end do
124
125#ifdef HAVE_MPI
126    call comm_allreduce(this%der%mesh%mpi_grp%comm, tmp)
127#endif
128
129    do ip = 1, this%der%mesh%np_global
130      if (part_v(ip) == this%der%mesh%vp%partno) then
131        pot(vec_global2local(this%der%mesh%vp, ip, this%der%mesh%vp%partno)) = tmp(ip)
132      end if
133    end do
134
135    SAFE_DEALLOCATE_A(pvec)
136    SAFE_DEALLOCATE_A(tmp)
137
138  else ! serial mode
139
140    do ip = 1, this%der%mesh%np - 4 + 1, 4
141
142      xx1(1:dim_ele) = this%der%mesh%x(ip    , 1:dim_ele)
143      xx2(1:dim_ele) = this%der%mesh%x(ip + 1, 1:dim_ele)
144      xx3(1:dim_ele) = this%der%mesh%x(ip + 2, 1:dim_ele)
145      xx4(1:dim_ele) = this%der%mesh%x(ip + 3, 1:dim_ele)
146
147      if (this%der%mesh%use_curvilinear) then
148
149        if (.not. include_diag) then
150          aa1 = prefactor*rho(ip    )*this%der%mesh%vol_pp(ip    )**(M_ONE - M_ONE/dim_ele)
151          aa2 = prefactor*rho(ip + 1)*this%der%mesh%vol_pp(ip + 1)**(M_ONE - M_ONE/dim_ele)
152          aa3 = prefactor*rho(ip + 2)*this%der%mesh%vol_pp(ip + 2)**(M_ONE - M_ONE/dim_ele)
153          aa4 = prefactor*rho(ip + 3)*this%der%mesh%vol_pp(ip + 3)**(M_ONE - M_ONE/dim_ele)
154        else
155          aa1 = M_ZERO
156          aa2 = M_ZERO
157          aa3 = M_ZERO
158          aa4 = M_ZERO
159        end if
160
161        !$omp parallel do reduction(+:aa1,aa2,aa3,aa4) firstprivate(yy)
162        do jp = 1, this%der%mesh%np
163          yy(1:dim_ele) = this%der%mesh%x(jp, 1:dim_ele)
164          if (ip     /= jp .or. include_diag) then
165            aa1 = aa1 + rho(jp)/sqrt(sum((xx1(1:dim_eff) - yy(1:dim_eff))**2)) * this%der%mesh%vol_pp(jp)
166          end if
167          if (ip + 1 /= jp .or. include_diag) then
168            aa2 = aa2 + rho(jp)/sqrt(sum((xx2(1:dim_eff) - yy(1:dim_eff))**2)) * this%der%mesh%vol_pp(jp)
169          end if
170          if (ip + 2 /= jp .or. include_diag) then
171            aa3 = aa3 + rho(jp)/sqrt(sum((xx3(1:dim_eff) - yy(1:dim_eff))**2)) * this%der%mesh%vol_pp(jp)
172          end if
173          if (ip + 3 /= jp .or. include_diag) then
174            aa4 = aa4 + rho(jp)/sqrt(sum((xx4(1:dim_eff) - yy(1:dim_eff))**2)) * this%der%mesh%vol_pp(jp)
175          end if
176        end do
177
178      else
179
180        if(.not. include_diag) then
181          aa1 = prefactor*rho(ip    )
182          aa2 = prefactor*rho(ip + 1)
183          aa3 = prefactor*rho(ip + 2)
184          aa4 = prefactor*rho(ip + 3)
185        else
186          aa1 = M_ZERO
187          aa2 = M_ZERO
188          aa3 = M_ZERO
189          aa4 = M_ZERO
190        end if
191
192        !$omp parallel do reduction(+:aa1,aa2,aa3,aa4) firstprivate(yy)
193        do jp = 1, this%der%mesh%np
194          yy(1:dim_ele) = this%der%mesh%x(jp, 1:dim_ele)
195          if (ip     /= jp .or. include_diag) aa1 = aa1 + rho(jp)/sqrt(sum((xx1(1:dim_eff) - yy(1:dim_eff))**2))
196          if (ip + 1 /= jp .or. include_diag) aa2 = aa2 + rho(jp)/sqrt(sum((xx2(1:dim_eff) - yy(1:dim_eff))**2))
197          if (ip + 2 /= jp .or. include_diag) aa3 = aa3 + rho(jp)/sqrt(sum((xx3(1:dim_eff) - yy(1:dim_eff))**2))
198          if (ip + 3 /= jp .or. include_diag) aa4 = aa4 + rho(jp)/sqrt(sum((xx4(1:dim_eff) - yy(1:dim_eff))**2))
199        end do
200
201      end if
202
203      pot(ip    ) = this%der%mesh%volume_element*aa1
204      pot(ip + 1) = this%der%mesh%volume_element*aa2
205      pot(ip + 2) = this%der%mesh%volume_element*aa3
206      pot(ip + 3) = this%der%mesh%volume_element*aa4
207
208    end do
209
210    do ip = ip, this%der%mesh%np
211
212      aa1 = CNST(0.0)
213
214      xx1(1:dim_ele) = this%der%mesh%x(ip,1:dim_ele)
215
216      if (this%der%mesh%use_curvilinear) then
217        do jp = 1, this%der%mesh%np
218          yy(1:dim_ele) = this%der%mesh%x(jp, 1:dim_ele)
219          if (ip == jp .and. .not. include_diag) then
220            aa1 = aa1 + prefactor*rho(ip)*this%der%mesh%vol_pp(jp)**(M_ONE - M_ONE/dim_ele)
221          else
222            aa1 = aa1 + rho(jp)/sqrt(sum((xx1(1:dim_eff) - yy(1:dim_eff))**2))*this%der%mesh%vol_pp(jp)
223          end if
224        end do
225      else
226        do jp = 1, this%der%mesh%np
227          yy(1:dim_ele) = this%der%mesh%x(jp, 1:dim_ele)
228          if (ip == jp .and. .not. include_diag) then
229            aa1 = aa1 + prefactor*rho(ip)
230          else
231            aa1 = aa1 + rho(jp)/sqrt(sum((xx1(1:dim_eff) - yy(1:dim_eff))**2))
232          end if
233        end do
234      end if
235
236      pot(ip) = this%der%mesh%volume_element*aa1
237
238    end do
239
240  end if
241
242  POP_SUB(poisson_solve_direct)
243end subroutine poisson_solve_direct
244
245!! Local Variables:
246!! mode: f90
247!! coding: utf-8
248!! End:
249