Line data Source code
1 : module dof_mod
2 : use shr_kind_mod, only: r8=>shr_kind_r8, i8=>shr_kind_i8
3 : use dimensions_mod, only: np, npsq, nelem, nelemd
4 : use quadrature_mod, only: quadrature_t
5 : use element_mod, only: element_t,index_t
6 : use spmd_utils, only: mpi_integer
7 : use parallel_mod, only: parallel_t
8 : use edge_mod, only: initedgebuffer,freeedgebuffer, &
9 : longedgevpack, longedgevunpackmin
10 : use edgetype_mod, only: longedgebuffer_t
11 : use bndry_mod, only: bndry_exchange
12 : implicit none
13 : private
14 : ! public data
15 : ! public subroutines
16 : public :: global_dof
17 : public :: UniquePoints
18 : public :: PutUniquePoints
19 : public :: UniqueNcolsP
20 : public :: UniqueCoords
21 : public :: CreateUniqueIndex
22 : public :: SetElemOffset
23 : public :: CreateMetaData
24 :
25 : interface UniquePoints
26 : module procedure UniquePoints2D
27 : module procedure UniquePoints3D
28 : module procedure UniquePoints4D
29 : end interface
30 : interface PutUniquePoints
31 : module procedure PutUniquePoints2D
32 : module procedure PutUniquePoints3D
33 : module procedure PutUniquePoints4D
34 : end interface
35 :
36 :
37 : contains
38 :
39 21600 : subroutine genLocalDof(ig,npts,ldof)
40 :
41 : integer, intent(in) :: ig
42 : integer, intent(in) :: npts
43 : integer, intent(inout) :: ldof(:,:)
44 :
45 : integer :: i,j,npts2
46 :
47 :
48 21600 : npts2=npts*npts
49 108000 : do j=1,npts
50 453600 : do i=1,npts
51 432000 : ldof(i,j) = (ig-1)*npts2 + (j-1)*npts + i
52 : enddo
53 : enddo
54 :
55 21600 : end subroutine genLocalDOF
56 :
57 : ! ===========================================
58 : ! global_dof
59 : !
60 : ! Compute the global degree of freedom for each element...
61 : ! ===========================================
62 :
63 1536 : subroutine global_dof(par,elem)
64 :
65 : type (parallel_t),intent(in) :: par
66 : type (element_t) :: elem(:)
67 :
68 : type (LongEdgeBuffer_t) :: edge
69 :
70 : real(kind=r8) da ! area element
71 :
72 : type (quadrature_t) :: gp
73 :
74 3072 : integer :: ldofP(np,np,nelemd)
75 :
76 : integer ii
77 : integer i,j,ig,ie
78 : integer kptr
79 : integer iptr
80 :
81 : ! ===================
82 : ! begin code
83 : ! ===================
84 1536 : call initEdgeBuffer(edge,1)
85 :
86 : ! =================================================
87 : ! mass matrix on the velocity grid
88 : ! =================================================
89 :
90 :
91 12336 : do ie=1,nelemd
92 10800 : ig = elem(ie)%vertex%number
93 10800 : call genLocalDOF(ig,np,ldofP(:,:,ie))
94 :
95 10800 : kptr=0
96 12336 : call LongEdgeVpack(edge,ldofP(:,:,ie),1,kptr,elem(ie)%desc)
97 : end do
98 :
99 : ! ==============================
100 : ! Insert boundary exchange here
101 : ! ==============================
102 :
103 1536 : call bndry_exchange(par,edge)
104 :
105 12336 : do ie=1,nelemd
106 : ! we should unpack directly into elem(ie)%gdofV, but we dont have
107 : ! a VunpackMIN that takes integer*8. gdofV integer*8 means
108 : ! more than 2G grid points.
109 10800 : kptr=0
110 10800 : call LongEdgeVunpackMIN(edge,ldofP(:,:,ie),1,kptr,elem(ie)%desc)
111 228336 : elem(ie)%gdofP(:,:)=ldofP(:,:,ie)
112 : end do
113 : !$OMP BARRIER
114 1536 : call FreeEdgeBuffer(edge)
115 :
116 1536 : end subroutine global_dof
117 :
118 :
119 0 : subroutine UniquePoints2D(idxUnique,src,dest)
120 : type (index_t) :: idxUnique
121 : real (kind=r8) :: src(:,:)
122 : real (kind=r8) :: dest(:)
123 :
124 : integer :: i,j,ii
125 :
126 :
127 0 : do ii=1,idxUnique%NumUniquePts
128 0 : i=idxUnique%ia(ii)
129 0 : j=idxUnique%ja(ii)
130 0 : dest(ii)=src(i,j)
131 : enddo
132 :
133 0 : end subroutine UniquePoints2D
134 :
135 : ! putUniquePoints first zeros out the destination array, then fills the unique points of the
136 : ! array with values from src. A boundary communication should then be called to fill in the
137 : ! redundent points of the array
138 :
139 0 : subroutine putUniquePoints2D(idxUnique,src,dest)
140 : type (index_t) :: idxUnique
141 : real (kind=r8),intent(in) :: src(:)
142 : real (kind=r8),intent(out) :: dest(:,:)
143 :
144 : integer :: i,j,ii
145 :
146 0 : dest=0.0D0
147 0 : do ii=1,idxUnique%NumUniquePts
148 0 : i=idxUnique%ia(ii)
149 0 : j=idxUnique%ja(ii)
150 0 : dest(i,j)=src(ii)
151 : enddo
152 :
153 0 : end subroutine putUniquePoints2D
154 :
155 0 : subroutine UniqueNcolsP(elem,idxUnique,cid)
156 : use element_mod, only : GetColumnIdP, element_t
157 : type (element_t), intent(in) :: elem
158 : type (index_t), intent(in) :: idxUnique
159 : integer,intent(out) :: cid(:)
160 : integer :: i,j,ii
161 :
162 :
163 0 : do ii=1,idxUnique%NumUniquePts
164 0 : i=idxUnique%ia(ii)
165 0 : j=idxUnique%ja(ii)
166 0 : cid(ii)=GetColumnIdP(elem,i,j)
167 : enddo
168 :
169 0 : end subroutine UniqueNcolsP
170 :
171 :
172 0 : subroutine UniqueCoords(idxUnique,src,lat,lon)
173 :
174 : use coordinate_systems_mod, only : spherical_polar_t
175 : type (index_t), intent(in) :: idxUnique
176 :
177 : type (spherical_polar_t) :: src(:,:)
178 : real (kind=r8), intent(out) :: lat(:)
179 : real (kind=r8), intent(out) :: lon(:)
180 :
181 : integer :: i,j,ii
182 :
183 0 : do ii=1,idxUnique%NumUniquePts
184 0 : i=idxUnique%ia(ii)
185 0 : j=idxUnique%ja(ii)
186 0 : lat(ii)=src(i,j)%lat
187 0 : lon(ii)=src(i,j)%lon
188 : enddo
189 :
190 0 : end subroutine UniqueCoords
191 :
192 0 : subroutine UniquePoints3D(idxUnique,nlyr,src,dest)
193 : type (index_t) :: idxUnique
194 : integer :: nlyr
195 : real (kind=r8) :: src(:,:,:)
196 : real (kind=r8) :: dest(:,:)
197 :
198 : integer :: i,j,k,ii
199 :
200 0 : do ii=1,idxUnique%NumUniquePts
201 0 : i=idxUnique%ia(ii)
202 0 : j=idxUnique%ja(ii)
203 0 : do k=1,nlyr
204 0 : dest(ii,k)=src(i,j,k)
205 : enddo
206 : enddo
207 :
208 0 : end subroutine UniquePoints3D
209 0 : subroutine UniquePoints4D(idxUnique,d3,d4,src,dest)
210 : type (index_t) :: idxUnique
211 : integer :: d3,d4
212 : real (kind=r8) :: src(:,:,:,:)
213 : real (kind=r8) :: dest(:,:,:)
214 :
215 : integer :: i,j,k,n,ii
216 :
217 0 : do n=1,d4
218 0 : do k=1,d3
219 0 : do ii=1,idxUnique%NumUniquePts
220 0 : i=idxUnique%ia(ii)
221 0 : j=idxUnique%ja(ii)
222 0 : dest(ii,k,n)=src(i,j,k,n)
223 : enddo
224 : end do
225 : enddo
226 :
227 0 : end subroutine UniquePoints4D
228 :
229 : ! putUniquePoints first zeros out the destination array, then fills the unique points of the
230 : ! array with values from src. A boundary communication should then be called to fill in the
231 : ! redundent points of the array
232 :
233 0 : subroutine putUniquePoints3D(idxUnique,nlyr,src,dest)
234 : type (index_t) :: idxUnique
235 : integer :: nlyr
236 : real (kind=r8),intent(in) :: src(:,:)
237 : real (kind=r8),intent(out) :: dest(:,:,:)
238 :
239 : integer :: i,j,k,ii
240 :
241 0 : dest=0.0D0
242 0 : do k=1,nlyr
243 0 : do ii=1,idxUnique%NumUniquePts
244 0 : i=idxUnique%ia(ii)
245 0 : j=idxUnique%ja(ii)
246 0 : dest(i,j,k)=src(ii,k)
247 : enddo
248 : enddo
249 :
250 0 : end subroutine putUniquePoints3D
251 :
252 0 : subroutine putUniquePoints4D(idxUnique,d3,d4,src,dest)
253 : type (index_t) :: idxUnique
254 : integer :: d3,d4
255 : real (kind=r8),intent(in) :: src(:,:,:)
256 : real (kind=r8),intent(out) :: dest(:,:,:,:)
257 :
258 : integer :: i,j,k,n,ii
259 :
260 0 : dest=0.0D0
261 0 : do n=1,d4
262 0 : do k=1,d3
263 0 : do ii=1,idxunique%NumUniquePts
264 0 : i=idxUnique%ia(ii)
265 0 : j=idxUnique%ja(ii)
266 0 : dest(i,j,k,n)=src(ii,k,n)
267 : enddo
268 : enddo
269 : end do
270 0 : end subroutine putUniquePoints4D
271 :
272 1536 : subroutine SetElemOffset(par,elem,GlobalUniqueColsP)
273 : use spmd_utils, only : mpi_sum
274 :
275 : type (parallel_t) :: par
276 : type (element_t) :: elem(:)
277 : integer, intent(out) :: GlobalUniqueColsP
278 :
279 1536 : integer, allocatable :: numElemP(:),numElem2P(:)
280 : integer, allocatable :: numElemV(:),numElem2V(:)
281 1536 : integer, allocatable :: gOffset(:)
282 :
283 : integer :: ie, ig, nprocs, ierr
284 : logical, parameter :: Debug = .FALSE.
285 :
286 1536 : nprocs = par%nprocs
287 4608 : allocate(numElemP(nelem))
288 3072 : allocate(numElem2P(nelem))
289 3072 : allocate(gOffset(nelem))
290 24887808 : numElemP=0;numElem2P=0;gOffset=0
291 :
292 12336 : do ie = 1, nelemd
293 10800 : ig = elem(ie)%GlobalId
294 12336 : numElemP(ig) = elem(ie)%idxP%NumUniquePts
295 : end do
296 1536 : call MPI_Allreduce(numElemP,numElem2P,nelem,MPI_INTEGER,MPI_SUM,par%comm,ierr)
297 :
298 1536 : gOffset(1)=1
299 8294400 : do ig = 2, nelem
300 8294400 : gOffset(ig) = gOffset(ig-1)+numElem2P(ig-1)
301 : end do
302 12336 : do ie = 1, nelemd
303 10800 : ig = elem(ie)%GlobalId
304 12336 : elem(ie)%idxP%UniquePtOffset=gOffset(ig)
305 : end do
306 1536 : GlobalUniqueColsP = gOffset(nelem)+numElem2P(nelem)-1
307 :
308 1536 : deallocate(numElemP)
309 1536 : deallocate(numElem2P)
310 1536 : deallocate(gOffset)
311 1536 : end subroutine SetElemOffset
312 :
313 10800 : subroutine CreateUniqueIndex(ig,gdof,idx)
314 :
315 : integer :: ig
316 : type (index_t) :: idx
317 : integer(i8) :: gdof(:,:)
318 :
319 10800 : integer, allocatable :: ldof(:,:)
320 : integer :: i,j,ii,npts
321 :
322 :
323 10800 : npts = size(gdof,dim=1)
324 43200 : allocate(ldof(npts,npts))
325 : ! ====================
326 : ! Form the local DOF
327 : ! ====================
328 10800 : call genLocalDOF(ig,npts,ldof)
329 :
330 10800 : ii=1
331 :
332 54000 : do j=1,npts
333 226800 : do i=1,npts
334 : ! ==========================
335 : ! check for point ownership
336 : ! ==========================
337 216000 : if(gdof(i,j) .eq. ldof(i,j)) then
338 97204 : idx%ia(ii) = i
339 97204 : idx%ja(ii) = j
340 97204 : ii=ii+1
341 : endif
342 : enddo
343 : enddo
344 :
345 10800 : idx%NumUniquePts=ii-1
346 10800 : deallocate(ldof)
347 :
348 10800 : end subroutine CreateUniqueIndex
349 :
350 :
351 0 : subroutine CreateMetaData(par,elem,subelement_corners, fdofp)
352 : type (parallel_t), intent(in) :: par
353 : type (element_t), target :: elem(:)
354 :
355 : integer, optional, intent(out) :: subelement_corners((np-1)*(np-1)*nelemd,4)
356 : integer, optional :: fdofp(np,np,nelemd)
357 :
358 : type (index_t), pointer :: idx
359 : type (LongEdgeBuffer_t) :: edge
360 : integer :: i, j, ii, ie, base
361 : integer(i8), pointer :: gdof(:,:)
362 0 : integer :: fdofp_local(np,np,nelemd)
363 :
364 0 : call initEdgeBuffer(edge,1)
365 0 : fdofp_local=0
366 :
367 0 : do ie=1,nelemd
368 0 : idx => elem(ie)%idxP
369 0 : do ii=1,idx%NumUniquePts
370 0 : i=idx%ia(ii)
371 0 : j=idx%ja(ii)
372 :
373 0 : fdofp_local(i,j,ie) = -(idx%UniquePtoffset+ii-1)
374 : end do
375 0 : call LongEdgeVpack(edge,fdofp_local(:,:,ie),1,0,elem(ie)%desc)
376 : end do
377 0 : call bndry_exchange(par,edge)
378 0 : do ie=1,nelemd
379 0 : base = (ie-1)*(np-1)*(np-1)
380 0 : call LongEdgeVunpackMIN(edge,fdofp_local(:,:,ie),1,0,elem(ie)%desc)
381 0 : if(present(subelement_corners)) then
382 : ii=0
383 0 : do j=1,np-1
384 0 : do i=1,np-1
385 0 : ii=ii+1
386 0 : subelement_corners(base+ii,1) = -fdofp_local(i,j,ie)
387 0 : subelement_corners(base+ii,2) = -fdofp_local(i,j+1,ie)
388 0 : subelement_corners(base+ii,3) = -fdofp_local(i+1,j+1,ie)
389 0 : subelement_corners(base+ii,4) = -fdofp_local(i+1,j,ie)
390 : end do
391 : end do
392 : end if
393 : end do
394 0 : if(present(fdofp)) then
395 0 : fdofp=-fdofp_local
396 : end if
397 :
398 :
399 :
400 0 : end subroutine CreateMetaData
401 :
402 : end module dof_mod
|
