xref: /phasta/phSolver/incompressible/forces.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      subroutine forces (y ,  ilwork)
c
c----------------------------------------------------------------------
c
c This subroutine calculates and updates the aerodynamic forces
c
c----------------------------------------------------------------------
c
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
c
      integer ilwork(nlwork)

      real*8  y(nshg,ndof)
      real*8  Forin(5), Forout(5),spmasss
      real*8  ftots(3,0:MAXSURF),ftot(3)

!MR CHANGE
        integer :: iv_rankpernode, iv_totnodes, iv_totcores
        integer :: iv_node, iv_core, iv_thread
!MR CHANGE

c
c  START OF DISTURBANCE BLOCK
c
      dimension diste(1), distesum(1)
      distcalc=0  ! we would want to activate this for T-S instability studies
      if(distcalc.eq.1) then
        if(iter.eq.nitr) then  ! we have completed the last iteration
          if (numpe > 1) then
          call MPI_REDUCE (dke,dkesum,1,MPI_DOUBLE_PRECISION,
     &                     MPI_SUM, master, MPI_COMM_WORLD,ierr)
          endif
          if (numpe.eq.1)
     &      write (76,1000) lstep+1, dke,dke

          if ((myrank .eq. master).and.(numpe > 1)) then
             write (76,1000) lstep+1, dkesum,dkesum
c
             call flush(76)
c
          endif
          dke=0.0               ! we must zero it back out for next step

       endif
       endif

C
C END OF DISTURBANCE BLOCK
C

