fix_conserve – AGILE

public subroutine fix_conserve(bgstep, idimmin, idimmax, nw0, nwfluxin)

Uses

- mod_comm_lib
- mod_global_parameters
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.
Correct the coarse cells abutting a refinement boundary: take out the coarse face flux the update used and put back the sum of the fine ones.
The stored flux is Cartesian qdtf/dx or curvilinear qdtfA, see the flux-fixing stores in mod_finite_volume. The Cartesian form is already a per-volume quantity but of the fine* cell, so it is rescaled by CoFiratio; the curvilinear form is extensive and is divided here by the coarse cell's own volume, with no ratio, because the four fine face areas already sum to the coarse face's.
Arguments

Type	Intent		Name
integer,	intent(in)	::	bgstep
integer,	intent(in)	::	idimmin
integer,	intent(in)	::	idimmax
integer,	intent(in)	::	nw0
integer,	intent(in)	::	nwfluxin
Calls

Help
Called by

Help
Source Code

   subroutine fix_conserve(bgstep,idimmin,idimmax,nw0,nwfluxin)
     use mod_global_parameters
     use mod_comm_lib, only: mpistop

     integer, intent(in) :: bgstep, idimmin,idimmax, nw0, nwfluxin

     integer :: iigrid, igrid, idims, iside, iotherside, i1,i2,i3, ic1,ic2,ic3,&
         inc1,inc2,inc3, ixmin1,ixmin2,ixmin3,ixmax1,ixmax2,ixmax3
     integer :: ix1, ix2, ix3 !JESSE ADDED
     integer :: nxCo1,nxCo2,nxCo3, iw, ix, ipe_neighbor, ineighbor, nbuf,&
         ibufnext, nw1
     double precision :: CoFiratio

     nw1=nw0-1+nwfluxin
     ! The flux is divided by volume of fine cell. We need, however,
     ! to divide by volume of coarse cell => muliply by volume ratio
     CoFiratio=one/dble(2**ndim)

     if (n_recv_pe>0) then
       call MPI_WAITALL(n_recv_pe,fc_recvreq,fc_recvstat,ierrmpi)
       ! With one aggregated message per peer, a disagreement about which
       ! chunks the run holds no longer shows up as an MPI mismatch, so check
       ! the arrival lengths. This catches any difference in the *set* of
       ! chunks the two ranks enumerated; a difference in their *order* cannot
       ! arise, since both sort the same run by the same key.
       do irecv=1,n_recv_pe
         call MPI_GET_COUNT(fc_recvstat(:,irecv),MPI_DOUBLE_PRECISION,nbuf,&
            ierrmpi)
         if (nbuf /= recv_pe_len(irecv)) call mpistop(&
            "fix_conserve: flux message length disagrees with the layout")
       end do
#ifdef NOGPUDIRECT
       ! Without GPU-direct the IRECVs landed in host memory; the unpack below
       ! runs on the device, so push the payload across.  Sliced, not
       ! whole-array: the buffer is grown and not shrunk.
       
