| 1 | c===============================================================================
|
|---|
| 2 | c pff@cfm.brown.edu 3/19/96
|
|---|
| 3 | c
|
|---|
| 4 | c
|
|---|
| 5 | c This is a suite of routines for solving overlapping subdomain
|
|---|
| 6 | c problems with finite elements on distributed memory machines.
|
|---|
| 7 | c
|
|---|
| 8 | c The overall hierarchy of data structures is as follows:
|
|---|
| 9 | c
|
|---|
| 10 | c System - index set denoted by _glob
|
|---|
| 11 | c
|
|---|
| 12 | c Processor - index set denoted by _pglob
|
|---|
| 13 | c
|
|---|
| 14 | c .Domain - index set denoted by _dloc (1,2,3,...,n_dloc)
|
|---|
| 15 | c
|
|---|
| 16 | c .Sp.Elem - index set denoted by _sloc (1,2,3,...,n_sloc)
|
|---|
| 17 | c
|
|---|
| 18 | c
|
|---|
| 19 | c A critical component of the parallel DD solver is the gather-scatter
|
|---|
| 20 | c operation. As there may be more than one domain on a processor, and
|
|---|
| 21 | c communication must be minimized, it is critical that communication be
|
|---|
| 22 | c processor oriented, not domain oriented. Hence domains will access data
|
|---|
| 23 | c via the dloc_to_pglob interface, and the pglob indexed variables will
|
|---|
| 24 | c be accessed via a gather-scatter which interacts via the pglob_glob
|
|---|
| 25 | c interface. Noticed that, in a uni-processor application, the pglob_glob
|
|---|
| 26 | c map will be simply the identity.
|
|---|
| 27 | c
|
|---|
| 28 | c===============================================================================
|
|---|
| 29 | c
|
|---|
| 30 | subroutine set_overlap
|
|---|
| 31 | c
|
|---|
| 32 | c Set up arrays for overlapping Schwartz algorithm *for pressure solver*
|
|---|
| 33 | c
|
|---|
| 34 | include 'SIZE'
|
|---|
| 35 | include 'DOMAIN'
|
|---|
| 36 | include 'ESOLV'
|
|---|
| 37 | include 'INPUT'
|
|---|
| 38 | include 'TSTEP'
|
|---|
| 39 | c
|
|---|
| 40 | REAL*8 dnekclock,t0
|
|---|
| 41 | c
|
|---|
| 42 | parameter ( n_tri = 7*ltotd )
|
|---|
| 43 | common /scrns/ tri (n_tri)
|
|---|
| 44 | integer tri,elem
|
|---|
| 45 | c
|
|---|
| 46 | common /screv/ x(2*ltotd)
|
|---|
| 47 | common /scrvh/ y(2*ltotd)
|
|---|
| 48 | common /scrch/ z(2*ltotd)
|
|---|
| 49 | c
|
|---|
| 50 | common /ctmp0/ nv_to_t(2*ltotd)
|
|---|
| 51 | c
|
|---|
| 52 | parameter (lia = ltotd - 2 - 2*lelt)
|
|---|
| 53 | common /scrcg/ ntri(lelt+1),nmask(lelt+1)
|
|---|
| 54 | $ , ia(lia)
|
|---|
| 55 | c
|
|---|
| 56 | common /scruz/ color (4*ltotd)
|
|---|
| 57 | common /scrmg/ ddmask (4*ltotd)
|
|---|
| 58 | common /ctmp1/ mask (4*ltotd)
|
|---|
| 59 |
|
|---|
| 60 | parameter(lxx=lx1*lx1, levb=lelv+lbelv)
|
|---|
| 61 | common /fastd/ df(lx1*ly1*lz1,levb)
|
|---|
| 62 | $ , sr(lxx*2,levb),ss(lxx*2,levb),st(lxx*2,levb)
|
|---|
| 63 |
|
|---|
| 64 |
|
|---|
| 65 | integer e
|
|---|
| 66 |
|
|---|
| 67 | if (lx1.eq.2) param(43)=1.
|
|---|
| 68 | if (lx1.eq.2.and.nid.eq.0) write(6,*) 'No mgrid for lx1=2!'
|
|---|
| 69 |
|
|---|
| 70 | if (ifaxis) ifmgrid = .false.
|
|---|
| 71 | if (param(43).ne.0) ifmgrid = .false.
|
|---|
| 72 |
|
|---|
| 73 | npass = 1
|
|---|
| 74 | if (ifmhd) npass = 2
|
|---|
| 75 | do ipass=1,npass
|
|---|
| 76 | ifield = 1
|
|---|
| 77 |
|
|---|
| 78 | if (ifsplit.and.ifmgrid) then
|
|---|
| 79 |
|
|---|
| 80 | if (ipass.gt.1) ifield = ifldmhd
|
|---|
| 81 |
|
|---|
| 82 | call swap_lengths
|
|---|
| 83 | call gen_fast_spacing(x,y,z)
|
|---|
| 84 |
|
|---|
| 85 | call hsmg_setup
|
|---|
| 86 | call h1mg_setup
|
|---|
| 87 |
|
|---|
| 88 | elseif (.not.ifsplit) then ! Pn-Pn-2
|
|---|
| 89 |
|
|---|
| 90 | if (ipass.gt.1) ifield = ifldmhd
|
|---|
| 91 |
|
|---|
| 92 | if (param(44).eq.1) then ! Set up local overlapping solves
|
|---|
| 93 | call set_fem_data_l2(nel_proc,ndom,n_o,x,y,z,tri)
|
|---|
| 94 | else
|
|---|
| 95 | call swap_lengths
|
|---|
| 96 | endif
|
|---|
| 97 |
|
|---|
| 98 | e = 1
|
|---|
| 99 | if (ifield.gt.1) e = nelv+1
|
|---|
| 100 |
|
|---|
| 101 | call gen_fast_spacing(x,y,z)
|
|---|
| 102 | call gen_fast(df(1,e),sr(1,e),ss(1,e),st(1,e),x,y,z)
|
|---|
| 103 |
|
|---|
| 104 | call init_weight_op
|
|---|
| 105 | if (param(43).eq.0) call hsmg_setup
|
|---|
| 106 | endif
|
|---|
| 107 |
|
|---|
| 108 | call set_up_h1_crs
|
|---|
| 109 |
|
|---|
| 110 | enddo
|
|---|
| 111 |
|
|---|
| 112 | return
|
|---|
| 113 | end
|
|---|
| 114 | c-----------------------------------------------------------------------
|
|---|
| 115 | subroutine overflow_ck(n_req,n_avail,signal)
|
|---|
| 116 | c
|
|---|
| 117 | c Check for buffer overflow
|
|---|
| 118 | c
|
|---|
| 119 | character*11 signal
|
|---|
| 120 | c
|
|---|
| 121 | nid = mynode()
|
|---|
| 122 | c
|
|---|
| 123 | if (n_req.gt.n_avail) then
|
|---|
| 124 | write(6,9) nid,n_req,n_avail,nid,signal
|
|---|
| 125 | 9 format(i7,' ERROR: requested array space (',i12
|
|---|
| 126 | $ ,') exceeds allocated amount (',i12,').'
|
|---|
| 127 | $ ,/,i12,' ABORTING.',3x,a11
|
|---|
| 128 | $ ,/,i12,' ABORTING.',3x,'from overflow_ck call.')
|
|---|
| 129 | call exitt
|
|---|
| 130 | endif
|
|---|
| 131 | return
|
|---|
| 132 | end
|
|---|
| 133 | c-----------------------------------------------------------------------
|
|---|
| 134 | subroutine iunswap(b,ind,n,temp)
|
|---|
| 135 | integer b(1),ind(1),temp(1)
|
|---|
| 136 | c
|
|---|
| 137 | c sort associated elements by putting item(jj)
|
|---|
| 138 | c into item(i), where jj=ind(i).
|
|---|
| 139 | c
|
|---|
| 140 | do 20 i=1,n
|
|---|
| 141 | jj=ind(i)
|
|---|
| 142 | temp(jj)=b(i)
|
|---|
| 143 | 20 continue
|
|---|
| 144 | do 30 i=1,n
|
|---|
| 145 | 30 b(i)=temp(i)
|
|---|
| 146 | return
|
|---|
| 147 | end
|
|---|
| 148 | c-----------------------------------------------------------------------
|
|---|
| 149 | subroutine set_fem_data_l2(nep,nd,no,x,y,z,p)
|
|---|
| 150 | include 'SIZE'
|
|---|
| 151 | include 'DOMAIN'
|
|---|
| 152 | include 'NONCON'
|
|---|
| 153 | include 'TOTAL'
|
|---|
| 154 | c
|
|---|
| 155 | real x(lx1,ly1,lz1,1),y(lx1,ly1,lz1,1),z(lx1,ly1,lz1,1)
|
|---|
| 156 | real p(lx1,ly1,lz1,1)
|
|---|
| 157 | integer ce,cf
|
|---|
| 158 | c
|
|---|
| 159 | common /ctmp0/ w1(lx1,ly1),w2(lx1,ly1)
|
|---|
| 160 | c
|
|---|
| 161 | c First, copy local geometry to temporary, expanded, arrays
|
|---|
| 162 | c
|
|---|
| 163 | nxy1 = lx1*ly1
|
|---|
| 164 | nxy2 = lx2*ly2
|
|---|
| 165 | nxyz1 = lx1*ly1*lz1
|
|---|
| 166 | nxyz2 = lx2*ly2*lz2
|
|---|
| 167 | ntot1 = nelv*nxyz1
|
|---|
| 168 | ntot2 = nelv*nxyz2
|
|---|
| 169 | c
|
|---|
| 170 | call rzero(x,ntot1)
|
|---|
| 171 | call rzero(y,ntot1)
|
|---|
| 172 | call rzero(z,ntot1)
|
|---|
| 173 | c
|
|---|
| 174 | c Take care of surface geometry
|
|---|
| 175 | c
|
|---|
| 176 | call map_face12(x,xm1,w1,w2)
|
|---|
| 177 | call map_face12(y,ym1,w1,w2)
|
|---|
| 178 | call map_face12(z,zm1,w1,w2)
|
|---|
| 179 | c
|
|---|
| 180 | c Take care of interior geometry
|
|---|
| 181 | c
|
|---|
| 182 | iz1 = 0
|
|---|
| 183 | if (if3d) iz1=1
|
|---|
| 184 | do ie=1,nelv
|
|---|
| 185 | do iz=1,lz2
|
|---|
| 186 | do iy=1,ly2
|
|---|
| 187 | do ix=1,lx2
|
|---|
| 188 | x(ix+1,iy+1,iz+iz1,ie) = xm2(ix,iy,iz,ie)
|
|---|
| 189 | y(ix+1,iy+1,iz+iz1,ie) = ym2(ix,iy,iz,ie)
|
|---|
| 190 | z(ix+1,iy+1,iz+iz1,ie) = zm2(ix,iy,iz,ie)
|
|---|
| 191 | enddo
|
|---|
| 192 | enddo
|
|---|
| 193 | enddo
|
|---|
| 194 | enddo
|
|---|
| 195 | c
|
|---|
| 196 | c
|
|---|
| 197 | c Compute absolute value of distance between pressure and vel. planes
|
|---|
| 198 | c
|
|---|
| 199 | call dface_add1sa(x)
|
|---|
| 200 | call dface_add1sa(y)
|
|---|
| 201 | call dface_add1sa(z)
|
|---|
| 202 | c
|
|---|
| 203 | c Sum this distance
|
|---|
| 204 | c
|
|---|
| 205 | ifielt = ifield
|
|---|
| 206 | ifield = 1
|
|---|
| 207 | c call dssum(x,lx1,ly1,lz1)
|
|---|
| 208 | c call dssum(y,lx1,ly1,lz1)
|
|---|
| 209 | c call dssum(z,lx1,ly1,lz1)
|
|---|
| 210 | c
|
|---|
| 211 | c Scale children values by 0.5. This is a hack, based on assumption
|
|---|
| 212 | c that number of children sharing a parent edge is 2
|
|---|
| 213 | c
|
|---|
| 214 | scale = 0.5
|
|---|
| 215 | if (if3d) scale = 0.25
|
|---|
| 216 | do im = 1 , mort_m
|
|---|
| 217 | cf = noncon_f(im)
|
|---|
| 218 | ce = noncon_e(im)
|
|---|
| 219 | call faces(x,scale,ce,cf,lx1,ly1,lz1)
|
|---|
| 220 | call faces(y,scale,ce,cf,lx1,ly1,lz1)
|
|---|
| 221 | if (if3d) call faces(z,scale,ce,cf,lx1,ly1,lz1)
|
|---|
| 222 | enddo
|
|---|
| 223 | call fgslib_gs_op_many(gsh_fld(ifield), x,y,z,x,x,x,ldim, 1,1,0)
|
|---|
| 224 | c
|
|---|
| 225 | ifield = ifielt
|
|---|
| 226 | c
|
|---|
| 227 | c Dirichlet (pmask=0) faces all set...
|
|---|
| 228 | c
|
|---|
| 229 | nel_proc = nelv
|
|---|
| 230 | ndom = nelv
|
|---|
| 231 | c ndom = 1
|
|---|
| 232 | c
|
|---|
| 233 | n_o = 1
|
|---|
| 234 | c
|
|---|
| 235 | c These are the return values, since the calling routine doesn't
|
|---|
| 236 | c know DOMAIN
|
|---|
| 237 | c
|
|---|
| 238 | nep = nel_proc
|
|---|
| 239 | nd = ndom
|
|---|
| 240 | no = n_o
|
|---|
| 241 | ifield = ifielt
|
|---|
| 242 | return
|
|---|
| 243 | end
|
|---|
| 244 | c-----------------------------------------------------------------------
|
|---|
| 245 | subroutine map_face12(x2,x1,w1,w2)
|
|---|
| 246 | include 'SIZE'
|
|---|
| 247 | include 'INPUT'
|
|---|
| 248 | include 'IXYZ'
|
|---|
| 249 | c
|
|---|
| 250 | c Interpolate iface of x1 (GLL pts) onto x2 (GL pts).
|
|---|
| 251 | c Work arrays should be of size lx1*lx1 each.
|
|---|
| 252 | c
|
|---|
| 253 | real x2(lx1*ly1*lz1,1)
|
|---|
| 254 | real x1(lx1*ly1*lz1,1)
|
|---|
| 255 | real w1(1),w2(1)
|
|---|
| 256 | c
|
|---|
| 257 | c
|
|---|
| 258 | do ie=1,nelv
|
|---|
| 259 | if (if3d) then
|
|---|
| 260 | call map_one_face12(x2(1,ie),x1(1,ie),1,ixm12,ixtm12,w1,w2)
|
|---|
| 261 | call map_one_face12(x2(1,ie),x1(1,ie),2,ixm12,ixtm12,w1,w2)
|
|---|
| 262 | call map_one_face12(x2(1,ie),x1(1,ie),3,ixm12,ixtm12,w1,w2)
|
|---|
| 263 | call map_one_face12(x2(1,ie),x1(1,ie),4,ixm12,ixtm12,w1,w2)
|
|---|
| 264 | call map_one_face12(x2(1,ie),x1(1,ie),5,ixm12,ixtm12,w1,w2)
|
|---|
| 265 | call map_one_face12(x2(1,ie),x1(1,ie),6,ixm12,ixtm12,w1,w2)
|
|---|
| 266 | else
|
|---|
| 267 | call map_one_face12(x2(1,ie),x1(1,ie),1,ixm12,ixtm12,w1,w2)
|
|---|
| 268 | call map_one_face12(x2(1,ie),x1(1,ie),2,ixm12,ixtm12,w1,w2)
|
|---|
| 269 | call map_one_face12(x2(1,ie),x1(1,ie),3,ixm12,ixtm12,w1,w2)
|
|---|
| 270 | call map_one_face12(x2(1,ie),x1(1,ie),4,ixm12,ixtm12,w1,w2)
|
|---|
| 271 | endif
|
|---|
| 272 | enddo
|
|---|
| 273 | c
|
|---|
| 274 | return
|
|---|
| 275 | end
|
|---|
| 276 | c-----------------------------------------------------------------------
|
|---|
| 277 | subroutine map_one_face12(x2,x1,iface,i12,i12t,w1,w2)
|
|---|
| 278 | include 'SIZE'
|
|---|
| 279 | include 'INPUT'
|
|---|
| 280 | c
|
|---|
| 281 | c Interpolate iface of x1 (GLL pts) onto interior of face of x2 (GL pts).
|
|---|
| 282 | c Work arrays should be of size lx1*lx1 each.
|
|---|
| 283 | c
|
|---|
| 284 | real x2(lx1,ly1,lz1)
|
|---|
| 285 | real x1(lx1,ly1,lz1)
|
|---|
| 286 | real w1(1),w2(1)
|
|---|
| 287 | real i12(1),i12t(1)
|
|---|
| 288 | c
|
|---|
| 289 | c
|
|---|
| 290 | c Extract surface data from x1
|
|---|
| 291 | c
|
|---|
| 292 | CALL FACIND (KX1,KX2,KY1,KY2,KZ1,KZ2,lx1,ly1,lz1,IFACE)
|
|---|
| 293 | i=0
|
|---|
| 294 | do iz=kz1,kz2
|
|---|
| 295 | do iy=ky1,ky2
|
|---|
| 296 | do ix=kx1,kx2
|
|---|
| 297 | i = i+1
|
|---|
| 298 | w1(i) = x1(ix,iy,iz)
|
|---|
| 299 | enddo
|
|---|
| 300 | enddo
|
|---|
| 301 | enddo
|
|---|
| 302 | c
|
|---|
| 303 | c Interpolate from mesh 1 to 2
|
|---|
| 304 | c
|
|---|
| 305 | if (if3d) then
|
|---|
| 306 | call mxm(i12 ,lx2,w1,lx1,w2,lx1)
|
|---|
| 307 | call mxm(w2,lx2,i12t,lx1,w1,lx2)
|
|---|
| 308 | else
|
|---|
| 309 | call mxm(i12 ,lx2,w1,lx1,w2, 1)
|
|---|
| 310 | call copy(w1,w2,lx2)
|
|---|
| 311 | endif
|
|---|
| 312 | c
|
|---|
| 313 | c Write to interior points on face of x2
|
|---|
| 314 | c
|
|---|
| 315 | kx1=min(kx1+1,lx1,kx2)
|
|---|
| 316 | kx2=max(kx2-1, 1,kx1)
|
|---|
| 317 | ky1=min(ky1+1,ly1,ky2)
|
|---|
| 318 | ky2=max(ky2-1, 1,ky1)
|
|---|
| 319 | kz1=min(kz1+1,lz1,kz2)
|
|---|
| 320 | kz2=max(kz2-1, 1,kz1)
|
|---|
| 321 | c
|
|---|
| 322 | i=0
|
|---|
| 323 | do iz=kz1,kz2
|
|---|
| 324 | do iy=ky1,ky2
|
|---|
| 325 | do ix=kx1,kx2
|
|---|
| 326 | i = i+1
|
|---|
| 327 | x2(ix,iy,iz) = w1(i)
|
|---|
| 328 | enddo
|
|---|
| 329 | enddo
|
|---|
| 330 | enddo
|
|---|
| 331 | c
|
|---|
| 332 | return
|
|---|
| 333 | end
|
|---|
| 334 | c-----------------------------------------------------------------------
|
|---|
| 335 | subroutine dface_add1sa(x)
|
|---|
| 336 | c Compute |face-interior|
|
|---|
| 337 | c
|
|---|
| 338 | include 'SIZE'
|
|---|
| 339 | include 'INPUT'
|
|---|
| 340 | real x(lx1,ly1,lz1,1)
|
|---|
| 341 | c
|
|---|
| 342 | do ie=1,nelv
|
|---|
| 343 | c
|
|---|
| 344 | if (if3d) then
|
|---|
| 345 | c
|
|---|
| 346 | do iz=2,lz1-1
|
|---|
| 347 | do ix=2,lx1-1
|
|---|
| 348 | x(ix,1 ,iz,ie)=abs(x(ix,1 ,iz,ie) - x(ix,2 ,iz,ie))
|
|---|
| 349 | x(ix,ly1,iz,ie)=abs(x(ix,ly1,iz,ie) - x(ix,ly1-1,iz,ie))
|
|---|
| 350 | enddo
|
|---|
| 351 | enddo
|
|---|
| 352 | c
|
|---|
| 353 | do iz=2,lz1-1
|
|---|
| 354 | do iy=2,ly1-1
|
|---|
| 355 | x(1 ,iy,iz,ie)=abs(x(1 ,iy,iz,ie) - x(2 ,iy,iz,ie))
|
|---|
| 356 | x(lx1,iy,iz,ie)=abs(x(lx1,iy,iz,ie) - x(lx1-1,iy,iz,ie))
|
|---|
| 357 | enddo
|
|---|
| 358 | enddo
|
|---|
| 359 | c
|
|---|
| 360 | do iy=2,ly1-1
|
|---|
| 361 | do ix=2,lx1-1
|
|---|
| 362 | x(ix,iy,1 ,ie)=abs(x(ix,iy,1 ,ie) - x(ix,iy,2 ,ie))
|
|---|
| 363 | x(ix,iy,lz1,ie)=abs(x(ix,iy,lz1,ie) - x(ix,iy,lz1-1,ie))
|
|---|
| 364 | enddo
|
|---|
| 365 | enddo
|
|---|
| 366 | c
|
|---|
| 367 | else
|
|---|
| 368 | c
|
|---|
| 369 | c 2D
|
|---|
| 370 | c
|
|---|
| 371 | do ix=2,lx1-1
|
|---|
| 372 | x(ix,1 ,1,ie)=abs(x(ix,1 ,1,ie) - x(ix,2 ,1,ie))
|
|---|
| 373 | x(ix,ly1,1,ie)=abs(x(ix,ly1,1,ie) - x(ix,ly1-1,1,ie))
|
|---|
| 374 | enddo
|
|---|
| 375 | do iy=2,ly1-1
|
|---|
| 376 | x(1 ,iy,1,ie)=abs(x(1 ,iy,1,ie) - x(2 ,iy,1,ie))
|
|---|
| 377 | x(lx1,iy,1,ie)=abs(x(lx1,iy,1,ie) - x(lx1-1,iy,1,ie))
|
|---|
| 378 | enddo
|
|---|
| 379 | c
|
|---|
| 380 | endif
|
|---|
| 381 | enddo
|
|---|
| 382 | return
|
|---|
| 383 | end
|
|---|
| 384 | c-----------------------------------------------------------------------
|
|---|
| 385 | subroutine faces(a,s,ie,iface,nx,ny,nz)
|
|---|
| 386 | C
|
|---|
| 387 | C Scale face(IFACE,IE) of array A by s.
|
|---|
| 388 | C IFACE is the input in the pre-processor ordering scheme.
|
|---|
| 389 | C
|
|---|
| 390 | INCLUDE 'SIZE'
|
|---|
| 391 | DIMENSION A(NX,NY,NZ,LELT)
|
|---|
| 392 | CALL FACIND (KX1,KX2,KY1,KY2,KZ1,KZ2,NX,NY,NZ,IFACE)
|
|---|
| 393 | DO 100 IZ=KZ1,KZ2
|
|---|
| 394 | DO 100 IY=KY1,KY2
|
|---|
| 395 | DO 100 IX=KX1,KX2
|
|---|
| 396 | A(IX,IY,IZ,IE)=S*A(IX,IY,IZ,IE)
|
|---|
| 397 | 100 CONTINUE
|
|---|
| 398 | RETURN
|
|---|
| 399 | END
|
|---|
| 400 | c-----------------------------------------------------------------------
|
|---|
