xref: /phasta/phSolver/common/proces.f (revision 53242657c9ca47d9327f9d6987e18a64b4277295)

        subroutine proces
c
c----------------------------------------------------------------------
c
c This subroutine generates the problem data and calls the solution
c  driver.
c
c
c Zdenek Johan, Winter 1991.  (Fortran 90)
c----------------------------------------------------------------------
c
        use readarrays          ! used to access x, iper, ilwork
        use turbsa          ! used to access d2wall
        use dtnmod
        use periodicity
        use pvsQbi
        include "common.h"
        include "mpif.h"
c
c arrays in the following 2 lines are now dimensioned in readnblk
c        dimension x(numnp,nsd)
c        dimension iper(nshg), ilwork(nlwork)
c
        dimension y(nshg,ndof),
     &            iBC(nshg),
     &            BC(nshg,ndofBC),
     &            ac(nshg,ndof)
c
c.... shape function declarations
c
        dimension shp(MAXTOP,maxsh,MAXQPT),
     &            shgl(MAXTOP,nsd,maxsh,MAXQPT),
     &            shpb(MAXTOP,maxsh,MAXQPT),
     &            shglb(MAXTOP,nsd,maxsh,MAXQPT)
c
c  stuff for dynamic model s.w.avg and wall model
c
c        dimension ifath(numnp),    velbar(nfath,nflow),
c     &            nsons(nfath)

        dimension velbar(nfath,nflow)
c
c stuff to interpolate profiles at inlet
c
        real*8 bcinterp(100,ndof+1),interp_mask(ndof)
        logical exlog

        !Duct
        real*8 c0, c1, c2, c3, x1, x2
        integer nn

c        if ((irscale .ge. 0) .and. (myrank.eq.master)) then
c           call setSPEBC(numnp, point2nfath, nsonmax)
c        endif
c
c.... generate the geometry and boundary conditions data
c
        call gendat (y,              ac,             point2x,
     &               iBC,            BC,
     &               point2iper,     point2ilwork,   shp,
     &               shgl,           shpb,           shglb,
     &               point2ifath,    velbar,         point2nsons )
        call setper(nshg)
        call perprep(iBC,point2iper,nshg)
        if (iLES/10 .eq. 1) then
        call keeplhsG ! Allocating space for the mass (Gram) matrix to be
                      ! used for projection filtering and reconstruction
                      ! of the strain rate tensor.

        call setrls   ! Allocating space for the resolved Leonard stresses
c                         See bardmc.f
        endif
c
c.... time averaged statistics
c
        if (ioform .eq. 2) then
           call initStats(point2x, iBC, point2iper, point2ilwork)
        endif
c
c.... RANS turbulence model
c
!        if (iRANS .lt. 0) then
           call initTurb( point2x )
!        endif
c
c.... p vs. Q boundary
c
           call initNABI( point2x, shpb )
c
c.... check for execution mode
c
        if (iexec .eq. 0) then
           lstep = 0
           call restar ('out ',  y  ,ac)
           return
        endif
c
c.... initialize AutoSponge
c
        if(matflg(5,1).ge.4) then ! cool case (sponge)
           call initSponge( y,point2x)
        endif
c
c
c.... adjust BC's to interpolate from file
c

        inquire(file="inlet.dat",exist=exlog)
        if(exlog) then
           open (unit=654,file="inlet.dat",status="old")
           read(654,*) ninterp,ibcmatch,(interp_mask(j),j=1,ndof)
! with no surfID control, this will also get applied to isothermal walls but! usually ok and can be extended if not.

           do i=1,ninterp
              read(654,*) (bcinterp(i,j),j=1,ndof+1) ! distance to wall+
                        ! ndof but note order of BC's
                        ! p,t,u,v,w,scalars. Also note that must be in
                        ! increasing distance from the wall order.

           enddo
           do i=1,nshg  ! only correct for linears at this time
             iupper=0
             do j=2,ninterp
               if(bcinterp(j,1).gt.d2wall(i)) then !bound found
                 xi=(d2wall(i)-bcinterp(j-1,1))/
     &           (bcinterp(j,1)-bcinterp(j-1,1))
                 iupper=j
                 exit
               endif
             enddo
             if(iupper.eq.0) then
               write(*,*) 'failure in finterp, ynint=',d2wall(i)
               stop
             endif
             if(lstep.eq.0) then  ! apply this inlet bl as an IC if on step 0
              y(i,1:3)=(xi*bcinterp(iupper,4:6)
     &                +(one-xi)*bcinterp(iupper-1,4:6))
              y(i,4)=(xi*bcinterp(iupper,2)
     &                +(one-xi)*bcinterp(iupper-1,2))
              y(i,5)=(xi*bcinterp(iupper,3)
     &                +(one-xi)*bcinterp(iupper-1,3))
             endif
             if(point2x(i,1).lt.1.0e-5 .and. mod(iBC(i),1024).eq.ibcmatch) then
               do j=1,ndof
                 if(interp_mask(j).ne.zero) then
                   BC(i,j)=(xi*bcinterp(iupper,j+1)
     &             +(one-xi)*bcinterp(iupper-1,j+1))
                 endif
               enddo
             endif
           enddo
        endif
