1!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira 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 mesh_oct_m 22 use basis_set_abst_oct_m 23 use curvilinear_oct_m 24 use geometry_oct_m 25 use global_oct_m 26 use hypercube_oct_m 27 use index_oct_m 28 use io_oct_m 29 use io_binary_oct_m 30 use mesh_cube_map_oct_m 31 use messages_oct_m 32 use mpi_oct_m 33 use namespace_oct_m 34 use par_vec_oct_m 35 use partition_oct_m 36 use parser_oct_m 37 use profiling_oct_m 38 use simul_box_oct_m 39 use symmetries_oct_m 40 use symm_op_oct_m 41 use species_oct_m 42 use unit_oct_m 43 use unit_system_oct_m 44 45 implicit none 46 47 private 48 public :: & 49 mesh_t, & 50 mesh_plane_t, & 51 mesh_line_t, & 52 mesh_dump, & 53 mesh_load, & 54 mesh_check_dump_compatibility, & 55 mesh_end, & 56 mesh_double_box, & 57 mesh_inborder, & 58 mesh_r, & 59 mesh_gcutoff, & 60 mesh_write_info, & 61 mesh_nearest_point, & 62 mesh_nearest_point_infos, & 63 mesh_periodic_point, & 64 mesh_global_memory, & 65 mesh_local_memory, & 66 mesh_x_global, & 67 mesh_write_fingerprint, & 68 mesh_read_fingerprint, & 69 mesh_compact_boundaries, & 70 mesh_check_symmetries 71 72 !> Describes mesh distribution to nodes. 73 !! 74 !! Some general things: 75 !! All members of type(mesh_t) are equal on all 76 !! nodes when running parallel except 77 !! - np, np_part 78 !! - x, vol_pp 79 !! These four are defined for all the points the node is responsible for. 80 type, extends(basis_set_abst_t) :: mesh_t 81 ! Components are public by default 82 type(simul_box_t), pointer :: sb !< simulation box 83 type(curvilinear_t), pointer :: cv 84 type(index_t) :: idx 85 logical :: use_curvilinear 86 87 FLOAT :: spacing(MAX_DIM) !< the (constant) spacing between the points 88 89 !> When running serially, the local number of points is 90 !! equal to the global number of points. 91 !! Otherwise, the next two are different on each node. 92 integer :: np !< Local number of points in mesh 93 integer :: np_part !< Local points plus ghost points plus boundary points. 94 integer :: np_global !< Global number of points in mesh. 95 integer :: np_part_global !< Global number of inner points and boundary points. 96 !> will I run parallel in domains? 97 !! yes or no?? 98 logical :: parallel_in_domains 99 type(mpi_grp_t) :: mpi_grp !< the mpi group describing parallelization in domains 100 type(pv_t) :: vp !< describes parallel vectors defined on the mesh. 101 type(partition_t) :: inner_partition !< describes how the inner points are assigned to the domains 102 type(partition_t) :: bndry_partition !< describes how the boundary points are assigned to the domains 103 104 FLOAT, allocatable :: x(:,:) !< The (local) \b points 105 integer, allocatable :: resolution(:, :, :) 106 FLOAT :: volume_element !< The global volume element. 107 FLOAT :: surface_element(MAX_DIM) 108 FLOAT, allocatable :: vol_pp(:) !< Element of volume for curvilinear coordinates. 109 110 type(mesh_cube_map_t) :: cube_map 111 112 logical :: masked_periodic_boundaries 113 character(len=256) :: periodic_boundary_mask 114 contains 115 procedure :: load => mesh_load 116 procedure :: dump => mesh_dump 117 procedure :: end => mesh_end 118 procedure :: init => mesh_init 119 procedure :: write_info => mesh_write_info 120 end type mesh_t 121 122 !> This data type defines a plane, and a regular grid defined on 123 !! this plane (or, rather, on a portion of this plane) 124 !! n should be a unit vector, that determines the normal of the plane. 125 !! Origin is a point belonging to the plane 126 !! u and v are unit orthogonal vectors belonging to the plane 127 !! The grid is generated by the vectors u and v: 128 !! x_{i,j} = origin + i*spacing*u + j*spacing*v, 129 !! for nu <= i <= mu and nv <= j <= mv 130 type mesh_plane_t 131 ! Components are public by default 132 FLOAT :: n(MAX_DIM) 133 FLOAT :: u(MAX_DIM), v(MAX_DIM) 134 FLOAT :: origin(MAX_DIM) 135 FLOAT :: spacing 136 integer :: nu, mu, nv, mv 137 end type mesh_plane_t 138 139 !> This data type defines a line, and a regular grid defined on this 140 !! line (or rather, on a portion of this line). 141 type mesh_line_t 142 ! Components are public by default 143 FLOAT :: n(MAX_DIM) 144 FLOAT :: u(MAX_DIM) 145 FLOAT :: origin(MAX_DIM) 146 FLOAT :: spacing 147 integer :: nu, mu 148 end type mesh_line_t 149 150 character(len=17), parameter :: dump_tag = '*** mesh_dump ***' 151 152contains 153 154 subroutine mesh_init(this) 155 class(mesh_t), intent(inout) :: this 156 157 PUSH_SUB(mesh_init) 158 159 call this%set_time_dependent(.false.) 160 161 POP_SUB(mesh_init) 162 end subroutine mesh_init 163 164! --------------------------------------------------------- 165 !> finds the dimension of a box doubled in the non-periodic dimensions 166 subroutine mesh_double_box(sb, mesh, alpha, db) 167 type(simul_box_t), intent(in) :: sb 168 type(mesh_t), intent(in) :: mesh 169 FLOAT, intent(in) :: alpha !< enlargement factor for double box 170 integer, intent(out) :: db(MAX_DIM) 171 172 integer :: idir 173 174 PUSH_SUB(mesh_double_box) 175 176 db = 1 177 178 ! double mesh with 2n points 179 do idir = 1, sb%periodic_dim 180 db(idir) = mesh%idx%ll(idir) 181 end do 182 do idir = sb%periodic_dim + 1, sb%dim 183 db(idir) = nint(alpha * (mesh%idx%ll(idir) - 1)) + 1 184 end do 185 186 POP_SUB(mesh_double_box) 187 end subroutine mesh_double_box 188 189 190 ! --------------------------------------------------------- 191 subroutine mesh_write_info(this, unit) 192 class(mesh_t), intent(in) :: this 193 integer, intent(in) :: unit 194 195 integer :: ii 196 FLOAT :: cutoff 197 198 if(.not.mpi_grp_is_root(mpi_world)) return 199 200 PUSH_SUB(mesh_write_info) 201 202 write(message(1),'(3a)') ' Spacing [', trim(units_abbrev(units_out%length)), '] = (' 203 do ii = 1, this%sb%dim 204 if(ii > 1) write(message(1), '(2a)') trim(message(1)), ',' 205 write(message(1), '(a,f6.3)') trim(message(1)), units_from_atomic(units_out%length, this%spacing(ii)) 206 end do 207 write(message(1), '(5a,f12.5)') trim(message(1)), ') ', & 208 ' volume/point [', trim(units_abbrev(units_out%length**this%sb%dim)), '] = ', & 209 units_from_atomic(units_out%length**this%sb%dim, this%vol_pp(1)) 210 211 write(message(2),'(a, i10)') ' # inner mesh = ', this%np_global 212 write(message(3),'(a, i10)') ' # total mesh = ', this%np_part_global 213 214 cutoff = mesh_gcutoff(this)**2 / M_TWO 215 write(message(4),'(3a,f12.6,a,f12.6)') ' Grid Cutoff [', trim(units_abbrev(units_out%energy)),'] = ', & 216 units_from_atomic(units_out%energy, cutoff), ' Grid Cutoff [Ry] = ', cutoff * M_TWO 217 call messages_info(4, unit) 218 219 POP_SUB(mesh_write_info) 220 end subroutine mesh_write_info 221 222 223 ! --------------------------------------------------------- 224 subroutine mesh_r(mesh, ip, rr, origin, coords) 225 type(mesh_t), intent(in) :: mesh 226 integer, intent(in) :: ip 227 FLOAT, intent(out) :: rr 228 FLOAT, intent(in), optional :: origin(:) !< origin(sb%dim) 229 FLOAT, intent(out), optional :: coords(:) !< coords(sb%dim) 230 231 FLOAT :: xx(MAX_DIM) 232 233 ! no push_sub because it is called too frequently 234 235 xx(1:mesh%sb%dim) = mesh%x(ip, 1:mesh%sb%dim) 236 if(present(origin)) xx(1:mesh%sb%dim) = xx(1:mesh%sb%dim) - origin(1:mesh%sb%dim) 237 rr = sqrt(dot_product(xx(1:mesh%sb%dim), xx(1:mesh%sb%dim))) 238 239 if(present(coords)) then 240 coords(1:MAX_DIM) = M_ZERO 241 coords(1:mesh%sb%dim) = xx(1:mesh%sb%dim) 242 end if 243 244 end subroutine mesh_r 245 246 247 !--------------------------------------------------------------------- 248 !> Finds out if a given point of a mesh belongs to the "border" of the 249 !! mesh. A point belongs to the border of the mesh if it is too close 250 !! to any of the walls of the mesh. The criterion is set by input 251 !! parameter "width". 252 !! 253 !! n : on output, the number (0<=n<=3) of "walls" of the mesh that 254 !! the point is too close to, in order to consider it belonging 255 !! to a mesh. 256 !! 257 !! So, if n > 0, the point is in the border. 258 ! ---------------------------------------------------------------------- 259 logical function mesh_inborder(mesh, geo, ip, dist, width) result(is_on_border) 260 type(mesh_t), intent(in) :: mesh !< the mesh 261 type(geometry_t), intent(in) :: geo 262 integer, intent(in) :: ip !< the point in the mesh 263 FLOAT, intent(in) :: width !< the width of the border 264 !> distance from border. The distances of the point to the walls, 265 !! for each of the walls that the point is too close to. 266 FLOAT, intent(out) :: dist 267 268 integer :: iatom, jatom, idir 269 FLOAT :: xx(MAX_DIM), rr, dd, radius 270 271 ! no PUSH SUB, called too often 272 273 is_on_border = .false. 274 dist = M_ZERO 275 276 select case(mesh%sb%box_shape) 277 case(SPHERE) 278 call mesh_r(mesh, ip, rr, coords=xx) 279 dd = rr - (mesh%sb%rsize - width) 280 if(dd > M_ZERO) then 281 is_on_border = .true. 282 dist = dd 283 end if 284 285 case(CYLINDER) 286 call mesh_r(mesh, ip, rr, coords=xx) 287 dd = sqrt(xx(2)**2 + xx(3)**2) - (mesh%sb%rsize - width) 288 if(dd > M_ZERO) then 289 is_on_border = .true. 290 dist = dd 291 end if 292 if(mesh%sb%periodic_dim == 0) then 293 dd = abs(xx(1)) - (mesh%sb%xsize - width) 294 if(dd > M_ZERO) then 295 is_on_border = .true. 296 dist = sqrt(dist*dist + dd*dd) 297 end if 298 end if 299 300 case(MINIMUM) 301 radius = mesh%sb%rsize 302 do iatom = 1, geo%natoms 303 call mesh_r(mesh, ip, rr, origin=geo%atom(iatom)%x, coords=xx) 304 if(mesh%sb%rsize < M_ZERO) radius = species_def_rsize(geo%atom(iatom)%species) 305 dd = rr - (radius - width) 306 ! check if the point is on the spherical shell of atom # iatom 307 if ((dd < M_ZERO) .or. (rr > radius)) cycle 308 309 ! make sure that the point is not inside some other atomic sphere 310 is_on_border = .true. 311 do jatom = 1, geo%natoms 312 if(jatom == iatom) cycle 313 call mesh_r(mesh, ip, rr, origin=geo%atom(jatom)%x) 314 if(mesh%sb%rsize < M_ZERO) radius = species_def_rsize(geo%atom(jatom)%species) 315 if(rr < radius - width) then 316 is_on_border = .false. 317 exit 318 end if 319 end do 320 321 if(is_on_border) dist = dd 322 end do 323 324 case(PARALLELEPIPED, HYPERCUBE) 325 call mesh_r(mesh, ip, rr, coords=xx) 326 do idir = mesh%sb%periodic_dim+1, mesh%sb%dim 327 dd = abs(xx(idir)) - (mesh%sb%lsize(idir) - width) 328 if(dd > M_ZERO) then 329 is_on_border = .true. 330 dist = dist + dd*dd 331 end if 332 end do 333 dist = sqrt(dist) 334 335 case(BOX_IMAGE, BOX_USDEF) 336 ! not implemented 337 dist = -1 338 339 end select 340 341 ! This may happen if the point is on more than one border at the same time. 342 if(dist > width) dist = width 343 344 end function mesh_inborder 345 346 347 !--------------------------------------------------------------------- 348 !> Returns the index of the point which is nearest to a given vector 349 !! position pos. Variable dmin will hold, on exit, the distance between 350 !! pos and this nearest mesh point. rankmin will be zero, if the mesh is 351 !! not partitioned, and the rank of the processor which holds the point 352 !! ind if the mesh is partitioned. 353 ! ---------------------------------------------------------------------- 354 integer function mesh_nearest_point(mesh, pos, dmin, rankmin) result(ind) 355 type(mesh_t), intent(in) :: mesh 356 FLOAT, intent(in) :: pos(MAX_DIM) 357 FLOAT, intent(out) :: dmin 358 integer, intent(out) :: rankmin 359 360 FLOAT :: dd 361 integer :: imin, ip 362#if defined(HAVE_MPI) 363 FLOAT :: min_loc_in(2), min_loc_out(2) 364#endif 365 366 PUSH_SUB(mesh_nearest_point) 367 368 !find the point of the grid that is closer to the atom 369 dmin = M_ZERO 370 do ip = 1, mesh%np 371 dd = sum((pos(1:mesh%sb%dim) - mesh%x(ip, 1:mesh%sb%dim))**2) 372 if((dd < dmin) .or. (ip == 1)) then 373 imin = ip 374 dmin = dd 375 end if 376 end do 377 378 rankmin = 0 379#if defined(HAVE_MPI) 380 if(mesh%parallel_in_domains) then 381 min_loc_in(1) = dmin 382 min_loc_in(2) = mesh%np_global * mesh%mpi_grp%rank + TOFLOAT(imin) 383 call MPI_Allreduce(min_loc_in, min_loc_out, 1, MPI_2FLOAT, & 384 MPI_MINLOC, mesh%mpi_grp%comm, mpi_err) 385 dmin = min_loc_out(1) 386 imin = mod(nint(min_loc_out(2)), mesh%np_global) 387 rankmin = nint(min_loc_out(2))/mesh%np_global 388 end if 389#endif 390 391 ind = imin 392 POP_SUB(mesh_nearest_point) 393 end function mesh_nearest_point 394 395 396 ! -------------------------------------------------------------- 397 subroutine mesh_nearest_point_infos(mesh, pos, dmin_global, rankmin, imin_local, imin_global) 398 type(mesh_t), intent(in) :: mesh 399 FLOAT, intent(in) :: pos(:) 400 FLOAT, intent(out) :: dmin_global 401 integer, intent(out) :: rankmin 402 integer, intent(out) :: imin_local 403 integer, intent(out) :: imin_global 404 405 integer :: ip, ip_global, idim, ipart 406 FLOAT :: dd, xx(3) 407 408 dmin_global = M_HUGE 409 if (mesh%parallel_in_domains) then 410 do ipart=1, mesh%vp%npart 411 do ip = 1, mesh%vp%np_local_vec(ipart) 412 ip_global = mesh%vp%local_vec(mesh%vp%xlocal_vec(ipart) + ip - 1) 413 do idim = 1, mesh%sb%dim 414 xx(idim) = mesh%idx%lxyz(ip_global,idim) * mesh%spacing(idim) 415 end do 416 dd = sqrt(sum((pos(1:3) - xx(1:3))**2)) 417 if (dd < dmin_global) then 418 imin_local = ip 419 rankmin = ipart-1 420 imin_global = ip_global 421 dmin_global = dd 422 end if 423 end do 424 end do 425 else 426 do ip = 1, mesh%np 427 do idim = 1, mesh%sb%dim 428 xx(idim) = mesh%idx%lxyz(ip,idim) * mesh%spacing(idim) 429 end do 430 dd = sqrt(sum((pos(1:3) - xx(1:3))**2)) 431 if (dd < dmin_global) then 432 imin_local = ip 433 rankmin = 0 434 imin_global = ip 435 dmin_global = dd 436 end if 437 end do 438 end if 439 440 end subroutine mesh_nearest_point_infos 441 442 443 ! -------------------------------------------------------------- 444 !> mesh_gcutoff returns the "natural" band limitation of the 445 !! grid mesh, in terms of the maximum G vector. For a cubic regular 446 !! grid, it is M_PI/spacing. 447 ! -------------------------------------------------------------- 448 FLOAT function mesh_gcutoff(mesh) result(gmax) 449 type(mesh_t), intent(in) :: mesh 450 451 PUSH_SUB(mesh_gcutoff) 452 gmax = M_PI / (maxval(mesh%spacing)) 453 454 POP_SUB(mesh_gcutoff) 455 end function mesh_gcutoff 456 457 458 ! -------------------------------------------------------------- 459 subroutine mesh_dump(this, dir, filename, mpi_grp, namespace, ierr) 460 class(mesh_t), intent(in) :: this 461 character(len=*), intent(in) :: dir 462 character(len=*), intent(in) :: filename 463 type(mpi_grp_t), intent(in) :: mpi_grp 464 type(namespace_t), intent(in) :: namespace 465 integer, intent(out) :: ierr 466 467 integer :: iunit, err 468 469 PUSH_SUB(mesh_dump) 470 471 ierr = 0 472 473 iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='write', & 474 position="append", die=.false., grp=mpi_grp) 475 if (iunit <= 0) then 476 ierr = ierr + 1 477 message(1) = "Unable to open file:" 478 message(2) = io_workpath(trim(dir)//"/"//trim(filename), namespace) 479 call messages_warning(2) 480 else 481 if (mpi_grp_is_root(mpi_grp)) then 482 write(iunit, '(a)') dump_tag 483 write(iunit, '(a20,1i10)') 'np_global= ', this%np_global 484 write(iunit, '(a20,1i10)') 'np_part_global= ', this%np_part_global 485 end if 486 call io_close(iunit, grp=mpi_grp) 487 end if 488 489 call index_dump(this%idx, dir, filename, mpi_grp, namespace, err) 490 if (err /= 0) ierr = ierr + 2 491 492 POP_SUB(mesh_dump) 493 end subroutine mesh_dump 494 495 496 ! -------------------------------------------------------------- 497 !> Read the mesh parameters from file that were written by mesh_dump. 498 subroutine mesh_load(this, dir, filename, mpi_grp, namespace, ierr) 499 class(mesh_t), intent(inout) :: this 500 character(len=*), intent(in) :: dir 501 character(len=*), intent(in) :: filename 502 type(mpi_grp_t), intent(in) :: mpi_grp 503 type(namespace_t), intent(in) :: namespace 504 integer, intent(out) :: ierr 505 506 integer :: iunit, err 507 character(len=20) :: str 508 character(len=100) :: lines(4) 509 510 PUSH_SUB(mesh_load) 511 512 ASSERT(this%sb%dim > 0 .and. this%sb%dim <= MAX_DIM) 513 514 ierr = 0 515 516 iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='read', & 517 status="old", die=.false., grp=mpi_grp) 518 if (iunit <= 0) then 519 ierr = ierr + 1 520 message(1) = "Unable to open file '"//trim(dir)//"/"//trim(filename)//"'." 521 call messages_warning(1) 522 else 523 ! Find the dump tag. 524 call iopar_find_line(mpi_grp, iunit, dump_tag, err) 525 if (err /= 0) ierr = ierr + 2 526 527 if (ierr == 0) then 528 call iopar_read(mpi_grp, iunit, lines, 2, err) 529 if (err /= 0) then 530 ierr = ierr + 4 531 else 532 read(lines(3), '(a20,1i10)') str, this%np_global 533 read(lines(4), '(a20,1i10)') str, this%np_part_global 534 this%parallel_in_domains = .false. 535 end if 536 end if 537 538 call io_close(iunit, grp=mpi_grp) 539 end if 540 541 call index_load(this%idx, dir, filename, mpi_grp, namespace, err) 542 if (err /= 0) ierr = ierr + 8 543 544 POP_SUB(mesh_load) 545 end subroutine mesh_load 546 547 548 ! -------------------------------------------------------------- 549 subroutine mesh_write_fingerprint(mesh, dir, filename, mpi_grp, namespace, ierr) 550 type(mesh_t), intent(in) :: mesh 551 character(len=*), intent(in) :: dir 552 character(len=*), intent(in) :: filename 553 type(mpi_grp_t), intent(in) :: mpi_grp 554 type(namespace_t),intent(in) :: namespace 555 integer, intent(out) :: ierr 556 557 integer :: iunit 558 559 PUSH_SUB(mesh_write_fingerprint) 560 561 ierr = 0 562 563 iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='write', & 564 die=.false., grp=mpi_grp) 565 if (iunit <= 0) then 566 message(1) = "Unable to open file '"//trim(dir)//"/"//trim(filename)//"'." 567 call messages_warning(1) 568 ierr = ierr + 1 569 else 570 if (mpi_grp_is_root(mpi_grp)) then 571 write(iunit, '(a20,i21)') 'box_shape = ', mesh%sb%box_shape 572 if(mesh%sb%box_shape /= HYPERCUBE) then 573 write(iunit, '(a20,i21)') 'np_part_global= ', mesh%np_part_global 574 write(iunit, '(a20,i21)') 'np_global= ', mesh%np_global 575 write(iunit, '(a20,i21)') 'algorithm= ', 1 576 write(iunit, '(a20,i21)') 'checksum= ', mesh%idx%checksum 577 end if 578 end if 579 call io_close(iunit, grp=mpi_grp) 580 end if 581 582 POP_SUB(mesh_write_fingerprint) 583 end subroutine mesh_write_fingerprint 584 585 586 ! ----------------------------------------------------------------------- 587 !> This function reads the fingerprint of a mesh written in 588 !! filename. If the meshes are equal (same fingerprint) return values 589 !! are 0, otherwise it returns the size of the mesh stored. 590 !! fingerprint cannot be read, it returns ierr /= 0. 591 subroutine mesh_read_fingerprint(mesh, dir, filename, mpi_grp, namespace, read_np_part, read_np, ierr) 592 type(mesh_t), intent(in) :: mesh 593 character(len=*), intent(in) :: dir 594 character(len=*), intent(in) :: filename 595 type(mpi_grp_t), intent(in) :: mpi_grp 596 type(namespace_t),intent(in) :: namespace 597 integer, intent(out) :: read_np_part 598 integer, intent(out) :: read_np 599 integer, intent(out) :: ierr 600 601 character(len=20) :: str 602 character(len=100) :: lines(4) 603 integer :: iunit, box_shape, algorithm, err 604 integer(8) :: checksum 605 606 PUSH_SUB(mesh_read_fingerprint) 607 608 ierr = 0 609 610 read_np_part = 0 611 read_np = 0 612 613 iunit = io_open(trim(dir)//"/"//trim(filename), namespace, action='read', & 614 status='old', die=.false., grp=mpi_grp) 615 if (iunit <= 0) then 616 ierr = ierr + 1 617 message(1) = "Unable to open file '"//trim(dir)//"/"//trim(filename)//"'." 618 call messages_warning(1) 619 else 620 box_shape = 0 621 call iopar_read(mpi_grp, iunit, lines, 1, err) 622 if (err /= 0) then 623 ierr = ierr + 2 624 else 625 read(lines(1), '(a20,i21)') str, box_shape 626 end if 627 628 if (box_shape == HYPERCUBE) then 629 ! We have a hypercube: we will assume everything is OK... 630 message(1) = "Simulation box is a hypercube: unable to check mesh compatibility." 631 call messages_warning(1) 632 633 else 634 call iopar_read(mpi_grp, iunit, lines, 4, err) 635 if (err /= 0) then 636 ierr = ierr + 4 637 else 638 read(lines(1), '(a20,i21)') str, read_np_part 639 read(lines(2), '(a20,i21)') str, read_np 640 read(lines(3), '(a20,i21)') str, algorithm 641 read(lines(4), '(a20,i21)') str, checksum 642 643 ASSERT(read_np_part >= read_np) 644 645 if (read_np_part == mesh%np_part_global & 646 .and. read_np == mesh%np_global & 647 .and. algorithm == 1 & 648 .and. checksum == mesh%idx%checksum) then 649 read_np_part = 0 650 read_np = 0 651 end if 652 end if 653 654 end if 655 656 call io_close(iunit, grp=mpi_grp) 657 end if 658 659 POP_SUB(mesh_read_fingerprint) 660 end subroutine mesh_read_fingerprint 661 662 ! --------------------------------------------------------- 663 subroutine mesh_check_dump_compatibility(mesh, dir, filename, namespace, mpi_grp, grid_changed, grid_reordered, map, ierr) 664 type(mesh_t), intent(in) :: mesh 665 character(len=*), intent(in) :: dir 666 character(len=*), intent(in) :: filename 667 type(namespace_t), intent(in) :: namespace 668 type(mpi_grp_t), intent(in) :: mpi_grp 669 logical, intent(out) :: grid_changed 670 logical, intent(out) :: grid_reordered 671 integer, pointer, intent(out) :: map(:) 672 integer, intent(out) :: ierr 673 674 integer :: ip, read_np_part, read_np, xx(MAX_DIM), err 675 integer, allocatable :: read_lxyz(:,:) 676 677 PUSH_SUB(mesh_check_dump_compatibility) 678 679 ierr = 0 680 681 nullify(map) 682 grid_changed = .false. 683 grid_reordered = .false. 684 685 ! Read the mesh fingerprint 686 call mesh_read_fingerprint(mesh, dir, filename, mpi_grp, namespace, read_np_part, read_np, err) 687 if (err /= 0) then 688 ierr = ierr + 1 689 message(1) = "Unable to read mesh fingerprint from '"//trim(dir)//"/"//trim(filename)//"'." 690 call messages_warning(1) 691 692 else if (read_np > 0) then 693 if (.not. associated(mesh%sb)) then 694 ! We can only check the compatibility of two meshes that have different fingerprints if we also 695 ! have the simulation box. In the case we do not, we will assume that the fingerprint is enough. 696 ierr = ierr + 2 697 else if (mesh%sb%box_shape /= HYPERCUBE) then 698 699 grid_changed = .true. 700 701 ! perhaps only the order of the points changed, this can only 702 ! happen if the number of points is the same and no points maps 703 ! to zero (this is checked below) 704 grid_reordered = (read_np == mesh%np_global) 705 706 ! the grid is different, so we read the coordinates. 707 SAFE_ALLOCATE(read_lxyz(1:read_np_part, 1:mesh%sb%dim)) 708 ASSERT(allocated(mesh%idx%lxyz)) 709 call io_binary_read(trim(io_workpath(dir, namespace))//'/lxyz.obf', read_np_part*mesh%sb%dim, read_lxyz, err) 710 if (err /= 0) then 711 ierr = ierr + 4 712 message(1) = "Unable to read index map from '"//trim(dir)//"'." 713 call messages_warning(1) 714 else 715 ! generate the map 716 SAFE_ALLOCATE(map(1:read_np)) 717 718 do ip = 1, read_np 719 xx = 0 720 xx(1:mesh%sb%dim) = read_lxyz(ip, 1:mesh%sb%dim) 721 if (any(xx(1:mesh%sb%dim) < mesh%idx%nr(1, 1:mesh%sb%dim)) .or. & 722 any(xx(1:mesh%sb%dim) > mesh%idx%nr(2, 1:mesh%sb%dim))) then 723 map(ip) = 0 724 grid_reordered = .false. 725 else 726 map(ip) = mesh%idx%lxyz_inv(xx(1), xx(2), xx(3)) 727 if(map(ip) > mesh%np_global) map(ip) = 0 728 end if 729 end do 730 end if 731 732 SAFE_DEALLOCATE_A(read_lxyz) 733 end if 734 end if 735 736 POP_SUB(mesh_check_dump_compatibility) 737 end subroutine mesh_check_dump_compatibility 738 739 740 ! -------------------------------------------------------------- 741 recursive subroutine mesh_end(this) 742 class(mesh_t), intent(inout) :: this 743 744 PUSH_SUB(mesh_end) 745 746 call mesh_cube_map_end(this%cube_map) 747 748 if(this%idx%is_hypercube) call hypercube_end(this%idx%hypercube) 749 750 SAFE_DEALLOCATE_A(this%resolution) 751 SAFE_DEALLOCATE_A(this%idx%lxyz) 752 SAFE_DEALLOCATE_A(this%idx%lxyz_inv) 753 SAFE_DEALLOCATE_A(this%x) 754 SAFE_DEALLOCATE_A(this%vol_pp) 755 756 if(this%parallel_in_domains) then 757#if defined(HAVE_MPI) 758 call vec_end(this%vp) 759 ! this is true if MeshUseTopology = false 760 if(this%mpi_grp%comm /= this%vp%comm) & 761 call MPI_Comm_free(this%vp%comm, mpi_err) 762 call partition_end(this%inner_partition) 763 call partition_end(this%bndry_partition) 764#endif 765 end if 766 767 POP_SUB(mesh_end) 768 end subroutine mesh_end 769 770 771 !> This function returns the point inside the grid corresponding to 772 !! a boundary point when PBCs are used. In case the point does not 773 !! have a correspondence (i.e. other BCs are used in that direction), 774 !! the same point is returned. Note that this function returns a 775 !! global point number when parallelization in domains is used. 776 ! --------------------------------------------------------- 777 integer function mesh_periodic_point(mesh, ip_global, ip_local) result(ipp) 778 type(mesh_t), intent(in) :: mesh 779 integer, intent(in) :: ip_global, ip_local 780 781 integer :: ix(MAX_DIM), nr(2, MAX_DIM), idim 782 FLOAT :: xx(MAX_DIM), rr, ufn_re, ufn_im 783 784 ! no push_sub, called too frequently 785 786 ix = mesh%idx%lxyz(ip_global, :) 787 nr(1, :) = mesh%idx%nr(1, :) + mesh%idx%enlarge(:) 788 nr(2, :) = mesh%idx%nr(2, :) - mesh%idx%enlarge(:) 789 790 do idim = 1, mesh%sb%periodic_dim 791 if(ix(idim) < nr(1, idim)) ix(idim) = ix(idim) + mesh%idx%ll(idim) 792 if(ix(idim) > nr(2, idim)) ix(idim) = ix(idim) - mesh%idx%ll(idim) 793 end do 794 795 ipp = mesh%idx%lxyz_inv(ix(1), ix(2), ix(3)) 796 797 if(mesh%masked_periodic_boundaries) then 798 call mesh_r(mesh, ip_local, rr, coords = xx) 799 call parse_expression(ufn_re, ufn_im, mesh%sb%dim, xx, rr, M_ZERO, mesh%periodic_boundary_mask) 800 if(int(ufn_re) == 0) ipp = ip_global ! Nothing will be done 801 end if 802 803 end function mesh_periodic_point 804 805 806 ! --------------------------------------------------------- 807 FLOAT pure function mesh_global_memory(mesh) result(memory) 808 type(mesh_t), intent(in) :: mesh 809 810 memory = M_ZERO 811 812 ! lxyz_inv 813 memory = memory + SIZEOF_UNSIGNED_INT * product(mesh%idx%nr(2, 1:mesh%sb%dim) - mesh%idx%nr(1, 1:mesh%sb%dim) + M_ONE) 814 ! resolution 815 if(mesh%sb%mr_flag) then 816 memory = memory + SIZEOF_UNSIGNED_INT * product(mesh%idx%nr(2, 1:mesh%sb%dim) - mesh%idx%nr(1, 1:mesh%sb%dim) + M_ONE) 817 end if 818 ! lxyz 819 memory = memory + SIZEOF_UNSIGNED_INT * TOFLOAT(mesh%np_part_global) * MAX_DIM 820 821 end function mesh_global_memory 822 823 824 ! --------------------------------------------------------- 825 FLOAT pure function mesh_local_memory(mesh) result(memory) 826 type(mesh_t), intent(in) :: mesh 827 828 memory = M_ZERO 829 830 ! x 831 memory = memory + REAL_PRECISION * TOFLOAT(mesh%np_part) * MAX_DIM 832 end function mesh_local_memory 833 834 835 ! --------------------------------------------------------- 836 function mesh_x_global(mesh, ip, force) result(xx) 837 type(mesh_t), intent(in) :: mesh 838 integer, intent(in) :: ip 839 logical, optional, intent(in) :: force 840 FLOAT :: xx(1:MAX_DIM) 841 842 FLOAT :: chi(1:MAX_DIM) 843 integer :: ix(1:MAX_DIM) 844 logical :: force_ 845 846! no push_sub because function is called too frequently 847 848 force_ = .false. 849 if (present(force)) force_ = force 850 851 if(mesh%parallel_in_domains .or. force_) then 852 call index_to_coords(mesh%idx, ip, ix) 853 chi(1:mesh%sb%dim) = ix(1:mesh%sb%dim) * mesh%spacing(1:mesh%sb%dim) 854 chi(mesh%sb%dim + 1:MAX_DIM) = M_ZERO 855 xx = M_ZERO ! this initialization is required by gfortran 4.4 or we get NaNs 856 call curvilinear_chi2x(mesh%sb, mesh%cv, chi, xx) 857 xx(mesh%sb%dim + 1:MAX_DIM) = M_ZERO 858 else 859 xx(1:MAX_DIM) = mesh%x(ip, 1:MAX_DIM) 860 end if 861 862 end function mesh_x_global 863 864 865 ! --------------------------------------------------------- 866 logical pure function mesh_compact_boundaries(mesh) result(cb) 867 type(mesh_t), intent(in) :: mesh 868 869 cb = .not. mesh%use_curvilinear .and. & 870 .not. mesh%parallel_in_domains .and. & 871 simul_box_has_zero_bc(mesh%sb) 872 873 end function mesh_compact_boundaries 874 875 876 ! --------------------------------------------------------- 877 subroutine mesh_check_symmetries(mesh, sb) 878 type(mesh_t), intent(in) :: mesh 879 type(simul_box_t), intent(in) :: sb 880 881 integer :: iop, ip, idim, nops 882 FLOAT :: destpoint(1:3), srcpoint(1:3), lsize(1:3), offset(1:3) 883 884 !If all the axis have the same spacing and the same length 885 !the grid is by obviously symmetric 886 !Indeed, reduced coordinates are proportional to the point index 887 !and the reduced rotation are integer matrices 888 !The result of the product is also proportional to an integer 889 !and therefore belong to the grid. 890 if(mesh%idx%ll(1) == mesh%idx%ll(2) .and. & 891 mesh%idx%ll(2) == mesh%idx%ll(3) .and. & 892 mesh%spacing(1) == mesh%spacing(2) .and. & 893 mesh%spacing(2) == mesh%spacing(3) ) return 894 895 PUSH_SUB(mesh_check_symmetries) 896 897 message(1) = "Checking if the real-space grid is symmetric"; 898 call messages_info(1) 899 900 lsize(1:3) = TOFLOAT(mesh%idx%ll(1:3)) 901 offset(1:3) = TOFLOAT(mesh%idx%nr(1, 1:3) + mesh%idx%enlarge(1:3)) 902 903 nops = symmetries_number(sb%symm) 904 905 do ip = 1, mesh%np 906 !We use floating point coordinates to check if the symmetric point 907 !belong to the grid. 908 !If yes, it should have integer reduced coordinates 909 if(mesh%parallel_in_domains) then 910 ! convert to global point 911 destpoint(1:3) = TOFLOAT(mesh%idx%lxyz(mesh%vp%local(mesh%vp%xlocal + ip - 1), 1:3)) - offset(1:3) 912 else 913 destpoint(1:3) = TOFLOAT(mesh%idx%lxyz(ip, 1:3)) - offset(1:3) 914 end if 915 ! offset moves corner of cell to origin, in integer mesh coordinates 916 ASSERT(all(destpoint >= 0)) 917 ASSERT(all(destpoint < lsize)) 918 919 ! move to center of cell in real coordinates 920 destpoint = destpoint - TOFLOAT(int(lsize)/2) 921 922 !convert to proper reduced coordinates 923 do idim = 1, 3 924 destpoint(idim) = destpoint(idim)/lsize(idim) 925 end do 926 927 ! iterate over all points that go to this point by a symmetry operation 928 do iop = 1, nops 929 srcpoint = symm_op_apply_red(sb%symm%ops(iop), destpoint) 930 931 !We now come back to what should be an integer, if the symmetric point beloings to the grid 932 do idim = 1, 3 933 srcpoint(idim) = srcpoint(idim)*lsize(idim) 934 end do 935 936 ! move back to reference to origin at corner of cell 937 srcpoint = srcpoint + TOFLOAT(int(lsize)/2) 938 939 ! apply periodic boundary conditions in periodic directions 940 do idim = 1, sb%periodic_dim 941 if(nint(srcpoint(idim)) < 0 .or. nint(srcpoint(idim)) >= lsize(idim)) then 942 srcpoint(idim) = modulo(srcpoint(idim)+M_HALF*SYMPREC, lsize(idim)) 943 end if 944 end do 945 ASSERT(all(srcpoint >= -SYMPREC)) 946 ASSERT(all(srcpoint < lsize)) 947 948 srcpoint(1:3) = srcpoint(1:3) + offset(1:3) 949 950 if(any(srcpoint-anint(srcpoint)> SYMPREC*M_TWO)) then 951 message(1) = "The real-space grid breaks at least one of the symmetries of the system." 952 message(2) = "Change your spacing or use SymmetrizeDensity=no." 953 call messages_fatal(2) 954 end if 955 end do 956 end do 957 958 POP_SUB(mesh_check_symmetries) 959 end subroutine 960 961end module mesh_oct_m 962 963 964!! Local Variables: 965!! mode: f90 966!! coding: utf-8 967!! End: 968