#endif
     end if

     nxCo1=(ixMhi1-ixMlo1+1)/2
     nxCo2=(ixMhi2-ixMlo2+1)/2
     nxCo3=(ixMhi3-ixMlo3+1)/2

     ! for all grids: perform flux update at Coarse-Fine interfaces
     ! Everything assigned per gang is named here rather than left to the
     ! implicit firstprivate OpenACC gives an unlisted scalar in a parallel
     ! region: the rule does cover them, but a half-populated list is
     ! indistinguishable from an oversight.
      
     do iigrid=1,igridstail
       igrid=igrids(iigrid)

       do idims=idimmin,idimmax
         select case (idims)
           case (1)
           do iside=1,2
             i1=kr(1,1)*(2*iside-3)
             i2=kr(2,1)*(2*iside-3)
             i3=kr(3,1)*(2*iside-3)

             if (neighbor_pole(i1,i2,i3,igrid)/=0) cycle

             if (neighbor_type(i1,i2,i3,igrid)/=4) cycle

 !opedit: skip over active/passive interface since flux for passive ones is
             ! not computed, keep the buffer counter up to date:
          !   if (.not.neighbor_active(i1,i2,i3,&
          !      igrid).or..not.neighbor_active(0,0,0,igrid) ) then
          !     do ic3=1+int((1-i3)/2),2-int((1+i3)/2)
          !     inc3=2*i3+ic3
          ! do ic2=1+int((1-i2)/2),2-int((1+i2)/2)
          !     inc2=2*i2+ic2
          ! do ic1=1+int((1-i1)/2),2-int((1+i1)/2)
          !     inc1=2*i1+ic1
          !     ipe_neighbor=neighbor_child(2,inc1,inc2,inc3,igrid)
          !     if (ipe_neighbor/=mype) then
          !       ibufnext=ibuf+isize(1)
          !       ibuf=ibufnext
          !     end if
          !     end do
          ! end do
          ! end do
          !     cycle
          !   end if
             !

             select case (iside)
             case (1)
               ix=ixMlo1
             case (2)
               ix=ixMhi1
             end select

             ! remove coarse flux
                
                do ix3=ixMlo3,ixMhi3
                  do ix2=ixMlo2,ixMhi2 
                    do iw=1,nwfluxin
                      bg(bgstep)%w(ix,ix2,ix3,nw0+iw-1,igrid) = &
                        bg(bgstep)%w(ix,ix2,ix3,nw0+iw-1,igrid) - &
                        pflux(iside,1)%flux(1,ix2-nghostcells,ix3-nghostcells,&
                                            iw,igrid)
                    end do
                  end do
                end do


             ! add fine flux

             do ic3=1+int((1-i3)/2),2-int((1+i3)/2)
               inc3=2*i3+ic3
             do ic2=1+int((1-i2)/2),2-int((1+i2)/2)
               inc2=2*i2+ic2
             do ic1=1+int((1-i1)/2),2-int((1+i1)/2)
               inc1=2*i1+ic1
               ineighbor=neighbor_child(1,inc1,inc2,inc3,igrid)
               ipe_neighbor=neighbor_child(2,inc1,inc2,inc3,igrid)
               ixmin1=ix
               ixmin2=ixMlo2+(ic2-1)*nxCo2
               ixmin3=ixMlo3+(ic3-1)*nxCo3
               ixmax1=ix
               ixmax2=ixmin2-1+nxCo2
               ixmax3=ixmin3-1+nxCo3
               if (ipe_neighbor==mype) then
                 iotherside=3-iside
                     
                     do ix3=1,nxCo3 
                        do ix2=1,nxCo2 
                          do iw=1,nwfluxin
                             bg(bgstep)%w(ix,ixmin2+ix2-1,ixmin3+ix3-1,nw0+iw-1,igrid) = &
                              bg(bgstep)%w(ix,ixmin2+ix2-1,ixmin3+ix3-1,nw0+iw-1,igrid) + &
                              pflux(iotherside,1)%flux(1,ix2,ix3,iw,&
                                ineighbor) * CoFiratio
                          end do
                        end do
                     end do
               !else
               !  if (slab_uniform) then
               !    ibufnext=ibuf+isize(1)
               !    if(stagger_grid) ibufnext=ibufnext-isize_stg(1)
               !    psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,ixmin3:ixmax3,&
               !       nw0:nw1) = psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !       ixmin3:ixmax3,nw0:nw1)+CoFiratio &
               !       *reshape(source=recvbuffer(ibuf:ibufnext-1),&
               !        shape=shape(psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !       ixmin3:ixmax3,nw0:nw1)))
               !    ibuf=ibuf+isize(1)
               !  else
               !    ibufnext=ibuf+isize(1)
               !    if(stagger_grid) then
               !      nbuf=(isize(1)-isize_stg(1))/nwfluxin
               !    else
               !      nbuf=isize(1)/nwfluxin
               !    end if
               !    do iw=nw0,nw1
               !      psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,ixmin3:ixmax3,&
               !         iw)=psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !         ixmin3:ixmax3,iw) &
               !         +reshape(source=recvbuffer(ibuf:ibufnext-1),&
               !          shape=shape(psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !         ixmin3:ixmax3,iw))) /ps(igrid)%dvolume(ixmin1:ixmax1,&
               !         ixmin2:ixmax2,ixmin3:ixmax3)
               !      ibuf=ibuf+nbuf
               !    end do
               !    ibuf=ibufnext
               !  end if
               else
                   ! Two transverse indices plus the variable index: every iteration
                   ! lands in a distinct cell of a distinct variable, and the buffer
                   ! offset is a pure function of the three, so all three collapse.
                   
                   do ix3=1,nxCo_fc(3,1)
                     do ix2=1,nxCo_fc(2,1)
                       do iw=1,nwfluxin
                         bg(bgstep)%w(ix,ixmin2+ix2-1,ixmin3+ix3-1,nw0+iw-1,igrid) = &
                           bg(bgstep)%w(ix,ixmin2+ix2-1,ixmin3+ix3-1,nw0+iw-1,igrid) + &
                           recvbuffer(ibuf_offset(4**3*(igrid-1)+inc1+4*inc2+16*inc3+1) &
                              +(ix2-1)+(ix3-1)*nxCo_fc(2,1)+(iw-1)*nxCo_fc(2,1)*nxCo_fc(3,1)) * CoFiratio
                       end do
                     end do
                   end do
               end if
            end do
           end do
           end do
           end do
           case (2)
           do iside=1,2
             i1=kr(1,2)*(2*iside-3)
             i2=kr(2,2)*(2*iside-3)
             i3=kr(3,2)*(2*iside-3)

             if (neighbor_pole(i1,i2,i3,igrid)/=0) cycle

             if (neighbor_type(i1,i2,i3,igrid)/=4) cycle

 !opedit: skip over active/passive interface since flux for passive ones is
             ! not computed, keep the buffer counter up to date:
            ! if (.not.neighbor_active(i1,i2,i3,&
            !    igrid).or..not.neighbor_active(0,0,0,igrid) ) then
            !   do ic3=1+int((1-i3)/2),2-int((1+i3)/2)
            !   inc3=2*i3+ic3
            !   do ic2=1+int((1-i2)/2),2-int((1+i2)/2)
            !   inc2=2*i2+ic2
            !   do ic1=1+int((1-i1)/2),2-int((1+i1)/2)
            !   inc1=2*i1+ic1
            !   ipe_neighbor=neighbor_child(2,inc1,inc2,inc3,igrid)
            !   if (ipe_neighbor/=mype) then
            !     ibufnext=ibuf+isize(2)
            !     ibuf=ibufnext
            !   end if
            !   end do
            !  end do
            !  end do
            !   cycle
            ! end if
             !

             select case (iside)
             case (1)
               ix=ixMlo2
             case (2)
               ix=ixMhi2
             end select

             ! remove coarse flux
                
                do ix3=ixMlo3,ixMhi3
                  do ix1=ixMlo1,ixMhi1 
                    do iw=1,nwfluxin
                      bg(bgstep)%w(ix1,ix,ix3,nw0+iw-1,igrid) = &
                       bg(bgstep)%w(ix1,ix,ix3,nw0+iw-1,igrid) - &
                       pflux(iside,2)%flux(ix1-nghostcells,1,ix3-nghostcells,&
                                    iw,igrid)
                    end do
                  end do
                end do


             ! add fine flux
             do ic3=1+int((1-i3)/2),2-int((1+i3)/2)
               inc3=2*i3+ic3
             do ic2=1+int((1-i2)/2),2-int((1+i2)/2)
               inc2=2*i2+ic2
             do ic1=1+int((1-i1)/2),2-int((1+i1)/2)
               inc1=2*i1+ic1
               ineighbor=neighbor_child(1,inc1,inc2,inc3,igrid)
               ipe_neighbor=neighbor_child(2,inc1,inc2,inc3,igrid)
               ixmin1=ixMlo1+(ic1-1)*nxCo1
               ixmin2=ix
               ixmin3=ixMlo3+(ic3-1)*nxCo3
               ixmax1=ixmin1-1+nxCo1
               ixmax2=ix
               ixmax3=ixmin3-1+nxCo3

               if (ipe_neighbor==mype) then
                 iotherside=3-iside

                   
                   do ix3=1,nxCo3 
                     do ix1=1,nxCo1 
                       do iw=1,nwfluxin
                         bg(bgstep)%w(ixmin1+ix1-1,ix,ixmin3+ix3-1,nw0+iw-1,igrid) = &
                           bg(bgstep)%w(ixmin1+ix1-1,ix,ixmin3+ix3-1,nw0+iw-1,igrid) + &
                           pflux(iotherside,2)%flux(ix1,1,ix3,&
                              iw,ineighbor) * CoFiratio
                       end do
                     end do
                   end do
               !else
               !  if (slab_uniform) then
               !    ibufnext=ibuf+isize(2)
               !    if(stagger_grid) ibufnext=ibufnext-isize_stg(2)
               !    psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,ixmin3:ixmax3,&
               !       nw0:nw1) = psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !       ixmin3:ixmax3,nw0:nw1)+CoFiratio &
               !       *reshape(source=recvbuffer(ibuf:ibufnext-1),&
               !        shape=shape(psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !       ixmin3:ixmax3,nw0:nw1)))
               !    ibuf=ibuf+isize(2)
               !  else
               !    ibufnext=ibuf+isize(2)
               !    if(stagger_grid) then
               !      nbuf=(isize(2)-isize_stg(2))/nwfluxin
               !    else
               !      nbuf=isize(2)/nwfluxin
               !    end if
               !    do iw=nw0,nw1
               !      psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,ixmin3:ixmax3,&
               !         iw)=psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !         ixmin3:ixmax3,iw) &
               !         +reshape(source=recvbuffer(ibuf:ibufnext-1),&
               !          shape=shape(psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !         ixmin3:ixmax3,iw))) /ps(igrid)%dvolume(ixmin1:ixmax1,&
               !         ixmin2:ixmax2,ixmin3:ixmax3)
               !      ibuf=ibuf+nbuf
               !    end do
               !    ibuf=ibufnext
               !  end if
               else
                   ! Two transverse indices plus the variable index: every iteration
                   ! lands in a distinct cell of a distinct variable, and the buffer
                   ! offset is a pure function of the three, so all three collapse.
                   
                   do ix3=1,nxCo_fc(3,2)
                     do ix1=1,nxCo_fc(1,2)
                       do iw=1,nwfluxin
                         bg(bgstep)%w(ixmin1+ix1-1,ix,ixmin3+ix3-1,nw0+iw-1,igrid) = &
                           bg(bgstep)%w(ixmin1+ix1-1,ix,ixmin3+ix3-1,nw0+iw-1,igrid) + &
                           recvbuffer(ibuf_offset(4**3*(igrid-1)+inc1+4*inc2+16*inc3+1) &
                              +(ix1-1)+(ix3-1)*nxCo_fc(1,2)+(iw-1)*nxCo_fc(1,2)*nxCo_fc(3,2)) * CoFiratio
                       end do
                     end do
                   end do
               end if
             end do
             end do
             end do
           end do

           case (3)
           do iside=1,2
             i1=kr(1,3)*(2*iside-3)
             i2=kr(2,3)*(2*iside-3)
             i3=kr(3,3)*(2*iside-3)

             if (neighbor_pole(i1,i2,i3,igrid)/=0) cycle

             if (neighbor_type(i1,i2,i3,igrid)/=4) cycle

 !opedit: skip over active/passive interface since flux for passive ones is
             ! not computed, keep the buffer counter up to date:
            !   if (.not.neighbor_active(i1,i2,i3,&
            !      igrid).or..not.neighbor_active(0,0,0,igrid) ) then
            !     do ic3=1+int((1-i3)/2),2-int((1+i3)/2)
            !     inc3=2*i3+ic3
            ! do ic2=1+int((1-i2)/2),2-int((1+i2)/2)
            !     inc2=2*i2+ic2
            ! do ic1=1+int((1-i1)/2),2-int((1+i1)/2)
            !     inc1=2*i1+ic1
            !     ipe_neighbor=neighbor_child(2,inc1,inc2,inc3,igrid)
            !     if (ipe_neighbor/=mype) then
            !       ibufnext=ibuf+isize(3)
            !       ibuf=ibufnext
            !     end if
            !     end do
            ! end do
            ! end do
            !     cycle
            !   end if

             select case (iside)
             case (1)
               ix=ixMlo3
             case (2)
               ix=ixMhi3
             end select

             ! remove coarse flux
               
               do ix2=ixMlo2,ixMhi2
                 do ix1=ixMlo1,ixMhi1 
                   do iw=1,nwfluxin
                     bg(bgstep)%w(ix1,ix2,ix,nw0+iw-1,igrid) = &
                       bg(bgstep)%w(ix1,ix2,ix,nw0+iw-1,igrid) - &
                       pflux(iside,3)%flux(ix1-nghostcells,ix2-nghostcells,&
                          1,iw,igrid)
                   end do
                 end do
               end do


             ! add fine flux
             do ic3=1+int((1-i3)/2),2-int((1+i3)/2)
               inc3=2*i3+ic3
             do ic2=1+int((1-i2)/2),2-int((1+i2)/2)
               inc2=2*i2+ic2
             do ic1=1+int((1-i1)/2),2-int((1+i1)/2)
               inc1=2*i1+ic1
               ineighbor=neighbor_child(1,inc1,inc2,inc3,igrid)
               ipe_neighbor=neighbor_child(2,inc1,inc2,inc3,igrid)
               ixmin1=ixMlo1+(ic1-1)*nxCo1
               ixmin2=ixMlo2+(ic2-1)*nxCo2
               ixmin3=ix
               ixmax1=ixmin1-1+nxCo1
               ixmax2=ixmin2-1+nxCo2
               ixmax3=ix
               if (ipe_neighbor==mype) then
                 iotherside=3-iside
                   
                   do ix2=1,nxCo2 
                     do ix1=1,nxCo1 
                       do iw=1,nwfluxin
                         bg(bgstep)%w(ixmin1+ix1-1,ixmin2+ix2-1,ix,nw0+iw-1,igrid) = &
                           bg(bgstep)%w(ixmin1+ix1-1,ixmin2+ix2-1,ix,nw0+iw-1,igrid) + &
                           pflux(iotherside,3)%flux(ix1,ix2,1,iw,&
                              ineighbor)* CoFiratio
                       end do
                     end do
                   end do
               !else
               !  if (slab_uniform) then
               !    ibufnext=ibuf+isize(3)
               !    if(stagger_grid) ibufnext=ibufnext-isize_stg(3)
               !    psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,ixmin3:ixmax3,&
               !       nw0:nw1) = psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !       ixmin3:ixmax3,nw0:nw1)+CoFiratio &
               !       *reshape(source=recvbuffer(ibuf:ibufnext-1),&
               !        shape=shape(psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !       ixmin3:ixmax3,nw0:nw1)))
               !    ibuf=ibuf+isize(3)
               !  else
               !    ibufnext=ibuf+isize(3)
               !    if(stagger_grid) then
               !      nbuf=(isize(3)-isize_stg(3))/nwfluxin
               !    else
               !      nbuf=isize(3)/nwfluxin
               !    end if
               !    do iw=nw0,nw1
               !      psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,ixmin3:ixmax3,&
               !         iw)=psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !         ixmin3:ixmax3,iw) &
               !         +reshape(source=recvbuffer(ibuf:ibufnext-1),&
               !          shape=shape(psb(igrid)%w(ixmin1:ixmax1,ixmin2:ixmax2,&
               !         ixmin3:ixmax3,iw))) /ps(igrid)%dvolume(ixmin1:ixmax1,&
               !         ixmin2:ixmax2,ixmin3:ixmax3)
               !      ibuf=ibuf+nbuf
               !    end do
               !    ibuf=ibufnext
               !  end if
               else
                   ! Two transverse indices plus the variable index: every iteration
                   ! lands in a distinct cell of a distinct variable, and the buffer
                   ! offset is a pure function of the three, so all three collapse.
                   
                   do ix2=1,nxCo_fc(2,3)
                     do ix1=1,nxCo_fc(1,3)
                       do iw=1,nwfluxin
                         bg(bgstep)%w(ixmin1+ix1-1,ixmin2+ix2-1,ix,nw0+iw-1,igrid) = &
                           bg(bgstep)%w(ixmin1+ix1-1,ixmin2+ix2-1,ix,nw0+iw-1,igrid) + &
                           recvbuffer(ibuf_offset(4**3*(igrid-1)+inc1+4*inc2+16*inc3+1) &
                              +(ix1-1)+(ix2-1)*nxCo_fc(1,3)+(iw-1)*nxCo_fc(1,3)*nxCo_fc(2,3)) * CoFiratio
                       end do
                     end do
                   end do
               end if
             end do
             end do
             end do
           end do
         end select
       end do
     end do

     if (n_send_pe>0) then
       call MPI_WAITALL(n_send_pe,fc_sendreq,fc_sendstat,ierrmpi)
     end if

   end subroutine fix_conserve
fix_conserve Subroutine

Contents

Source Code

public subroutine fix_conserve(bgstep, idimmin, idimmax, nw0, nwfluxin)

Uses

Arguments

Calls

Called by

Source Code
