1!! Copyright (C) 2005-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch, X. Andrade 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_multigrid_oct_m 22 use boundaries_oct_m 23 use derivatives_oct_m 24 use global_oct_m 25 use lalg_basic_oct_m 26 use mesh_oct_m 27 use mesh_function_oct_m 28 use messages_oct_m 29 use multigrid_oct_m 30 use namespace_oct_m 31 use operate_f_oct_m 32 use parser_oct_m 33 use par_vec_oct_m 34 use poisson_corrections_oct_m 35 use profiling_oct_m 36 use varinfo_oct_m 37 38 implicit none 39 40 integer, parameter :: & 41 GAUSS_SEIDEL = 1, & 42 GAUSS_JACOBI = 2, & 43 GAUSS_JACOBI2 = 3 44 45 private 46 public :: & 47 poisson_multigrid_solver, & 48 poisson_multigrid_init, & 49 poisson_multigrid_end, & 50 mg_solver_t 51 52 type mg_solver_t 53 private 54 55 FLOAT :: & 56 threshold, & 57 relax_factor 58 59 integer :: & 60 maxcycles, & 61 presteps, & 62 poststeps, & 63 restriction_method, & 64 relaxation_method 65 66 type(poisson_corr_t) :: corrector 67 68 end type mg_solver_t 69 70contains 71 72 ! --------------------------------------------------------- 73 subroutine poisson_multigrid_init(this, namespace, mesh, ml, thr) 74 type(mg_solver_t), intent(out) :: this 75 type(namespace_t), intent(in) :: namespace 76 type(mesh_t), intent(inout) :: mesh 77 integer, intent(in) :: ml 78 FLOAT, intent(in) :: thr 79 80 PUSH_SUB(poisson_multigrid_init) 81 82 call poisson_corrections_init(this%corrector, namespace, ml, mesh) 83 84 this%threshold = thr 85 86 !%Variable PoissonSolverMGPresmoothingSteps 87 !%Type integer 88 !%Default 1 89 !%Section Hamiltonian::Poisson::Multigrid 90 !%Description 91 !% Number of Gauss-Seidel smoothing steps before coarse-level 92 !% correction in the multigrid Poisson solver. 93 !%End 94 call parse_variable(namespace, 'PoissonSolverMGPresmoothingSteps', 1, this%presteps) 95 96 !%Variable PoissonSolverMGPostsmoothingSteps 97 !%Type integer 98 !%Default 4 99 !%Section Hamiltonian::Poisson::Multigrid 100 !%Description 101 !% Number of Gauss-Seidel smoothing steps after coarse-level 102 !% correction in the multigrid Poisson solver. 103 !%End 104 call parse_variable(namespace, 'PoissonSolverMGPostsmoothingSteps', 4, this%poststeps) 105 106 !%Variable PoissonSolverMGMaxCycles 107 !%Type integer 108 !%Default 60 109 !%Section Hamiltonian::Poisson::Multigrid 110 !%Description 111 !% Maximum number of multigrid cycles that are performed if 112 !% convergence is not achieved. 113 !%End 114 call parse_variable(namespace, 'PoissonSolverMGMaxCycles', 50, this%maxcycles) 115 116 !%Variable PoissonSolverMGRestrictionMethod 117 !%Type integer 118 !%Default fullweight 119 !%Section Hamiltonian::Poisson::Multigrid 120 !%Description 121 !% Method used from fine-to-coarse grid transfer. 122 !%Option injection 1 123 !% Injection 124 !%Option fullweight 2 125 !% Fullweight restriction 126 !%End 127 call parse_variable(namespace, 'PoissonSolverMGRestrictionMethod', 2, this%restriction_method) 128 if(.not.varinfo_valid_option('PoissonSolverMGRestrictionMethod', this%restriction_method)) then 129 call messages_input_error(namespace, 'PoissonSolverMGRestrictionMethod') 130 end if 131 call messages_print_var_option(stdout, "PoissonSolverMGRestrictionMethod", this%restriction_method) 132 133 !%Variable PoissonSolverMGRelaxationMethod 134 !%Type integer 135 !%Section Hamiltonian::Poisson::Multigrid 136 !%Description 137 !% Method used to solve the linear system approximately in each grid for the 138 !% multigrid procedure that solves Poisson equation. Default is <tt>gauss_seidel</tt>, 139 !% unless curvilinear coordinates are used, in which case the default is <tt>gauss_jacobi</tt>. 140 !%Option gauss_seidel 1 141 !% Gauss-Seidel. 142 !%Option gauss_jacobi 2 143 !% Gauss-Jacobi. 144 !%Option gauss_jacobi2 3 145 !% Alternative implementation of Gauss-Jacobi. 146 !%End 147 if ( mesh%use_curvilinear ) then 148 call parse_variable(namespace, 'PoissonSolverMGRelaxationMethod', GAUSS_JACOBI, this%relaxation_method) 149 else 150 call parse_variable(namespace, 'PoissonSolverMGRelaxationMethod', GAUSS_SEIDEL, this%relaxation_method) 151 end if 152 153 if(.not.varinfo_valid_option('PoissonSolverMGRelaxationMethod', this%relaxation_method)) then 154 call messages_input_error(namespace, 'PoissonSolverMGRelaxationMethod') 155 end if 156 call messages_print_var_option(stdout, "PoissonSolverMGRelaxationMethod", this%relaxation_method) 157 158 !%Variable PoissonSolverMGRelaxationFactor 159 !%Type float 160 !%Section Hamiltonian::Poisson::Multigrid 161 !%Description 162 !% Relaxation factor of the relaxation operator used for the 163 !% multigrid method. This is mainly for debugging, 164 !% since overrelaxation does not help in a multigrid scheme. 165 !% The default is 1.0, except 0.6666 for the <tt>gauss_jacobi</tt> method. 166 !%End 167 if ( this%relaxation_method == GAUSS_JACOBI) then 168 call parse_variable(namespace, 'PoissonSolverMGRelaxationFactor', CNST(0.6666), this%relax_factor ) 169 else 170 call parse_variable(namespace, 'PoissonSolverMGRelaxationFactor', M_ONE, this%relax_factor) 171 end if 172 173 POP_SUB(poisson_multigrid_init) 174 end subroutine poisson_multigrid_init 175 176 177 ! --------------------------------------------------------- 178 subroutine poisson_multigrid_end(this) 179 type(mg_solver_t), intent(inout) :: this 180 181 PUSH_SUB(poisson_multigrid_end) 182 183 call poisson_corrections_end(this%corrector) 184 185 POP_SUB(poisson_multigrid_end) 186 end subroutine poisson_multigrid_end 187 188 189 ! --------------------------------------------------------- 190 subroutine poisson_multigrid_solver(this, der, pot, rho) 191 type(mg_solver_t), intent(in) :: this 192 type(derivatives_t), intent(in) :: der 193 FLOAT, intent(inout) :: pot(:) 194 FLOAT, intent(in) :: rho(:) 195 196 integer :: iter, ip 197 FLOAT :: resnorm 198 FLOAT, allocatable :: vh_correction(:), res(:), cor(:), err(:) 199 200 PUSH_SUB(poisson_multigrid_solver) 201 202 ! correction for treating boundaries 203 SAFE_ALLOCATE(vh_correction(1:der%mesh%np_part)) 204 SAFE_ALLOCATE(res(1:der%mesh%np)) 205 SAFE_ALLOCATE(cor(1:der%mesh%np_part)) 206 SAFE_ALLOCATE(err(1:der%mesh%np)) 207 208 call correct_rho(this%corrector, der, rho, res, vh_correction) 209 call lalg_scal(der%mesh%np, -M_FOUR*M_PI, res) 210 211 do ip = 1, der%mesh%np 212 cor(ip) = pot(ip) - vh_correction(ip) 213 end do 214 215 do iter = 1, this%maxcycles 216 217 call poisson_multigrid_cycle(this, der, cor, res) 218 call dderivatives_lapl(der, cor, err) 219 do ip = 1, der%mesh%np 220 err(ip) = res(ip) - err(ip) 221 end do 222 resnorm = dmf_nrm2(der%mesh, err) 223 224 if(resnorm < this%threshold) exit 225 226 if(debug%info) then 227 write(message(1), '(a,i5,a,e13.6)') "Multigrid: base level: iter ", iter, " res ", resnorm 228 call messages_info(1) 229 end if 230 231 end do 232 233 if(resnorm >= this%threshold) then 234 message(1) = 'Multigrid Poisson solver did not converge.' 235 write(message(2), '(a,e14.6)') ' Res = ', resnorm 236 call messages_warning(2) 237 end if 238 239 do ip = 1, der%mesh%np 240 pot(ip) = cor(ip) + vh_correction(ip) 241 end do 242 243 SAFE_DEALLOCATE_A(vh_correction) 244 SAFE_DEALLOCATE_A(res) 245 SAFE_DEALLOCATE_A(cor) 246 SAFE_DEALLOCATE_A(err) 247 248 POP_SUB(poisson_multigrid_solver) 249 end subroutine poisson_multigrid_solver 250 251 ! --------------------------------------------------------- 252 253 recursive subroutine poisson_multigrid_cycle(this, der, pot, rho) 254 type(mg_solver_t), intent(in) :: this 255 type(derivatives_t), intent(in) :: der 256 FLOAT, intent(inout) :: pot(:) 257 FLOAT, intent(in) :: rho(:) 258 259 integer :: ip, iter 260 FLOAT :: resnorm 261 FLOAT, allocatable :: res(:), cres(:), cor(:), ccor(:) 262 263 PUSH_SUB(poisson_multigrid_cycle) 264 265 SAFE_ALLOCATE(res(1:der%mesh%np_part)) 266 267 if(associated(der%coarser)) then 268 SAFE_ALLOCATE(cor(1:der%mesh%np_part)) 269 SAFE_ALLOCATE(cres(1:der%coarser%mesh%np_part)) 270 SAFE_ALLOCATE(ccor(1:der%coarser%mesh%np_part)) 271 272 call multigrid_relax(this, der%mesh, der, pot, rho, this%presteps) 273 274 call dderivatives_lapl(der, pot, res) 275 do ip = 1, der%mesh%np 276 res(ip) = rho(ip) - res(ip) 277 end do 278 279 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, res, cres, this%restriction_method) 280 281 ccor = M_ZERO 282 call poisson_multigrid_cycle(this, der%coarser, ccor, cres) 283 284 cor = M_ZERO 285 call dmultigrid_coarse2fine(der%to_coarser, der%coarser, der%mesh, ccor, cor) 286 287 do ip = 1, der%mesh%np 288 pot(ip) = pot(ip) + cor(ip) 289 end do 290 291 SAFE_DEALLOCATE_A(cor) 292 SAFE_DEALLOCATE_A(cres) 293 SAFE_DEALLOCATE_A(ccor) 294 295 call multigrid_relax(this, der%mesh, der, pot, rho, this%poststeps) 296 297 else 298 299 do iter = 1, this%maxcycles 300 call multigrid_relax(this, der%mesh, der, pot, rho, this%presteps + this%poststeps) 301 call dderivatives_lapl(der, pot, res) 302 do ip = 1, der%mesh%np 303 res(ip) = rho(ip) - res(ip) 304 end do 305 resnorm = dmf_nrm2(der%mesh, res) 306 if(resnorm < this%threshold) exit 307 end do 308 309 end if 310 311 SAFE_DEALLOCATE_A(res) 312 313 POP_SUB(poisson_multigrid_cycle) 314 315 end subroutine poisson_multigrid_cycle 316 317 ! --------------------------------------------------------- 318 319 subroutine multigrid_relax(this, mesh, der, pot, rho, steps) 320 type(mg_solver_t), intent(in) :: this 321 type(mesh_t), intent(in) :: mesh 322 type(derivatives_t), intent(in) :: der 323 FLOAT, intent(inout) :: pot(:) 324 FLOAT, intent(in) :: rho(:) 325 integer, intent(in) :: steps 326 327 integer :: istep 328 integer :: ip, nn 329 FLOAT :: point_lap, factor 330 FLOAT, allocatable :: lpot(:), ldiag(:) 331 type(profile_t), save :: prof 332 333 PUSH_SUB(multigrid_relax) 334 call profiling_in(prof, "MG_GAUSS_SEIDEL") 335 336 select case(this%relaxation_method) 337 338 case(GAUSS_SEIDEL) 339 340 factor = CNST(-1.0)/der%lapl%w(der%lapl%stencil%center, 1)*this%relax_factor 341 342 do istep = 1, steps 343 344 call boundaries_set(der%boundaries, pot) 345 346#ifdef HAVE_MPI 347 if(mesh%parallel_in_domains) call dvec_ghost_update(mesh%vp, pot) 348#endif 349 350 nn = der%lapl%stencil%size 351 352 if(der%lapl%const_w) then 353 call dgauss_seidel(der%lapl%stencil%size, der%lapl%w(1, 1), der%lapl%nri, & 354 der%lapl%ri(1, 1), der%lapl%rimap_inv(1), der%lapl%rimap_inv(2), & 355 factor, pot(1), rho(1)) 356 else 357 do ip = 1, mesh%np 358 point_lap = sum(der%lapl%w(1:nn, ip)*pot(der%lapl%index(1:nn, ip))) 359 pot(ip) = pot(ip) - CNST(0.7)/der%lapl%w(der%lapl%stencil%center, ip)*(point_lap-rho(ip)) 360 end do 361 end if 362 end do 363 call profiling_count_operations(mesh%np*(steps + 1)*(2*nn + 3)) 364 365 case(GAUSS_JACOBI) 366 367 SAFE_ALLOCATE( lpot(1:mesh%np)) 368 SAFE_ALLOCATE(ldiag(1:mesh%np)) 369 370 call derivatives_lapl_diag(der, ldiag) 371 372 do istep = 1, steps 373 call dderivatives_lapl(der, pot, lpot) 374 pot(1:mesh%np) = pot(1:mesh%np) - this%relax_factor/ldiag(1:mesh%np)*(lpot(1:mesh%np) - rho(1:mesh%np)) 375 end do 376 377 SAFE_DEALLOCATE_A(ldiag) 378 SAFE_DEALLOCATE_A(lpot) 379 380 end select 381 382 call profiling_out(prof) 383 POP_SUB(multigrid_relax) 384 385 end subroutine multigrid_relax 386 387end module poisson_multigrid_oct_m 388 389!! Local Variables: 390!! mode: f90 391!! coding: utf-8 392!! End: 393