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