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