73 use,
intrinsic :: iso_fortran_env
97 type(grid_t),
pointer :: sys_grid
98 logical,
public :: compute_surface_correction
99 integer,
allocatable :: surface_points(:)
100 integer,
allocatable :: inner_stencil(:)
102 logical :: enforce_coulomb_gauge
104 real(real64) :: coulomb_gauge_tolerance
108 real(real64) :: poisson_prefactor
109 type(poisson_t) :: poisson_solver
115 procedure :: dapply_inner_stencil_mask, zapply_inner_stencil_mask
116 generic :: apply_inner_stencil_mask => dapply_inner_stencil_mask, zapply_inner_stencil_mask
118 procedure :: dget_vector_potential, zget_vector_potential
119 generic :: get_vector_potential => dget_vector_potential, zget_vector_potential
121 procedure :: dget_scalar_potential, zget_scalar_potential
122 generic :: get_scalar_potential => dget_scalar_potential, zget_scalar_potential
124 procedure :: dget_trans_field, zget_trans_field
125 generic :: get_trans_field => dget_trans_field, zget_trans_field
127 procedure :: dget_long_field, zget_long_field
128 generic :: get_long_field => dget_long_field, zget_long_field
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
143 logical,
allocatable :: mask(:)
144 logical :: visualize_boxes
145 integer :: number_of_points
150 this%sys_grid => sys_grid
153 call poisson_init(this%poisson_solver, namespace, space, sys_grid%der, system_mc, sys_grid%stencil, label =
"Helmholtz")
157 safe_allocate(mask(1:sys_grid%np))
160 number_of_points = count(mask)
162 safe_allocate(this%inner_stencil(1:number_of_points))
174 call parse_variable(namespace,
'HelmholtzEnforceCoulombGauge', .
true., this%enforce_coulomb_gauge)
175 if (this%enforce_coulomb_gauge)
then
184 call parse_variable(namespace,
'HelmholtzCoulombGaugeTolerance', 1e-5_real64, this%coulomb_gauge_tolerance)
196 call parse_variable(namespace,
'SurfaceCorrection', .false., this%compute_surface_correction)
197 if (this%compute_surface_correction)
then
201 if (this%compute_surface_correction)
then
203 mask = sys_grid%box%get_surface_points(namespace, sys_grid%spacing, sys_grid%np, sys_grid%x)
204 number_of_points = count(mask)
207 safe_allocate(this%surface_points(1:number_of_points))
211 select case (sys_grid%box%dim)
219 message(1) =
"Internal error: surface correction for helmholtz decomposition can only be called for 1D, 2D or 3D."
223 safe_deallocate_a(mask)
238 call parse_variable(namespace,
'HelmholtzVisualizeBoxes', .false., visualize_boxes)
239 if (visualize_boxes)
then
247 integer,
intent(in) :: np
248 logical,
intent(in) :: mask(:)
249 integer,
intent(out) :: indices(:)
271 safe_deallocate_a(this%inner_stencil)
272 safe_deallocate_a(this%surface_points)
281 class(
space_t),
intent(in) :: space
283 real(real64) :: coords(space%dim), normal_vector(space%dim), surface_element, ra
284 integer :: ip, ii, iunit
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(:)
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(:)
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
324#include "complex.F90"
325#include "helmholtz_decomposition_inc.F90"
329#include "helmholtz_decomposition_inc.F90"
This module contains interfaces for BLAS routines You should not use these routines directly....
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
real(real64), parameter, public m_four
real(real64), parameter, public m_pi
some mathematical constants
real(real64), parameter, public m_twothird
real(real64), parameter, public m_one
real(real64), parameter, public m_three
This module implements the underlying real-space grid.
The Helmholtz decomposition is intended to contain "only mathematical" functions and procedures to co...
subroutine helmholtz_finalize(this)
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.
subroutine, public io_close(iunit, grp)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
This module is intended to contain "only mathematical" functions and procedures.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
subroutine, public messages_not_implemented(feature, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
This module handles the communicators for the various parallelization strategies.
subroutine, public poisson_init(this, namespace, space, der, mc, stencil, qtot, label, solver, verbose, force_serial, force_cmplx)
subroutine, public poisson_end(this)
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
This module defines the unit system, used for input and output.