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> 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