1!! Copyright (C) 2002-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
21module grid_oct_m
22  use cube_oct_m
23  use curvilinear_oct_m
24  use derivatives_oct_m
25  use double_grid_oct_m
26  use geometry_oct_m
27  use global_oct_m
28  use mesh_oct_m
29  use mesh_init_oct_m
30  use messages_oct_m
31  use mpi_oct_m
32  use multicomm_oct_m
33  use multigrid_oct_m
34  use namespace_oct_m
35  use nl_operator_oct_m
36  use parser_oct_m
37  use profiling_oct_m
38  use space_oct_m
39  use simul_box_oct_m
40  use stencil_oct_m
41  use stencil_cube_oct_m
42  use unit_oct_m
43  use unit_system_oct_m
44
45  implicit none
46
47  private
48  public ::                &
49    grid_t,                &
50    grid_init_stage_0,     &
51    grid_init_stage_1,     &
52    grid_init_stage_2,     &
53    grid_end,              &
54    grid_write_info,       &
55    grid_create_multigrid
56
57  type grid_t
58    ! Components are public by default
59    type(simul_box_t)           :: sb
60    type(mesh_t)                :: mesh
61    type(multigrid_level_t)     :: fine
62    type(derivatives_t)         :: der
63    type(curvilinear_t)         :: cv
64    type(multigrid_t), pointer  :: mgrid
65    type(double_grid_t)         :: dgrid
66    logical                     :: have_fine_mesh
67    type(stencil_t)             :: stencil
68  end type grid_t
69
70
71contains
72
73  !-------------------------------------------------------------------
74  !>
75  !! "Zero-th" stage of grid initialization. It initializes the simulation box.
76  subroutine grid_init_stage_0(gr, namespace, geo, space)
77    type(grid_t),          intent(inout) :: gr
78    type(namespace_t),     intent(in)    :: namespace
79    type(geometry_t),      intent(inout) :: geo
80    type(space_t),         intent(in)    :: space
81
82    PUSH_SUB(grid_init_stage_0)
83
84    call simul_box_init(gr%sb, namespace, geo, space)
85
86    POP_SUB(grid_init_stage_0)
87  end subroutine grid_init_stage_0
88
89
90  !-------------------------------------------------------------------
91  subroutine grid_init_stage_1(gr, namespace, geo)
92    type(grid_t),      intent(inout) :: gr
93    type(namespace_t), intent(in)    :: namespace
94    type(geometry_t),  intent(in)    :: geo
95
96    type(stencil_t) :: cube
97    integer :: enlarge(1:MAX_DIM)
98    type(block_t) :: blk
99    integer :: idir
100    FLOAT :: def_h, def_rsize
101    FLOAT :: grid_spacing(1:MAX_DIM)
102
103    PUSH_SUB(grid_init_stage_1)
104
105    !%Variable UseFineMesh
106    !%Type logical
107    !%Default no
108    !%Section Mesh
109    !%Description
110    !% If enabled, <tt>Octopus</tt> will use a finer mesh for the calculation
111    !% of the forces or other sensitive quantities.
112    !% Experimental, and incompatible with domain-parallelization.
113    !%End
114    if (gr%sb%dim == 3) then
115      call parse_variable(namespace, 'UseFineMesh', .false., gr%have_fine_mesh)
116    else
117      gr%have_fine_mesh = .false.
118    end if
119
120    if(gr%have_fine_mesh) call messages_experimental("UseFineMesh")
121
122    call geometry_grid_defaults(geo, def_h, def_rsize)
123
124    ! initialize to -1
125    grid_spacing = -M_ONE
126
127    !%Variable Spacing
128    !%Type float
129    !%Section Mesh
130    !%Description
131    !% The spacing between the points in the mesh. This controls the
132    !% quality of the discretization: smaller spacing gives more
133    !% precise results but increased computational cost.
134    !%
135    !% When using curvilinear coordinates, this is a canonical spacing
136    !% that will be changed locally by the transformation. In periodic
137    !% directions, your spacing may be slightly different than what
138    !% you request here, since the box size must be an integer
139    !% multiple of the spacing.
140    !%
141    !% The default value is defined by the species if only default pseudopotentials are used
142    !% or by the image resolution if <tt>BoxShape = box_image</tt>. Otherwise, there is
143    !% no default.
144    !%
145    !% It is possible to have a different spacing in each one of the Cartesian directions
146    !% if we define <tt>Spacing</tt> as block of the form
147    !%
148    !% <tt>%Spacing
149    !% <br>&nbsp;&nbsp;spacing_x | spacing_y | spacing_z
150    !% <br>%</tt>
151    !%End
152
153    if(parse_block(namespace, 'Spacing', blk) == 0) then
154      if(parse_block_cols(blk,0) < gr%sb%dim) call messages_input_error(namespace, 'Spacing')
155      do idir = 1, gr%sb%dim
156        call parse_block_float(blk, 0, idir - 1, grid_spacing(idir), units_inp%length)
157        if(def_h > M_ZERO) call messages_check_def(grid_spacing(idir), .true., def_h, 'Spacing', units_out%length)
158      end do
159      call parse_block_end(blk)
160    else
161      call parse_variable(namespace, 'Spacing', -M_ONE, grid_spacing(1), units_inp%length)
162      grid_spacing(1:gr%sb%dim) = grid_spacing(1)
163      if(def_h > M_ZERO) call messages_check_def(grid_spacing(1), .true., def_h, 'Spacing', units_out%length)
164    end if
165
166#if defined(HAVE_GDLIB)
167    if(gr%sb%box_shape == BOX_IMAGE) then
168      do idir = 1, gr%sb%dim
169        ! default grid_spacing is determined from lsize and the size of the image
170        if(grid_spacing(idir) < M_ZERO) then
171          grid_spacing(idir) = M_TWO*gr%sb%lsize(idir)/TOFLOAT(gr%sb%image_size(idir))
172        end if
173      end do
174    end if
175#endif
176
177    if (any(grid_spacing(1:gr%sb%dim) < M_EPSILON)) then
178      if (def_h > M_ZERO .and. def_h < huge(def_h)) then
179        call geometry_grid_defaults_info(geo)
180        do idir = 1, gr%sb%dim
181          grid_spacing(idir) = def_h
182          write(message(1), '(a,i1,3a,f6.3)') "Info: Using default spacing(", idir, &
183            ") [", trim(units_abbrev(units_out%length)), "] = ",                        &
184            units_from_atomic(units_out%length, grid_spacing(idir))
185          call messages_info(1)
186        end do
187        ! Note: the default automatically matches the 'recommended' value compared by messages_check_def above.
188      else
189        message(1) = 'Either:'
190        message(2) = "   *) variable 'Spacing' is not defined and"
191        message(3) = "      I can't find a suitable default"
192        message(4) = "   *) your input for 'Spacing' is negative or zero"
193        call messages_fatal(4)
194      end if
195    end if
196
197    !%Variable PeriodicBoundaryMask
198    !%Type block
199    !%Section Mesh
200    !%Description
201    !% (Experimental) Defines a mask for which periodic boundaries are replaced by zero boundary conditions.
202    !%End
203    if(parse_block(namespace, 'PeriodicBoundaryMask', blk) < 0) then
204      gr%mesh%masked_periodic_boundaries = .false.
205    else
206      gr%mesh%masked_periodic_boundaries = .true.
207      call parse_block_string(blk, 0, 0, gr%mesh%periodic_boundary_mask)
208      call messages_experimental('PeriodicBoundaryMask')
209    end if
210
211    ! initialize curvilinear coordinates
212    call curvilinear_init(gr%cv, namespace, gr%sb, geo, grid_spacing)
213
214    ! initialize derivatives
215    call derivatives_nullify(gr%der)
216    call derivatives_init(gr%der, namespace, gr%sb, gr%cv%method /= CURV_METHOD_UNIFORM)
217
218    call double_grid_init(gr%dgrid, namespace, gr%sb)
219
220    enlarge = 0
221    enlarge(1:gr%sb%dim) = 2
222    enlarge = max(enlarge, double_grid_enlarge(gr%dgrid))
223    enlarge = max(enlarge, gr%der%n_ghost)
224
225    ! now we generate the mesh and the derivatives
226    call mesh_init_stage_1(gr%mesh, gr%sb, gr%cv, grid_spacing, enlarge)
227
228    ! the stencil used to generate the grid is a union of a cube (for
229    ! multigrid) and the Laplacian.
230    call stencil_cube_get_lapl(cube, gr%sb%dim, order = 2)
231    call stencil_union(gr%sb%dim, cube, gr%der%lapl%stencil, gr%stencil)
232    call stencil_end(cube)
233
234    call mesh_init_stage_2(gr%mesh, gr%sb, geo, gr%cv, gr%stencil, namespace)
235
236    POP_SUB(grid_init_stage_1)
237
238  end subroutine grid_init_stage_1
239
240
241  !-------------------------------------------------------------------
242  subroutine grid_init_stage_2(gr, namespace, mc, geo)
243    type(grid_t), target, intent(inout) :: gr
244    type(namespace_t),    intent(in)    :: namespace
245    type(multicomm_t),    intent(in)    :: mc
246    type(geometry_t),     intent(in)    :: geo
247
248    PUSH_SUB(grid_init_stage_2)
249
250    call mesh_init_stage_3(gr%mesh, namespace, gr%stencil, mc)
251
252    call nl_operator_global_init(namespace)
253    if(gr%have_fine_mesh) then
254      message(1) = "Info: coarse mesh"
255      call messages_info(1)
256    end if
257    call derivatives_build(gr%der, namespace, gr%mesh)
258
259    ! initialize a finer mesh to hold the density, for this we use the
260    ! multigrid routines
261
262    if(gr%have_fine_mesh) then
263
264      if(gr%mesh%parallel_in_domains) then
265        message(1) = 'UseFineMesh does not work with domain parallelization.'
266        call messages_fatal(1)
267      end if
268
269      SAFE_ALLOCATE(gr%fine%mesh)
270      SAFE_ALLOCATE(gr%fine%der)
271
272      call multigrid_mesh_double(geo, gr%cv, gr%mesh, gr%fine%mesh, gr%stencil, namespace)
273
274      call derivatives_nullify(gr%fine%der)
275      call derivatives_init(gr%fine%der, namespace, gr%mesh%sb, gr%cv%method /= CURV_METHOD_UNIFORM)
276
277      call mesh_init_stage_3(gr%fine%mesh, namespace, gr%stencil, mc)
278
279      call multigrid_get_transfer_tables(gr%fine%tt, gr%fine%mesh, gr%mesh)
280
281      message(1) = "Info: fine mesh"
282      call messages_info(1)
283      call derivatives_build(gr%fine%der, namespace, gr%fine%mesh)
284
285      gr%fine%der%coarser => gr%der
286      gr%der%finer =>  gr%fine%der
287      gr%fine%der%to_coarser => gr%fine%tt
288      gr%der%to_finer => gr%fine%tt
289
290    else
291      gr%fine%mesh => gr%mesh
292      gr%fine%der => gr%der
293    end if
294
295    ! multigrids are not initialized by default
296    nullify(gr%mgrid)
297
298    ! print info concerning the grid
299    call grid_write_info(gr, geo, stdout)
300
301    POP_SUB(grid_init_stage_2)
302  end subroutine grid_init_stage_2
303
304
305  !-------------------------------------------------------------------
306  subroutine grid_end(gr)
307    type(grid_t), intent(inout) :: gr
308
309    PUSH_SUB(grid_end)
310
311    call nl_operator_global_end()
312
313    if(gr%have_fine_mesh) then
314      call derivatives_end(gr%fine%der)
315      call mesh_end(gr%fine%mesh)
316      SAFE_DEALLOCATE_P(gr%fine%mesh)
317      SAFE_DEALLOCATE_P(gr%fine%der)
318      SAFE_DEALLOCATE_P(gr%fine%tt%to_coarse)
319      SAFE_DEALLOCATE_P(gr%fine%tt%to_fine1)
320      SAFE_DEALLOCATE_P(gr%fine%tt%to_fine2)
321      SAFE_DEALLOCATE_P(gr%fine%tt%to_fine4)
322      SAFE_DEALLOCATE_P(gr%fine%tt%to_fine8)
323      SAFE_DEALLOCATE_P(gr%fine%tt%fine_i)
324    end if
325
326    call double_grid_end(gr%dgrid)
327
328    call derivatives_end(gr%der)
329    call curvilinear_end(gr%cv)
330    call mesh_end(gr%mesh)
331
332    if(associated(gr%mgrid)) then
333      call multigrid_end(gr%mgrid)
334      SAFE_DEALLOCATE_P(gr%mgrid)
335    end if
336
337    call stencil_end(gr%stencil)
338
339    POP_SUB(grid_end)
340  end subroutine grid_end
341
342
343  !-------------------------------------------------------------------
344  subroutine grid_write_info(gr, geo, iunit)
345    type(grid_t),     intent(in) :: gr
346    type(geometry_t), intent(in) :: geo
347    integer,          intent(in) :: iunit
348
349    PUSH_SUB(grid_write_info)
350
351    if(.not.mpi_grp_is_root(mpi_world)) then
352      if(debug%info) call messages_debug_newlines(6)
353      POP_SUB(grid_write_info)
354      return
355    end if
356
357    call messages_print_stress(iunit, "Grid")
358    call simul_box_write_info(gr%sb, geo, iunit)
359
360    if(gr%have_fine_mesh) then
361      message(1) = "Wave-functions mesh:"
362      call messages_info(1, iunit)
363      call mesh_write_info(gr%mesh, iunit)
364      message(1) = "Density mesh:"
365    else
366      message(1) = "Main mesh:"
367    end if
368    call messages_info(1, iunit)
369    call mesh_write_info(gr%fine%mesh, iunit)
370
371    if (gr%mesh%use_curvilinear) then
372      call curvilinear_write_info(gr%cv, iunit)
373    end if
374
375    call messages_print_stress(iunit)
376
377    POP_SUB(grid_write_info)
378  end subroutine grid_write_info
379
380
381  !-------------------------------------------------------------------
382  subroutine grid_create_multigrid(gr, namespace, geo, mc)
383    type(grid_t),      intent(inout) :: gr
384    type(namespace_t), intent(in)    :: namespace
385    type(geometry_t),  intent(in)    :: geo
386    type(multicomm_t), intent(in)    :: mc
387
388    PUSH_SUB(grid_create_multigrid)
389
390    SAFE_ALLOCATE(gr%mgrid)
391    call multigrid_init(gr%mgrid, namespace, geo, gr%cv, gr%mesh, gr%der, gr%stencil, mc)
392
393    POP_SUB(grid_create_multigrid)
394  end subroutine grid_create_multigrid
395
396end module grid_oct_m
397
398!! Local Variables:
399!! mode: f90
400!! coding: utf-8
401!! End:
402