c
c.... -------------------->  Aerodynamic Forces  <----------------------
c
c.... output the forces and the heat flux
c
      if (numpe > 1) then
         Forin = (/ Force(1), Force(2), Force(3), HFlux,
     &              entrop /)
         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE (Forin(1), Forout(1),5,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
         Force = Forout(1:3)
         HFlux = Forout(4)
         entrop= Forout(5)
      endif

      if (numpe > 1) then
         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE (flxID(1,isrfIM), Atots,1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE (flxID(2,isrfIM), spmasss,1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE (flxID(3,:), Ftots(1,:),MAXSURF+1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE (flxID(4,:), Ftots(2,:),MAXSURF+1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
         call MPI_BARRIER(MPI_COMM_WORLD, ierr)
         call MPI_ALLREDUCE (flxID(5,:), Ftots(3,:),MAXSURF+1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
      else
         Ftots=flxID(3:5,:)
         Atots=flxID(1,isrfIM)
         spmasss=flxID(2,isrfIM)
      endif
      call sforce()

      ftot(1)=sum(Ftots(1,0:MAXSURF))
      ftot(2)=sum(Ftots(2,0:MAXSURF))
      ftot(3)=sum(Ftots(3,0:MAXSURF))

      if (myrank .eq. master) then
         write (iforce,1000) lstep, (Force(i), i=1,nsd),
     &        HFlux,spmasss,(ftot(i),i=1,nsd)
         call flush(iforce)
      endif

      if((spmasss.ne.0) .and. (matflg(5,1).ne.0))then ! we are doing
                                                      ! force adjustment
         tmp=Dtgl*0.025         ! .025= 1/2/tmp1 when tmp1=20

         select case ( matflg(5,1) )
         case ( 1 )             ! standard linear body force
            fwall=Force(1)
            vlngth=xlngth
         case ( 2 )
         case ( 3 )
            fwall=Force(3)
            vlngth=zlngth
         end select

         write(222,*) spmasss
         fnew=(vel-spmasss/(ro*Atots))*tmp +
     &        fwall/(vlngth*Atots)
         if(myrank.eq.0) then
            write(880,*) datmat(1,5,1),fnew
            call flush(880)
         endif
         datmat(1,5,1)=fnew     !* (one - tmp2) + datmat(1,5,1) * tmp2
      endif

      if (pzero .eq. 1) then
c

c  also adjust the pressure here so that the mean stays at zero (prevent drift)
c
c  This is only valid if there are NO Dirichlet Boundary conditions on p!!!
c
c  (method below counts partition boundary nodes twice. When we get ntopsh
c  into the method it will be easy to fix this.  Currently it would require
c  some work and it is probably not worth it
c
c  switched to avg of the on-processor averages to get around hierarchic
c  modes problem
c
c              pave=sum(y(1:numnp,1))
c              call MPI_ALLREDUCE (pave, pavet, 1,
c     &             MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
c     pavet=pavet/nshgt
         pave=sum(y(1:numnp,4))
         if (numpe .gt. 1) then
            xnpts=numnp
           call MPI_BARRIER(MPI_COMM_WORLD, ierr)
            call MPI_ALLREDUCE (pave, pavet, 1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
            call MPI_BARRIER(MPI_COMM_WORLD, ierr)
            call MPI_ALLREDUCE (xnpts, xnptst, 1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
            pavet=pavet/xnptst
         else
            pavet = pave/numnp
         endif

         y(1:numnp,4)=y(1:numnp,4)-pavet
         pave=sum(y(1:numnp,4))
c         if(myrank.eq.0) then
c            write(90+myrank,*) pavet,pave
c            call flush(900)
c         endif
c
      endif

      return
c
! 1000 format(1p,i6,8e13.5)
 1000 format(1p,i6,8e20.12)
c
      end

      subroutine sforce()

      include "common.h"
      include "mpif.h"

      real*8  SFlux(10),SFluxg(10)
      character*20 fname1
      character*5  cname
      character*10  cname2

      integer icalled
      data icalled /0/
      save icalled

      if(icalled.eq.0) then
        icalled = 1

!MR CHANGE
        iv_rankpernode = iv_rankpercore*iv_corepernode
        iv_totnodes = numpe/iv_rankpernode
        iv_totcores = iv_corepernode*iv_totnodes
        irankfilesforce(:) = -1
!MR CHANGE

        do isrf = 0,MAXSURF

!MR CHANGE
          ! Compute the adequate rank which will take care of isrf
          iv_node = (iv_totnodes-1)-mod(isrf,iv_totnodes)
          iv_core = (iv_corepernode-1) - mod((isrf -
     &              mod(isrf,iv_totnodes))/iv_totnodes,iv_corepernode)
          iv_thread = (iv_rankpercore-1) - mod((isrf-
     &              (mod(isrf,iv_totcores)))/iv_totcores,iv_rankpercore)
          irankfilesforce(isrf) = iv_node*iv_rankpernode
     &                     + iv_core*iv_rankpercore
     &                     + iv_thread

          if(myrank == 0) then
              write(*,*) '  sforce ', isrf, 'if any handled by rank',
     &                      irankfilesforce(isrf), ' on node', iv_node
          endif

          ! Verification just in case
          if(irankfilesforce(isrf) .lt.0 .or.
     &                 irankfilesforce(isrf) .ge. numpe) then
            write(*,*) 'WARNING: irankfilesforce(',isrf,') is ',
     &                                  irankfilesforce(isrf),
     &                      ' and reset to numpe-1'
            irankfilesforce(isrf) = numpe-1
          endif
!MR CHANGE

!          if ( nsrflist(isrf).ne.0 .and. myrank.eq.master) then
          if ( nsrflist(isrf).ne.0 .and.
     &                     myrank.eq.irankfilesforce(isrf)) then
            iunit=60+isrf
            fname1 = 'forces_s'//trim(cname(isrf))// cname2(lskeep+1)
            open(unit=iunit,file=trim(fname1),status='unknown')
          endif
        enddo
      endif

      do isrf = 0,MAXSURF
        if ( nsrflist(isrf).ne.0 ) then
          SFlux(:) = flxID(:,isrf) ! Area, Mass Flux, Force 1, Force 2, Force 3, Scalars
          if ( numpe > 1 ) then
            call MPI_BARRIER(MPI_COMM_WORLD, ierr)
            call MPI_ALLREDUCE (SFlux, SFluxg,10,
     .        MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr)
            SFlux = SFluxg
          endif
!          if(myrank.eq.master) then
          if(myrank.eq.irankfilesforce(isrf)) then
            iunit=60+isrf
!            write(iunit,"(i7,1p10e14.5)") lstep, SFlux
            write(iunit,"(i7,1p10e20.12)") lstep, SFlux
            call flush(iunit)
! done in itrdrv.f
!            if(istep.eq.nstep(1)) then !end step then close file
!              close(iunit)
!            endif
          endif
        endif
      enddo

      return
      end