Octopus
grid.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2021 S. Ohlmann
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
23
24module grid_oct_m
25 use accel_oct_m
27 use box_oct_m
32 use cube_oct_m
37 use global_oct_m
39 use mesh_oct_m
42 use mpi_oct_m
46 use parser_oct_m
48 use space_oct_m
53 use unit_oct_m
56
57 implicit none
58
59 private
60 public :: &
61 grid_t, &
65 grid_end, &
74
76 type, extends(mesh_t) :: grid_t
77 ! Components are public by default
78 type(derivatives_t) :: der
79 type(stencil_t) :: stencil
80 type(symmetries_t) :: symm
81 type(symmetrizer_t) :: symmetrizer
82 end type grid_t
83
84 integer, parameter :: &
85 CURV_AFFINE = 1, &
86 curv_gygi = 2, &
87 curv_briggs = 3, &
88 curv_modine = 4
89
90contains
91
92 !-------------------------------------------------------------------
101 subroutine grid_init_stage_1(gr, namespace, space, symm, latt, n_sites, site_position)
102 type(grid_t), intent(inout) :: gr
103 type(namespace_t), intent(in) :: namespace
104 class(space_t), intent(in) :: space
105 type(symmetries_t), optional, intent(in) :: symm
106 type(lattice_vectors_t), optional, intent(in) :: latt
107 integer, optional, intent(in) :: n_sites
108 real(real64), optional, intent(in) :: site_position(:,:)
109
110 type(stencil_t) :: cube
111 integer :: enlarge(space%dim)
112 type(block_t) :: blk
113 integer :: idir
114 real(real64) :: grid_spacing(space%dim)
115
116 push_sub(grid_init_stage_1)
118 gr%box => box_factory_create(namespace, space, latt=latt, n_sites=n_sites, site_position=site_position)
119
120 if (present(symm)) then
121 gr%symm = symm
122 end if
123
124 !%Variable Spacing
125 !%Type float
126 !%Section Mesh
127 !%Description
128 !% The spacing between the points in the mesh. This controls the
129 !% quality of the discretization: smaller spacing gives more
130 !% precise results but increased computational cost.
131 !%
132 !% When using curvilinear coordinates, this is a canonical spacing
133 !% that will be changed locally by the transformation. In periodic
134 !% directions, your spacing may be slightly different than what
135 !% you request here, since the box size must be an integer
136 !% multiple of the spacing.
137 !%
138 !% The default value is defined by the image resolution if <tt>BoxShape =
139 !% box_image</tt>. Othewise here is no default otherwise.
140 !%
141 !% It is possible to have a different spacing in each one of the Cartesian directions
142 !% if we define <tt>Spacing</tt> as block of the form
143 !%
144 !% <tt>%Spacing
145 !% <br>&nbsp;&nbsp;spacing_x | spacing_y | spacing_z
146 !% <br>%</tt>
147 !%End
148
149 if(parse_block(namespace, 'Spacing', blk) == 0) then
150 if(parse_block_cols(blk,0) < space%dim) call messages_input_error(namespace, 'Spacing')
151 do idir = 1, space%dim
152 call parse_block_float(blk, 0, idir - 1, grid_spacing(idir), units_inp%length)
153 end do
154 call parse_block_end(blk)
155 else
156 call parse_variable(namespace, 'Spacing', -m_one, grid_spacing(1), units_inp%length)
157 grid_spacing = grid_spacing(1)
158 end if
159
160 if (associated(gr%box)) then
161 select type (box => gr%box)
162 type is (box_image_t)
163 do idir = 1, space%dim
164 ! default grid_spacing is determined from the pixel size such that one grid point = one pixel.
165 if (grid_spacing(idir) < m_zero) then
166 grid_spacing(idir) = box%pixel_size(idir)
167 end if
168 end do
169 end select
170 end if
172 if (any(grid_spacing(1:space%dim) < m_epsilon)) then
173 message(1) = " Your input for 'Spacing' is negative or zero."
174 call messages_fatal(1, namespace=namespace)
175 end if
176
177 !%Variable PeriodicBoundaryMask
178 !%Type block
179 !%Section Mesh
180 !%Description
181 !% (Experimental) Defines a mask for which periodic boundaries are replaced by zero boundary conditions.
182 !%End
183 if (parse_block(namespace, 'PeriodicBoundaryMask', blk) < 0) then
184 gr%masked_periodic_boundaries = .false.
185 else
186 gr%masked_periodic_boundaries = .true.
187 call parse_block_string(blk, 0, 0, gr%periodic_boundary_mask)
188 call messages_experimental('PeriodicBoundaryMask')
189 call parse_block_end(blk)
190 end if
191
192 ! Initialize coordinate system
193 call initialize_coordinate_system(gr, namespace, space, latt, n_sites, site_position)
195 ! initialize derivatives
196 call derivatives_init(gr%der, namespace, space, gr%coord_system)
197 ! the stencil used to generate the grid is a union of a cube (for
198 ! multigrid) and the Laplacian
199 call stencil_cube_get_lapl(cube, space%dim, order = 2)
200 call stencil_union(cube, gr%der%lapl%stencil, gr%stencil)
201 call stencil_end(cube)
202
203 enlarge = 2
204 enlarge = max(enlarge, gr%der%n_ghost)
205
206 call mesh_init_stage_1(gr, namespace, space, gr%box, gr%coord_system, grid_spacing, enlarge)
207 call mesh_init_stage_2(gr, namespace, space, gr%box, gr%stencil)
208
209 pop_sub(grid_init_stage_1)
210 end subroutine grid_init_stage_1
211
212 !--------------------------------------------------------------------
214 subroutine grid_init_from_grid_stage_1(gr, original_gr, namespace, space, box, spacing_prefactor, latt, n_sites, site_position)
215 type(grid_t), intent(inout) :: gr
216 type(grid_t), intent(in) :: original_gr
217 type(namespace_t), intent(in) :: namespace
218 class(space_t), intent(in) :: space
219 class(box_t), target, optional, intent(in) :: box
220 real(real64), optional, intent(in) :: spacing_prefactor(:)
221 type(lattice_vectors_t), optional, intent(in) :: latt
222 integer, optional, intent(in) :: n_sites
223 real(real64), optional, intent(in) :: site_position(:,:)
224
225 type(stencil_t) :: cube
226 integer :: enlarge(space%dim)
227
229
230 ! First of all, create and associate the box that contains the grid
231 if (.not. present(box)) then
232 gr%box => box_factory_create(namespace, space)
233 else
234 gr%box => box
235 end if
236
237 !!! Copy all data !!!
238 ! Symmetries
239 gr%symm = original_gr%symm
240 ! Spacing
241 if (present(spacing_prefactor)) then
242 gr%spacing = spacing_prefactor * original_gr%spacing
243 if (any(abs(gr%spacing - m_one) > m_epsilon)) then
244 write(message(1), '(a)') "The two grids have different spacings, it will not be possible to establish a map between them"
245 call messages_warning(1)
246 end if
247 else
248 gr%spacing = original_gr%spacing
249 end if
250 ! Periodic boundaries
251 gr%masked_periodic_boundaries = original_gr%masked_periodic_boundaries
252 if (gr%masked_periodic_boundaries) gr%periodic_boundary_mask = original_gr%periodic_boundary_mask
253 ! Coordinate System
254 call initialize_coordinate_system(gr, namespace, space, latt, n_sites, site_position)
255
256 ! initialize derivatives
257 call derivatives_init(gr%der, namespace, space, gr%coord_system)
258 ! the stencil used to generate the grid is a union of a cube (for
259 ! multigrid) and the Laplacian
260 call stencil_cube_get_lapl(cube, space%dim, order = 2)
261 call stencil_union(cube, gr%der%lapl%stencil, gr%stencil)
262 call stencil_end(cube)
263
264 enlarge = 2
265 enlarge = max(enlarge, gr%der%n_ghost)
266
267 call mesh_init_stage_1(gr, namespace, space, gr%box, gr%coord_system, gr%spacing, enlarge)
268 call mesh_init_stage_2(gr, namespace, space, gr%box, gr%stencil)
269
271 end subroutine grid_init_from_grid_stage_1
272
273
274 !--------------------------------------------------------------------
278 subroutine initialize_coordinate_system(gr, namespace, space, latt, n_sites, site_position)
279 type(grid_t), intent(inout) :: gr
280 type(namespace_t), intent(in) :: namespace
281 class(space_t), intent(in) :: space
282 type(lattice_vectors_t), optional, intent(in) :: latt
283 integer, optional, intent(in) :: n_sites
284 real(real64), optional, intent(in) :: site_position(:,:)
285
286 integer :: cv_method
287
289
290 !%Variable CurvMethod
291 !%Type integer
292 !%Default curv_uniform
293 !%Section Mesh::Curvilinear
294 !%Description
295 !% The relevant functions in octopus are represented on a mesh in real space.
296 !% This mesh may be an evenly spaced regular rectangular grid (standard mode),
297 !% or else an adaptive or curvilinear grid. We have implemented
298 !% three kinds of adaptive meshes, although only one is currently working,
299 !% the one invented by F. Gygi (<tt>curv_gygi</tt>). The code will stop if any of
300 !% the other two is invoked. All are experimental with domain parallelization.
301 !%Option curv_affine 1
302 !% Regular, uniform rectangular grid.
303 !%Option curv_gygi 2
304 !% The deformation of the grid is done according to the scheme described by
305 !% F. Gygi [F. Gygi and G. Galli, <i>Phys. Rev. B</i> <b>52</b>, R2229 (1995)].
306 !%Option curv_briggs 3
307 !% The deformation of the grid is done according to the scheme described by
308 !% Briggs [E.L. Briggs, D.J. Sullivan, and J. Bernholc, <i>Phys. Rev. B</i> <b>54</b> 14362 (1996)]
309 !% (NOT WORKING).
310 !%Option curv_modine 4
311 !% The deformation of the grid is done according to the scheme described by
312 !% Modine [N.A. Modine, G. Zumbach and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289 (1997)]
313 !% (NOT WORKING).
314 !%End
315 call parse_variable(namespace, 'CurvMethod', curv_affine, cv_method)
316 if (.not. varinfo_valid_option('CurvMethod', cv_method)) call messages_input_error(namespace, 'CurvMethod')
317 if (cv_method /= curv_affine .and. accel_is_enabled()) then
318 call messages_write('Curvilinear coordinates on GPUs is not implemented')
319 call messages_fatal(namespace=namespace)
320 end if
321 if (present(latt)) then
322 if (cv_method /= curv_affine .and. latt%nonorthogonal) then
323 call messages_input_error(namespace, 'CurvMethod', 'Curvilinear coordinates with non-orthogonal cells are not implemented')
324 end if
325 end if
326 call messages_print_var_option("CurvMethod", cv_method, namespace=namespace)
327
328 ! FIXME: The other two methods are apparently not working
329 if (cv_method > curv_gygi) then
330 call messages_experimental('Selected curvilinear coordinates method')
331 end if
332
333 select case (cv_method)
334 case (curv_briggs)
335 gr%coord_system => curv_briggs_t(namespace, space%dim, &
336 gr%box%bounding_box_l(1:space%dim), gr%spacing(1:space%dim))
337 case (curv_gygi)
338 if (present(site_position) .and. present(n_sites)) then
339 gr%coord_system => curv_gygi_t(namespace, space%dim, n_sites, site_position)
340 else
341 message(1) = "Option CurvMethod = curv_gygi is not currently implemented without ions present."
342 call messages_fatal(1, namespace=namespace)
343 end if
344 case (curv_modine)
345 if (present(site_position) .and. present(n_sites)) then
346 gr%coord_system => curv_modine_t(namespace, space%dim, n_sites, site_position, &
347 gr%box%bounding_box_l(1:space%dim), gr%spacing(1:space%dim))
348 else
349 message(1) = "Option CurvMethod = curv_modine is not currently implemented without ions present."
350 call messages_fatal(1, namespace=namespace)
351 end if
352 case (curv_affine)
353 if (present(latt)) then
354 if (latt%nonorthogonal) then
355 gr%coord_system => affine_coordinates_t(namespace, space%dim, latt%rlattice_primitive)
356 else
357 gr%coord_system => cartesian_t(namespace, space%dim)
358 end if
359 else
360 gr%coord_system => cartesian_t(namespace, space%dim)
361 end if
362 end select
363
365 end subroutine initialize_coordinate_system
366
367
368 !-------------------------------------------------------------------
376 subroutine grid_init_stage_2(gr, namespace, space, mc, qvector)
377 type(grid_t), target, intent(inout) :: gr
378 type(namespace_t), intent(in) :: namespace
379 class(space_t), intent(in) :: space
380 type(multicomm_t), intent(in) :: mc
381 real(real64), optional, intent(in) :: qvector(:)
382
383 push_sub(grid_init_stage_2)
384
385 call mesh_init_stage_3(gr, namespace, space, gr%stencil, mc)
386
387 call nl_operator_global_init(namespace)
388 call derivatives_build(gr%der, namespace, space, gr, qvector)
389
390 ! print info concerning the grid
391 call grid_write_info(gr, namespace=namespace)
392
393 if(space%dim == 3) then
394 call symmetrizer_init(gr%symmetrizer, gr, gr%symm)
395 end if
396
397 pop_sub(grid_init_stage_2)
398 end subroutine grid_init_stage_2
399
400
401 !-------------------------------------------------------------------
404 subroutine grid_end(gr)
405 type(grid_t), intent(inout) :: gr
406
407 class(coordinate_system_t), pointer :: coord_system
408 class(box_t), pointer :: box
409
410 push_sub(grid_end)
411
413
414 call derivatives_end(gr%der)
415
416 ! We need to take a pointer here, otherwise we run into a gfortran bug.
417 coord_system => gr%coord_system
418 safe_deallocate_p(coord_system)
419 box => gr%box
420 safe_deallocate_p(box)
421
422 call mesh_end(gr)
423
424 call stencil_end(gr%stencil)
425
426 call symmetrizer_end(gr%symmetrizer)
427
428 pop_sub(grid_end)
429 end subroutine grid_end
430
431
432 !-------------------------------------------------------------------
433 subroutine grid_write_info(gr, iunit, namespace)
434 type(grid_t), intent(in) :: gr
435 integer, optional, intent(in) :: iunit
436 type(namespace_t), optional, intent(in) :: namespace
437
438 push_sub(grid_write_info)
439
440 call messages_print_with_emphasis(msg="Grid", iunit=iunit, namespace=namespace)
441
442 message(1) = 'Simulation Box:'
443 call messages_info(1, iunit=iunit, namespace=namespace)
444 call gr%box%write_info(iunit, namespace)
445
446 message(1) = "Main mesh:"
447 call messages_info(1, iunit=iunit, namespace=namespace)
448 call mesh_write_info(gr, iunit, namespace)
449
450 if (gr%use_curvilinear) then
451 call gr%coord_system%write_info(iunit, namespace)
452 end if
453
454 call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
455
456 pop_sub(grid_write_info)
457 end subroutine grid_write_info
458
459 !-------------------------------------------------------------------
463 subroutine grid_lattice_vectors_update(gr, space, namespace, mc, new_latt)
464 type(grid_t), intent(inout) :: gr
465 type(space_t), intent(in) :: space
466 type(namespace_t), intent(in) :: namespace
467 type(multicomm_t), intent(in) :: mc
468 type(lattice_vectors_t), intent(in) :: new_latt
470 real(real64) :: new_box_bounds(2, space%dim)
471
473
474 call new_latt%write_info(namespace)
475
476 ! Regenerate the coordinate system
477 select type (coord_system=>gr%coord_system)
478 type is (affine_coordinates_t)
479 if (new_latt%nonorthogonal) then
480 deallocate(gr%coord_system)
481 gr%coord_system => affine_coordinates_t(namespace, space%dim, new_latt%rlattice_primitive)
482 end if
483 end select
484
485 ! Regenerate the grid keeping the same number of points
486 select type (coord_system=>gr%coord_system)
487 class is (affine_coordinates_t)
488 new_box_bounds = gr%box%bounds(coord_system%basis)
489 gr%spacing = (new_box_bounds(2,:)-new_box_bounds(1,:))/gr%idx%ll(:)
490 write(message(1), '(a)') trim(gr%box%short_info(unit_angstrom))
491 call messages_info(1, namespace=namespace)
492 call messages_new_line()
493 end select
494
495 safe_deallocate_a(gr%x)
496 safe_deallocate_a(gr%vol_pp)
497 call mesh_init_stage_3(gr, namespace, space, gr%stencil, mc, regenerate=.true.)
498
499 call derivatives_build(gr%der, namespace, space, gr, regenerate=.true., verbose=.false.)
500
502 end subroutine grid_lattice_vectors_update
503
504#include "undef.F90"
505#include "real.F90"
506#include "grid_inc.F90"
507
508#include "undef.F90"
509#include "complex.F90"
510#include "grid_inc.F90"
511
512end module grid_oct_m
513
514!! Local Variables:
515!! mode: f90
516!! coding: utf-8
517!! End:
pure logical function, public accel_is_enabled()
Definition: accel.F90:402
Module, implementing a factory for boxes.
class(box_t) function, pointer, public box_factory_create(namespace, space, latt, n_sites, site_position)
initialize a box of any type
This module implements the curvilinear coordinates given in E.L. Briggs, D.J. Sullivan,...
This module implements the curvilinear coordinates given in F. Gygi and G. Galli, PRB 52 R2229 (1996)...
Definition: curv_gygi.F90:117
This module implements the curvilinear coordinates given in N. A. Modine, G. Zumbach,...
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public derivatives_build(der, namespace, space, mesh, qvector, regenerate, verbose)
build the derivatives object:
subroutine, public derivatives_init(der, namespace, space, coord_system, order)
subroutine, public derivatives_end(der)
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public grid_init_stage_1(gr, namespace, space, symm, latt, n_sites, site_position)
First stage of the grid initialization.
Definition: grid.F90:195
integer, parameter curv_gygi
Definition: grid.F90:177
integer, parameter curv_briggs
Definition: grid.F90:177
subroutine, public zgrid_symmetrize_scalar_field(gr, field, suppress_warning)
Definition: grid.F90:810
subroutine initialize_coordinate_system(gr, namespace, space, latt, n_sites, site_position)
this subroutine initializes the coordinate system
Definition: grid.F90:372
subroutine, public grid_init_from_grid_stage_1(gr, original_gr, namespace, space, box, spacing_prefactor, latt, n_sites, site_position)
this subroutine allows to create a grid from an existing grid
Definition: grid.F90:308
subroutine, public zgrid_symmetrize_vector_field(gr, field, suppress_warning)
Definition: grid.F90:836
subroutine, public grid_init_stage_2(gr, namespace, space, mc, qvector)
Second stage of the grid initialization.
Definition: grid.F90:470
subroutine, public dgrid_symmetrize_vector_field(gr, field, suppress_warning)
Definition: grid.F90:692
subroutine, public grid_write_info(gr, iunit, namespace)
Definition: grid.F90:527
subroutine, public dgrid_symmetrize_single(gr, iop, field, symm_field, suppress_warning)
Definition: grid.F90:720
subroutine, public grid_lattice_vectors_update(gr, space, namespace, mc, new_latt)
Regenerate the grid information after update of the lattice vectors.
Definition: grid.F90:557
subroutine, public zgrid_symmetrize_single(gr, iop, field, symm_field, suppress_warning)
Definition: grid.F90:864
subroutine, public grid_end(gr)
finalize a grid object
Definition: grid.F90:498
subroutine, public dgrid_symmetrize_scalar_field(gr, field, suppress_warning)
Definition: grid.F90:666
integer, parameter curv_modine
Definition: grid.F90:177
This module contains subroutines, related to the initialization of the mesh.
Definition: mesh_init.F90:117
subroutine, public mesh_init_stage_3(mesh, namespace, space, stencil, mc, parent, regenerate)
When running parallel in domains, stencil and np_stencil are needed to compute the ghost points....
Definition: mesh_init.F90:538
subroutine, public mesh_init_stage_1(mesh, namespace, space, box, coord_system, spacing, enlarge)
First stage mesh initialization.
Definition: mesh_init.F90:170
subroutine, public mesh_init_stage_2(mesh, namespace, space, box, stencil, regenerate)
This subroutine creates the global array of spatial indices and the inverse mapping.
Definition: mesh_init.F90:290
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public mesh_write_info(this, iunit, namespace)
Definition: mesh.F90:304
recursive subroutine, public mesh_end(this)
Definition: mesh.F90:680
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
character(len=512), private msg
Definition: messages.F90:165
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_new_line()
Definition: messages.F90:1134
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
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
This module defines non-local operators.
subroutine, public nl_operator_global_init(namespace)
initialize global settings for non-local operators
subroutine, public nl_operator_global_end()
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
This module defines routines, generating operators for a cubic stencil.
subroutine, public stencil_cube_get_lapl(this, dim, order)
This module defines stencils used in Octopus.
Definition: stencil.F90:135
subroutine, public stencil_end(this)
Definition: stencil.F90:214
subroutine, public stencil_union(st1, st2, stu)
Definition: stencil.F90:226
subroutine, public symmetrizer_end(this)
subroutine, public symmetrizer_init(this, mesh, symm)
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.
type(unit_t), public unit_angstrom
For XYZ files.
type(unit_system_t), public units_inp
the units systems for reading and writing
Class implementing a box generated from a 2D image.
Definition: box_image.F90:133
class to tell whether a point is inside or outside
Definition: box.F90:141
abstract class to describe coordinate systems
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:169
Describes mesh distribution to nodes.
Definition: mesh.F90:186
Stores all communicators and groups.
Definition: multicomm.F90:206
The class representing the stencil, which is used for non-local mesh operations.
Definition: stencil.F90:163
int true(void)