1! readnblk.f (pronounce "Reed and Block Dot Eff") contains: 2! 3! module readarrays ("Red Arrays") -- contains the arrays that 4! are read in from binary files but not immediately blocked 5! through pointers. 6! 7! subroutine readnblk ("Reed and Block") -- allocates space for 8! and reads data to be contained in module readarrays. Reads 9! all remaining data and blocks them with pointers. 10! 11 12 13 module readarrays 14 15 real*8, allocatable :: point2x(:,:) 16 real*8, allocatable :: qold(:,:) 17 real*8, allocatable :: dwal(:) 18 real*8, allocatable :: errors(:,:) 19 real*8, allocatable :: ybar(:,:) 20 real*8, allocatable :: yphbar(:,:,:) 21 real*8, allocatable :: vort(:,:) 22 real*8, allocatable :: uold(:,:) 23 real*8, allocatable :: acold(:,:) 24 integer, allocatable :: iBCtmp(:) 25 real*8, allocatable :: BCinp(:,:) 26 27 integer, allocatable :: point2ilwork(:) 28 integer, allocatable :: nBC(:) 29 integer, allocatable :: point2iper(:) 30 integer, allocatable :: point2ifath(:) 31 integer, allocatable :: point2nsons(:) 32 33 end module 34 35 36 37 subroutine readnblk 38! 39 use readarrays 40 include "commonM2N.h" 41 include "mpif.h" 42 43! 44 real*8, allocatable :: xread(:,:), qread(:,:), acread(:,:) 45 real*8, allocatable :: qreadN(:), qoldN(:,:) 46 real*8, allocatable :: pID(:), vID(:), dwalN(:), ybarN(:,:) 47 real*8, allocatable :: errorsN(:,:), yphbarN(:,:,:), vortN(:,:) 48 real*8, allocatable :: qTota2as(:,:), qTota2ar(:,:) 49 real*8, allocatable :: qrecv(:), qsend(:) 50 real*8 :: qreadN1 51 integer, allocatable :: iperread(:), iBCtmpread(:) 52 integer, allocatable :: ilworkread(:), nBCread(:) 53 integer, allocatable :: indexpart(:) 54 integer, allocatable :: getinfo(:) 55 character*10 cname2 56 character*8 mach2 57 character*30 fmt1 58 character*255 fname1,fnamer, fnamelr 59 character*255 warning 60 integer :: igeomBAK, ibndc, irstin, ierr 61 integer :: icountN(numpe) ! integers read from headers 62 integer :: intfromfile(50) ! integers read from headers 63 integer :: maxvID(numpe),maxvIDout(numpe) 64 integer :: nmap, ndim 65 integer :: maxicountN, maxicountNglob, numvar, sumvar 66 integer :: ncount, ifill, tag 67 integer :: stat (MPI_STATUS_SIZE) 68 integer :: ndofdwal,ndofybar,ndoferrors,ndofyphbar,ndofvort 69 logical :: exmap, exinput 70 integer :: my_local_comm, my_local_size, my_local_rank 71!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 72!cccccccccccccccccccccccccccccccccccccccc New PhastaIO Definition Part ccccccccccccccccccccccccccccccccccccccccc 73 74 integer :: descriptor, descriptorG, GPID, color, nfiles, nfields 75 integer :: numparts, nppf, islesseqN 76 integer :: ierr_io, numprocs, itmp, itmp2 77 integer :: i, nonzero 78 integer :: ireducemethod, isbinary, irstart, irstartmap, iybar 79 integer :: ierror, numphavg, ivort, idwal, idebug, iphavg, irankN 80 character*255 fname2, temp2 81 character*64 temp1 82 83!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 84!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 85 86! 87! Set some parameters for the reduce operation 88! 89! 0 for communication based approach with mpi_send and mpi_recv (recommended), 1 for file based approach (not maintained any more) 90 ireducemethod = 0 91 92! 1 to write the restartMap* files in binary format, 0 in ascii (in conjunction with ireducemethod = 1) 93 isbinary = 1 94! 95! Read the input parameter for the the M2N reduction 96! 97! open(unit=72,file='numstart.dat',status='old') 98! read(72,*) irstart 99! close(72) 100 101 if(myrank == 0) then 102 write(*,*) 'Reading M2N_input.dat' 103 fnamer='M2N_input.dat' 104 inquire(file=fnamer,exist=exinput) 105 if(exinput) then 106 open(unit=72,file=fnamer,status='old') 107 read(72,*) irstart 108 read(72,*) irstartmap 109 read(72,*) iybar 110 read(72,*) ierror 111 read(72,*) numphavg 112 read(72,*) ivort 113 read(72,*) idwal 114 read(72,*) idebug 115 close(72) 116 else 117 write(*,*) 'ERROR: Input file ', 118 & trim(fnamer),' does not exist!' 119 endif 120 121 endif 122 123 call mpi_barrier(mpi_comm_world, ierr) 124 call mpi_bcast(exinput,1,MPI_LOGICAL,0,mpi_comm_world,ierr) 125 if(.not. exinput) then ! M2N_input.dat does not exist. Quit 126 return 127 else ! broadcast the information read by rank 0 128 call mpi_bcast(irstart,1,MPI_INTEGER,0,mpi_comm_world,ierr) 129 call mpi_bcast(irstartmap,1,MPI_INTEGER,0,mpi_comm_world,ierr) 130 call mpi_bcast(iybar,1,MPI_INTEGER,0,mpi_comm_world,ierr) 131 call mpi_bcast(ierror,1,MPI_INTEGER,0,mpi_comm_world,ierr) 132 call mpi_bcast(numphavg,1,MPI_INTEGER,0,mpi_comm_world,ierr) 133 call mpi_bcast(ivort,1,MPI_INTEGER,0,mpi_comm_world,ierr) 134 call mpi_bcast(idwal,1,MPI_INTEGER,0,mpi_comm_world,ierr) 135 call mpi_bcast(idebug,1,MPI_INTEGER,0,mpi_comm_world,ierr) 136 endif 137 138 if(myrank == 0 ) then 139 ! Print some info 140 write(*,*) 'The solution field is reduced by default' 141 142 if(iybar .gt. 0) then 143 write(*,*) 'The ybar field will also be reduced' 144 iybar = 1 ! security 145 else 146 write(*,*) 'The ybar field will NOT be reduced' 147 endif 148 149 if(ierror .gt. 0) then 150 write(*,*) 'The error field will also be reduced' 151 ierror = 1 ! security 152 else 153 write(*,*) 'The error field will NOT be reduced' 154 endif 155 156 if(numphavg .gt. 0) then 157 write(*,*) 'The phase average fields (',numphavg, 158 & ') will also be reduced' 159 else 160 write(*,*) 'The phase average fields will NOT be reduced' 161 endif 162 163 if(ivort .gt. 0) then 164 write(*,*) 'The vorticity field will also be reduced' 165 ivort = 1 ! security 166 else 167 write(*,*) 'The vorticity field will NOT be reduced' 168 endif 169 170 if(idwal .gt. 0) then 171 write(*,*) 'The dwal field will also be reduced' 172 idwal = 1 ! security 173 else 174 write(*,*) 'The dwal field will NOT be reduced' 175 endif 176 write(*,*) '' 177 endif 178 179 call mpi_barrier(mpi_comm_world, ierr) 180 if(myrank == 0) then 181 write(*,*) 'Reading the geombc-dat files' 182 endif 183 184 lstep=irstart ! in case restart files have no fields 185 186 nfiles = nsynciofiles 187 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 188 189 color = int(myrank/(numparts/nfiles)) !Should call the color routine in SyncIO here 190 itmp2 = int(log10(float(color+1)))+1 191 write (temp2,"('(''geombc-dat.'',i',i1,')')") itmp2 192 write (fnamer,temp2) (color+1) 193 fnamer = trim(fnamer)//char(0) 194 195 itwo=2 196 ione=1 197 ieleven=11 198 itmp = int(log10(float(myrank+1)))+1 199 200 call queryphmpiio(fnamer, nfields, nppf); 201 if (myrank == 0) then 202 write(*,*) 'Number of fields in geombc-dat: ',nfields 203 write(*,*) 'Number of parts per file geombc-dat: ',nppf 204 endif 205 call initphmpiio( nfields, nppf, nfiles, igeom, 206 & 'read'//char(0)) 207 call openfile( fnamer, 'read'//char(0), igeom ) 208 209 write (temp1,"('(''number of nodes@'',i',i1,',A1)')") itmp 210 write (fname2,temp1) (myrank+1),'?' 211 call readheader(igeom,fname2//char(0),numnp,ione, 212 & 'integer'//char(0), iotype) 213 214 write (temp1,"('(''number of modes@'',i',i1,',A1)')") itmp 215 write (fname2,temp1) (myrank+1),'?' 216 call readheader(igeom,fname2//char(0),nshg,ione, 217 & 'integer'//char(0), iotype) 218 219 write (temp1,"('(''number of interior elements@'',i',i1,',A1)')") 220 & itmp 221 write (fname2,temp1) (myrank+1),'?' 222 call readheader(igeom,fname2//char(0),numel,ione, 223 & 'integer'//char(0), iotype) 224 225 write (temp1,"('(''number of boundary elements@'',i',i1,',A1)')") 226 & itmp 227 write (fname2,temp1) (myrank+1),'?' 228 call readheader(igeom,fname2//char(0),numelb,ione, 229 & 'integer'//char(0),iotype) 230 231 write (temp1, 232 & "('(''maximum number of element nodes@'',i',i1,',A1)')") itmp 233 write (fname2,temp1) (myrank+1),'?' 234 call readheader(igeom,fname2//char(0),nen,ione, 235 &'integer'//char(0),iotype) 236 237 write (temp1,"('(''number of interior tpblocks@'',i',i1,',A1)')") 238 & itmp 239 write (fname2,temp1) (myrank+1),'?' 240 call readheader(igeom,fname2//char(0),nelblk,ione, 241 & 'integer'//char(0) ,iotype) 242 243 write (temp1,"('(''number of boundary tpblocks@'',i',i1,',A1)')") 244 & itmp 245 write (fname2,temp1) (myrank+1),'?' 246 call readheader(igeom,fname2//char(0),nelblb,ione, 247 & 'integer'//char(0), iotype) 248 249 write (temp1, 250 & "('(''number of nodes with Dirichlet BCs@'',i',i1,',A1)')") itmp 251 write (fname2,temp1) (myrank+1),'?' 252 call readheader(igeom,fname2//char(0),numpbc,ione, 253 & 'integer'//char(0),iotype) 254 255 write (temp1,"('(''number of shape functions@'',i',i1,',A1)')") 256 & itmp 257 write (fname2,temp1) (myrank+1),'?' 258 call readheader(igeom,fname2//char(0),ntopsh,ione, 259 & 'integer'//char(0),iotype) 260 261!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 262! 263!.... calculate the maximum number of boundary element nodes 264! 265 nenb = 0 266 do i = 1, melCat 267 if (nen .eq. nenCat(i,nsd)) nenb = max(nenCat(i,nsd-1), nenb) 268 enddo 269! 270! 271!.... setup some useful constants 272! 273 I3nsd = nsd / 3 ! nsd=3 integer flag 274 E3nsd = float(I3nsd) ! nsd=3 real flag 275! 276 if(matflg(1,1).lt.0) then 277 nflow = nsd + 1 278 else 279 nflow = nsd + 2 280 endif 281 ndof = nsd + 2 282 nsclr=impl(1)/100 283 ndof=ndof+nsclr ! number of sclr transport equations to solve 284 285 ndofBC = ndof + I3nsd ! dimension of BC array 286 ndiBCB = 2 ! dimension of iBCB array 287 ndBCB = ndof + 1 ! dimension of BCB array 288! 289 nsymdf = (ndof*(ndof + 1)) / 2 ! symm. d.o.f.'s 290! 291!.... ----------------------> Communication tasks <-------------------- 292! 293 call closefile( igeom, "read"//char(0) ) 294 call finalizephmpiio( igeom ) 295 296! 297!.... Read restart files for the solution, dwal, error and ybar 298! 299 call mpi_barrier(mpi_comm_world, ierr) 300 if(myrank == 0) then 301 write(*,*)'Reading the restart-dat files for the s-d-e-y fields' 302 endif 303 304 itmp=1 305 if (irstart .gt. 0) itmp = int(log10(float(irstart+1)))+1 306 307 write (fmt1,"('(''restart-dat.'',i',i1,',1x)')") itmp 308 309 write (fnamer,fmt1) irstart 310 fnamer = trim(fnamer) // cname2(color+1) 311 312 call queryphmpiio(fnamer//char(0), nfields, nppf); 313 if (myrank == 0) then 314 write(*,*) 'Number of fields in restart-dat: ',nfields 315 write(*,*) 'Number of parts per file restart-dat: ',nppf 316 endif 317 call initphmpiio(nfields,nppf,nfiles,descriptor, 318 & 'read'//char(0)) 319 call openfile( fnamer//char(0) , 320 & 'read'//char(0), descriptor ) 321 322 ithree=3 323! call creadlist(irstin,ithree,nshg2,ndof2,lstep) 324 325 itmp = int(log10(float(myrank+1)))+1 326 write (temp1,"('(''solution@'',i',i1,',A1)')") itmp 327 write (fname1,temp1) (myrank+1),'?' 328 fname1 = trim(fname1) 329 330! print *, "Solution is : ", fname1 331 332 intfromfile=0 333 call readheader(descriptor,fname1//char(0) ,intfromfile, 334 & ithree,'integer'//char(0), iotype) 335! 336!.... read the values of primitive variables into q 337! 338 339! print *, "intfromfile(1) is ", intfromfile(1) 340! print *, "intfromfile(2) is ", intfromfile(2) 341! print *, "intfromfile(3) is ", intfromfile(3) 342 343 if(intfromfile(1).ne.0) then 344 nshg2=intfromfile(1) 345 ndof2=intfromfile(2) 346 lstep=intfromfile(3) 347 if(nshg2.ne.nshg) then 348 write(*,*) 'ERROR: mistmatch between nshg '// 349 & 'from the geombc and restart files',nshg,nshg2 350 endif 351 352 ndof=ndof2 353 allocate( qold(nshg,ndof2) ) 354 allocate( qread(nshg,ndof2) ) 355 iqsiz=nshg*ndof2 356 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 357 & 'double'//char(0),iotype) 358 qold(:,1:ndof)=qread(:,1:ndof) 359 deallocate(qread) 360 else 361 if (myrank.eq.master) then 362 if (matflg(1,1).eq.0) then ! compressible 363 warning='WARNING: Solution is set to zero '// 364 & '(with p and T to one)' 365 else 366 warning='WARNING: Solution is set to zero' 367 endif 368 write(*,*) warning 369 endif 370 qold=zero 371 if (matflg(1,1).eq.0) then ! compressible 372 qold(:,1)=one ! avoid zero pressure 373 qold(:,nflow)=one ! avoid zero temperature 374 endif 375 endif 376 377! 378! Above this line is the usual loading of fields in readnblk. Quite a bit of stuff has been removed. 379! 380 381 382! 383! follow the usual convention for loading the ybar 384! 385 ndofybar=0 386 if(iybar == 1) then 387 itmp = int(log10(float(myrank+1)))+1 388 write (temp1,"('(''ybar@'',i',i1,',a1)')") itmp 389 write (fname1,temp1) (myrank+1),'?' 390 fname1 = trim(fname1) 391 intfromfile=0 392 call readheader(descriptor,fname1//char(0),intfromfile, 393 & ithree,'integer'//char(0),iotype) 394 395 !nshg2=intfromfile(1) 396 ndofybar=intfromfile(2) 397 !lstep=intfromfile(3) 398 399 if(ndofybar .ne. 0) then 400 allocate( ybar(nshg,ndofybar) ) 401 iqsiz=nshg*ndofybar 402 call readdatablock(descriptor,fname1//char(0) ,ybar,iqsiz, 403 & 'double'//char(0),iotype) 404 else 405 if(myrank==0) then 406 write(*,*) 'WARNING: ybar is missing in the restart files' 407 endif 408 iybar = 0 409 endif 410 endif 411 412! 413! follow the usual convention for loading the error 414! 415 ndoferrors=0 416 if(ierror == 1) then 417 itmp = int(log10(float(myrank+1)))+1 418 write (temp1,"('(''errors@'',i',i1,',a1)')") itmp 419 write (fname1,temp1) (myrank+1),'?' 420 fname1 = trim(fname1) 421 intfromfile=0 422 call readheader(descriptor,fname1//char(0),intfromfile, 423 & ithree,'integer'//char(0),iotype) 424 425 !nshg2=intfromfile(1) 426 ndoferrors=intfromfile(2) 427 !lstep=intfromfile(3) 428 429 if(ndoferrors .ne. 0) then 430 allocate( errors(nshg,ndoferrors) ) 431 iqsiz=nshg*ndoferrors 432 call readdatablock(descriptor,fname1//char(0),errors,iqsiz, 433 & 'double'//char(0),iotype) 434 else 435 if(myrank==0) then 436 write(*,*) 'WARNING: errors is missing in the restart files' 437 endif 438 ierror = 0 439 endif 440 endif 441 442! 443! Read the phase_average fields 444! 445 ndofyphbar=0 446 if(numphavg .gt. 0) then 447 do iphavg = 1,numphavg 448 itmp = int(log10(float(myrank+1)))+1 449 itmp2 = int(log10(float(iphavg)))+1 450 write (temp1, 451 & "('(''phase_average'',i',i1,',''@'',i',i1,',A1)')") 452 & itmp2, itmp 453 write (fname1,temp1) iphavg,(myrank+1),'?' 454 fname1 = trim(fname1) 455 intfromfile=0 456 call readheader(descriptor,fname1//char(0),intfromfile, 457 & ithree,'integer'//char(0),iotype) 458 459 !nshg=intfromfile(1) !Do not use nshg and ndof from common.h here! 460 ndofyphbar=intfromfile(2) 461 !lstep=intfromfile(3) 462 463 if(ndofyphbar.ne.0) then 464 ! Allocate some memory for the first ts only 465 if(iphavg==1) then 466 allocate( yphbar(nshg,ndofyphbar,numphavg) ) 467 endif 468 469 allocate( qread(nshg,ndofyphbar) ) ; qread = 0.d0 470 iqsiz = nshg*ndofyphbar 471 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 472 & 'double'//char(0),iotype) 473 yphbar(:,:,iphavg) = qread(:,:) 474 deallocate(qread) 475 else 476 if(myrank==0) then 477 write(*,*)'WARNING: phase_average is missing '// 478 & 'in the restart files' 479 endif 480 numphavg = 0 481 if(iphavg > 1) then 482 deallocate(yphbar) 483 endif 484 exit 485 endif 486 enddo 487 endif 488 489! 490! follow the usual convention for loading the vorticity field 491! 492 ndofvort=0 493 if(ivort == 1) then 494 itmp = int(log10(float(myrank+1)))+1 495 write (temp1,"('(''vorticity@'',i',i1,',a1)')") itmp 496 write (fname1,temp1) (myrank+1),'?' 497 fname1 = trim(fname1) 498 intfromfile=0 499 call readheader(descriptor,fname1//char(0),intfromfile, 500 & ithree,'integer'//char(0),iotype) 501 502 !nshg2=intfromfile(1) 503 ndofvort=intfromfile(2) 504 !lstep=intfromfile(3) 505 506 if(ndofvort .ne. 0) then 507 allocate( vort(nshg,ndofvort) ) 508 iqsiz=nshg*ndofvort 509 call readdatablock(descriptor,fname1//char(0) ,vort ,iqsiz, 510 & 'double'//char(0),iotype) 511 else 512 if(myrank==0) then 513 write(*,*) 'WARNING: vorticity is missing '// 514 & 'in the restart files' 515 endif 516 ivort = 0 517 endif 518 endif 519 520 521! 522! follow the usual convention for loading the d2wal 523! 524 ndofdwal=0 525 if(idwal == 1 ) then 526 write (temp1,"('(''dwal@'',i',i1,',a1)')") 527 & itmp 528 write (fname1,temp1) (myrank+1),'?' 529 fname1 = trim(fname1) 530 intfromfile=0 531 call readheader(descriptor,fname1//char(0),intfromfile, 532 & ithree,'integer'//char(0),iotype) 533 534 !nshg2=intfromfile(1) 535 ndofdwal=intfromfile(2) 536 !lstep=intfromfile(3) 537 538 if(ndofdwal .ne. 0) then 539 if(ndofdwal.ne.1) then 540 warning='WARNING: ndofdwal not equal 1' 541 write(*,*) warning, ndofdwal 542 endif 543 allocate( dwal(nshg) ) 544 iqsiz=nshg 545 call readdatablock(descriptor,fname1//char(0), dwal, iqsiz, 546 & 'double'//char(0),iotype) 547 else 548 if(myrank==0) then 549 write(*,*) 'WARNING: dwal is missing in the restart files' 550 endif 551 idwal = 0 552 endif 553 endif 554 555! 556!.... close c-binary files 557! 558 call closefile( descriptor, "read"//char(0) ) 559 call finalizephmpiio( descriptor ) 560 561! 562!.... Read restart files for the vtx mapping between M and N parts 563! 564 call mpi_barrier(mpi_comm_world, ierr) 565 if(myrank == 0) then 566 write(*,*) 'Reading the restart-dat files for vtx mapping' 567 endif 568 569 itmp=1 570 if (irstartmap .gt. 0) itmp = int(log10(float(irstartmap+1)))+1 571 572 write (fmt1,"('(''restart-dat.'',i',i1,',1x)')") itmp 573 574 write (fnamer,fmt1) irstartmap 575 fnamer = trim(fnamer) // cname2(color+1) 576 577 call queryphmpiio(fnamer//char(0), nfields, nppf); 578 if (myrank == 0) then 579 write(*,*) 'Number of fields in restart-dat: ',nfields 580 write(*,*) 'Number of parts per file restart-dat: ',nppf 581 endif 582 call initphmpiio(nfields,nppf,nfiles,descriptor, 583 & 'read'//char(0)) 584 call openfile( fnamer//char(0) , 585 & 'read'//char(0), descriptor ) 586 587! 588! follow the usual convention for loading the partid 589! 590 itmp = int(log10(float(myrank+1)))+1 591 write (temp1,"('(''mapping_partid@'',i',i1,',A1)')") 592 & itmp 593 write (fname1,temp1) (myrank+1),'?' 594 fname1 = trim(fname1) 595 intfromfile=0 596 call readheader(descriptor,fname1//char(0),intfromfile, 597 & ithree,'integer'//char(0),iotype) 598 599 nshg2=intfromfile(1) 600 nmap=intfromfile(2) 601 !lstep=intfromfile(3) !lstep comes from another restart file 602 if(nmap.ne.1) then 603 warning='WARNING mapping_partid not equal 1' 604 write(*,*) warning , nmap, intfromfile(1), intfromfile(2) 605 endif 606 607 allocate( pID(nshg) ) 608 iqsiz=nshg 609 call readdatablock(descriptor,fname1//char(0) ,pID,iqsiz, 610 & 'double'//char(0),iotype) 611 612! 613! follow the usual convention for loading the vtxid 614! 615 write (temp1,"('(''mapping_vtxid@'',i',i1,',a1)')") 616 & itmp 617 write (fname1,temp1) (myrank+1),'?' 618 fname1 = trim(fname1) 619 intfromfile=0 620 call readheader(descriptor,fname1//char(0),intfromfile, 621 & ithree,'integer'//char(0),iotype) 622 623 nshg2=intfromfile(1) 624 nmap=intfromfile(2) 625 !lstep=intfromfile(3) !lstep comes from another restart file 626 if(nmap.ne.1) then 627 warning='warning mapping_partid not equal 1' 628 write(*,*) warning , nmap 629 endif 630 631 allocate( vid(nshg) ) 632 iqsiz=nshg 633 call readdatablock(descriptor,fname1//char(0) ,vid,iqsiz, 634 & 'double'//char(0),iotype) 635 636! 637!.... close c-binary files 638! 639 call closefile( descriptor, "read"//char(0) ) 640 call finalizephmpiio( descriptor ) 641 642! 643! Debugging information 644! 645! minvid=minval(vid) 646! write(*,*) minvid,myrank 647! do i=1,nshg 648! write(800+myrank,*) pID(i),vID(i),qold(i,1), qold(i,2) 649! ivID=vID(i)+0.1 650! ipID=pID(i)+0.1 651! write(900+myrank,*) ipID,ivID,qold(i,1), qold(i,2) 652! enddo 653 654! 655! here starts the mapping code 656! 657 658! first, on each of the M parts, count the number of vertices on each of the 659! N parts and put the result in icountN for each of the M parts 660! 661 call mpi_barrier(mpi_comm_world, ierr) 662 if(myrank == 0) then 663 write(*,*) 'Computing the size of the mapped partition' 664 endif 665 666 maxPID=0 667 icountN=0 668 do i=1,nshg 669 ifileN=pID(i)+1.1 670 icountN(ifileN)=icountN(ifileN)+1 671 maxPID=max(maxPID,ifileN) 672 enddo 673! 674! ALLREDUCE so that every processor knows N the size of the mapped partition 675! 676 call MPI_ALLREDUCE (maxPID, maxPIDout, 1, 677 & MPI_INT, MPI_MAX, mpi_comm_world, ierr) 678 679 irankN=maxPIDout 680 681 call mpi_barrier(mpi_comm_world, ierr) 682 if(myrank.eq.0) then 683 write (*,*) 'M2N reduction with M =', 684 & numpe, 'and N =',irankN 685 endif 686! 687! DEBUG echo to the screen the number of number of vtx that will be mapped to 688! each of the N for each of the M parts 689! 690 if (idebug == 1) then 691 call mpi_barrier(mpi_comm_world, ierr) 692 do i=1,numpe 693 if(myrank == i) then 694 write(*,926) myrank, (icountN(j),j=1,irankN) 695 endif 696 call mpi_barrier(mpi_comm_world, ierr) 697 enddo 698 endif 699 700! 701! Compute the max vertex ID for every N part 702! 703 call mpi_barrier(mpi_comm_world, ierr) 704 if(myrank == 0) then 705 write(*,*) 'Computing the max vertex ID for every N part' 706 endif 707 708 maxvID(:)=0 !only maxvID(1:irankN) should have non-zero entries at the end of the program 709 do i=1,nshg 710 ifileN=pID(i)+1.1 !shift pID to start from 1 711 intvID=(vID(i)+1.1) !shift vID to start from 1 712 maxvID(ifileN)=max(maxvID(ifileN),intvID) 713 enddo 714 715! 716! Sanity check for maxvID 717! 718 do i=irankN+1,numpe 719 if (maxvID(i) .ne. 0) then 720 write(*,*) 'ERROR: maxvID(',i,') not equal to 0 on rank', 721 & myrank 722 endif 723 enddo 724! 725! ALLREDUCE so that every processor knows the number of vtx on each of the N 726! parts in the reduced partition. No need to reduce maxvID(irankN+1:numpe), as it should be 0 727! 728 call MPI_ALLREDUCE (maxvID, maxvIDout, irankN, 729 & MPI_INT, MPI_MAX, mpi_comm_world, ierr) 730 731! 732! echo debug information to the screen 733! 734 if (idebug == 1) then 735 if(myrank.eq.0) write(*,*) 'Writing maxvID on each of N ranks' 736 if(myrank.eq.0) write(*,926) (maxvIDout(j),j=1,irankN) 737 endif 738 739! Beware! maxvIDout does not represent the largest vID present in the parts 740! but the largest vID referenced in the parts through the mapping. 741! On part boundaries, it is possible to have unreferenced vID that will need 742! to be fixed through P2P communication. 743 744! 745! Allocate the memory for qoldN on N parts handled by ranks 1..N. 746! Memory available per comupte node could be improved by chosing others ranks 747! 748 call mpi_barrier(mpi_comm_world, ierr) 749 if(myrank == 0) then 750 write(*,*) 'Allocating memory for the reduced s-d-e-y vectors' 751 endif 752 753 if(myrank.lt.irankN) then 754 nshgN=maxvIDout(myrank+1) ! number of vtx on my part computed above 755 756 allocate(qoldN(nshgN,ndof))! size of my final mapped field array 757 qoldN=-9.87654321e32 ! set to absurd, easily recognized value to find holes 758 759 if(iybar == 1) then 760 allocate(ybarN(nshgN,ndofybar)) ! size of my final mapped field array 761 ybarN=-9.87654321e32 ! set to absurd, easily recognized value to find holes 762 endif 763 764 if(ierror == 1) then 765 allocate(errorsN(nshgN,ndoferrors)) ! size of my final mapped field array 766 errorsN=-9.87654321e32 ! set to absurd, easily recognized value to find holes 767 endif 768 769 if(numphavg .gt. 0) then 770 allocate(yphbarN(nshgN,ndofyphbar,numphavg)) ! size of my final mapped field array 771 yphbarN=-9.87654321e32 ! set to absurd, easily recognized value to find holes 772 endif 773 774 if(ivort == 1) then 775 allocate( vortN(nshgN,ndofvort)) 776 vortN=-9.87654321e32 777 endif 778 779 if(idwal == 1) then 780 allocate(dwalN(nshgN)) ! size of my final mapped field array 781 dwalN=-9.87654321e32 ! set to absurd, easily recognized value to find holes 782 endif 783 784 785 else ! Not used but sane, as passed to write_M2N 786 allocate(qoldN(1,1)) 787 if(iybar == 1) then 788 allocate(ybarN(1,1)) 789 endif 790 if(ierror == 1) then 791 allocate(errorsN(1,1)) 792 endif 793 if(numphavg .gt. 0) then 794 allocate(yphbarN(1,1,1)) 795 endif 796 if(ivort == 1) then 797 allocate(vortN(1,1)) 798 endif 799 if(idwal == 1) then 800 allocate(dwalN(1)) 801 endif 802 endif 803 804! 805! Most important routine: build the qoldN array for the reduction of the solution from M to N 806! 807 call mpi_barrier(mpi_comm_world, ierr) 808 if(myrank == 0) then 809 write(*,*) 'Building the reduced field vectors' 810 endif 811 812 if (ireducemethod == 0) then 813! 814! Communication approach through mpi_isend and mpi_irecv to limit memory consumption 815! and avoid disk access 816! 817 numvar = 1+ndof+ndofybar+ndoferrors 818 & +(numphavg*ndofyphbar)+ndofvort+ndofdwal 819! write(*,*) numvar,ndof,ndofybar,ndoferrors,numphavg,ndofyphbar, 820! & ndofvort 821 call mpi_barrier(mpi_comm_world, ierr) 822 !1 for vtxid, ndof for solution, ndofdwal for dwal, ndoferrors for error, ndofybar for ybar 823 824 do irN = 0,irankN-1 ! every irN rank 1..N receives information from 825 826 827 ! Get which of the numpe-1 ranks have some information to communicate to rank irN 828 if (myrank == irN) then 829 allocate(getinfo(numpe)) 830 getinfo(:) = 0 831 endif 832 833 ifill = icountN(irN+1) 834 call mpi_gather(ifill, 1, MPI_INT, getinfo, 1, MPI_INT, 835 & irN, mpi_comm_world, ierr) 836 837! if (myrank == irN) then 838! write(*,*) getinfo(:) 839! endif 840 841 do irM = 0,numpe-1 ! all the irM ranks 1..M 842 843! DEBUG 844! if (myrank == 0) then 845! write(*,*) 'NEW data exchange between irN:', 846! & irN,'and irM:',irM 847! endif 848! call mpi_barrier(mpi_comm_world, ierr) 849!DEBUG 850 851 if (myrank == irN .and. myrank == irM) then ! no need to send the data. Already there so get them 852! write(*,*) ' myrank=irN=irM=',irN,' so treat data locally' 853 do i=1,nshg 854 ifileN = pID(i)+1.1 !shift to start from 1 855 if(ifileN == irN+1) then ! these data need to be packed because rank irN is expecting them 856 intvID = vID(i)+1.1 ! shift to start from 1 857 ! Solution qoldN 858 do idof=1,ndof 859 qoldN(intvID,idof) = qold(i,idof) 860 enddo 861 ! ybarN 862 if(iybar==1) then 863 do idof=1,ndofybar 864 ybarN(intvID,idof) = ybar(i,idof) 865 enddo 866 endif 867 ! errorsN 868 if(ierror==1) then 869 do idof=1,ndoferrors 870 errorsN(intvID,idof) = errors(i,idof) 871 enddo 872 endif 873 ! phase average 874 if(numphavg .gt. 0) then 875 do iphavg = 1,numphavg 876 do idof=1,ndofyphbar 877 yphbarN(intvID,idof,iphavg) = 878 & yphbar(i,idof,iphavg) 879 enddo 880 enddo 881 endif 882 ! vortN 883 if(ivort==1) then 884 do idof=1,ndofvort 885 vortN(intvID,idof) = vort(i,idof) 886 enddo 887 endif 888 ! Distance to the wall dwalN 889 if(idwal==1) then 890 dwalN(intvID) = dwal(i) 891 endif 892 endif 893 enddo 894 895 else !myrank ne to irN and irM 896 tag=irN 897 if(myrank == irN) then 898 ! rank irN receives data from rank j and update qoldN. 899! call mpi_recv(ifill, 1, MPI_INTEGER, irM, tag, 900! & mpi_comm_world, stat, ierr) 901! write(*,*) ' ', irN,'receives',ifill,getinfo(irM+1), 902! & 'data from',irM 903 ifill = getinfo(irM+1) 904 if(ifill .gt. 0) then 905 ! there are really data to receive 906! write(*,*) ' ', irN,'receives',ifill,'data from',irM 907 allocate(qrecv(ifill*numvar)) 908 qrecv = -9.87654321e32 909 call mpi_recv(qrecv, ifill*numvar, MPI_DOUBLE, irM, 910 & tag, mpi_comm_world, stat, ierr) 911 do i=1,ifill 912 intvID = qrecv( (i-1)*numvar+1 )+1.1 ! shift to start from 1 913 ! Solution qold 914 do idof=1,ndof 915 sumvar = 1+idof 916 qoldN(intvID,idof) = qrecv( (i-1)*numvar+sumvar ) 917 enddo 918 ! ybarN 919 if(iybar==1) then 920 do idof=1,ndofybar 921 sumvar = 1+ndof+idof 922 ybarN(intvID,idof) = qrecv( (i-1)*numvar+sumvar) 923 enddo 924 endif 925 ! errorN 926 if(ierror==1) then 927 do idof=1,ndoferrors 928 sumvar = 1+ndof+ndofybar+idof 929 errorsN(intvID,idof) = 930 & qrecv( (i-1)*numvar+sumvar) 931 enddo 932 endif 933 ! yphbarN 934 if(numphavg .gt. 0) then 935 do iphavg = 1,numphavg 936 do idof=1,ndofyphbar 937 sumvar = 1+ndof+ndofybar+ndoferrors+ 938 & (iphavg-1)*ndofyphbar+idof 939 yphbarN(intvID,idof,iphavg) = 940 & qrecv( (i-1)*numvar+sumvar) 941 enddo 942 enddo 943 endif 944 ! vortN 945 if(ivort==1) then 946 do idof=1,ndofvort 947 sumvar = 1+ndof+ndofybar+ndoferrors+ 948 & numphavg*ndofyphbar+idof 949 vortN(intvID,idof) = 950 & qrecv( (i-1)*numvar+sumvar) 951 enddo 952 endif 953 ! Distance to the wall dwalN 954 if(idwal==1) then 955 sumvar = 1+ndof+ndofybar+ndoferrors+ 956 & numphavg*ndofyphbar+ndofvort+1 957 dwalN(intvID) = qrecv ( (i-1)*numvar+sumvar) 958 endif 959 enddo 960 deallocate(qrecv) 961 endif 962 963 elseif(myrank == irM) then ! rank irM sends data to rank irN 964 ! First build the data to send to rank i for qoldN 965 ifill = 0 966 if (icountN(irN+1) .gt. 0) then !there are data to transfer 967 allocate(qsend(icountN(irN+1)*numvar)) 968 qsend(:) = -9.87654321e32 969 do i=1,nshg 970 ifileN = pID(i)+1.1 !shift to start from 1 971 if(ifileN == irN+1) then 972 ! these data need to be packed because rank irN is expecting them 973 ifill = ifill+1 974 ! vtx ID 975 qsend( (ifill-1)*numvar+1) = vID(i) ! no shift to start from 1 yet 976 ! Solution qold 977 do idof=1,ndof 978 sumvar = 1+idof 979 qsend( (ifill-1)*numvar+sumvar) = qold(i,idof) 980 enddo 981 ! ybar 982 if(iybar==1) then 983 do idof=1,ndofybar 984 sumvar = 1+ndof+idof 985 qsend( (ifill-1)*numvar+sumvar) = ybar(i,idof) 986 enddo 987 endif 988 ! errors 989 if(ierror==1) then 990 do idof=1,ndoferrors 991 sumvar = 1+ndof+ndofybar+idof 992 qsend( (ifill-1)*numvar+sumvar) = 993 & errors(i,idof) 994 enddo 995 endif 996 ! yphbarN 997 if(numphavg .gt. 0) then 998 do iphavg = 1,numphavg 999 do idof=1,ndofyphbar 1000 sumvar = 1+ndof+ndofybar+ndoferrors+ 1001 & (iphavg-1)*ndofyphbar+idof 1002 qsend( (ifill-1)*numvar+sumvar) = 1003 & yphbar(i,idof,iphavg) 1004 enddo 1005 enddo 1006 endif 1007 ! vortN 1008 if(ivort==1) then 1009 do idof=1,ndofvort 1010 sumvar = 1+ndof+ndofybar+ndoferrors+ 1011 & numphavg*ndofyphbar+idof 1012 qsend( (ifill-1)*numvar+sumvar) = 1013 & vort(i,idof) 1014 enddo 1015 endif 1016 ! Distance to the wall dwal 1017 if(idwal==1) then 1018 sumvar = 1+ndof+ndofybar+ndoferrors+ 1019 & numphavg*ndofyphbar+ndofvort+1 1020 qsend( (ifill-1)*numvar+sumvar) = dwal(i) 1021 endif 1022 endif 1023 enddo 1024 endif 1025 if(ifill .ne. icountN(irN+1)) then 1026 write(*,*) 'ERROR with data set:', irM, irN, 1027 & ifill, icountN(irN+1) 1028 endif 1029! write(*,*) ' ', irM,'sends', ifill,'data to',irN 1030! call mpi_send(ifill, 1 ,MPI_INTEGER, irN, tag, 1031! & mpi_comm_world, ierr) ! Send the size of the data 1032 if(ifill .gt. 0) then !if the size of the data is >0, send it for good 1033! write(*,*) ' ', irM,'sends', ifill,'data to',irN 1034 call mpi_send(qsend(1), ifill*numvar, MPI_DOUBLE, irN, 1035 & tag, mpi_comm_world, ierr) 1036 deallocate(qsend) 1037 endif 1038 1039 endif !if(myrank==irN) 1040 endif !if (myrank == irN = irM) 1041 1042 enddo !loop over irM 1043 1044 call mpi_barrier(mpi_comm_world,ierr) 1045 if (myrank == irN) then 1046 nonzero = 0 1047 do i = 1,numpe 1048 if (getinfo(i) .gt. 0) then 1049 nonzero = nonzero +1 1050 endif 1051 enddo 1052 write(*,*)"Part ",irN,"out of N = ",irankN," updated from", 1053 & nonzero, "parts" 1054 deallocate(getinfo) 1055 endif 1056 enddo !loop over irN 1057 1058! 1059! Do it without files through mpialltoall 1060! 1061! maxicountN=maxval(icountN(:)) 1062! call MPI_ALLREDUCE (maxicountN, maxicountNglob, 1, 1063! & MPI_INT, MPI_MAX, mpi_comm_world, ierr) 1064! 1065! numvar = 1+ndof ! vID + ndof from solution 1066! write(*,*) 'maxicountN:',myrank, maxicountN, maxicountNglob, nshg 1067! call mpi_barrier(mpi_comm_world, ierr) 1068! allocate(qTota2as(numvar*maxicountNglob,numpe)) 1069! qTota2as(:,:) = -9.87654321e32 1070 1071! allocate(indexpart(numpe)); indexpart(:)=0 1072! do i=1,nshg 1073! ifileN = pID(i)+1.1 !shift pID to start from 1 1074! indexpart(ifileN) = indexpart(ifileN)+1 1075! qTota2as( (indexpart(ifileN)-1)*numvar+1 , ifileN ) = vID(i) ! vID not yet shifted to start from 1 here 1076! do j=1,ndof 1077! 1078! if((indexpart(ifileN)-1)*numvar+1+j>numvar*maxicountNglob) then 1079! write(*,*) 'ERROR DIMENSION1:',myrank, 1080! & ifileN, indexpart(ifileN), 1081! & (indexpart(ifileN)-1)*numvar+1+j, numvar*maxicountNglob 1082! endif 1083! if (ifileN > numpe ) then 1084! write(*,*) 'ERROR DIMENSION2:',myrank, ifileN, numpe 1085! endif 1086! qTota2as( (indexpart(ifileN)-1)*numvar+1+j, ifileN )=qold(i,j) 1087! enddo 1088! enddo 1089 1090! allocate(qTota2ar(numvar*maxicountNglob,numpe)) 1091! qTota2ar(:,:) = -9.87654321e32 1092 1093! ncount = numvar*maxicountNglob 1094! call mpi_alltoall(qTota2as(1,1), ncount, MPI_DOUBLE, 1095! & qTota2ar(1,1), ncount, MPI_DOUBLE, 1096! & mpi_comm_world, ierr) 1097! if (ierr .ne. 0) then 1098! write(*,*) 'Error with mpi_alltoall:', myrank, ierr 1099! endif 1100! deallocate(qTota2as) 1101 1102 1103 elseif (ireducemethod == 1) then ! file based approach 1104! 1105! Now for each process on M parts push the fields into an i,j pair named 1106! ascii file restartMap.i.j 1<=i<=N 1<=j<=M. This is embarassingly 1107! parallel with no data collision but writes a lot of files 1108! 1109 do l=1,irankN 1110 if(icountN(l).gt.0) then ! I have data to write for rank l 1111! 1112! generate the filename since this is a file we have to write 1113! note this avoids opening files we don't need but each rank does open 1114! all the files it will need to write to....lots of files=weakness 1115! 1116! 1117! write the number of vtx being mapped from part j to part i 1118! note ndof will need to change to the total size of vector field being mapped 1119! and this also should be computed in some routine above for allocation 1120! purposes. currently it is at risk of being confused with ndof 1121! 1122 itmp=1 1123 if (l .gt. 0) itmp = int(log10(float(l)))+1 1124 write (fmt1,"('(''restartMap.'',i',i1,',1x)')") itmp 1125 write (fnamer,fmt1) l 1126 fnamer = trim(fnamer) // cname2(myrank+1) 1127 if (isbinary == 1) then 1128 open(unit=700+l,file=fnamer,status='unknown', 1129 & form='unformatted') 1130 else 1131 open(unit=700+l,file=fnamer,status='unknown') 1132 endif 1133 endif 1134 enddo 1135! 1136! so that we can read directly into the mapped array later, it is helpful to 1137! find the number of verts on each part in the N partition which is easily 1138! obtained by looking at the max of each vtx id going to each N rank 1139! and then doing an allReduce to find global maximum 1140! 1141! Write solution 1142! 1143 do l=1,irankN 1144 if(icountN(l).gt.0) then ! I have data to write for rank l 1145 if (isbinary == 1) then 1146 write(700+l) icountN(l), ndof 1147 else 1148 write(700+l,*) icountN(l), ndof 1149 endif 1150 endif 1151 enddo 1152 1153 do i=1,nshg 1154 ifileN=pID(i)+1.1 !shift pID to start from 1 1155 intvID=(vID(i)+1.1) !shift vID to start from 1 1156 if (isbinary == 1) then 1157 write(700+ifileN) intvID,(qold(i,j),j=1,ndof) 1158 else 1159 write(700+ifileN,925) intvID,(qold(i,j),j=1,ndof) 1160 endif 1161! the above "dealing" of records to files based on pID is safe because each 1162! process 1..M has a unique set of files to write to 1163 enddo 1164 1165! 1166! Write dwal 1167! 1168 ione = 1 1169 do l=1,irankN 1170 if(icountN(l).gt.0) then ! I have data to write for rank l 1171 if (isbinary == 1) then 1172 write(700+l) icountN(l), ione 1173 else 1174 write(700+l,*) icountN(l), ione 1175 endif 1176 endif 1177 enddo 1178 1179 do i=1,nshg 1180 ifileN=pID(i)+1.1 !shift pID to start from 1 1181 intvID=(vID(i)+1.1) !shift vID to start from 1 1182 if (isbinary == 1) then 1183 write(700+ifileN) intvID, dwal(i) 1184 else 1185 write(700+ifileN,925) intvID, dwal(i) 1186 endif 1187! the above "dealing" of records to files based on pID is safe because each 1188! process 1..M has a unique set of files to write to 1189 enddo 1190! 1191! close all the files on each of the M ranks and deallocate as we are done 1192! with the ascii file creation step 1193! 1194 do l=1,irankN 1195 if(icountN(l).gt.0) close(700+l) 1196 enddo 1197 1198 endif ! ireducemethod 1199 1200! 1201! Deallocate some memory 1202! 1203 call mpi_barrier(mpi_comm_world, ierr) 1204 if(myrank == 0) then 1205 write(*,*) 'Deallocating some memory' 1206 endif 1207 deallocate(vID) 1208 deallocate(pID) 1209 deallocate(qold) 1210 if(iybar == 1) then 1211 deallocate(ybar) 1212 endif 1213 if(ierror == 1) then 1214 deallocate(errors) 1215 endif 1216 if(numphavg .gt. 0) then 1217 deallocate(yphbar) 1218 endif 1219 if(ivort == 1) then 1220 deallocate(vort) 1221 endif 1222 if(idwal == 1) then 1223 deallocate(dwal) 1224 endif 1225! 1226! Make sure all the ranks are done with their restartMap files before 1227! trying to read them for ireducemethod = 1 1228! 1229 call mpi_barrier(mpi_comm_world, ierr) 1230 1231! 1232! Create a subcommunicator for the first N ranks 1233! 1234 1235! islesseqN = 0 1236! if(myrank.lt.irankN) then 1237! islesseqN = 1 1238! endif 1239! call mpi_comm_split(mpi_comm_world, islesseqN, myrank, 1240! & my_local_comm, ierr) 1241 1242! my_local_size = -1 1243! call mpi_comm_size(my_local_comm, my_local_size, ierr) 1244! my_local_rank = -1 1245! call mpi_comm_rank(my_local_comm, my_local_rank, ierr) 1246! write(*,*) "RANKS:", myrank, my_local_rank, my_local_size, irankN 1247 1248! 1249! Now for all N ranks, load the partial files written above and merge them 1250! 1251! if(myrank.lt.irankN) then ! only parts 1..N do work in this phase 1252 1253 if(ireducemethod == 1) then ! Read the files back for the file based approach. 1254 ! Beware, this method is not maintained any more 1255 1256 l=myrank+1 ! help me find my mapped files 1257 do j=1,numpe ! potentially all M have a map for me 1258! 1259! create the filename for all M of my potential Maps 1260! 1261 itmp=1 1262 if (l .gt. 0) itmp = int(log10(float(l)))+1 1263 write (fmt1,"('(''restartMap.'',i',i1,',1x)')") itmp 1264 write (fnamer,fmt1) l 1265 fnamer = trim(fnamer) // cname2(j) 1266! 1267! most won't exist so inquire for existance 1268! 1269 inquire(file=fnamer,exist=exmap) 1270 if(exmap) then ! so this map does exist so we open the Map files 1271 1272 if(isbinary == 1) then 1273 open(unit=700+l,file=fnamer,status='unknown', 1274 & form='unformatted') 1275 else 1276 open(unit=700+l,file=fnamer,status='unknown') 1277 endif 1278! 1279! Read solution first 1280! 1281 if(isbinary == 1) then 1282 read(700+l) icountRead, ndofRead ! read it's header 1283 else 1284 read(700+l,*) icountRead, ndofRead ! read it's header 1285 endif 1286 allocate(qreadN(ndofRead)) 1287 do i=1, icountRead ! loop over its contents 1288 if(isbinary == 1) then 1289 read(700+l) intvID,(qreadN(k),k=1,ndofRead) ! read a line in binary format 1290 else 1291 read(700+l,925) intvID,(qreadN(k),k=1,ndofRead) ! read a line in ascii format 1292 endif 1293! DEBUG 1294! if((i.eq.1).and.(myrank.eq.0)) 1295! & write(*,*) vidRead, ndofRead, (qreadN(k),k=1,ndofRead) 1296! DEBUG 1297 qoldN(intvID,1:ndof)=qreadN(1:ndof) ! fill directly into qoldN 1298! 1299! the above avoids sorting and sifting as duplicated vtx values just overwrite 1300! and you get the last one as a final value. Assuming the field was commu-ed 1301! all duplicat values are the same so this is safe. For error fields this 1302! might not be true but it should be o.k. 1303! 1304! DEBUG 1305! if((i.eq.1).and.(myrank.eq.0)) 1306! & write(*,*) intvID,(qoldN(intvID,k),k=1,ndofRead) 1307! DEBUG 1308 enddo ! end of this map file so close the file and look for more 1309 deallocate(qreadN) 1310! 1311! Read dwalN second 1312! 1313 if(isbinary == 1) then 1314 read(700+l) icountRead, ndofRead ! read it's header 1315 else 1316 read(700+l,*) icountRead, ndofRead ! read it's header 1317 endif 1318 do i=1, icountRead 1319 if(isbinary == 1) then 1320 read(700+l) intvID,qreadN1 1321 else 1322 read(700+l,925) intvID,qreadN1 1323 endif 1324 dwalN(intvID)=qreadN1 ! fill directly into dwalN 1325 enddo ! end of this map file so close the file and look for more 1326 1327 close(700+l) 1328 endif ! map existed 1329 enddo ! loop over all M ranks 1330 1331 endif ! ireducemethod = 1 1332! 1333! At this point I should have a complete field for my part so I write the 1334! mapped file 1335! 1336 1337! DEBUG 1338! itmp=1 1339! if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1 1340! write (fmt1,"('(''restartMapped.'',i',i1,',1x)')") itmp 1341! write (fnamer,fmt1) lstep 1342! fnamer = trim(fnamer) // cname2(myrank+1) 1343! open(unit=700,file=fnamer,status='unknown') !, 1344! write(700,*) maxvIDout(myrank+1),ndof 1345! do i=1,maxVIDout(myrank+1)+1 1346! write(700,927) i, (qoldN(i,j),j=1,ndof) 1347! enddo 1348! close(700) 1349! DEBUG 1350 1351! 1352! Write the N posix files. Note that the nodes on the part boundaries still need to be updated 1353! 1354! Posix format 1355! call write_restartonly(myrank, lstep, nshgN, ndof, qoldN) 1356! call write_field (myrank,'a','dwal',4,dwalN,'d',nshgN,1,lstep) 1357! call write_field (myrank,'a','errors',6,errorsN,'d', 1358! & nshgN,ndoferrors,lstep) 1359! call write_field (myrank,'a','ybar',4,ybarN,'d', 1360! & nshgN,ndofybar,lstep) 1361 1362! SyncIO format 1363 1364 1365! 1366! Test the result from mpi_alltoall 1367! 1368! qoldN=-9.87654321e32 1369! l=myrank+1 1370! do j=1,numpe ! should really go to irankN 1371! do i=1,maxicountNglob 1372! intvID = qTota2ar( (i-1)*numvar+1 , j) + 1.1 1373! if (intvID > 0) then 1374! do idof=1,ndof 1375! qoldN(intvID,idof) = qTota2ar( (i-1)*numvar+1+idof, j) 1376! enddo 1377! endif 1378! enddo 1379! enddo 1380 1381! itmp=1 1382! if (lstep .gt. 0) itmp = int(log10(float(lstep)))+1 1383! write (fmt1,"('(''restartMapped2.'',i',i1,',1x)')") itmp 1384! write (fnamer,fmt1) lstep 1385! fnamer = trim(fnamer) // cname2(myrank+1) 1386! open(unit=700,file=fnamer,status='unknown') !, 1387! write(700,*) maxvIDout(myrank+1),ndof 1388! do i=1,maxVIDout(myrank+1)+1 1389! write(700,927) i, (qoldN(i,j),j=1,ndof) 1390! enddo 1391! close(700) 1392 1393! endif ! my rank is within range 1..N 1394 1395 call mpi_barrier(mpi_comm_world, ierr) 1396 if(myrank == 0) then 1397 write(*,*) 'Writing the reduced restartRedTmp-dat files' 1398 endif 1399 1400! call Write_M2N(myrank, irankN, lstep, nshgN, 1401! & ndof, ndofybar, ndoferrors, 1402! & qoldN, ybarN, errorsN, dwalN) 1403 1404 1405 nsynciofieldswriterestart = 1+iybar+ierror+numphavg+ivort+idwal 1406 1407! We check here how many geombcRed files have been red. 1408! If 0, then two possibilities: 1409! - Reduction to 1 part so geombcRed are not required. Set nsynciofilesred to 1 1410! - User forgot to create the link to geombcRed. 1411! Then, assumes the number of geombc and geombcRed files ia the same 1412 1413 if(nsynciofilesred == 0) then 1414 if(irankN == 1) then ! 1 part reduction -> set to 1 1415 nsynciofilesred = 1 1416 else ! Links are missing 1417 nsynciofilesred = nsynciofiles 1418 endif 1419 endif 1420 1421 call Write_M2N_SolOnly(myrank,irankN,lstep,nshgN,ndof,qoldN) 1422 deallocate(qoldN) 1423 1424 if(iybar == 1) then 1425 call Write_M2N_Field(myrank,irankN,'a','ybar',4,ybarN,'d', 1426 & nshgN,ndofybar,lstep) 1427 deallocate(ybarN) 1428 endif 1429 1430 if(ierror == 1) then 1431 call Write_M2N_Field(myrank,irankN,'a','errors',6,errorsN,'d', 1432 & nshgN,ndoferrors,lstep) 1433 deallocate(errorsN) 1434 endif 1435 1436 if(numphavg .gt. 0) then 1437 do iphavg=1,numphavg 1438 call write_M2N_phavg2(myrank,irankN,'a','phase_average',13, 1439 & iphavg,numphavg,yphbarN(:,:,iphavg),'d', 1440 & nshgN,ndofyphbar,lstep) 1441 enddo 1442 deallocate(yphbarN) 1443 endif 1444 1445 if(ivort == 1) then 1446 call Write_M2N_Field(myrank,irankN,'a','vorticity',9, 1447 & vortN,'d',nshgN,ndofvort,lstep) 1448 deallocate(vortN) 1449 endif 1450 1451 if(idwal == 1) then 1452 call Write_M2N_Field(myrank,irankN,'a','dwal',4,dwalN,'d', 1453 & nshgN,1,lstep) 1454 deallocate(dwalN) 1455 endif 1456 1457 1458! 1459! Below is debug stuff trying to figure out why the code crashes on exit 1460! 1461 1462! call mpi_barrier(mpi_comm_world, ierr) 1463! write(*,*) myrank,'leaving readnblk' 1464 1465!925 format(i5,2x,6(e14.7,2x)) 1466925 format(i9,1x,6(e24.17,1x)) ! used for file based approach 1467 1468926 format(9(i5,2x)) !used for debugging 1469927 format(i5,2x,6(e14.7,2x)) !used for debugging 1470 1471 return 1472! 1473! 1474 end 1475 1476