Octopus
submesh.F90
Go to the documentation of this file.
1!! Copyright (C) 2007 X. Andrade
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 submesh_oct_m
22 use accel_oct_m
23 use batch_oct_m
25 use box_oct_m
26 use comm_oct_m
27 use debug_oct_m
28 use global_oct_m
29 use index_oct_m
30 use, intrinsic :: iso_fortran_env
34 use sort_oct_m
35 use mesh_oct_m
37 use mpi_oct_m
40 use space_oct_m
41 use types_oct_m
42
43 implicit none
44 private
45
46 public :: &
47 submesh_t, &
60 dsm_nrm2, &
61 zsm_nrm2, &
79
82 type submesh_t
83 ! Components are public by default
84 real(real64), allocatable :: center(:)
85 real(real64) :: radius = m_zero
86 class(box_t), pointer :: box => null()
87 integer :: np = -1
88 integer, allocatable :: map(:)
89 integer :: num_regions
90 integer, allocatable :: regions(:)
91 type(accel_mem_t) :: buff_map
92 real(real64), allocatable :: rel_x(:,:)
93 real(real64), allocatable :: r(:)
94 type(mesh_t), pointer :: mesh => null()
95 logical :: overlap
97 integer :: np_global = -1
98 real(real64), allocatable :: rel_x_global(:,:)
99 integer, allocatable :: part_v(:)
100 integer, allocatable :: global2local(:)
101
102 type(mesh_cube_map_t) :: cube_map
103 end type submesh_t
104
105 interface submesh_add_to_mesh
107 end interface submesh_add_to_mesh
108
109 interface submesh_to_mesh_dotp
111 end interface submesh_to_mesh_dotp
112
113contains
115 ! -------------------------------------------------------------
116 ! Multipliers for recursive formulation of n-ellipsoid volume
117 ! simplifying the Gamma function
118 ! f(n) = 2f(n-2)/n, f(0)=1, f(1)=2
119 recursive real(real64) function f_n(dims) result(fn)
120 integer :: dims
121
122 if (dims == 0) then
123 fn = m_one
124 else if (dims == 1) then
125 fn = m_two
126 else
127 fn = m_two * f_n(dims - 2) / dims
128 end if
129
130 end function f_n
131
133 subroutine submesh_init_box(this, space, mesh, box, center)
134 type(submesh_t), intent(inout) :: this
135 class(space_t), intent(in) :: space
136 class(mesh_t), target, intent(in) :: mesh
137 class(box_t), target, intent(in) :: box
138 real(real64), intent(in) :: center(:)
139
140 integer :: original_mesh_index, submesh_index
141 logical :: submesh_points(1:mesh%np)
142 real(real64), allocatable :: xtmp(:, :), rtmp(:)
143
144 push_sub(submesh_init_box)
145
146 assert(.not. space%is_periodic())
147
148 ! Assign the mesh and the box
149 this%mesh => mesh
150 this%box => box
151
152 ! Get the total number of points inside the submesh. First of all get the number of points inside the submesh box
153 submesh_points = box%contains_points(mesh%np, mesh%x)
154 this%np = count(submesh_points)
155 ! Verify that the box of the submesh is contained into the box of the mesh
156 assert(this%np <= mesh%np)
157
158 safe_allocate(this%map(1:this%np))
159 safe_allocate(xtmp(1:space%dim, 1:this%np))
160 safe_allocate(rtmp(1:this%np))
161 xtmp = m_zero
162 rtmp = m_zero
163
164 ! Now we generate a map between the points of the submesh and the points of the original mesh
165 submesh_index = 0
166 do original_mesh_index = 1, mesh%np
167 ! check that the point is contained in the new box
168 if (submesh_points(original_mesh_index)) then
169 submesh_index = submesh_index + 1
170 this%map(submesh_index) = original_mesh_index
171 ! Assign the coordinates to this new point
172 xtmp(:, submesh_index) = mesh%x(original_mesh_index, :) - center(:)
173 rtmp(submesh_index) = norm2(xtmp(:, submesh_index))
174 end if
175 end do
176
177 call submesh_reorder_points(this, space, xtmp, rtmp)
179 safe_deallocate_a(xtmp)
180 safe_deallocate_a(rtmp)
183 end subroutine submesh_init_box
185 ! -------------------------------------------------------------
186 subroutine submesh_init(this, space, mesh, latt, center, rc)
187 type(submesh_t), intent(inout) :: this
188 class(space_t), intent(in) :: space
189 class(mesh_t), target, intent(in) :: mesh
190 type(lattice_vectors_t), intent(in) :: latt
191 real(real64), intent(in) :: center(1:space%dim)
192 real(real64), intent(in) :: rc
194 real(real64) :: r2, rc2, xx(space%dim), rc_norm_n
195 real(real64), allocatable :: center_copies(:,:), xtmp(:, :), rtmp(:), x(:,:)
196 integer :: icell, is, ip, ix, iy, iz
197 integer(int64) :: max_elements_count
198 type(lattice_iterator_t) :: latt_iter
199 integer, allocatable :: map_inv(:), map_temp(:)
200 integer :: nmax(3), nmin(3)
201 real(real64), parameter :: tol = 1e-13_real64
203
204 push_sub(submesh_init)
205 call profiling_in("SUBMESH_INIT")
206
207 assert(space%dim <= 3)
208
209 this%mesh => mesh
210
211 safe_allocate(this%center(1:space%dim))
212 this%center = center
213
214 this%radius = rc
215 ! We add a small number of avoid instabilities due to rounding errors
216 rc2 = rc**2 + tol
217
218 ! The spheres are generated differently for periodic coordinates,
219 ! mainly for performance reasons.
220 if (.not. space%is_periodic()) then
221
222 call profiling_in("SUBMESH_INIT_MAP_INV")
223 safe_allocate(map_inv(0:this%mesh%np))
224 map_inv(0:this%mesh%np) = 0
225
226 nmin = 0
227 nmax = 0
228
229 ! get a cube of points that contains the sphere
230 nmin(1:space%dim) = int((center(1:space%dim) - abs(rc))/mesh%spacing(1:space%dim)) - 1
231 nmax(1:space%dim) = int((center(1:space%dim) + abs(rc))/mesh%spacing(1:space%dim)) + 1
232
233 ! make sure that the cube is inside the grid
234 ! parts of the cube which would fall outside the simulation box are chopped off.
235 nmin(1:space%dim) = max(mesh%idx%nr(1, 1:space%dim), nmin(1:space%dim))
236 nmax(1:space%dim) = min(mesh%idx%nr(2, 1:space%dim), nmax(1:space%dim))
237
238 ! Get the total number of points inside the sphere
239 is = 0 ! this index counts inner points
240 do iz = nmin(3), nmax(3)
241 do iy = nmin(2), nmax(2)
242 do ix = nmin(1), nmax(1)
243 ip = mesh_local_index_from_coords(mesh, [ix, iy, iz])
244 if (ip == 0 .or. ip > mesh%np) cycle
245 r2 = sum((mesh%x(ip, :) - center)**2)
246 if (r2 <= rc2) then
247 is = is + 1
248 map_inv(ip) = is
249 end if
250 end do
251 end do
252 end do
253 this%np = is
254 call profiling_out("SUBMESH_INIT_MAP_INV")
255
256 call profiling_in("SUBMESH_INIT_RTMP")
257 safe_allocate(this%map(1:this%np))
258 safe_allocate(xtmp(1:space%dim, 1:this%np))
259 safe_allocate(rtmp(1:this%np))
260
261 ! Generate the table and the positions
262 do iz = nmin(3), nmax(3)
263 do iy = nmin(2), nmax(2)
264 do ix = nmin(1), nmax(1)
265 ip = mesh_local_index_from_coords(mesh, [ix, iy, iz])
266 if (ip == 0 .or. ip > mesh%np) cycle
267 is = map_inv(ip)
268 if (is == 0) cycle
269 this%map(is) = ip
270 xtmp(:, is) = mesh%x(ip, :) - center
271 rtmp(is) = norm2(xtmp(:,is))
272 end do
273 end do
274 end do
275
276 safe_deallocate_a(map_inv)
277 call profiling_out("SUBMESH_INIT_RTMP")
278
279 ! This is the case for a periodic system
280 else
281
282 ! Get the total number of points inside the sphere considering
283 ! replicas along PBCs
284
285 ! this requires some optimization
286
287 latt_iter = lattice_iterator_t(latt, rc)
288
289 safe_allocate(center_copies(1:space%dim, 1:latt_iter%n_cells))
290 do icell = 1, latt_iter%n_cells
291 center_copies(:, icell) = center + latt_iter%get(icell)
292 end do
293
294 !Recursive formulation for the volume of n-ellipsoid
295 !Garry Tee, NZ J. Mathematics Vol. 34 (2005) p. 165 eqs. 53,55
296 rc_norm_n = product(ceiling(rc / mesh%spacing(1:space%dim), int64) + m_one)
297 if (mesh%use_curvilinear) rc_norm_n = rc_norm_n / mesh%coord_system%min_mesh_scaling_product
298 max_elements_count = 3**space%dim * int(m_pi**floor(0.5 * space%dim) * rc_norm_n * f_n(space%dim), int64)
299
300 safe_allocate(x(1:space%dim,1:mesh%np))
301 do ip = 1, mesh%np
302 x(:,ip) = mesh%x(ip,:)
303 end do
304 call profiling_in("SUBMESH_INIT_PERIODIC_R")
305 safe_allocate(map_temp(1:max_elements_count))
306 safe_allocate(xtmp(1:space%dim, 1:max_elements_count))
307 safe_allocate(rtmp(1:max_elements_count))
308
309 is = 0
310 do ip = 1, mesh%np
311 do icell = 1, latt_iter%n_cells
312 xx = x(:,ip) - center_copies(:, icell)
313 if(any(abs(xx)>rc+tol)) cycle
314 r2 = sum(xx**2)
315 if (r2 > rc2) cycle
316 is = is + 1
317 map_temp(is) = ip
318 rtmp(is) = sqrt(r2)
319 xtmp(:, is) = xx
320 ! Note that xx can be outside the unit cell
321 end do
322 end do
323 assert(is < huge(is))
324 this%np = is
325
326 safe_allocate(this%map(1:this%np))
327 this%map(1:this%np) = map_temp(1:this%np)
328 call profiling_out("SUBMESH_INIT_PERIODIC_R")
329
330 safe_deallocate_a(map_temp)
331 safe_deallocate_a(center_copies)
332 safe_deallocate_a(x)
333
334 end if
335
336 call submesh_reorder_points(this, space, xtmp, rtmp)
337 call profiling_out("SUBMESH_INIT")
338
339 safe_deallocate_a(xtmp)
340 safe_deallocate_a(rtmp)
341
342 pop_sub(submesh_init)
343 end subroutine submesh_init
344
345 subroutine submesh_reorder_points(this, space, xtmp, rtmp)
346 type(submesh_t), intent(inout) :: this
347 class(space_t), intent(in) :: space
348 real(real64), intent(in) :: xtmp(:, :), rtmp(:)
349
350 integer :: ip, i_region, offset
351 integer, allocatable :: order(:), order_new(:)
352 integer, allocatable :: map_new(:)
353 integer, allocatable :: np_region(:), tmp_array(:)
354
355 push_sub(submesh_reorder_points)
356
357 ! now order points for better locality
358
359 call profiling_in("SUBMESH_INIT_ORDER")
360 safe_allocate(order(1:this%np))
361 safe_allocate(this%rel_x(1:space%dim, 1:this%np))
362 safe_allocate(this%r(1:this%np))
363
364 do ip = 1, this%np
365 order(ip) = ip
366 end do
367
368 ! First we just reorder in order to determine overlap:
369 call sort(this%map, order)
370
371 !check whether points overlap (i.e. whether a submesh contains the same point more than once)
372 this%overlap = .false.
373 do ip = 1, this%np - 1
374 if (this%map(ip) == this%map(ip + 1)) then
375 ! this simplified test works, as the points are ordered.
376 this%overlap = .true.
377 exit
378 end if
379 end do
380
381 this%num_regions = 1
382 call profiling_out("SUBMESH_INIT_ORDER")
383
384 if(this%overlap) then
385 call profiling_in("SUBMESH_INIT_OVERLAP")
386 !disentangle the map into injective regions
387
388 safe_allocate(tmp_array(1:this%np))
389 safe_allocate(order_new(1:this%np))
390 safe_allocate(np_region(1:this%np))
391 safe_allocate(map_new( 1:this%np))
392
393 np_region(1) = 1
394 tmp_array(1) = 1
395 i_region = 1
396
397 do ip = 2, this%np
398 if (this%map(ip) == this%map(ip - 1)) then
399 i_region = i_region + 1
400 if (i_region > this%num_regions) then
401 this%num_regions = i_region
402 np_region(i_region) = 0
403 end if
404 else
405 i_region = 1
406 end if
407 tmp_array(ip) = i_region ! which region does ip belong to
408 np_region(i_region) = np_region(i_region) + 1 ! increase number of points in i_region
409 end do
410
411 assert( .not. allocated(this%regions))
412
413 ! construct array of offsets
414 safe_allocate(this%regions(1:this%num_regions+1))
415
416 this%regions(1) = 1
417
418 if(this%num_regions > 1) then
419 do i_region = 1, this%num_regions
420 this%regions(i_region + 1) = this%regions(i_region) + np_region(i_region)
421 end do
422 else
423 this%regions(2) = this%np + 1
424 end if
425
426 np_region(1:this%np) = 0
427 order_new(1:this%np) = -1
428 map_new(1:this%np) = -1
429
430
431 !reassemble regions into global map array
432 do ip = 1, this%np
433 i_region = tmp_array(ip)
434 np_region(i_region) = np_region(i_region) + 1
435 offset = this%regions(i_region) - 1
436 map_new( offset + np_region(i_region) ) = this%map(ip)
437 order_new( offset + np_region(i_region) ) = order(ip)
438 end do
439
440 order(1:this%np) = order_new(1:this%np)
441 this%map(1:this%np) = map_new(1:this%np)
442
443 safe_deallocate_a(tmp_array)
444 safe_deallocate_a(order_new)
445 safe_deallocate_a(np_region)
446 safe_deallocate_a(map_new)
447 call profiling_out("SUBMESH_INIT_OVERLAP")
448
449 else
450 this%num_regions = 1
451 safe_allocate(this%regions(1:2))
452 this%regions(1) = 1
453 this%regions(2) = this%np + 1
454 end if
455
456 ! Lastly, reorder the points according to the new scheme
457 do ip = 1, this%np
458 this%rel_x(:, ip) = xtmp(:, order(ip))
459 this%r(ip) = rtmp(order(ip))
460 end do
461
462 safe_deallocate_a(order)
463
465 end subroutine submesh_reorder_points
466
467
468 ! --------------------------------------------------------------
469 !This routine takes two submeshes and merge them into a bigger submesh
470 !The grid is centered on the first center
471 subroutine submesh_merge(this, space, mesh, sm1, sm2, shift)
472 type(submesh_t), intent(inout) :: this
473 class(space_t), intent(in) :: space
474 class(mesh_t), target, intent(in) :: mesh
475 type(submesh_t), intent(in) :: sm1
476 type(submesh_t), intent(in) :: sm2
477 real(real64), optional, intent(in) :: shift(:)
478
479 real(real64) :: r2
480 integer :: ip, is
481 real(real64) :: xx(space%dim), diff_centers(space%dim)
482
483 push_sub(submesh_merge)
484 call profiling_in("SUBMESH_MERGE")
485
486 this%mesh => mesh
487
488 safe_allocate(this%center(1:space%dim))
489 this%center = sm1%center
490 this%radius = sm1%radius
491
492 ! This is a quick fix to prevent uninitialized variables. To properly check the self-overlap,
493 ! a similar approach as in submesh_init should be taken with respect to the merged map.
494 this%overlap = sm1%overlap .or. sm2%overlap
495
496 diff_centers = sm1%center - sm2%center
497 if (present(shift)) diff_centers = diff_centers - shift
498
499 !As we take the union of the two submeshes, we know that we have all the points from the first one included.
500 !The extra points from the second submesh are those which are not included in the first one
501 is = sm1%np
502 do ip = 1, sm2%np
503 !sm2%x contains points coordinates defined with respect to sm2%center
504 xx = sm2%rel_x(:, ip) - diff_centers
505 !If the point is not in sm1, we add it
506 if (sum(xx**2) > sm1%radius**2) is = is + 1
507 end do
508
509 this%np = is
510
511 safe_allocate(this%map(1:this%np))
512 safe_allocate(this%rel_x(1:space%dim, 1:this%np))
513 safe_allocate(this%r(1:this%np))
514 this%map(1:sm1%np) = sm1%map(1:sm1%np)
515 this%rel_x(:, 1:sm1%np) = sm1%rel_x(:, 1:sm1%np)
516 this%r(1:sm1%np) = sm1%r(1:sm1%np)
517
518 !iterate again to fill the tables
519 is = sm1%np
520 do ip = 1, sm2%np
521 xx = sm2%rel_x(:, ip) - diff_centers
522 r2 = sum(xx**2)
523 if (r2 > sm1%radius**2) then
524 is = is + 1
525 this%map(is) = sm2%map(ip)
526 this%r(is) = sqrt(r2)
527 this%rel_x(:, is) = xx
528 end if
529 end do
530
531 call profiling_out("SUBMESH_MERGE")
532 pop_sub(submesh_merge)
533 end subroutine submesh_merge
534
535 ! --------------------------------------------------------------
536 !This routine shifts the center of a submesh, without changing the grid points
537 subroutine submesh_shift_center(this, space, newcenter)
538 type(submesh_t), intent(inout) :: this
539 class(space_t), intent(in) :: space
540 real(real64), intent(in) :: newcenter(:)
541
542 real(real64) :: xx(space%dim), diff_centers(space%dim), oldcenter(space%dim)
543 integer :: ip
544
545 push_sub(submesh_shift_center)
546 call profiling_in("SUBMESH_SHIFT")
547
548 oldcenter = this%center
549 this%center = newcenter
550
551 diff_centers = newcenter - oldcenter
552
553 do ip = 1, this%np
554 xx = this%rel_x(:, ip) - diff_centers
555 this%r(ip) = norm2(xx)
556 this%rel_x(:, ip) = xx
557 end do
558
559 call profiling_out("SUBMESH_SHIFT")
560 pop_sub(submesh_shift_center)
561 end subroutine submesh_shift_center
562
563 ! --------------------------------------------------------------
564 subroutine submesh_broadcast(this, space, mesh, center, radius, root, mpi_grp)
565 type(submesh_t), intent(inout) :: this
566 class(space_t), intent(in) :: space
567 type(mesh_t), target, intent(in) :: mesh
568 real(real64), intent(in) :: center(1:space%dim)
569 real(real64), intent(in) :: radius
570 integer, intent(in) :: root
571 type(mpi_grp_t), intent(in) :: mpi_grp
572
573 integer :: nparray(1:3)
574
575 push_sub(submesh_broadcast)
576 call profiling_in('SUBMESH_BCAST')
577
578 if (root /= mpi_grp%rank) then
579 this%mesh => mesh
580 this%center = center
581 this%radius = radius
582 end if
583
584 if (mpi_grp%size > 1) then
585
586 if (root == mpi_grp%rank) then
587 nparray(1) = this%np
588 nparray(2) = this%num_regions
589 if (this%overlap) then
590 nparray(3) = 1
591 else
592 nparray(3) = 0
593 end if
594 end if
595
596 call mpi_grp%bcast(nparray, 2, mpi_integer, root)
597 this%np = nparray(1)
598 this%num_regions = nparray(2)
599 this%overlap = (nparray(3) == 1)
600
601 if (root /= mpi_grp%rank) then
602 safe_allocate(this%map(1:this%np))
603 safe_allocate(this%rel_x(1:space%dim, 1:this%np))
604 safe_allocate(this%r(1:this%np))
605 safe_allocate(this%regions(1:this%num_regions+1))
606 end if
607
608 call mpi_grp%bcast(this%regions(1), this%num_regions+1, mpi_integer, root)
609
610 if (this%np > 0) then
611 call mpi_grp%bcast(this%map(1), this%np, mpi_integer, root)
612 call mpi_grp%bcast(this%rel_x(1, 1), this%np*space%dim, mpi_double_precision, root)
613 call mpi_grp%bcast(this%r(1), this%np, mpi_double_precision, root)
614 end if
615
616 end if
617
618 call profiling_out('SUBMESH_BCAST')
619 pop_sub(submesh_broadcast)
620 end subroutine submesh_broadcast
621
622 ! --------------------------------------------------------------
623 logical function submesh_compatible(this, radius, center, dx) result(compatible)
624 type(submesh_t), intent(in) :: this
625 real(real64), intent(in) :: radius
626 real(real64), intent(in) :: center(:)
627 real(real64), intent(in) :: dx
628
629 compatible =.false.
630 if (allocated(this%center)) then
633 if (radius <= this%radius+dx*1e-6_real64 .and. all(abs(this%center - center) < 0.25*dx)) then
634 compatible = .true.
635 end if
636 end if
637
638 end function submesh_compatible
639
640 ! --------------------------------------------------------------
641 subroutine submesh_end(this)
642 type(submesh_t), intent(inout) :: this
643
644 push_sub(submesh_end)
645
646 if (this%np /= -1) then
647 nullify(this%mesh)
648 this%np = -1
649 safe_deallocate_a(this%center)
650 safe_deallocate_a(this%map)
651 safe_deallocate_a(this%rel_x)
652 safe_deallocate_a(this%r)
653 safe_deallocate_a(this%regions)
654 end if
655
656 if (accel_is_enabled()) then
657 call accel_release_buffer(this%buff_map)
658 end if
659
660 pop_sub(submesh_end)
661 end subroutine submesh_end
662
663 ! --------------------------------------------------------------
664
665 subroutine submesh_get_inv(this, map_inv)
666 type(submesh_t), intent(in) :: this
667 integer, intent(out) :: map_inv(:)
668
669 integer :: is
670
671 push_sub(submesh_get_inv)
672
673 map_inv(1:this%mesh%np_part) = 0
674 do is = 1, this%np
675 map_inv(this%map(is)) = is
676 end do
677
678 pop_sub(submesh_get_inv)
679 end subroutine submesh_get_inv
680
681 ! --------------------------------------------------------------
682 logical function submesh_overlap(sm1, sm2, space) result(overlap)
683 type(submesh_t), intent(in) :: sm1
684 type(submesh_t), intent(in) :: sm2
685 class(space_t), intent(in) :: space
686
687 integer :: ii, jj, dd
688 real(real64) :: distance
689
690 !no PUSH_SUB, called too often
691
692 if (.not. space%is_periodic()) then
693 !first check the distance
694 distance = sum((sm1%center - sm2%center)**2)
695 overlap = distance <= (1.5_real64*(sm1%radius + sm2%radius))**2
696
697 ! if they are very far, no need to check in detail
698 if (.not. overlap) return
699 end if
700
701 ! Otherwise check whether they have the some point in common. We
702 ! can make the comparison faster using that the arrays are sorted.
703 overlap = .false.
704 ii = 1
705 jj = 1
706 do while(ii <= sm1%np .and. jj <= sm2%np)
707 dd = sm1%map(ii) - sm2%map(jj)
708 if (dd < 0) then
709 ii = ii + 1
710 else if (dd > 0) then
711 jj = jj + 1
712 else
713 overlap = .true.
714 exit
715 end if
716 end do
717
718 if (sm1%mesh%parallel_in_domains) then
719 call sm1%mesh%mpi_grp%allreduce_inplace(overlap, 1, mpi_logical, mpi_lor)
720 end if
721
722 end function submesh_overlap
723
724 ! -------------------------------------------------------------
725 subroutine submesh_build_global(this, space)
726 type(submesh_t), intent(inout) :: this
727 class(space_t), intent(in) :: space
728
729 integer, allocatable :: part_np(:)
730 integer :: ipart, ind, ip
731
732 push_sub(submesh_build_global)
733
734 if (.not. this%mesh%parallel_in_domains) then
735 this%np_global = this%np
736 pop_sub(submesh_build_global)
737 return
738 end if
739
740 safe_allocate(part_np(this%mesh%pv%npart))
741 part_np = 0
742 part_np(this%mesh%pv%partno) = this%np
743
744 call this%mesh%allreduce(part_np)
745 this%np_global = sum(part_np)
746
747 safe_allocate(this%rel_x_global(1:space%dim, 1:this%np_global))
748 safe_allocate(this%part_v(1:this%np_global))
749 safe_allocate(this%global2local(1:this%np_global))
750 this%rel_x_global(1:space%dim, 1:this%np_global) = m_zero
751 this%part_v(1:this%np_global) = 0
752 this%global2local(1:this%np_global) = 0
753
754 ind = 0
755 do ipart = 1, this%mesh%pv%npart
756 if (ipart == this%mesh%pv%partno) then
757 do ip = 1, this%np
758 this%rel_x_global(:, ind + ip) = this%rel_x(:, ip)
759 this%part_v(ind + ip) = this%mesh%pv%partno
760 this%global2local(ind + ip) = ip
761 end do
762 end if
763 ind = ind + part_np(ipart)
764 end do
765
766 call this%mesh%allreduce(this%rel_x_global)
767 call this%mesh%allreduce(this%part_v)
768 call this%mesh%allreduce(this%global2local)
769
770 safe_deallocate_a(part_np)
771
772 pop_sub(submesh_build_global)
773 end subroutine submesh_build_global
774
775 ! -----------------------------------------------------------
776 subroutine submesh_end_global(this)
777 type(submesh_t), intent(inout) :: this
778
779 push_sub(submesh_end_global)
780
781 safe_deallocate_a(this%rel_x_global)
782 this%np_global = -1
783 safe_deallocate_a(this%part_v)
784 safe_deallocate_a(this%global2local)
785
786 pop_sub(submesh_end_global)
787 end subroutine submesh_end_global
788
789
790 ! -----------------------------------------------------------
791 subroutine zzsubmesh_add_to_mesh(this, sphi, phi, factor)
792 type(submesh_t), intent(in) :: this
793 complex(real64), intent(in) :: sphi(:)
794 complex(real64), intent(inout) :: phi(:)
795 complex(real64), optional, intent(in) :: factor
796
797 integer :: ip, m
798
799 push_sub(zzdsubmesh_add_to_mesh)
800
801 if (present(factor)) then
802 !Loop unrolling inspired by BLAS axpy routine
803 m = mod(this%np,4)
804 do ip = 1, m
805 phi(this%map(ip)) = phi(this%map(ip)) + factor*sphi(ip)
806 end do
807 if (this%np.ge.4) then
808 do ip = m+1, this%np, 4
809 phi(this%map(ip)) = phi(this%map(ip)) + factor*sphi(ip)
810 phi(this%map(ip+1)) = phi(this%map(ip+1)) + factor*sphi(ip+1)
811 phi(this%map(ip+2)) = phi(this%map(ip+2)) + factor*sphi(ip+2)
812 phi(this%map(ip+3)) = phi(this%map(ip+3)) + factor*sphi(ip+3)
813 end do
814 end if
815 else
816 m = mod(this%np,4)
817 do ip = 1, m
818 phi(this%map(ip)) = phi(this%map(ip)) + sphi(ip)
819 end do
820 if (this%np.ge.4) then
821 do ip = m+1, this%np, 4
822 phi(this%map(ip)) = phi(this%map(ip)) + sphi(ip)
823 phi(this%map(ip+1)) = phi(this%map(ip+1)) + sphi(ip+1)
824 phi(this%map(ip+2)) = phi(this%map(ip+2)) + sphi(ip+2)
825 phi(this%map(ip+3)) = phi(this%map(ip+3)) + sphi(ip+3)
826 end do
827 end if
828 end if
829
830 pop_sub(zzdsubmesh_add_to_mesh)
831 end subroutine zzsubmesh_add_to_mesh
832
833 !------------------------------------------------------------
834 complex(real64) function zzsubmesh_to_mesh_dotp(this, sphi, phi, reduce) result(dotp)
835 type(submesh_t), intent(in) :: this
836 complex(real64), intent(in) :: sphi(:)
837 complex(real64), intent(in) :: phi(:)
838 logical, optional, intent(in) :: reduce
839
840 integer :: is, m, ip
841
842 push_sub(zzsubmesh_to_mesh_dotp)
843
844 dotp = m_z0
845
846 if (this%mesh%use_curvilinear) then
847 do is = 1, this%np
848 dotp = dotp + this%mesh%vol_pp(this%map(is))*phi(this%map(is))*conjg(sphi(is))
849 end do
850 else
851 m = mod(this%np,4)
852 do ip = 1, m
853 dotp = dotp + phi(this%map(ip))*conjg(sphi(ip))
854 end do
855 if (this%np.ge.4) then
856 do ip = m+1, this%np, 4
857 dotp = dotp + phi(this%map(ip))*conjg(sphi(ip)) &
858 + phi(this%map(ip+1))*conjg(sphi(ip+1)) &
859 + phi(this%map(ip+2))*conjg(sphi(ip+2)) &
860 + phi(this%map(ip+3))*conjg(sphi(ip+3))
861 end do
862 end if
863 dotp = dotp*this%mesh%vol_pp(1)
864 end if
865
866 if (optional_default(reduce, .true.) .and. this%mesh%parallel_in_domains) then
867 call profiling_in("SM_REDUCE_DOTP")
868 call this%mesh%allreduce(dotp)
869 call profiling_out("SM_REDUCE_DOTP")
870 end if
871
873 end function zzsubmesh_to_mesh_dotp
874
875 !------------------------------------------------------------
877 subroutine submesh_get_cube_dim(sm, space, db)
878 type(submesh_t), target, intent(in) :: sm
879 class(space_t), intent(in) :: space
880 integer, intent(out) :: db(1:space%dim)
881
882 integer :: ip, idir
883 real(real64) :: chi(space%dim), max_chi(space%dim)
884 integer :: db_red(1:space%dim)
885 real(real64), parameter :: tol=1.0e-10_real64
886
887 push_sub(submesh_get_cube_dim)
888
889 max_chi = abs(sm%mesh%coord_system%from_cartesian(sm%rel_x(:, 1)))/sm%mesh%spacing
890 do ip = 2, sm%np
891 chi = sm%mesh%coord_system%from_cartesian(sm%rel_x(:, ip))
892 do idir = 1, space%dim
893 max_chi(idir) = max(max_chi(idir), abs(chi(idir))/sm%mesh%spacing(idir))
894 end do
895 end do
896
897 do idir = 1, space%dim
898 db(idir) = nint(max_chi(idir)-tol)
899 end do
900
901 if(sm%mesh%parallel_in_domains) then
902 call sm%mesh%mpi_grp%allreduce(db(1), db_red(1), space%dim, mpi_integer, mpi_max)
903 db(1:space%dim) = db_red(1:space%dim)
904 end if
905
906 db = 2 * db + 1
907
908 pop_sub(submesh_get_cube_dim)
909 end subroutine submesh_get_cube_dim
910
911 !------------------------------------------------------------
912 subroutine submesh_init_cube_map(sm, space)
913 type(submesh_t), target, intent(inout) :: sm
914 class(space_t), intent(in) :: space
915
916 integer(int64) :: ip
917 integer ::idir
918 real(real64) :: chi(space%dim), shift(space%dim)
919
920 push_sub(submesh_init_cube_map)
921
922 sm%cube_map%nmap = sm%np
923
924 safe_allocate(sm%cube_map%map(1:space%dim, 1:sm%cube_map%nmap))
925
926 !The center of the submesh does not belong to the mesh
927 !So we first need to find the closest grid point, and center the cube to it
928 chi = sm%mesh%coord_system%from_cartesian(sm%center)
929 do idir = 1, space%dim
930 shift(idir) = nint(chi(idir)/sm%mesh%spacing(idir))*sm%mesh%spacing(idir)
931 end do
932 shift = sm%mesh%coord_system%to_cartesian(shift)
933 shift = shift - sm%center
934
935 do ip = 1, sm%cube_map%nmap
936 chi = sm%mesh%coord_system%from_cartesian(sm%rel_x(:,ip) - shift)
937 do idir = 1, space%dim
938 sm%cube_map%map(idir, ip) = nint(chi(idir)/sm%mesh%spacing(idir))
939 end do
940 end do
941
942 if (accel_is_enabled()) then
943 call accel_create_buffer(sm%cube_map%map_buffer, accel_mem_read_only, type_integer, sm%cube_map%nmap*space%dim)
944 call accel_write_buffer(sm%cube_map%map_buffer, sm%cube_map%nmap*space%dim, sm%cube_map%map)
945 end if
946
947 pop_sub(submesh_init_cube_map)
948 end subroutine submesh_init_cube_map
949
950 !------------------------------------------------------------
951 subroutine submesh_end_cube_map(sm)
952 type(submesh_t), intent(inout) :: sm
953
954 push_sub(submesh_end_cube_map)
955
956 call mesh_cube_map_end(sm%cube_map)
957
958 pop_sub(submesh_end_cube_map)
959 end subroutine submesh_end_cube_map
960
967 subroutine dzsubmesh_batch_add(this, ss, mm)
968 type(submesh_t), intent(in) :: this
969 class(batch_t), intent(in) :: ss
970 class(batch_t), intent(inout) :: mm
971
972 integer :: ist, idim, jdim, is
973
974 push_sub(dzsubmesh_batch_add)
975
976 assert(.not. mm%is_packed())
977 assert(ss%type() == type_float)
978 assert(mm%type() == type_cmplx)
979 assert(ss%nst_linear == mm%nst_linear)
980 assert(ss%status() == mm%status())
981 assert(ss%dim == mm%dim)
982
983 assert(mm%nst == ss%nst)
984
985 !$omp parallel do private(ist, idim, jdim, is) if(.not. this%overlap)
986 do ist = 1, mm%nst
987 do idim = 1, mm%dim
988 jdim = min(idim, ss%dim)
989
990 do is = 1, this%np
991 mm%zff(this%map(is), idim, ist) = &
992 mm%zff(this%map(is), idim, ist) + ss%dff(is, jdim, ist)
993 end do
994
995 end do
996 end do
997 !$omp end parallel do
998
999 pop_sub(dzsubmesh_batch_add)
1000 end subroutine dzsubmesh_batch_add
1001
1002
1003#include "undef.F90"
1004#include "real.F90"
1005#include "submesh_inc.F90"
1006
1007#include "undef.F90"
1008#include "complex.F90"
1009#include "submesh_inc.F90"
1010
1011end module submesh_oct_m
1012
1013!! Local Variables:
1014!! mode: f90
1015!! coding: utf-8
1016!! End:
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:149
double sqrt(double __x) __attribute__((__nothrow__
double floor(double __x) __attribute__((__nothrow__
subroutine, public accel_release_buffer(this)
Definition: accel.F90:1246
pure logical function, public accel_is_enabled()
Definition: accel.F90:400
integer, parameter, public accel_mem_read_only
Definition: accel.F90:183
This module implements batches of mesh functions.
Definition: batch.F90:133
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:122
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
complex(real64), parameter, public m_z0
Definition: global.F90:198
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the index, used for the mesh points.
Definition: index.F90:122
subroutine, public mesh_cube_map_end(this)
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
integer function, public mesh_local_index_from_coords(mesh, ix)
This function returns the local index of the point for a given vector of integer coordinates.
Definition: mesh.F90:935
Some general things and nomenclature:
Definition: par_vec.F90:171
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
real(real64) function, public dsm_integrate_frommesh(mesh, sm, ff, reduce)
Definition: submesh.F90:1205
subroutine zdsubmesh_add_to_mesh(this, sphi, phi, factor)
Definition: submesh.F90:1775
subroutine, public submesh_init_cube_map(sm, space)
Definition: submesh.F90:1006
logical function, public submesh_overlap(sm1, sm2, space)
Definition: submesh.F90:776
subroutine, public zsubmesh_batch_dotp_matrix(this, mm, ss, dot, reduce)
Definition: submesh.F90:2080
real(real64) function, public dsm_nrm2(sm, ff, reduce)
this function returns the the norm of a vector
Definition: submesh.F90:1316
subroutine zzsubmesh_add_to_mesh(this, sphi, phi, factor)
Definition: submesh.F90:885
subroutine, public submesh_end_global(this)
Definition: submesh.F90:870
subroutine, public zsubmesh_copy_from_mesh(this, phi, sphi, conjugate)
Definition: submesh.F90:1823
subroutine submesh_reorder_points(this, space, xtmp, rtmp)
Definition: submesh.F90:439
complex(real64) function, public zsm_integrate(mesh, sm, ff, reduce)
Definition: submesh.F90:1702
subroutine, public submesh_shift_center(this, space, newcenter)
Definition: submesh.F90:631
subroutine, public submesh_init_box(this, space, mesh, box, center)
This subroutine creates a submesh which has the same box of the one that the user passes.
Definition: submesh.F90:227
complex(real64) function zzsubmesh_to_mesh_dotp(this, sphi, phi, reduce)
Definition: submesh.F90:928
real(real64) function, public dsm_integrate(mesh, sm, ff, reduce)
Definition: submesh.F90:1165
subroutine, public dsubmesh_batch_add_matrix(this, factor, ss, mm)
The following functions takes a batch of functions defined in submesh (ss) and adds all of them to ea...
Definition: submesh.F90:1398
recursive real(real64) function f_n(dims)
Definition: submesh.F90:213
subroutine, public dsubmesh_copy_from_mesh(this, phi, sphi, conjugate)
Definition: submesh.F90:1286
real(real64) function ddsubmesh_to_mesh_dotp(this, sphi, phi, reduce)
Definition: submesh.F90:1353
subroutine, public submesh_merge(this, space, mesh, sm1, sm2, shift)
Definition: submesh.F90:565
subroutine, public dzsubmesh_batch_add(this, ss, mm)
The following function takes a batch of functions defined in submesh (ss) and adds one of them to eac...
Definition: submesh.F90:1061
subroutine, public zsubmesh_batch_add(this, ss, mm)
The following function takes a batch of functions defined in submesh (ss) and adds one of them to eac...
Definition: submesh.F90:2028
real(real64) function, public zsm_nrm2(sm, ff, reduce)
this function returns the the norm of a vector
Definition: submesh.F90:1853
subroutine, public zsubmesh_batch_add_matrix(this, factor, ss, mm)
The following functions takes a batch of functions defined in submesh (ss) and adds all of them to ea...
Definition: submesh.F90:1935
subroutine, public submesh_end_cube_map(sm)
Definition: submesh.F90:1045
subroutine, public submesh_end(this)
Definition: submesh.F90:735
subroutine ddsubmesh_add_to_mesh(this, sphi, phi, factor)
Definition: submesh.F90:1238
subroutine, public submesh_get_cube_dim(sm, space, db)
finds the dimension of a box containing the submesh
Definition: submesh.F90:971
subroutine, public dsubmesh_batch_dotp_matrix(this, mm, ss, dot, reduce)
Definition: submesh.F90:1543
complex(real64) function, public zsm_integrate_frommesh(mesh, sm, ff, reduce)
Definition: submesh.F90:1742
logical function, public submesh_compatible(this, radius, center, dx)
Definition: submesh.F90:717
subroutine, public submesh_build_global(this, space)
Definition: submesh.F90:819
subroutine, public submesh_broadcast(this, space, mesh, center, radius, root, mpi_grp)
Definition: submesh.F90:658
subroutine, public submesh_get_inv(this, map_inv)
Definition: submesh.F90:759
subroutine, public dsubmesh_batch_add(this, ss, mm)
The following function takes a batch of functions defined in submesh (ss) and adds one of them to eac...
Definition: submesh.F90:1491
subroutine, public submesh_init(this, space, mesh, latt, center, rc)
Definition: submesh.F90:280
complex(real64) function zdsubmesh_to_mesh_dotp(this, sphi, phi, reduce)
Definition: submesh.F90:1890
type(type_t), public type_float
Definition: types.F90:133
type(type_t), public type_cmplx
Definition: types.F90:134
type(type_t), public type_integer
Definition: types.F90:135
Class defining batches of mesh functions.
Definition: batch.F90:159
The following class implements a lattice iterator. It allows one to loop over all cells that are with...
Describes mesh distribution to nodes.
Definition: mesh.F90:186
This is defined even when running serial.
Definition: mpi.F90:142
A submesh is a type of mesh, used for the projectors in the pseudopotentials It contains points on a ...
Definition: submesh.F90:175
int true(void)
void distance(const int iatom, const int jatom, const double coordinates[], double *rr, double *rr2, double *rr6, double *rr7)
Definition: vdw_ts_low.c:2059