| 1 | #ifndef NOMPIIO
|
|---|
| 2 |
|
|---|
| 3 | subroutine gfldr(sourcefld)
|
|---|
| 4 | c
|
|---|
| 5 | c generic field file reader
|
|---|
| 6 | c reads sourcefld and interpolates all avaiable fields
|
|---|
| 7 | c onto current mesh
|
|---|
| 8 | c
|
|---|
| 9 | c memory requirement:
|
|---|
| 10 | c nelgs*nxs**ldim < np*(4*lelt*lx1**ldim)
|
|---|
| 11 | c
|
|---|
| 12 | include 'SIZE'
|
|---|
| 13 | include 'TOTAL'
|
|---|
| 14 | include 'RESTART'
|
|---|
| 15 | include 'GFLDR'
|
|---|
| 16 |
|
|---|
| 17 | character sourcefld*(*)
|
|---|
| 18 |
|
|---|
| 19 | common /scrcg/ pm1(lx1*ly1*lz1,lelv)
|
|---|
| 20 | common /nekmpi/ nidd,npp,nekcomm,nekgroup,nekreal
|
|---|
| 21 |
|
|---|
| 22 | character*1 hdr(iHeaderSize)
|
|---|
| 23 |
|
|---|
| 24 | integer*8 dtmp8
|
|---|
| 25 |
|
|---|
| 26 | logical if_byte_swap_test
|
|---|
| 27 | real*4 bytetest
|
|---|
| 28 |
|
|---|
| 29 | etime_t = dnekclock_sync()
|
|---|
| 30 | if(nio.eq.0) write(6,*) 'call gfldr ',trim(sourcefld)
|
|---|
| 31 |
|
|---|
| 32 | ! open source field file
|
|---|
| 33 | ierr = 0
|
|---|
| 34 | if(nid.eq.0) then
|
|---|
| 35 | open (90,file=sourcefld,status='old',err=100)
|
|---|
| 36 | close(90)
|
|---|
| 37 | goto 101
|
|---|
| 38 | 100 ierr = 1
|
|---|
| 39 | 101 endif
|
|---|
| 40 | call err_chk(ierr,' Cannot open source fld file!$')
|
|---|
| 41 | call byte_open_mpi(sourcefld,fldh_gfldr,.true.,ierr)
|
|---|
| 42 |
|
|---|
| 43 | ! read and parse header
|
|---|
| 44 | call byte_read_mpi(hdr,iHeaderSize/4,0,fldh_gfldr,ierr)
|
|---|
| 45 | call byte_read_mpi(bytetest,1,0,fldh_gfldr,ierr)
|
|---|
| 46 |
|
|---|
| 47 | call mfi_parse_hdr(hdr,ierr)
|
|---|
| 48 | call err_chk(ierr,' Invalid header!$')
|
|---|
| 49 | ifbswp = if_byte_swap_test(bytetest,ierr)
|
|---|
| 50 | call err_chk(ierr,' Invalid endian tag!$')
|
|---|
| 51 |
|
|---|
| 52 | nelgs = nelgr
|
|---|
| 53 | nxs = nxr
|
|---|
| 54 | nys = nyr
|
|---|
| 55 | nzs = nzr
|
|---|
| 56 | if(nzs.gt.1) then
|
|---|
| 57 | ldims = 3
|
|---|
| 58 | else
|
|---|
| 59 | ldims = 2
|
|---|
| 60 | endif
|
|---|
| 61 | if (ifgtim) time = timer
|
|---|
| 62 |
|
|---|
| 63 | ! distribute elements across all ranks
|
|---|
| 64 | nels = nelgs/np
|
|---|
| 65 | do i = 0,mod(nelgs,np)-1
|
|---|
| 66 | if(i.eq.nid) nels = nels + 1
|
|---|
| 67 | enddo
|
|---|
| 68 | nxyzs = nxs*nys*nzs
|
|---|
| 69 | dtmp8 = nels
|
|---|
| 70 | ntots_b = dtmp8*nxyzs*wdsizr
|
|---|
| 71 | rankoff_b = igl_running_sum(nels) - dtmp8
|
|---|
| 72 | rankoff_b = rankoff_b*nxyzs*wdsizr
|
|---|
| 73 | dtmp8 = nelgs
|
|---|
| 74 | nSizeFld_b = dtmp8*nxyzs*wdsizr
|
|---|
| 75 | noff0_b = iHeaderSize + iSize + iSize*dtmp8
|
|---|
| 76 |
|
|---|
| 77 | ! do some checks
|
|---|
| 78 | if(ldims.ne.ldim)
|
|---|
| 79 | $ call exitti('ldim of source does not match target!$',0)
|
|---|
| 80 | if(ntots_b/wdsize .gt. ltots) then
|
|---|
| 81 | dtmp8 = nelgs
|
|---|
| 82 | lelt_req = dtmp8*nxs*nys*nzs / (np*ltots/lelt)
|
|---|
| 83 | lelt_req = lelt_req + 1
|
|---|
| 84 | if(nio.eq.0) write(6,*)
|
|---|
| 85 | $ 'ABORT: buffer too small, increase lelt > ', lelt_req
|
|---|
| 86 | call exitt
|
|---|
| 87 | endif
|
|---|
| 88 |
|
|---|
| 89 | ifldpos = 0
|
|---|
| 90 | if(ifgetxr) then
|
|---|
| 91 | ! read source mesh coordinates
|
|---|
| 92 | call gfldr_getxyz(xm1s,ym1s,zm1s)
|
|---|
| 93 | ifldpos = ldim
|
|---|
| 94 | else
|
|---|
| 95 | call exitti('source does not contain a mesh!$',0)
|
|---|
| 96 | endif
|
|---|
| 97 |
|
|---|
| 98 | if(if_full_pres) then
|
|---|
| 99 | call exitti('no support for if_full_pres!$',0)
|
|---|
| 100 | endif
|
|---|
| 101 |
|
|---|
| 102 | ! initialize interpolation tool using source mesh
|
|---|
| 103 | nxf = 2*nxs
|
|---|
| 104 | nyf = 2*nys
|
|---|
| 105 | nzf = 2*nzs
|
|---|
| 106 | nhash = nels*nxs*nys*nzs
|
|---|
| 107 | nmax = 128
|
|---|
| 108 |
|
|---|
| 109 | call fgslib_findpts_setup(inth_gfldr,nekcomm,np,ldim,
|
|---|
| 110 | & xm1s,ym1s,zm1s,nxs,nys,nzs,
|
|---|
| 111 | & nels,nxf,nyf,nzf,bb_t,
|
|---|
| 112 | & nhash,nhash,nmax,tol)
|
|---|
| 113 |
|
|---|
| 114 | ! read source fields and interpolate
|
|---|
| 115 | if(ifgetur) then
|
|---|
| 116 | if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'reading vel'
|
|---|
| 117 | ntot = nx1*ny1*nz1*nelv
|
|---|
| 118 | call gfldr_getfld(vx,vy,vz,ntot,ldim,ifldpos+1)
|
|---|
| 119 | ifldpos = ifldpos + ldim
|
|---|
| 120 | endif
|
|---|
| 121 | if(ifgetpr) then
|
|---|
| 122 | if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'reading pr'
|
|---|
| 123 | ntot = nx1*ny1*nz1*nelv
|
|---|
| 124 | call gfldr_getfld(pm1,dum,dum,ntot,1,ifldpos+1)
|
|---|
| 125 | ifldpos = ifldpos + 1
|
|---|
| 126 | if (ifaxis) call axis_interp_ic(pm1)
|
|---|
| 127 | call map_pm1_to_pr(pm1,1)
|
|---|
| 128 | endif
|
|---|
| 129 | if(ifgettr .and. ifheat) then
|
|---|
| 130 | if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'reading temp'
|
|---|
| 131 | ntot = nx1*ny1*nz1*nelfld(2)
|
|---|
| 132 | call gfldr_getfld(t(1,1,1,1,1),dum,dum,ntot,1,ifldpos+1)
|
|---|
| 133 | ifldpos = ifldpos + 1
|
|---|
| 134 | endif
|
|---|
| 135 | do i = 1,ldimt-1
|
|---|
| 136 | if(ifgtpsr(i)) then
|
|---|
| 137 | if(nid.eq.0 .and. loglevel.gt.2)
|
|---|
| 138 | $ write(6,*) 'reading scalar',i
|
|---|
| 139 | ntot = nx1*ny1*nz1*nelfld(i+2)
|
|---|
| 140 | call gfldr_getfld(t(1,1,1,1,i+1),dum,dum,ntot,1,ifldpos+1)
|
|---|
| 141 | ifldpos = ifldpos + 1
|
|---|
| 142 | endif
|
|---|
| 143 | enddo
|
|---|
| 144 |
|
|---|
| 145 | call byte_close_mpi(fldh_gfldr,ierr)
|
|---|
| 146 | etime_t = dnekclock_sync() - etime_t
|
|---|
| 147 | call fgslib_findpts_free(inth_gfldr)
|
|---|
| 148 | if(nio.eq.0) write(6,'(A,1(1g9.2),A)')
|
|---|
| 149 | & ' done :: gfldr ', etime_t, ' sec'
|
|---|
| 150 |
|
|---|
| 151 | return
|
|---|
| 152 | end
|
|---|
| 153 | c-----------------------------------------------------------------------
|
|---|
| 154 | subroutine gfldr_getxyz(xout,yout,zout)
|
|---|
| 155 |
|
|---|
| 156 | include 'SIZE'
|
|---|
| 157 | include 'GFLDR'
|
|---|
| 158 | include 'RESTART'
|
|---|
| 159 |
|
|---|
| 160 | real xout(*)
|
|---|
| 161 | real yout(*)
|
|---|
| 162 | real zout(*)
|
|---|
| 163 |
|
|---|
| 164 | integer*8 ioff_b
|
|---|
| 165 |
|
|---|
| 166 |
|
|---|
| 167 | ioff_b = noff0_b + ldim*rankoff_b
|
|---|
| 168 | call byte_set_view(ioff_b,fldh_gfldr)
|
|---|
| 169 |
|
|---|
| 170 | nread = ldim*ntots_b/4
|
|---|
| 171 | call byte_read_mpi(bufr,nread,-1,fldh_gfldr,ierr)
|
|---|
| 172 | if(ifbswp) then
|
|---|
| 173 | if(wdsizr.eq.4) call byte_reverse (bufr,nread,ierr)
|
|---|
| 174 | if(wdsizr.eq.8) call byte_reverse8(bufr,nread,ierr)
|
|---|
| 175 | endif
|
|---|
| 176 |
|
|---|
| 177 | call gfldr_buf2vi (xout,1,bufr,ldim,wdsizr,nels,nxyzs)
|
|---|
| 178 | call gfldr_buf2vi (yout,2,bufr,ldim,wdsizr,nels,nxyzs)
|
|---|
| 179 | if(ldim.eq.3)
|
|---|
| 180 | $ call gfldr_buf2vi(zout,3,bufr,ldim,wdsizr,nels,nxyzs)
|
|---|
| 181 |
|
|---|
| 182 | return
|
|---|
| 183 | end
|
|---|
| 184 | c-----------------------------------------------------------------------
|
|---|
| 185 | subroutine gfldr_getfld(out1,out2,out3,nout,nldim,ifldpos)
|
|---|
| 186 |
|
|---|
| 187 | include 'SIZE'
|
|---|
| 188 | include 'GEOM'
|
|---|
| 189 | include 'GFLDR'
|
|---|
| 190 | include 'RESTART'
|
|---|
| 191 |
|
|---|
| 192 | real out1(*)
|
|---|
| 193 | real out2(*)
|
|---|
| 194 | real out3(*)
|
|---|
| 195 |
|
|---|
| 196 | integer*8 ioff_b
|
|---|
| 197 |
|
|---|
| 198 | logical ifpts
|
|---|
| 199 |
|
|---|
| 200 | integer icalld
|
|---|
| 201 | save icalld
|
|---|
| 202 | data icalld /0/
|
|---|
| 203 |
|
|---|
| 204 |
|
|---|
| 205 | ifpts = .false.
|
|---|
| 206 | if(icalld.eq.0) then
|
|---|
| 207 | ifpts = .true. ! find points
|
|---|
| 208 | icalld = 1
|
|---|
| 209 | endif
|
|---|
| 210 |
|
|---|
| 211 | ! read field data from source fld file
|
|---|
| 212 | ioff_b = noff0_b + (ifldpos-1)*nSizeFld_b
|
|---|
| 213 | ioff_b = ioff_b + nldim*rankoff_b
|
|---|
| 214 | call byte_set_view(ioff_b,fldh_gfldr)
|
|---|
| 215 | nread = nldim*ntots_b/4
|
|---|
| 216 | call byte_read_mpi(bufr,nread,-1,fldh_gfldr,ierr)
|
|---|
| 217 | if(ifbswp) then
|
|---|
| 218 | if(wdsizr.eq.4) call byte_reverse (bufr,nread,ierr)
|
|---|
| 219 | if(wdsizr.eq.8) call byte_reverse8(bufr,nread,ierr)
|
|---|
| 220 | endif
|
|---|
| 221 |
|
|---|
| 222 | ! interpolate onto current mesh
|
|---|
| 223 | call gfldr_buf2vi (buffld,1,bufr,nldim,wdsizr,nels,nxyzs)
|
|---|
| 224 | call gfldr_intp (out1,nout,buffld,ifpts)
|
|---|
| 225 | if(nldim.eq.1) return
|
|---|
| 226 |
|
|---|
| 227 | call gfldr_buf2vi (buffld,2,bufr,nldim,wdsizr,nels,nxyzs)
|
|---|
| 228 | call gfldr_intp (out2,nout,buffld,.false.)
|
|---|
| 229 | if(nldim.eq.2) return
|
|---|
| 230 |
|
|---|
| 231 | if(nldim.eq.3) then
|
|---|
| 232 | call gfldr_buf2vi(buffld,3,bufr,nldim,wdsizr,nels,nxyzs)
|
|---|
| 233 | call gfldr_intp (out3,nout,buffld,.false.)
|
|---|
| 234 | endif
|
|---|
| 235 |
|
|---|
| 236 | return
|
|---|
| 237 | end
|
|---|
| 238 | c-----------------------------------------------------------------------
|
|---|
| 239 | subroutine gfldr_buf2vi(vi,index,buf,ldim,wds,nel,nxyz)
|
|---|
| 240 |
|
|---|
| 241 | real vi(*)
|
|---|
| 242 | real*4 buf(*)
|
|---|
| 243 | integer wds
|
|---|
| 244 |
|
|---|
| 245 |
|
|---|
| 246 | do iel = 1,nel
|
|---|
| 247 | j = (iel-1)*nxyz
|
|---|
| 248 | k = (iel-1)*ldim*nxyz
|
|---|
| 249 |
|
|---|
| 250 | if(index.eq.2) k = k+nxyz
|
|---|
| 251 | if(index.eq.3) k = k+2*nxyz
|
|---|
| 252 |
|
|---|
| 253 | if(wds.eq.4) call copy4r(vi(j+1),buf(k+1) ,nxyz)
|
|---|
| 254 | if(wds.eq.8) call copy (vi(j+1),buf(2*k+1),nxyz)
|
|---|
| 255 | enddo
|
|---|
| 256 |
|
|---|
| 257 | return
|
|---|
| 258 | end
|
|---|
| 259 | c-----------------------------------------------------------------------
|
|---|
| 260 | subroutine gfldr_intp(fieldout,nout,fieldin,iffpts)
|
|---|
| 261 |
|
|---|
| 262 | include 'SIZE'
|
|---|
| 263 | include 'RESTART'
|
|---|
| 264 | include 'GEOM'
|
|---|
| 265 | include 'GFLDR'
|
|---|
| 266 |
|
|---|
| 267 | real fieldout(nout)
|
|---|
| 268 | real fieldin (*)
|
|---|
| 269 | logical iffpts
|
|---|
| 270 |
|
|---|
| 271 | integer*8 i8glsum,nfail,nfail_sum
|
|---|
| 272 |
|
|---|
| 273 | if(iffpts) then ! locate points (iel,iproc,r,s,t)
|
|---|
| 274 | nfail = 0
|
|---|
| 275 | toldist = 5e-6
|
|---|
| 276 | if(wdsizr.eq.8) toldist = 5e-14
|
|---|
| 277 |
|
|---|
| 278 | ntot = lx1*ly1*lz1*nelt
|
|---|
| 279 | call fgslib_findpts(inth_gfldr,
|
|---|
| 280 | & grcode,1,
|
|---|
| 281 | & gproc,1,
|
|---|
| 282 | & gelid,1,
|
|---|
| 283 | & grst,ldim,
|
|---|
| 284 | & gdist,1,
|
|---|
| 285 | & xm1,1,
|
|---|
| 286 | & ym1,1,
|
|---|
| 287 | & zm1,1,ntot)
|
|---|
| 288 |
|
|---|
| 289 | do i=1,ntot
|
|---|
| 290 | if(grcode(i).eq.1 .and. sqrt(gdist(i)).gt.toldist)
|
|---|
| 291 | & nfail = nfail + 1
|
|---|
| 292 | if(grcode(i).eq.2) nfail = nfail + 1
|
|---|
| 293 | enddo
|
|---|
| 294 |
|
|---|
| 295 | nfail_sum = i8glsum(nfail,1)
|
|---|
| 296 | if(nfail_sum.gt.0) then
|
|---|
| 297 | if(nio.eq.0) write(6,*)
|
|---|
| 298 | & ' WARNING: Unable to find all mesh points in source fld ',
|
|---|
| 299 | & nfail_sum
|
|---|
| 300 | endif
|
|---|
| 301 | endif
|
|---|
| 302 |
|
|---|
| 303 | ! evaluate inut field at given points
|
|---|
| 304 | npt = nout
|
|---|
| 305 | call fgslib_findpts_eval(inth_gfldr,
|
|---|
| 306 | & fieldout,1,
|
|---|
| 307 | & grcode,1,
|
|---|
| 308 | & gproc,1,
|
|---|
| 309 | & gelid,1,
|
|---|
| 310 | & grst,ldim,npt,
|
|---|
| 311 | & fieldin)
|
|---|
| 312 |
|
|---|
| 313 | return
|
|---|
| 314 | end
|
|---|
| 315 |
|
|---|
| 316 | #endif
|
|---|
