Octopus
helmholtz_decomposition.F90
Go to the documentation of this file.
1!! Copyright (C) 2019 R. Jestaedt, H. Appel, F. Bonafe, M. Oliveira, N. Tancogne-Dejean
2!! Copyright (C) 2022-2023 F. Troisi
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
61
62! TODO: Issue 705 (ftroisi): Add a test for non Coulomb Gauge
64 use blas_oct_m
65 use box_oct_m
67 use debug_oct_m
69 use global_oct_m
70 use grid_oct_m
71 use io_oct_m
73 use, intrinsic :: iso_fortran_env
75 use loct_oct_m
76 use math_oct_m
77 use mesh_oct_m
82 use parser_oct_m
85 use space_oct_m
86 use unit_oct_m
88
89 implicit none
90
91 private
93
95 private
96 ! Helmholtz solution handling
97 type(grid_t), pointer :: sys_grid
98 logical, public :: compute_surface_correction
99 integer, allocatable :: surface_points(:)
100 integer, allocatable :: inner_stencil(:)
101 ! Coulomb Gauge
102 logical :: enforce_coulomb_gauge
104 real(real64) :: coulomb_gauge_tolerance
105 ! !! The code will evaluate the following expression:
106 ! !! $spacing * \frac{\nabla \cdot A}{|A|}$
107 ! Poisson Solver
108 real(real64) :: poisson_prefactor
109 type(poisson_t) :: poisson_solver
110
111 contains
112 ! Init
113 procedure :: init => helmholtz_decomposition_init
114 ! Apply mask
115 procedure :: dapply_inner_stencil_mask, zapply_inner_stencil_mask
116 generic :: apply_inner_stencil_mask => dapply_inner_stencil_mask, zapply_inner_stencil_mask
117 ! Vector potential
118 procedure :: dget_vector_potential, zget_vector_potential
119 generic :: get_vector_potential => dget_vector_potential, zget_vector_potential
120 ! Scalar potential
121 procedure :: dget_scalar_potential, zget_scalar_potential
122 generic :: get_scalar_potential => dget_scalar_potential, zget_scalar_potential
123 ! Transverse field
124 procedure :: dget_trans_field, zget_trans_field
125 generic :: get_trans_field => dget_trans_field, zget_trans_field
126 ! Longitudinal field
127 procedure :: dget_long_field, zget_long_field
128 generic :: get_long_field => dget_long_field, zget_long_field
129 ! Finalize
130 final :: helmholtz_finalize
132
133contains
134
136 subroutine helmholtz_decomposition_init(this, namespace, sys_grid, system_mc, space)
137 class(helmholtz_decomposition_t), intent(inout) :: this
138 type(namespace_t), intent(in) :: namespace
139 type(grid_t), target, intent(in) :: sys_grid
140 type(multicomm_t), intent(in) :: system_mc
141 class(space_t), intent(in) :: space
142
143 logical, allocatable :: mask(:)
144 logical :: visualize_boxes
145 integer :: number_of_points
146
148
149 ! Allocate pointer to system grid
150 this%sys_grid => sys_grid
151
152 ! 1. Initialize poisson solver
153 call poisson_init(this%poisson_solver, namespace, space, sys_grid%der, system_mc, sys_grid%stencil, label = "Helmholtz")
154
155 ! 2. Create a mask for the points contained in a layer of the width of the stencil between the innermost layer of points and
156 ! last row of points in the grid.
157 safe_allocate(mask(1:sys_grid%np))
158 ! Get the mask
159 mask = derivatives_get_inner_boundary_mask(sys_grid%der)
160 number_of_points = count(mask)
161 ! Store the indices of the mask points
162 safe_allocate(this%inner_stencil(1:number_of_points))
163 call get_indices_from_mask(sys_grid%np, mask, this%inner_stencil)
164
165 ! 3. Coulomb Gauge
166 !%Variable HelmholtzEnforceCoulombGauge
167 !%Type logical
168 !%Default yes
169 !%Section Calculation Modes::Test
170 !%Description
171 !% If true, the Vector Potential is enforced to be in Coulomb Gauge.
172 !% See proof of eq 30 of: https://www.scirp.org/pdf/jmp_2016053115275279.pdf
173 !%End
174 call parse_variable(namespace, 'HelmholtzEnforceCoulombGauge', .true., this%enforce_coulomb_gauge)
175 if (this%enforce_coulomb_gauge) then
176 !%Variable HelmholtzCoulombGaugeTolerance
177 !%Type float
178 !%Default 1e-5
179 !%Section Calculation Modes::Test
180 !%Description
181 !% If the vector potential should enforced to be in Coulomb Gauge, this variable defines the tolerance. The code will
182 !% evaluate the following expression: $spacing * \frac{\nabla \cdot A}{|A|}$
183 !%End
184 call parse_variable(namespace, 'HelmholtzCoulombGaugeTolerance', 1e-5_real64, this%coulomb_gauge_tolerance)
185 end if
186
187 ! 3. Surface correction TODO: See Issue 705 (@ftroisi)
188
189 !%Variable SurfaceCorrection
190 !%Type logical
191 !%Default no
192 !%Section Calculation Modes::Test
193 !%Description
194 !% Compute the surface correction for Helmholtz decomposition?
195 !%End
196 call parse_variable(namespace, 'SurfaceCorrection', .false., this%compute_surface_correction)
197 if (this%compute_surface_correction) then
198 call messages_not_implemented("Surface correction for Helmholtz decomposition")
199 end if
200
201 if (this%compute_surface_correction) then
202 ! First of all we have to get the surface points
203 mask = sys_grid%box%get_surface_points(namespace, sys_grid%spacing, sys_grid%np, sys_grid%x)
204 number_of_points = count(mask)
205
206 ! Allocate array surface_points and assign those points
207 safe_allocate(this%surface_points(1:number_of_points))
208 call get_indices_from_mask(sys_grid%np, mask, this%surface_points)
210 ! Initialize the poisson prefactor to treat the singular points (for the surface correction)
211 select case (sys_grid%box%dim)
212 case (3)
213 this%poisson_prefactor = m_two * m_pi * (m_three / (m_pi * m_four))**(m_twothird)
214 case (2)
215 this%poisson_prefactor = m_two * sqrt(m_pi)
216 case (1)
217 this%poisson_prefactor = m_one
218 case default
219 message(1) = "Internal error: surface correction for helmholtz decomposition can only be called for 1D, 2D or 3D."
220 call messages_fatal(1, namespace = namespace)
221 end select
222 end if
223 safe_deallocate_a(mask)
224
225 ! 4. Optionally, the user may want to visualize the relevant regions used to compute the Helmholtz decomposition
226 !%Variable HelmholtzVisualizeBoxes
227 !%Type logical
228 !%Default no
229 !%Section Calculation Modes::Test
230 !%Description
231 !% If true, output the volume points for the three boxes of the Helmholtz surface correction.
232 !% 1) The volume points of the system box
233 !% 2) The inner mask for the system box. This region has the thickness of the stencil and it is used to set to zero
234 !% the longitudinal or transverse field after computing the final divergence or curl (to avoid spikes)
235 !% 3) The surface points of the system box
236 !%
237 !%End
238 call parse_variable(namespace, 'HelmholtzVisualizeBoxes', .false., visualize_boxes)
239 if (visualize_boxes) then
240 call helmholtz_visualize_boxes(this, namespace, space)
241 end if
242
244 end subroutine helmholtz_decomposition_init
245
246 subroutine get_indices_from_mask(np, mask, indices)
247 integer, intent(in) :: np
248 logical, intent(in) :: mask(:)
249 integer, intent(out) :: indices(:)
250
251 integer :: ip, j
252 push_sub(get_indices_from_mask)
253
254 j = 1
255 do ip = 1, np
256 if (mask(ip)) then
257 indices(j) = ip
258 j = j + 1
259 end if
260 end do
261
262 pop_sub(get_indices_from_mask)
263 end subroutine get_indices_from_mask
264
265 subroutine helmholtz_finalize(this)
266 type(helmholtz_decomposition_t), intent(inout) :: this
267 push_sub(helmholtz_finalize)
268
269 call poisson_end(this%poisson_solver)
270
271 safe_deallocate_a(this%inner_stencil)
272 safe_deallocate_a(this%surface_points)
273
274 pop_sub(helmholtz_finalize)
275 end subroutine helmholtz_finalize
276
278 subroutine helmholtz_visualize_boxes(this, namespace, space)
279 type(helmholtz_decomposition_t), intent(inout) :: this
280 type(namespace_t), intent(in) :: namespace
281 class(space_t), intent(in) :: space
282
283 real(real64) :: coords(space%dim), normal_vector(space%dim), surface_element, ra
284 integer :: ip, ii, iunit
285
287
288 ! SYSTEM BOX - Volume points
289 iunit = io_open("system_box_volume_points", namespace, action='write')
290 write(iunit, '(a6, 5X, a1, 12X, a1, 12X, a1)') "#point", "x", "y", "z"
291 do ip = 1, this%sys_grid%np
292 call mesh_r(this%sys_grid, ip, ra, coords = coords)
293 write(iunit, '(i7, 3f12.7)') ip, coords(:)
294 end do
295 call io_close(iunit)
296
297 ! SYSTEM BOX - Mask for inner stencil
298 iunit = io_open("system_box_inner_stencil", namespace, action='write')
299 write(iunit, '(a6, 5X, a1, 12X, a1, 12X, a1)') "#point", "x", "y", "z"
300 do ip = 1, SIZE(this%inner_stencil, dim = 1)
301 ii = this%inner_stencil(ip)
302 call mesh_r(this%sys_grid, ii, ra, coords = coords)
303 write(iunit, '(i7, 3f12.7)') ip, coords(:)
304 end do
305 call io_close(iunit)
306
307 ! SYSTEM BOX - Surface Points
308 if (this%compute_surface_correction) then
309 iunit = io_open("system_box_surface_points", namespace, action='write')
310 write(iunit, '(a,12X,a,12X,a,12X,a,12X,a,12X,a,12X,a)') "x", "y", "z", "norm_x", "norm_y", "norm_z", "surf_element"
311 do ip = 1, SIZE(this%surface_points, dim = 1)
312 ii = this%surface_points(ip)
313 call mesh_r(this%sys_grid, ii, ra, coords = coords)
314 call this%sys_grid%box%get_surface_point_info(coords, this%sys_grid%spacing, normal_vector, surface_element)
315 write(iunit, '(7f12.7)') coords(:), normal_vector(:), surface_element
316 end do
317 call io_close(iunit)
318 end if
319
321 end subroutine helmholtz_visualize_boxes
322
323#include "undef.F90"
324#include "complex.F90"
325#include "helmholtz_decomposition_inc.F90"
326
327#include "undef.F90"
328#include "real.F90"
329#include "helmholtz_decomposition_inc.F90"
330
332
333!! Local Variables:
334!! mode: f90
335!! coding: utf-8
336!! End:
This module contains interfaces for BLAS routines You should not use these routines directly....
Definition: blas.F90:118
Module, implementing a factory for boxes.
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
logical function, dimension(1:this%mesh%np), public derivatives_get_inner_boundary_mask(this)
This function tells whether a point in the grid is contained in a layer of the width of the stencil b...
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
real(real64), parameter, public m_twothird
Definition: global.F90:196
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:117
The Helmholtz decomposition is intended to contain "only mathematical" functions and procedures to co...
subroutine get_indices_from_mask(np, mask, indices)
subroutine helmholtz_decomposition_init(this, namespace, sys_grid, system_mc, space)
Initialize Helmholtz decomposition object.
subroutine helmholtz_visualize_boxes(this, namespace, space)
Visualise boxes for use in Helmholtz Decomposition.
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:336
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
subroutine, public poisson_init(this, namespace, space, der, mc, stencil, qtot, label, solver, verbose, force_serial, force_cmplx)
Definition: poisson.F90:233
subroutine, public poisson_end(this)
Definition: poisson.F90:709
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
This module defines the unit system, used for input and output.
int true(void)