Octopus
poisson_cg.F90
Go to the documentation of this file.
1!! Cropyright (C) 2004-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
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
22 use debug_oct_m
24 use global_oct_m
25 use, intrinsic :: iso_c_binding, only: c_ptr
26 use, intrinsic :: iso_fortran_env
28 use mesh_oct_m
35
36 implicit none
37
38 private
39 public :: &
44
45 real(real64), public :: threshold
46 integer, public :: maxiter
47
48contains
49
50
51 ! ---------------------------------------------------------
52 subroutine poisson_cg_init(thr, itr)
53 integer, intent(in) :: itr
54 real(real64), intent(in) :: thr
55
56 push_sub(poisson_cg_init)
57 threshold = thr
58 maxiter = itr
59 pop_sub(poisson_cg_init)
60 end subroutine poisson_cg_init
61
62
63 ! ---------------------------------------------------------
64 subroutine poisson_cg_end
65
66 end subroutine poisson_cg_end
67
68
69 ! ---------------------------------------------------------
70 subroutine poisson_cg1(namespace, der, corrector, pot, rho)
71 type(namespace_t), intent(in) :: namespace
72 type(derivatives_t), target, intent(in) :: der
73 type(poisson_corr_t), intent(in) :: corrector
74 real(real64), intent(inout) :: pot(:)
75 real(real64), intent(in) :: rho(:)
76
77 integer :: iter
78 real(real64) :: res
79 real(real64), allocatable :: wk(:), lwk(:), zk(:), pk(:)
80 type(c_ptr) :: userdata(0) ! workaround, literal [c_ptr::] crashes ifort, eventually encapsulate der_pointer and mesh_pointer
81
82 push_sub(poisson_cg1)
83
84 safe_allocate( wk(1:der%mesh%np_part))
85 safe_allocate(lwk(1:der%mesh%np_part))
86 safe_allocate( zk(1:der%mesh%np_part))
87 safe_allocate( pk(1:der%mesh%np_part))
88
89 ! build initial guess for the potential
90 wk(1:der%mesh%np) = pot(1:der%mesh%np)
91 call poisson_boundary_conditions(corrector, der%mesh, rho, wk)
92 call dderivatives_lapl(der, wk, lwk, ghost_update = .true., set_bc = .false.)
93
94 ! Solving -\Delta = zk, as CG require SDP operator
95 zk(1:der%mesh%np) = (m_four*m_pi*rho(1:der%mesh%np) + lwk(1:der%mesh%np))
96 safe_deallocate_a(wk)
97 safe_deallocate_a(lwk) ! they are no longer needed
98
99 der_pointer => der
100 mesh_pointer => der%mesh
101 call lalg_copy(der%mesh%np, zk, pk)
102 pk(der%mesh%np + 1:) = m_zero
103 iter = maxiter
104 call dconjugate_gradients(der%mesh%np, pk, zk, internal_laplacian_op, internal_dotp, iter, res, threshold, userdata)
105 if (res >= threshold) then
106 message(1) = 'Conjugate-gradients Poisson solver did not converge.'
107 write(message(2), '(a,i8)') ' Iter = ',iter
108 write(message(3), '(a,e14.6)') ' Res = ', res
109 call messages_warning(3, namespace=namespace)
110 end if
111 nullify(der_pointer, mesh_pointer)
112 pot(1:der%mesh%np) = pot(1:der%mesh%np) + pk(1:der%mesh%np)
113
114 safe_deallocate_a(zk)
115 safe_deallocate_a(pk)
116 pop_sub(poisson_cg1)
117 end subroutine poisson_cg1
118
119
120 ! ---------------------------------------------------------
121 subroutine poisson_cg2(namespace, der, pot, rho)
122 type(namespace_t), intent(in) :: namespace
123 type(derivatives_t), target, intent(in) :: der
124 real(real64), contiguous, intent(inout) :: pot(:)
125 real(real64), intent(in) :: rho(:)
126
127 integer :: iter, ip
128 real(real64), allocatable :: potc(:), rhs(:)
129 real(real64) :: res
130 type(c_ptr) :: userdata(0) ! workaround, literal [c_ptr::] crashes ifort, eventually encapsulate der_pointer and mesh_pointer
131
132 push_sub(poisson_cg2)
133
134 iter = maxiter
135 der_pointer => der
136 mesh_pointer => der%mesh
137
138 safe_allocate(rhs(1:der%mesh%np))
139 safe_allocate(potc(1:der%mesh%np_part))
141 ! Solving -\Delta = rhs, as CG require SDP operator
142 do ip = 1, der%mesh%np
143 rhs(ip) = 4.0_real64*m_pi*rho(ip)
144 end do
145 call lalg_copy(der%mesh%np, pot, potc)
146
147 call dconjugate_gradients(der%mesh%np, potc, rhs, internal_laplacian_op, internal_dotp, iter, res, threshold, userdata)
148
149 if (res >= threshold) then
150 message(1) = 'Conjugate-gradients Poisson solver did not converge.'
151 write(message(2), '(a,i8)') ' Iter = ', iter
152 write(message(3), '(a,e14.6)') ' Res = ', res
153 call messages_warning(3, namespace=namespace)
154 end if
155
156 call lalg_copy(der%mesh%np, potc, pot)
157
158 nullify(der_pointer, mesh_pointer)
160 safe_deallocate_a(rhs)
161 safe_deallocate_a(potc)
162
163 pop_sub(poisson_cg2)
164 end subroutine poisson_cg2
166end module poisson_cg_oct_m
167
168!! Local Variables:
169!! mode: f90
170!! coding: utf-8
171!! End:
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
This module calculates the derivatives (gradients, Laplacians, etc.) of a 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:190
real(real64), parameter, public m_four
Definition: global.F90:194
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:188
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:531
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public poisson_cg2(namespace, der, pot, rho)
Definition: poisson_cg.F90:217
subroutine, public poisson_cg1(namespace, der, corrector, pot, rho)
Definition: poisson_cg.F90:166
subroutine, public poisson_cg_init(thr, itr)
Definition: poisson_cg.F90:148
subroutine, public poisson_cg_end
Definition: poisson_cg.F90:160
subroutine, public internal_laplacian_op(xx, yy, userdata)
real(real64) function, public internal_dotp(xx, yy)
subroutine, public poisson_boundary_conditions(this, mesh, rho, pot)
type(mesh_t), pointer, public mesh_pointer
type(derivatives_t), pointer, public der_pointer
This module is intended to contain "only mathematical" functions and procedures.
Definition: solvers.F90:117
int true(void)