mpi_amr_checkpoint_mpiio.f90

做网格的好程序

        icheckp_on_ec(:,i) = 0
        if(checkp_on_ec(1,i)) then
           nvar_chk_ec = nvar_chk_ec + 1
           icheckp_on_ec(:,i) = 1
        end if
      enddo
      nvar_chk_nc =  0
      do i=1,nvarcorn
        icheckp_on_nc(i) = 0
        if(checkp_on_nc(i)) then
           nvar_chk_nc = nvar_chk_nc + 1
           icheckp_on_nc(i) = 1
        end if
      enddo

      ngid = nfaces + 1 + nchild

      do i = 1,4
         iorder(i) = size(ordering,dim=i)
      end do

#ifdef REAL8

      call read_blocks_mpiio_r8(filename, & 
     &                         tot_blocks, & 
     &                         lnblocks, & 
     &                         mdim, & 
     &                         ndim, & 
     &                         ngid, & 
     &                         mflags, & 
     &                         lrefine, & 
     &                         nodetype, & 
     &                         which_child, & 
     &                         gid, & 
     &                         bflags, & 
     &                         coord, & 
     &                         bnd_box, & 
     &                         work_block, & 
     &                         unk, & 
     &                         nvar, & 
     &                         nvar_chk_cc, & 
     &                         icheckp_on_cc, & 
     &                         facevarx, facevary, facevarz, & 
     &                         nbndvar, & 
     &                         nvar_chk_fc, & 
     &                         icheckp_on_fc, & 
     &                         unk_e_x, unk_e_y, unk_e_z, & 
     &                         nbndvare, & 
     &                         nvar_chk_ec, & 
     &                         icheckp_on_ec, & 
     &                         unk_n, & 
     &                         nbndvarc, & 
     &                         nvar_chk_nc, & 
     &                         icheckp_on_nc, & 
     &                         iu_bnd, ju_bnd, ku_bnd, & 
     &                         il0, iu0, & 
     &                         jl0, ju0, & 
     &                         kl0, ku0, & 
     &                         user_attr_1_value, & 
     &                         user_attr_2_value, & 
     &                         user_attr_3_value, & 
     &                         user_attr_4_value, & 
     &                         user_attr_5_value, &
     &                         iorder)

#else

      call read_blocks_mpiio_r4(filename, & 
     &                         tot_blocks, & 
     &                         lnblocks, & 
     &                         mdim, & 
     &                         ndim, & 
     &                         ngid, & 
     &                         mflags, & 
     &                         lrefine, & 
     &                         nodetype, & 
     &                         which_child, & 
     &                         gid, & 
     &                         bflags, & 
     &                         coord, & 
     &                         bnd_box, & 
     &                         work_block, & 
     &                         unk, & 
     &                         nvar, & 
     &                         nvar_chk_cc, & 
     &                         icheckp_on_cc, & 
     &                         facevarx, facevary, facevarz, & 
     &                         nbndvar, & 
     &                         nvar_chk_fc, & 
     &                         icheckp_on_fc, & 
     &                         unk_e_x, unk_e_y, unk_e_z, & 
     &                         nbndvare, & 
     &                         nvar_chk_ec, & 
     &                         icheckp_on_ec, & 
     &                         unk_n, & 
     &                         nbndvarc, & 
     &                         nvar_chk_nc, & 
     &                         icheckp_on_nc, & 
     &                         iu_bnd, ju_bnd, ku_bnd, & 
     &                         il0, iu0, & 
     &                         jl0, ju0, & 
     &                         kl0, ku0, & 
     &                         user_attr_1_value, & 
     &                         user_attr_2_value, & 
     &                         user_attr_3_value, & 
     &                         user_attr_4_value, & 
     &                         user_attr_5_value, &
     &                         iorder)

#endif

      if(ndim.lt.3) then
         bnd_box(1,3,:) = 0.
         bnd_box(2,3,:) = 1.
         coord(3,:) = .5*(bnd_box(2,3,:)+bnd_box(1,3,:))
      endif
      if(ndim.lt.2) then
        bnd_box(1,2,:) = 0.
        bnd_box(2,2,:) = 1.
        coord(2,:) = .5*(bnd_box(2,2,:)+bnd_box(1,2,:))
      endif

      do block_no = 1,lnblocks
         bsize(:,block_no) = bnd_box(2,:,block_no)- & 
     &        bnd_box(1,:,block_no)
      enddo

! COMPUTE TREE DATA FROM gid

      if(allocated(n_to_left)) deallocate( n_to_left )
      allocate ( n_to_left(0:nprocs-1) )

      proc = 0
      icount = 0
      n_to_left(:) = 0
      do while (icount < tot_blocks)
         if (proc < mype) then
            n_to_left(mype) = n_to_left(mype) + 1
         end if
         proc = proc + 1
         if (proc > nprocs-1) proc = 0
         icount = icount + 1
      end do

      call MPI_Allgather(n_to_left(mype), 1,MPI_INTEGER, & 
     &                   n_to_left,1,MPI_INTEGER, & 
     &                   MPI_COMM_WORLD,ierr)

      do block_no = 1,lnblocks

! neighbor data
         ngid = 0
         do j = 1,nfaces
            ngid = ngid + 1
            if (gid(ngid,block_no).gt.0) then
               proc = 0
               do while (n_to_left(proc) < gid(ngid,block_no))
                  proc = proc + 1
                  if (proc >= nprocs) exit
               end do
               proc = proc - 1
               neigh(2,j,block_no) = proc 
               neigh(1,j,block_no) = gid(ngid,block_no) - & 
     &                               n_to_left(proc)
            else
               neigh(1,j,block_no) = gid(ngid,block_no)
               neigh(2,j,block_no) = gid(ngid,block_no)
            end if
         end do


