Octopus
multigrid_solver.F90
Go to the documentation of this file.
1!! Copyright (C) 2005-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch, X. Andrade
2!! Copyright (C) 2024 N. Tancogne-Dejean
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
25 use debug_oct_m
27 use global_oct_m
28 use, intrinsic :: iso_fortran_env
30 use mesh_oct_m
37 use parser_oct_m
40 use space_oct_m
42
43 implicit none
44
45 integer, parameter :: &
46 GAUSS_SEIDEL = 1, &
48
49 private
50 public :: &
57
58 type mg_solver_t
59 private
60
61 real(real64), public :: threshold
62 real(real64) :: relax_factor
63
64 integer, public :: maxcycles = 0
65 integer :: presteps
66 integer :: poststeps
67 integer :: restriction_method
68 integer :: relaxation_method
69 end type mg_solver_t
70
71 integer, parameter, public :: &
72 MG_V_SHAPE = 1, &
73 mg_w_shape = 2, &
74 mg_fmg = 3
75
76contains
77
78 ! ---------------------------------------------------------
79 subroutine multigrid_solver_init(this, namespace, mesh, thr)
80 type(mg_solver_t), intent(out) :: this
81 type(namespace_t), intent(in) :: namespace
82 type(mesh_t), intent(inout) :: mesh
83 real(real64), intent(in) :: thr
84
85 push_sub(multigrid_solver_init)
86
87 this%threshold = thr
88
89 !%Variable MultigridPresmoothingSteps
90 !%Type integer
91 !%Default 1
92 !%Section Hamiltonian::Poisson::Multigrid
93 !%Description
94 !% Number of Gauss-Seidel smoothing steps before coarse-level
95 !% correction in the multigrid solver.
96 !%End
97 call parse_variable(namespace, 'MultigridPresmoothingSteps', 1, this%presteps)
98
99 !%Variable MultigridPostsmoothingSteps
100 !%Type integer
101 !%Default 4
102 !%Section Hamiltonian::Poisson::Multigrid
103 !%Description
104 !% Number of Gauss-Seidel smoothing steps after coarse-level
105 !% correction in the multigrid solver.
106 !%End
107 call parse_variable(namespace, 'MultigridPostsmoothingSteps', 4, this%poststeps)
108
109 !%Variable MultigridMaxCycles
110 !%Type integer
111 !%Default 100
112 !%Section Hamiltonian::Poisson::Multigrid
113 !%Description
114 !% Maximum number of multigrid cycles that are performed if
115 !% convergence is not achieved.
116 !%End
117 call parse_variable(namespace, 'MultigridMaxCycles', 100, this%maxcycles)
119 !%Variable MultigridRestrictionMethod
120 !%Type integer
121 !%Default fullweight
122 !%Section Hamiltonian::Poisson::Multigrid
123 !%Description
124 !% Method used from fine-to-coarse grid transfer.
125 !%Option injection 1
126 !% Injection
127 !%Option fullweight 2
128 !% Fullweight restriction
129 !%End
130 call parse_variable(namespace, 'MultigridRestrictionMethod', 2, this%restriction_method)
131 if (.not. varinfo_valid_option('MultigridRestrictionMethod', this%restriction_method)) then
132 call messages_input_error(namespace, 'MultigridRestrictionMethod')
133 end if
134 call messages_print_var_option("MultigridRestrictionMethod", this%restriction_method, namespace=namespace)
135
136 !%Variable MultigridRelaxationMethod
137 !%Type integer
138 !%Section Hamiltonian::Poisson::Multigrid
139 !%Description
140 !% Method used to solve the linear system approximately in each grid for the
141 !% multigrid procedure that solves a linear equation like the Poisson equation. Default is <tt>gauss_seidel</tt>.
142 !%Option gauss_seidel 1
143 !% Gauss-Seidel.
144 !%Option weighted_jacobi 2
145 !% Jacobi relaxation with a weight. The weight is determined by by MultigridRelaxationFactor.
146 !%End
147 call parse_variable(namespace, 'MultigridRelaxationMethod', gauss_seidel, this%relaxation_method)
148
149 if (.not. varinfo_valid_option('MultigridRelaxationMethod', this%relaxation_method)) then
150 call messages_input_error(namespace, 'MultigridRelaxationMethod')
151 end if
152 call messages_print_var_option("MultigridRelaxationMethod", this%relaxation_method, namespace=namespace)
154 if (this%relaxation_method == weighted_jacobi) then
155 !%Variable MultigridRelaxationFactor
156 !%Type float
157 !%Section Hamiltonian::Poisson::Multigrid
158 !%Description
159 !% Relaxation factor of the relaxation operator used for the
160 !% multigrid method. Only used for the <tt>gauss_jacobi</tt> method.
161 !% The default is 0.6666 for the <tt>gauss_jacobi</tt> method.
162 !%End
163 call parse_variable(namespace, 'MultigridRelaxationFactor', 0.6666_real64, this%relax_factor)
164 end if
165
167 end subroutine multigrid_solver_init
168
169
170 ! ---------------------------------------------------------
174 recursive subroutine multigrid_solver_v_cycle(this, der, op, sol, rhs)
175 type(mg_solver_t), intent(in) :: this
176 type(derivatives_t), intent(in) :: der
177 type(nl_operator_t), intent(in) :: op
178 real(real64), contiguous, intent(inout) :: sol(:)
179 real(real64), contiguous, intent(in) :: rhs(:)
180
181 real(real64), allocatable :: residue(:), coarse_residue(:), correction(:), coarse_correction(:)
182
184
185 safe_allocate(residue(1:der%mesh%np_part))
186
187 if (associated(der%coarser)) then
188 assert(associated(op%coarser))
189
190 safe_allocate(correction(1:der%mesh%np_part))
191 safe_allocate(coarse_residue(1:der%coarser%mesh%np_part))
192 safe_allocate(coarse_correction(1:der%coarser%mesh%np_part))
193
194 ! Pre-Smoothing
195 call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%presteps)
196
197 ! Compute the residual error
198 call get_residual(op, der, sol, rhs, residue)
199
200 ! Restriction of the residual is the next r.h.s
201 message(1) = "Debug: Multigrid restriction"
202 call messages_info(1, debug_only=.true.)
203
204 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, residue, coarse_residue, this%restriction_method)
205
206 ! Recursive call for the coarse-grid correction
207 coarse_correction = m_zero
208 call multigrid_solver_v_cycle(this, der%coarser, op%coarser, coarse_correction, coarse_residue)
209
210 !Prolongation
211 message(1) = "Debug: Multigrid prolongation"
212 call messages_info(1, debug_only=.true.)
213
214 correction = m_zero
215 call dmultigrid_coarse2fine(der%to_coarser, der%coarser, der%mesh, coarse_correction, correction)
216
217 ! Correction
218 call lalg_axpy(der%mesh%np, m_one, correction, sol)
219
220 safe_deallocate_a(correction)
221 safe_deallocate_a(coarse_residue)
222 safe_deallocate_a(coarse_correction)
223
224 ! Post-Smoothing
225 call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%poststeps)
226
227 else ! Coarsest grid
228
229 call multigrid_solver_solve_coarsest(this, der, op, sol, rhs, residue)
230
231 end if
232
233 safe_deallocate_a(residue)
234
236 end subroutine multigrid_solver_v_cycle
237
238 ! ---------------------------------------------------------
242 recursive subroutine multigrid_solver_w_cycle(this, der, op, sol, rhs)
243 type(mg_solver_t), intent(in) :: this
244 type(derivatives_t), intent(in) :: der
245 type(nl_operator_t), intent(in) :: op
246 real(real64), contiguous, intent(inout) :: sol(:)
247 real(real64), contiguous, intent(in) :: rhs(:)
248
249 real(real64), allocatable :: residue(:), coarse_residue(:), correction(:), coarse_correction(:)
250
252
253 safe_allocate(residue(1:der%mesh%np_part))
254
255 if (associated(der%coarser)) then
256 assert(associated(op%coarser))
257
258 safe_allocate(correction(1:der%mesh%np_part))
259 safe_allocate(coarse_residue(1:der%coarser%mesh%np_part))
260 safe_allocate(coarse_correction(1:der%coarser%mesh%np_part))
261
262 ! Pre-Smoothing
263 call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%presteps)
264
265 ! Compute the residual error
266 call get_residual(op, der, sol, rhs, residue)
267
268 ! Restriction of the residual is the next r.h.s
269 message(1) = "Debug: Multigrid restriction"
270 call messages_info(1, debug_only=.true.)
271
272 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, residue, coarse_residue, this%restriction_method)
273
274 ! Recursive call for the coarse-grid correction
275 coarse_correction = m_zero
276 call multigrid_solver_w_cycle(this, der%coarser, op%coarser, coarse_correction, coarse_residue)
277
278 !Prolongation
279 message(1) = "Debug: Multigrid prolongation"
280 call messages_info(1, debug_only=.true.)
281
282 correction = m_zero
283 call dmultigrid_coarse2fine(der%to_coarser, der%coarser, der%mesh, coarse_correction, correction)
284
285 ! Correction
286 call lalg_axpy(der%mesh%np, m_one, correction, sol)
287
288 ! Re-Smoothing
289 call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%presteps)
290
291 ! Compute the residual error
292 call get_residual(op, der, sol, rhs, residue)
293
294 ! Restriction of the residual is the next r.h.s
295 message(1) = "Debug: Multigrid restriction"
296 call messages_info(1, debug_only=.true.)
297
298 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, residue, coarse_residue, this%restriction_method)
299
300 ! Recursive call for the coarse-grid correction
301 coarse_correction = m_zero
302 call multigrid_solver_w_cycle(this, der%coarser, op%coarser, coarse_correction, coarse_residue)
303
304 !Prolongation
305 message(1) = "Debug: Multigrid prolongation"
306 call messages_info(1, debug_only=.true.)
307
308 correction = m_zero
309 call dmultigrid_coarse2fine(der%to_coarser, der%coarser, der%mesh, coarse_correction, correction)
310
311 ! Correction
312 call lalg_axpy(der%mesh%np, m_one, correction, sol)
313
314 safe_deallocate_a(correction)
315 safe_deallocate_a(coarse_residue)
316 safe_deallocate_a(coarse_correction)
317
318 ! Post-Smoothing
319 call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%poststeps)
320
321 else ! Coarsest grid
322
323 call multigrid_solver_solve_coarsest(this, der, op, sol, rhs, residue)
324
325 end if
326
327 safe_deallocate_a(residue)
328
330
331 end subroutine multigrid_solver_w_cycle
332
333 ! ---------------------------------------------------------
337 subroutine multigrid_iterative_solver(this, namespace, der, op, sol, rhs, multigrid_shape)
338 type(mg_solver_t), intent(in) :: this
339 type(namespace_t), intent(in) :: namespace
340 type(derivatives_t), intent(in) :: der
341 type(nl_operator_t), intent(in) :: op
342 real(real64), contiguous, intent(inout) :: sol(:)
343 real(real64), contiguous, intent(inout) :: rhs(:)
344 integer, intent(in) :: multigrid_shape
345
346 integer :: iter
347 real(real64) :: resnorm
348 real(real64), allocatable :: err(:)
349
351
352 safe_allocate(err(1:der%mesh%np))
353
354 do iter = 1, this%maxcycles
355
356 select case (multigrid_shape)
357 case(mg_v_shape)
358 call multigrid_solver_v_cycle(this, der, op, sol, rhs)
359 case(mg_w_shape)
360 call multigrid_solver_w_cycle(this, der, op, sol, rhs)
361 case default
362 assert(.false.)
363 end select
364 ! Compute the residual
365 call dderivatives_lapl(der, sol, err)
366 call lalg_axpy(der%mesh%np, -m_one, rhs, err)
367 resnorm = dmf_nrm2(der%mesh, err)
368
369 if (resnorm < this%threshold) exit
370
371 write(message(1), '(a,i5,a,e13.6)') "Multigrid: base level: iter ", iter, " res ", resnorm
372 call messages_info(1, namespace=namespace, debug_only=.true.)
373
374 end do
375
376 if (resnorm >= this%threshold) then
377 message(1) = 'Multigrid Poisson solver did not converge.'
378 write(message(2), '(a,e14.6)') ' Abs. norm of the residue = ', resnorm
379 call messages_warning(2, namespace=namespace)
380 else
381 write(message(1), '(a,i4,a)') "Multigrid Poisson solver converged in ", iter, " iterations."
382 write(message(2), '(a,e14.6)') ' Abs. norm of the residue = ', resnorm
383 call messages_info(2, namespace=namespace, debug_only=.true.)
384 end if
385
386 safe_deallocate_a(err)
387
389 end subroutine multigrid_iterative_solver
390
391 ! ---------------------------------------------------------
396 recursive subroutine multigrid_fmg_solver(this, namespace, der, op, sol, rhs)
397 type(mg_solver_t), intent(in) :: this
398 type(namespace_t), intent(in) :: namespace
399 type(derivatives_t), intent(in) :: der
400 type(nl_operator_t), intent(in) :: op
401 real(real64), contiguous, intent(inout) :: sol(:)
402 real(real64), contiguous, intent(inout) :: rhs(:)
403
404 real(real64), allocatable :: coarse_solution(:), coarse_rhs(:), residue(:)
405
406 push_sub(multigrid_fmg_solver)
407
408 if (associated(der%coarser)) then
409 assert(associated(op%coarser))
410
411 safe_allocate(coarse_rhs(1:der%coarser%mesh%np_part))
412 safe_allocate(coarse_solution(1:der%coarser%mesh%np_part))
413
414 ! Restriction of the r.h.s
415 message(1) = "Debug: Full Multigrid restriction"
416 call messages_info(1, debug_only=.true.)
417
418 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, rhs, coarse_rhs, this%restriction_method)
419
420 ! Recursive call for the coarse-grid correction
421 coarse_solution = m_zero
422 call multigrid_fmg_solver(this, namespace, der%coarser, op%coarser, coarse_solution, coarse_rhs)
423
424 !Prolongation
425 message(1) = "Debug: Full Multigrid prolongation"
426 call messages_info(1, debug_only=.true.)
427
428 sol = m_zero
429 call dmultigrid_coarse2fine(der%to_coarser, der%coarser, der%mesh, coarse_solution, sol)
430
431 ! Perform N times a V cycle, up to convergence at each step
432 call multigrid_iterative_solver(this, namespace, der, op, sol, rhs, mg_v_shape)
433
434 safe_deallocate_a(coarse_solution)
435 safe_deallocate_a(coarse_rhs)
436
437 else ! Coarsest grid - solve the problem
438
439 safe_allocate(residue(1:der%mesh%np))
440 call multigrid_solver_solve_coarsest(this, der, op, sol, rhs, residue)
441 safe_deallocate_a(residue)
442 end if
443
444 pop_sub(multigrid_fmg_solver)
445 end subroutine multigrid_fmg_solver
446
447
448 ! ---------------------------------------------------------
450 subroutine multigrid_solver_solve_coarsest(this, der, op, sol, rhs, residue)
451 type(mg_solver_t), intent(in) :: this
452 type(derivatives_t), intent(in) :: der
453 type(nl_operator_t), intent(in) :: op
454 real(real64), contiguous, intent(inout) :: sol(:)
455 real(real64), contiguous, intent(in) :: rhs(:)
456 real(real64), contiguous, intent(inout) :: residue(:)
457
458 integer :: iter
459 real(real64) :: resnorm
460
462
463 ! Solution of the problem, i.e., multiple call to multigrid_relax up to convergence
464 do iter = 1, this%maxcycles
465
466 call multigrid_relax(this, der%mesh, der, op, sol, rhs, 1)
467
468 call get_residual(op, der, sol, rhs, residue)
469 resnorm = dmf_nrm2(der%mesh, residue)
470 if (resnorm < this%threshold) exit
471
472 end do
473
474 write(message(1), '(a,i4,a)') "Debug: Multigrid coarsest grid solver converged in ", iter, " iterations."
475 write(message(2), '(a,es18.6)') " Residue norm is ", resnorm
476 call messages_info(2, debug_only=.true.)
477
479 end subroutine
480
481 ! ---------------------------------------------------------
483 subroutine get_residual(op, der, sol, rhs, residue)
484 type(nl_operator_t), intent(in) :: op
485 type(derivatives_t), intent(in) :: der
486 real(real64), contiguous, intent(inout) :: sol(:)
487 real(real64), contiguous, intent(in) :: rhs(:)
488 real(real64), contiguous, intent(inout) :: residue(:)
489
490 integer :: ip
492 ! Compute the residue
493 call dderivatives_perform(op, der, sol, residue)
494 !$omp parallel do
495 do ip = 1, der%mesh%np
496 residue(ip) = rhs(ip) - residue(ip)
497 end do
498 end subroutine get_residual
499
500
501 ! ---------------------------------------------------------
505 subroutine multigrid_relax(this, mesh, der, op, sol, rhs, steps)
506 type(mg_solver_t), intent(in) :: this
507 type(mesh_t), intent(in) :: mesh
508 type(derivatives_t), intent(in) :: der
509 type(nl_operator_t), intent(in) :: op
510 real(real64), contiguous, intent(inout) :: sol(:)
511 real(real64), contiguous, intent(in) :: rhs(:)
512 integer, intent(in) :: steps
513
514 integer :: istep, index
515 integer :: ip, nn, is
516 real(real64) :: point, factor
517 real(real64), allocatable :: op_sol(:), diag(:)
518
519 push_sub(multigrid_relax)
520 call profiling_in("MG_GAUSS_SEIDEL")
521
522 select case (this%relaxation_method)
523
524 case (gauss_seidel)
525
526 do istep = 1, steps
527
528 call boundaries_set(der%boundaries, der%mesh, sol)
529
530 if (mesh%parallel_in_domains) then
531 call dpar_vec_ghost_update(mesh%pv, sol)
532 end if
533
534 nn = op%stencil%size
535
536 if (op%const_w) then
537 factor = -m_one/op%w(op%stencil%center, 1)
538 call dgauss_seidel(op%stencil%size, op%w(1, 1), op%nri, &
539 op%ri(1, 1), op%rimap_inv(1), op%rimap_inv(2), factor, sol(1), rhs(1))
540 else
541 !$omp parallel do private(point, index)
542 do ip = 1, mesh%np
543 point = m_zero
544 do is = 1, nn
545 index = nl_operator_get_index(op, is, ip)
546 point = point + op%w(is, ip)*sol(index)
547 end do
548 sol(ip) = sol(ip) - (point-rhs(ip))/op%w(op%stencil%center, ip)
549 end do
550 end if
551
552 end do
553 call profiling_count_operations(mesh%np*(steps + 1)*(2*nn + 3))
554
555 case (weighted_jacobi)
556
557 safe_allocate(op_sol(1:mesh%np))
558 safe_allocate(diag(1:mesh%np))
559
560 call dnl_operator_operate_diag(op, diag)
561 !$omp parallel do
562 do ip = 1, mesh%np
563 diag(ip) = this%relax_factor/diag(ip)
564 end do
565
566 do istep = 1, steps
567 call dderivatives_perform(op, der, sol, op_sol)
568 !$omp parallel do
569 do ip = 1, mesh%np
570 sol(ip) = sol(ip) - diag(ip)*(op_sol(ip) - rhs(ip))
571 end do
572 end do
573
574 safe_deallocate_a(diag)
575 safe_deallocate_a(op_sol)
576
577 end select
579 call profiling_out("MG_GAUSS_SEIDEL")
580 pop_sub(multigrid_relax)
581
582 end subroutine multigrid_relax
583
584end module multigrid_solver_oct_m
585
586!! Local Variables:
587!! mode: f90
588!! coding: utf-8
589!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:124
subroutine, public dpar_vec_ghost_update(pv, v_local)
Updates ghost points of every node.
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_perform(op, der, ff, op_ff, ghost_update, set_bc, factor)
apply a nl_operator to a mesh function
subroutine, public dderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_one
Definition: global.F90:192
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
subroutine, public dmultigrid_coarse2fine(tt, coarse_der, fine_mesh, f_coarse, f_fine, set_bc)
Definition: multigrid.F90:688
subroutine, public dmultigrid_fine2coarse(tt, fine_der, coarse_mesh, f_fine, f_coarse, method_p)
Definition: multigrid.F90:767
This modules provides the routines for solving Ax=b using the V-shaped multigrid method.
recursive subroutine, public multigrid_fmg_solver(this, namespace, der, op, sol, rhs)
Full multigrid (FMG) solver.
recursive subroutine, public multigrid_solver_v_cycle(this, der, op, sol, rhs)
Performs one cycle of a V-shaped multigrid solver.
subroutine, public multigrid_solver_init(this, namespace, mesh, thr)
recursive subroutine, public multigrid_solver_w_cycle(this, der, op, sol, rhs)
Performs one cycle of a W-shaped multigrid solver.
subroutine multigrid_solver_solve_coarsest(this, der, op, sol, rhs, residue)
Computes the solution on the coarsest grid.
subroutine, public multigrid_iterative_solver(this, namespace, der, op, sol, rhs, multigrid_shape)
An iterative multigrid solver.
integer, parameter weighted_jacobi
integer, parameter, public mg_fmg
subroutine get_residual(op, der, sol, rhs, residue)
Computes the residual.
subroutine multigrid_relax(this, mesh, der, op, sol, rhs, steps)
Given a nonlocal operator op, perform the relaxation operator.
integer, parameter, public mg_w_shape
This module defines non-local operators.
subroutine, public dnl_operator_operate_diag(op, fo)
integer pure function, public nl_operator_get_index(op, is, ip)
This module contains interfaces for routines in operate.c.
Definition: operate_f.F90:119
Some general things and nomenclature:
Definition: par_vec.F90:173
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
class representing derivatives
Describes mesh distribution to nodes.
Definition: mesh.F90:187
data type for non local operators
int true(void)