28 use,
intrinsic :: iso_fortran_env
45 integer,
parameter :: &
61 real(real64),
public :: threshold
62 real(real64) :: relax_factor
64 integer,
public :: maxcycles = 0
67 integer :: restriction_method
68 integer :: relaxation_method
71 integer,
parameter,
public :: &
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
97 call parse_variable(namespace,
'MultigridPresmoothingSteps', 1, this%presteps)
107 call parse_variable(namespace,
'MultigridPostsmoothingSteps', 4, this%poststeps)
117 call parse_variable(namespace,
'MultigridMaxCycles', 100, this%maxcycles)
130 call parse_variable(namespace,
'MultigridRestrictionMethod', 2, this%restriction_method)
147 call parse_variable(namespace,
'MultigridRelaxationMethod', gauss_seidel, this%relaxation_method)
163 call parse_variable(namespace,
'MultigridRelaxationFactor', 0.6666_real64, this%relax_factor)
178 real(real64),
contiguous,
intent(inout) :: sol(:)
179 real(real64),
contiguous,
intent(in) :: rhs(:)
181 real(real64),
allocatable :: residue(:), coarse_residue(:), correction(:), coarse_correction(:)
185 safe_allocate(residue(1:der%mesh%np_part))
187 if (
associated(der%coarser))
then
188 assert(
associated(op%coarser))
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))
201 message(1) =
"Debug: Multigrid restriction"
204 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, residue, coarse_residue, this%restriction_method)
207 coarse_correction =
m_zero
211 message(1) =
"Debug: Multigrid prolongation"
220 safe_deallocate_a(correction)
221 safe_deallocate_a(coarse_residue)
222 safe_deallocate_a(coarse_correction)
225 call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%poststeps)
233 safe_deallocate_a(residue)
246 real(real64),
contiguous,
intent(inout) :: sol(:)
247 real(real64),
contiguous,
intent(in) :: rhs(:)
249 real(real64),
allocatable :: residue(:), coarse_residue(:), correction(:), coarse_correction(:)
253 safe_allocate(residue(1:der%mesh%np_part))
255 if (
associated(der%coarser))
then
256 assert(
associated(op%coarser))
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))
272 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, residue, coarse_residue, this%restriction_method)
275 coarse_correction =
m_zero
279 message(1) =
"Debug: Multigrid prolongation"
295 message(1) =
"Debug: Multigrid restriction"
298 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, residue, coarse_residue, this%restriction_method)
301 coarse_correction =
m_zero
305 message(1) =
"Debug: Multigrid prolongation"
314 safe_deallocate_a(correction)
315 safe_deallocate_a(coarse_residue)
316 safe_deallocate_a(coarse_correction)
319 call multigrid_relax(this, der%mesh, der, op, sol, rhs, this%poststeps)
327 safe_deallocate_a(residue)
342 real(real64),
contiguous,
intent(inout) :: sol(:)
343 real(real64),
contiguous,
intent(inout) :: rhs(:)
344 integer,
intent(in) :: multigrid_shape
347 real(real64) :: resnorm
348 real(real64),
allocatable :: err(:)
352 safe_allocate(err(1:der%mesh%np))
354 do iter = 1, this%maxcycles
356 select case (multigrid_shape)
369 if (resnorm < this%threshold)
exit
371 write(
message(1),
'(a,i5,a,e13.6)')
"Multigrid: base level: iter ", iter,
" res ", resnorm
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
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
386 safe_deallocate_a(err)
401 real(real64),
contiguous,
intent(inout) :: sol(:)
402 real(real64),
contiguous,
intent(inout) :: rhs(:)
404 real(real64),
allocatable :: coarse_solution(:), coarse_rhs(:), residue(:)
408 if (
associated(der%coarser))
then
409 assert(
associated(op%coarser))
411 safe_allocate(coarse_rhs(1:der%coarser%mesh%np_part))
412 safe_allocate(coarse_solution(1:der%coarser%mesh%np_part))
415 message(1) =
"Debug: Full Multigrid restriction"
418 call dmultigrid_fine2coarse(der%to_coarser, der, der%coarser%mesh, rhs, coarse_rhs, this%restriction_method)
422 call multigrid_fmg_solver(this, namespace, der%coarser, op%coarser, coarse_solution, coarse_rhs)
425 message(1) =
"Debug: Full Multigrid prolongation"
434 safe_deallocate_a(coarse_solution)
435 safe_deallocate_a(coarse_rhs)
439 safe_allocate(residue(1:der%mesh%np))
441 safe_deallocate_a(residue)
454 real(real64),
contiguous,
intent(inout) :: sol(:)
455 real(real64),
contiguous,
intent(in) :: rhs(:)
456 real(real64),
contiguous,
intent(inout) :: residue(:)
459 real(real64) :: resnorm
464 do iter = 1, this%maxcycles
469 resnorm =
dmf_nrm2(der%mesh, residue)
470 if (resnorm < this%threshold)
exit
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
486 real(real64),
contiguous,
intent(inout) :: sol(:)
487 real(real64),
contiguous,
intent(in) :: rhs(:)
488 real(real64),
contiguous,
intent(inout) :: residue(:)
495 do ip = 1, der%mesh%np
496 residue(ip) = rhs(ip) - residue(ip)
507 type(
mesh_t),
intent(in) :: mesh
510 real(real64),
contiguous,
intent(inout) :: sol(:)
511 real(real64),
contiguous,
intent(in) :: rhs(:)
512 integer,
intent(in) :: steps
514 integer :: istep, index
515 integer :: ip, nn, is
516 real(real64) :: point, factor
517 real(real64),
allocatable :: op_sol(:), diag(:)
522 select case (this%relaxation_method)
530 if (mesh%parallel_in_domains)
then
537 factor = -
m_one/op%w(op%stencil%center, 1)
539 op%ri(1, 1), op%rimap_inv(1), op%rimap_inv(2), factor, sol(1), rhs(1))
546 point = point + op%w(is, ip)*sol(index)
548 sol(ip) = sol(ip) - (point-rhs(ip))/op%w(op%stencil%center, ip)
557 safe_allocate(op_sol(1:mesh%np))
558 safe_allocate(diag(1:mesh%np))
563 diag(ip) = this%relax_factor/diag(ip)
570 sol(ip) = sol(ip) - diag(ip)*(op_sol(ip) - rhs(ip))
574 safe_deallocate_a(diag)
575 safe_deallocate_a(op_sol)
constant times a vector plus a vector
Module implementing boundary conditions in Octopus.
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
real(real64), parameter, public m_one
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_warning(no_lines, all_nodes, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
subroutine, public dmultigrid_coarse2fine(tt, coarse_der, fine_mesh, f_coarse, f_fine, set_bc)
subroutine, public dmultigrid_fine2coarse(tt, fine_der, coarse_mesh, f_fine, f_coarse, method_p)
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.
Some general things and nomenclature:
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
class representing derivatives
Describes mesh distribution to nodes.
data type for non local operators