c$$$$$$$$$$$$$$$$$$$$

!======================================================================
!Modifications for Duct. Need to encapsulate in a function call.
        !specify an initial eddy viscosity ramp
        if(isetEVramp .gt. 0) then
          if(myrank .eq. 0) then
            write(*,*) "Setting eddy viscosity ramp with:"
            write(*,*) "  - ramp X min = ", EVrampXmin
            write(*,*) "  - ramp X max = ", EVrampXmax
            write(*,*) "  - EV min = ", EVrampMin
            write(*,*) "  - EV max = ", EVrampMax
          endif

          x1 = EVrampXmin  !stuff in a shorter variable name to
          x2 = EVrampXmax  !make the formulas cleaner
          !Newton Divided differences to generate a polynomial of
          !the form P(x) = c0 + x*(c1 + x*(c2 + (x - (x2 - x1))*c3))
          !satisfying P(x1) = EVrampMin, P(x2) = EVrampMax,
          ! P'(x1) = 0, and P'(x2) = 0

          c0 = EVrampMin
          c1 = 0            !zero derivative
          c2 = (EVrampMax - EVrampMin)/(x2 - x1)
          c3 = 0            !zero derivative
          c3 = (c3 - c2)/(x2 - x1)
          c2 = (c2 - c1)/(x2 - x1)
          c3 = (c3 - c2)/(x2 - x1)

          do nn = 1,nshg
            if(y(nn,6) .eq. 0) cycle  !don't change wall boundary conditions, should be iTurbWall == 1

            if(point2x(nn,1) .gt. EVrampXmax) then !downstream of the ramp
              y(nn,6) = EVrampMax
            elseif(point2x(nn,1) .gt. EVrampXmin) then !and x(:,1) <= EVrampXmax

              !P(x) = c0 + x*(c1 + x*(c2 + (x - (x2 - x1))*c3))
              !     = c0 + x*(c1 + x*(c2 - (x2 - x1)*c3 + x*c3
              y(nn,6) = c0                 + (point2x(nn,1) - x1)*(
     &                  c1                 + (point2x(nn,1) - x1)*(
     &                 (c2 - (x2 - x1)*c3) + (point2x(nn,1) - x1)*c3))
            else
              y(nn,6) = EVrampMin
            endif
          enddo
        endif
!End modifications for Duct
!======================================================================
c
c
c.... call the semi-discrete predictor multi-corrector iterative driver
c
        call itrdrv (y,              ac,
     &               uold,           point2x,
     &               iBC,            BC,
     &               point2iper,     point2ilwork,   shp,
     &               shgl,           shpb,           shglb,
     &               point2ifath,    velbar,         point2nsons )
c
c.... return
c
c
c.... stop CPU-timer
c
CAD        call timer ('End     ')
c
c.... close echo file
c

!MR CHANGE
      if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
      if(myrank.eq.0)  then
          write(*,*) 'process - before closing iecho'
      endif
!MR CHANGE END

        close (iecho)

!MR CHANGE
      if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr)
      if(myrank.eq.0)  then
          write(*,*) 'process - after closing iecho'
      endif
!MR CHANGE END


c
c.... end of the program
c
CAD        write(6,*) 'Life: ', second(0) - ttim(100)
        deallocate(point2iper)
        if(numpe.gt.1) then
          call Dctypes(point2ilwork(1))
        endif
        deallocate(point2ilwork)
        deallocate(point2x)
        deallocate(point2nsons)
        deallocate(point2ifath)
        deallocate(uold)
        deallocate(wnrm)
        deallocate(otwn)
        call finalizeDtN
        call clearper
        call finalizeNABI

        return
        end