Line data Source code
1 : module GridGraph_mod
2 : !-------------------------
3 : use shr_kind_mod, only: r8=>shr_kind_r8
4 : !-------------------------------
5 : use dimensions_mod, only: max_neigh_edges
6 : !-------------------------
7 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
8 : !-----
9 : use cam_logfile, only: iulog
10 : !-----
11 : implicit none
12 :
13 :
14 : private
15 :
16 : integer, public, parameter :: num_neighbors=8 ! for north, south, east, west, neast, nwest, seast, swest
17 :
18 :
19 : type, public :: GridVertex_t
20 :
21 : integer, pointer :: nbrs(:) => null() ! The numbers of the neighbor elements
22 : integer, pointer :: nbrs_face(:) => null() ! The cube face number of the neighbor element (nbrs array)
23 : integer, pointer :: nbrs_wgt(:) => null() ! The weights for edges defined by nbrs array
24 : integer, pointer :: nbrs_wgt_ghost(:) => null() ! The weights for edges defined by nbrs array
25 : integer :: nbrs_ptr(num_neighbors + 1) !index into the nbrs array for each neighbor direction
26 :
27 : integer :: face_number ! which face of the cube this vertex is on
28 : integer :: number ! element number
29 : integer :: processor_number ! processor number
30 : integer :: SpaceCurve ! index in Space-Filling curve
31 : end type GridVertex_t
32 :
33 : type, public :: GridEdge_t
34 : integer :: head_face ! needed if head vertex has shape (i.e. square)
35 : integer :: tail_face ! needed if tail vertex has shape (i.e. square)
36 : integer :: head_dir !which of 8 neighbor directions is the head
37 : integer :: tail_dir !which of 8 neighbor directions is the tail
38 : integer :: wgtP, wgtS
39 : type (GridVertex_t),pointer :: head => null() ! edge head vertex
40 : type (GridVertex_t),pointer :: tail => null() ! edge tail vertex
41 : logical :: reverse
42 :
43 : end type GridEdge_t
44 :
45 : ! ==========================================
46 : ! Public Interfaces
47 : ! ==========================================
48 :
49 : public :: set_GridVertex_number
50 : public :: PrintGridVertex
51 :
52 : public :: allocate_gridvertex_nbrs
53 : public :: deallocate_gridvertex_nbrs
54 : public :: initgridedge
55 : public :: gridedge_search
56 : public :: gridedge_type
57 : public :: grid_edge_uses_vertex
58 : public :: PrintGridEdge
59 : public :: CheckGridNeighbors
60 : public :: PrintChecksum
61 :
62 : public :: CreateSubGridGraph
63 :
64 : public :: assignment ( = )
65 :
66 : interface assignment ( = )
67 : module procedure copy_gridedge
68 : module procedure copy_gridvertex
69 : end interface
70 :
71 : contains
72 :
73 : !======================================================================
74 :
75 16621200 : subroutine allocate_gridvertex_nbrs(vertex, dim)
76 :
77 : type (GridVertex_t), intent(inout) :: vertex
78 : integer, optional, intent(in) :: dim
79 : integer :: num
80 :
81 16621200 : if (present(dim)) then
82 0 : num = dim
83 : else
84 16621200 : num = max_neigh_edges
85 : end if
86 :
87 49863600 : allocate(vertex%nbrs(num))
88 33242400 : allocate(vertex%nbrs_face(num))
89 33242400 : allocate(vertex%nbrs_wgt(num))
90 33242400 : allocate(vertex%nbrs_wgt_ghost(num))
91 :
92 :
93 16621200 : end subroutine allocate_gridvertex_nbrs
94 : !======================================================================
95 :
96 16599600 : subroutine deallocate_gridvertex_nbrs(vertex)
97 :
98 : type (GridVertex_t), intent(inout) :: vertex
99 :
100 16599600 : deallocate(vertex%nbrs)
101 16599600 : deallocate(vertex%nbrs_face)
102 16599600 : deallocate(vertex%nbrs_wgt)
103 16599600 : deallocate(vertex%nbrs_wgt_ghost)
104 :
105 16599600 : end subroutine deallocate_gridvertex_nbrs
106 :
107 : !======================================================================
108 :
109 : ! =====================================
110 : ! copy edge:
111 : ! copy device for overloading = sign.
112 : ! =====================================
113 :
114 :
115 134168 : recursive subroutine copy_gridedge(edge2, edge1)
116 :
117 : type (GridEdge_t), intent(out) :: edge2
118 : type (GridEdge_t), intent(in) :: edge1
119 :
120 134168 : edge2%tail_face = edge1%tail_face
121 134168 : edge2%head_face = edge1%head_face
122 134168 : edge2%tail_dir = edge1%tail_dir
123 134168 : edge2%head_dir = edge1%head_dir
124 134168 : edge2%reverse = edge1%reverse
125 134168 : edge2%wgtP = edge1%wgtP
126 134168 : edge2%wgtS = edge1%wgtS
127 :
128 :
129 134168 : if (associated(edge1%tail)) then
130 134168 : edge2%tail=>edge1%tail
131 : end if
132 134168 : if (associated(edge1%head)) then
133 134168 : edge2%head=>edge1%head
134 : end if
135 :
136 134168 : end subroutine copy_gridedge
137 :
138 : !======================================================================
139 :
140 21600 : recursive subroutine copy_gridvertex(vertex2, vertex1)
141 :
142 : implicit none
143 :
144 : type (GridVertex_t), intent(out) :: vertex2
145 : type (GridVertex_t), intent(in) :: vertex1
146 :
147 : integer :: i,j,n
148 :
149 21600 : n = SIZE(vertex1%nbrs)
150 :
151 : if (associated(vertex2%nbrs)) then
152 : nullify(vertex2%nbrs)
153 : end if
154 : if (associated(vertex2%nbrs_face)) then
155 : nullify(vertex2%nbrs_face)
156 : end if
157 : if (associated(vertex2%nbrs_wgt)) then
158 : nullify(vertex2%nbrs_wgt)
159 : end if
160 : if (associated(vertex2%nbrs_wgt_ghost)) then
161 : nullify(vertex2%nbrs_wgt_ghost)
162 : end if
163 :
164 21600 : call allocate_gridvertex_nbrs(vertex2)
165 :
166 194400 : do i=1,n
167 172800 : vertex2%nbrs(i) = vertex1%nbrs(i)
168 172800 : vertex2%nbrs_face(i) = vertex1%nbrs_face(i)
169 172800 : vertex2%nbrs_wgt(i) = vertex1%nbrs_wgt(i)
170 194400 : vertex2%nbrs_wgt_ghost(i) = vertex1%nbrs_wgt_ghost(i)
171 : enddo
172 :
173 216000 : do i=1, num_neighbors+1
174 216000 : vertex2%nbrs_ptr(i) = vertex1%nbrs_ptr(i)
175 : enddo
176 :
177 21600 : vertex2%face_number = vertex1%face_number
178 21600 : vertex2%number = vertex1%number
179 21600 : vertex2%processor_number = vertex1%processor_number
180 21600 : vertex2%SpaceCurve = vertex1%SpaceCurve
181 :
182 21600 : end subroutine copy_gridvertex
183 :
184 : !===========================
185 : ! search edge list for match
186 : !===========================
187 :
188 0 : function gridedge_search(nvert1, nvert2, edge) result(number)
189 :
190 : integer, intent(in) :: nvert1
191 : integer, intent(in) :: nvert2
192 : type(GridEdge_t), intent(in) :: edge(:)
193 : integer :: number
194 :
195 : integer :: tmp
196 : integer :: head
197 : integer :: tail
198 :
199 : integer :: nedge
200 : integer :: i
201 :
202 0 : nedge=SIZE(edge)
203 :
204 0 : tail=nvert1
205 0 : head=nvert2
206 :
207 0 : if (tail > head) then
208 0 : tmp = tail
209 0 : tail = head
210 0 : head = tmp
211 : end if
212 :
213 0 : do i=1,nedge
214 0 : if (edge(i)%tail%number==tail .and. edge(i)%head%number==head)then
215 0 : number=i
216 : end if
217 : end do
218 :
219 0 : end function gridedge_search
220 :
221 : !======================================================================
222 :
223 21320 : function gridedge_type(edge) result(type)
224 :
225 : use params_mod, only : INTERNAL_EDGE, EXTERNAL_EDGE
226 : type (GridEdge_t), intent(in) :: edge
227 : integer :: type
228 :
229 21320 : if (edge%head%processor_number==edge%tail%processor_number) then
230 : type=INTERNAL_EDGE
231 : else
232 19784 : type=EXTERNAL_EDGE
233 : endif
234 :
235 21320 : end function gridedge_type
236 :
237 : !======================================================================
238 :
239 :
240 :
241 0 : function grid_edge_uses_vertex(Vertex,Edge) result(log)
242 :
243 : type(GridVertex_t), intent(in) :: Vertex
244 : type(GridEdge_t), intent(in) :: Edge
245 : logical :: log
246 : integer :: number
247 :
248 0 : number = Vertex%number
249 0 : if(number == Edge%head%number .or. number == Edge%tail%number) then
250 : log = .TRUE.
251 : else
252 0 : log = .FALSE.
253 : endif
254 :
255 0 : end function grid_edge_uses_vertex
256 :
257 : !======================================================================
258 :
259 0 : subroutine PrintChecksum(TestPattern,Checksum)
260 :
261 : use dimensions_mod, only : nlev, nelemd, np
262 :
263 : implicit none
264 :
265 : real(kind=r8), target,intent(in) :: TestPattern(:,:,:,:)
266 : real(kind=r8), target,intent(in) :: Checksum(:,:,:,:)
267 :
268 : integer :: i,k,ix,iy
269 :
270 0 : print *
271 0 : write (iulog,*) 'checksums:'
272 0 : do i=1,nelemd
273 : ! Lets start out only looking at the first element
274 0 : write(iulog,*)
275 0 : do k=1,nlev
276 0 : do iy=1,np
277 0 : do ix=1,np
278 0 : write(iulog,*)INT(TestPattern(ix,iy,k,i))," checksum = ",INT(Checksum(ix,iy,k,i))
279 : enddo
280 : enddo
281 : enddo
282 : enddo
283 :
284 :
285 0 : end subroutine PrintChecksum
286 :
287 : !======================================================================
288 :
289 0 : subroutine CreateSubGridGraph(Vertex, SVertex, local2global)
290 :
291 : implicit none
292 :
293 : type (GridVertex_t),intent(in) :: Vertex(:)
294 : type (GridVertex_t),intent(inout) :: SVertex(:)
295 : integer,intent(in) :: local2global(:)
296 :
297 : integer :: nelem,nelem_s,n,ncount,cnt,pos, orig_start
298 : integer :: inbr,i,ig,j,k, new_pos
299 :
300 0 : integer,allocatable :: global2local(:)
301 :
302 0 : nelem = SIZE(Vertex)
303 0 : nelem_s = SiZE(SVertex)
304 :
305 0 : allocate(global2local(nelem))
306 :
307 0 : global2local(:) = 0
308 0 : do i=1,nelem_s
309 0 : ig = local2global(i)
310 0 : global2local(ig) = i
311 : enddo
312 :
313 0 : do i=1,nelem_s
314 0 : ig = local2global(i)
315 :
316 0 : call copy_gridvertex(SVertex(i),Vertex(ig)) !svertex(i) = vertex(ig)
317 :
318 0 : n = SIZE(SVertex(i)%nbrs(:))
319 : ! ==============================================
320 : ! Apply the correction to the neighbors list to
321 : ! reflect new subgraph numbers
322 : ! ==============================================
323 :
324 0 : orig_start = 1
325 :
326 0 : do j=1,num_neighbors
327 :
328 0 : cnt = Svertex(i)%nbrs_ptr(j+1) - orig_start !number of neighbors for this direction
329 0 : ncount = 0
330 0 : do k = 1, cnt
331 0 : pos = orig_start + k-1
332 0 : inbr = global2local(Svertex(i)%nbrs(pos))
333 :
334 0 : if(inbr .gt. 0) then
335 0 : new_pos = Svertex(i)%nbrs_ptr(j) + ncount
336 :
337 0 : Svertex(i)%nbrs(new_pos) = inbr
338 0 : Svertex(i)%nbrs_face(new_pos) = Svertex(i)%nbrs_face(pos)
339 0 : Svertex(i)%nbrs_wgt(new_pos) = Svertex(i)%nbrs_wgt(pos)
340 0 : Svertex(i)%nbrs_wgt_ghost(new_pos) = Svertex(i)%nbrs_wgt_ghost(pos)
341 0 : ncount = ncount+1
342 : endif
343 : enddo
344 : !set neighbors ptr
345 0 : orig_start = Svertex(i)%nbrs_ptr(j+1);
346 0 : Svertex(i)%nbrs_ptr(j+1) = Svertex(i)%nbrs_ptr(j) + ncount
347 :
348 :
349 : enddo !num_neighbors loop
350 :
351 :
352 0 : Svertex(i)%number = i
353 : enddo !nelem_s loop
354 0 : deallocate(global2local)
355 :
356 0 : end subroutine CreateSubGridGraph
357 :
358 : !======================================================================
359 :
360 0 : subroutine PrintGridEdge(Edge)
361 :
362 : implicit none
363 : type (GridEdge_t), intent(in) :: Edge(:)
364 :
365 : integer :: i,nedge,ii,wgtP
366 :
367 0 : nedge = SIZE(Edge)
368 :
369 0 : write(iulog,95)
370 0 : do i=1,nedge
371 0 : ii=Edge(i)%tail_face
372 :
373 : !map to correct location - for now all on same nbr side have same wgt, so take the first one
374 0 : ii = Edge(i)%tail%nbrs_ptr(ii)
375 :
376 0 : wgtP=Edge(i)%tail%nbrs_wgt(ii)
377 0 : write(iulog,100) i, &
378 0 : Edge(i)%tail%number,Edge(i)%tail_face, wgtP, &
379 0 : Edge(i)%head%number,Edge(i)%head_face, gridedge_type(Edge(i))
380 : enddo
381 : 95 format(5x,'GRIDEDGE #',3x,'Tail (face)',5x,'Head (face)',3x,'Type')
382 : 100 format(10x,I6,8x,I4,1x,'(',I1,') --',I2,'--> ',I6,1x,'(',I1,')',5x,'[',I1,']')
383 :
384 0 : end subroutine PrintGridEdge
385 :
386 : !======================================================================
387 : ! ==========================================
388 : ! set_GridVertex_neighbors:
389 : !
390 : ! Set global element number for element elem
391 : ! ==========================================
392 :
393 0 : subroutine set_GridVertex_number(elem,number)
394 :
395 : type(GridVertex_t) :: elem
396 : integer :: number
397 :
398 0 : elem%number=number
399 :
400 0 : end subroutine set_GridVertex_number
401 :
402 : !======================================================================
403 0 : subroutine PrintGridVertex(Vertex)
404 :
405 : implicit none
406 : type (GridVertex_t), intent(in),target :: Vertex(:)
407 :
408 : integer :: i,nvert
409 : integer ::n_west, n_east, n_south, n_north, n_swest, n_seast, n_nwest, n_neast
410 : integer ::w_west, w_east, w_south, w_north, w_swest, w_seast, w_nwest, w_neast
411 : integer ::n, print_buf(90), nbr(8), j, k, start, cnt, nbrs_cnt(8)
412 :
413 0 : nbr = (/ west, east, south, north, swest, seast, nwest, neast/)
414 :
415 0 : nvert = SIZE(Vertex)
416 :
417 0 : write(iulog,98)
418 0 : do i=1,nvert
419 :
420 0 : print_buf(:) = 0
421 0 : nbrs_cnt(:) = 0
422 0 : cnt = 1
423 0 : do j = 1,num_neighbors
424 0 : n = Vertex(i)%nbrs_ptr(nbr(j)+1) - Vertex(i)%nbrs_ptr(nbr(j)) !num neigbors in that directions
425 0 : start = Vertex(i)%nbrs_ptr(nbr(j)) !start in array
426 0 : nbrs_cnt(j) = n
427 0 : do k = 1, n
428 0 : print_buf(cnt) = Vertex(i)%nbrs(start+k-1)
429 0 : print_buf(cnt+1) = Vertex(i)%nbrs_wgt(start+k-1)
430 0 : print_buf(cnt+2) = Vertex(i)%nbrs_face(start+k-1)
431 0 : cnt = cnt + 3
432 : end do
433 : enddo
434 :
435 0 : write(iulog,991) Vertex(i)%number, Vertex(i)%processor_number, &
436 0 : Vertex(i)%face_number, &
437 0 : print_buf(1:cnt-1)
438 :
439 0 : write(iulog,992) nbrs_cnt(1:8)
440 :
441 :
442 : enddo
443 : 98 format(5x,'GRIDVERTEX #',2x,'PART',2x,'DEG',4x,'W',9x,'E',9x, &
444 : 'S',9x,'N',9x,'SW',9x,'SE',9x,'NW',9x,'NE')
445 :
446 : 991 format(10x,I3,8x,I4,8x,I4,2x,30(1x,I4,1x,'(',I2,I2,')'))
447 : 992 format(30x,'nbrs_cnt:', 2x,8(1x,I4))
448 :
449 0 : end subroutine PrintGridVertex
450 :
451 :
452 : !======================================================================
453 :
454 0 : subroutine CheckGridNeighbors(Vertex)
455 :
456 : implicit none
457 : type (GridVertex_t), intent(in) :: Vertex(:)
458 :
459 : integer :: i,j,k,l,m,nnbrs,inbrs,nvert
460 0 : nvert = SIZE(Vertex)
461 :
462 0 : do i=1,nvert
463 0 : nnbrs = SIZE(Vertex(i)%nbrs)
464 0 : do j=1,nnbrs
465 0 : inbrs = Vertex(i)%nbrs(j)
466 0 : if(inbrs > 0) then
467 0 : do k=1,nnbrs
468 0 : if( inbrs .eq. Vertex(i)%nbrs(k) .and. (j/=k) ) &
469 0 : write(iulog,*)'CheckGridNeighbors: ERROR identical neighbors detected for Vertex ',i
470 :
471 : enddo
472 : endif
473 : enddo
474 : enddo
475 :
476 0 : end subroutine CheckGridNeighbors
477 :
478 : !======================================================================
479 1536 : subroutine initgridedge(GridEdge,GridVertex)
480 : use cam_abortutils, only : endrun
481 : use dimensions_mod, only : max_corner_elem
482 :
483 : type (GridEdge_t), intent(inout) :: GridEdge(:)
484 : type (GridVertex_t), intent(in),target :: GridVertex(:)
485 :
486 : integer :: i,j,k,iptr,m,n,wgtV,wgtP
487 : integer :: nelem,nelem_edge,inbr
488 : logical :: Verbose=.FALSE.
489 : integer :: mynbr_cnt, cnt, mystart, start
490 :
491 1536 : nelem = SIZE(GridVertex)
492 1536 : nelem_edge = SIZE(GridEdge)
493 :
494 66319872 : GridEdge(:)%reverse=.FALSE.
495 66319872 : GridEdge(:)%wgtP=-1
496 66319872 : GridEdge(:)%wgtS=-1
497 :
498 : iptr=1
499 8295936 : do j=1,nelem
500 74651136 : do i=1,num_neighbors
501 66355200 : mynbr_cnt = GridVertex(j)%nbrs_ptr(i+1) - GridVertex(j)%nbrs_ptr(i) !length of neighbor location
502 : mystart = GridVertex(j)%nbrs_ptr(i)
503 140967936 : do m=0,mynbr_cnt-1
504 132673536 : if((GridVertex(j)%nbrs_wgt(mystart + m) .gt. 0)) then ! Do this only if has a non-zero weight
505 66318336 : if (nelem_edge<iptr) call endrun('Error in initgridedge: Number of edges greater than expected.')
506 66318336 : GridEdge(iptr)%tail => GridVertex(j)
507 66318336 : GridEdge(iptr)%tail_face = mystart + m ! needs to be mystart + m (location in array)
508 66318336 : GridEdge(iptr)%tail_dir = i*max_corner_elem + m !conversion needed for setcycle
509 66318336 : inbr = GridVertex(j)%nbrs(mystart+m)
510 66318336 : GridEdge(iptr)%head => GridVertex(inbr)
511 :
512 : ! ===========================================
513 : ! Need this awful piece of code to determine
514 : ! which "face" of the neighbor element the
515 : ! edge links (i.e. the "head_face")
516 : ! ===========================================
517 596865024 : do k=1,num_neighbors
518 530546688 : cnt = GridVertex(inbr)%nbrs_ptr(k+1) -GridVertex(inbr)%nbrs_ptr(k)
519 : start = GridVertex(inbr)%nbrs_ptr(k)
520 1127153664 : do n = 0, cnt-1
521 1060835328 : if(GridVertex(inbr)%nbrs(start+n) == GridVertex(j)%number) then
522 66318336 : GridEdge(iptr)%head_face=start+n !needs to be start + n (location in array)
523 66318336 : GridEdge(iptr)%head_dir=k*max_corner_elem+n !conversion (un-done in setcycle)
524 : endif
525 : enddo
526 : enddo
527 66318336 : GridEdge(iptr)%wgtP = GridVertex(j)%nbrs_wgt(mystart+m)
528 66318336 : GridEdge(iptr)%wgtS = 1
529 66318336 : iptr=iptr+1
530 : end if
531 : end do ! m loop
532 : end do !end i loop
533 : end do !end j loop
534 1536 : if (nelem_edge+1 /= iptr) then
535 0 : call endrun('Error in initgridedge: Number of edges less than expected.')
536 : end if
537 1536 : if (Verbose) then
538 :
539 0 : print *
540 0 : write(iulog,*)"element edge tail,head list: (TEST)"
541 0 : do i=1,nelem_edge
542 0 : write(iulog,*)GridEdge(i)%tail%number,GridEdge(i)%head%number
543 : end do
544 :
545 0 : print *
546 0 : write(iulog,*)"element edge tail_face, head_face list: (TEST)"
547 0 : do i=1,nelem_edge
548 0 : write(iulog,*)GridEdge(i)%tail_face,GridEdge(i)%head_face
549 : end do
550 : end if
551 :
552 1536 : end subroutine initgridedge
553 : !======================================================================
554 :
555 0 : end module GridGraph_mod
|