! parent data
         ngid = ngid + 1
         if (gid(ngid,block_no).gt.0) then
            proc = 0
            do while (n_to_left(proc) < gid(ngid,block_no))
               proc = proc + 1
               if (proc >= nprocs) exit
            end do
            proc = proc - 1
            parent(2,block_no) = proc
            parent(1,block_no) = gid(ngid,block_no) - & 
     &                           n_to_left(proc)
         else
            parent(1,block_no) = gid(ngid,block_no)
            parent(2,block_no) = gid(ngid,block_no)
         end if

! children data
         do j = 1,nchild
            ngid = ngid + 1
            if (gid(ngid,block_no).gt.0) then
               proc = 0
               do while (n_to_left(proc) < gid(ngid,block_no))
                  proc = proc + 1
                  if (proc >= nprocs) exit
               end do
               proc = proc - 1
               child(2,j,block_no) = proc
               child(1,j,block_no) = gid(ngid,block_no) - & 
     &                               n_to_left(proc)
            else
               child(1,j,block_no) = gid(ngid,block_no)
               child(2,j,block_no) = gid(ngid,block_no)
            end if
         end do

      end do

! Now reorder blocks such that they are better balanced
! NOTE: this assumes that the total number of blocks is > nprocs

! NOTE: We cannot do a morton ordering here if l_with_guardcells is defined since
!       amr_redist_blks do not move the guardcells.

      if (.not.l_with_guardcells2) then  

      lnblocks_old = lnblocks
      l_move_solution = .true.
      call amr_morton_order (lnblocks_old,nprocs,mype, & 
     &                       l_move_solution)

      write(*,*) 'after amr_morton_order : pe lnb ',mype,lnblocks

      end if

!---------------------------------------------
! compute grid_xmax, etc
      if (timing_mpi) then
         time1 = mpi_wtime()
      endif
      call mpi_amr_global_domain_limits
      if (timing_mpi) then
      timer_amr_global_domain_limits = & 
     &     timer_amr_global_domain_limits + mpi_wtime() - time1
      endif

#ifdef SAVE_MORTS
! Compute xmin,ymin,zmin,xmax,ymax,zmax or get them from storage
      xmin = grid_xmin
      ymin = grid_ymin
      zmin = grid_zmin
      xmax = grid_xmax
      ymax = grid_ymax
      zmax = grid_zmax

      write(*,*) 'checkre xmin etc ',xmin,ymin,zmin,xmax,ymax,zmax
      write(*,*) 'checkre lperiodicxtc ', & 
     &                        lperiodicx,lperiodicy,lperiodicz
      do lb = 1,lnblocks
        call morton_neighbors(xmin,ymin,zmin,xmax,ymax,zmax, & 
     &                        lperiodicx,lperiodicy,lperiodicz, & 
     &                        coord(:,lb),bsize(:,lb),ndim, & 
     &                        lrefine(lb),lrefine_max,mort_neigh)
        surr_morts(:,:,:,:,lb) = & 
     &                    mort_neigh(:,:,2-k2d:2+k2d,2-k3d:2+k3d)
      write(*,*) 'check_re set up surr_morts for initial blk ',lb,mype, & 
     &            ' surr_morts(6,:,:,1,lb) ',surr_morts(6,:,:,1,lb)
      enddo
#endif /* SAVE_MORTS */

      call amr_morton_process()

!
! St up an array of cell sizes for each grid refinement level.
! These can be used to minimize variation due to roundoff, but
! should ONLY be used with a uniformly spaced grid.
      level_cell_sizes = 0.
      level_cell_sizes(1,1) = (grid_xmax-grid_xmin)/real(nxb)
      if(ndim.gt.1) & 
     &  level_cell_sizes(2,1) = (grid_ymax-grid_ymin)/real(nyb)
      if(ndim.eq.3) & 
     &  level_cell_sizes(3,1) = (grid_zmax-grid_zmin)/real(nzb)
      do i=2,lrefine_max
        level_cell_sizes(1:ndim,i) = .5*level_cell_sizes(1:ndim,i-1)
      enddo
!---------------------------------------------

!
! mark grid as changed
      grid_changed = 1
      grid_analysed_mpi = 1


! Now make sure guardcell information is up to date

      if (.not.l_with_guardcells2)  & 
     &     call amr_guardcell(mype,1,nguard)

      call mpi_amr_boundary_block_info(mype,nprocs)

      deallocate(n_to_left)

      deallocate(icheckp_on_cc)
      deallocate(icheckp_on_fc)
      deallocate(icheckp_on_ec)
      deallocate(icheckp_on_nc)

! handle user attributes

      if (present(user_attr_1)) then
         user_attr_1 = user_attr_1_value
      end if

      if (present(user_attr_2)) then
         user_attr_2 = user_attr_2_value
      end if
         
      if (present(user_attr_3)) then
         user_attr_3 = user_attr_3_value
      end if
         
      if (present(user_attr_4)) then
         user_attr_4 = user_attr_4_value
      end if
         
      if (present(user_attr_5)) then
         user_attr_5 = user_attr_5_value
      end if

      return
      end subroutine amr_checkpoint_re_mpiio