xref: /phasta/phSolver/common/INIprofile.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)

        subroutine INIprofile(BC,y,x)

         include "common.h"
         real*8 BC(nshg,ndofBC)
         real*8 y(nshg,ndof)
         real*8 x(nshg,nsd)

         if(iI2Binlet .gt. 0)then ! useless for geometry with contraction
              call TakeBCfromIC_ScaleInlVelX(BC,y,x)
         endif

         if(isetInitial .gt.0)then ! obscure
              call setInitial(x,y)
         endif

        return
        end
c..............................................................................
        subroutine TakeBCfromIC_ScaleInlVelX(BC,y,x)

        use BCsfIDmap
        include "common.h"
        include "mpif.h"
        include "auxmpi.h"

        integer inlf_tmp(nshg)
        real*8 BC(nshg,ndofBC)
        real*8 y(nshg,ndof)
        real*8 x(nshg,nsd)
        real*8 u_max,u_max_all

         call sfID2np(iI2Binlet,ninlet,inlf_tmp)
         if(ninlet.gt.0) then  ! only for process including inlet nodes
            allocate(inlf(ninlet))
            inlf(1:ninlet)=inlf_tmp(1:ninlet) ! generate map mapping inlet nodes to nshg nodes
            do i=1,ninlet
               nn=inlf(i)
               if(y(nn,1).gt.u_max)u_max=y(nn,1)
            enddo
         endif
         call MPI_BARRIER(MPI_COMM_WORLD,ierr)
         call MPI_ALLREDUCE( u_max, u_max_all, 1,
     &        MPI_REAL8, MPI_MAX, MPI_COMM_WORLD,ierr)
         if(myrank.eq.0)then
           write(*,*)''
           write(*,556)'u_max_all=',u_max_all ! find the largest x velocity at inlet
         endif
c apply IC to BC, must be run before applying BC to IC later in genini
         if(ninlet.gt.0) then
            BC(inlf(:),2) = y(inlf(:),5) ! set temp
            BC(inlf(:),3) = y(inlf(:),1)/u_max_all*inletVelX  ! set and scale x velocity
            BC(inlf(:),4) = 0             ! set y velocity
            BC(inlf(:),5) = 0             ! set z velocity
            if(nsclr.eq.1)then
              BC(inlf(:),7) = y(inlf(:),6)  ! set scalar_1
            endif
         endif
555          format(I3,F10.3)
556          format(A,F10.3)

        return
        end

c..............................................................................
c... this subroutine gives the formular of initial conditions, user may want to
c... change this file and recompile the code, just like other commercial code which
c... supports User Defined Function, you need to recompile and relink the code
        subroutine setInitial(x,y)

        use turbSA   ! this gives us d2wall(1:nshg) nodal distance to the closest wall
        use BCsfIDmap
        include "common.h"
        include "mpif.h"
        include "auxmpi.h"

        real*8 x(nshg,nsd)
        real*8 y(nshg,ndof)
        if(isetinitial.eq.1) then

c... user gives the formular of initial conditions from here
           y(:,1)= xvel_ini
           y(:,4)= pres_ini
           y(:,2)=yvel_ini
           y(:,3)=zvel_ini
           y(:,5)=temp_ini
           if(nsclr.eq.1)then
              y(:,6)=evis_ini
           endif
!moving previously coded spatially varying pressure to 10 so that 1 is simple
        elseif (isetInitial.eq.10) then
        if(ninlet.gt.0)then
           xinlet=x(inlf(1),1)
        else
           xinlet=1e18
        endif

        if(noutlet.gt.0)then
           xoutlet=x(outf(1),1)
        else
           xoutlet=-1e18
        endif

        call MPI_BARRIER(MPI_COMM_WORLD,ierr)
        call MPI_ALLREDUCE(xinlet, xinlet_all, 1,
     &       MPI_REAL8, MPI_MIN, MPI_COMM_WORLD,ierr)
        call MPI_ALLREDUCE(xoutlet, xoutlet_all, 1,
     &       MPI_REAL8, MPI_MAX, MPI_COMM_WORLD,ierr)
        call MPI_BARRIER(MPI_COMM_WORLD,ierr)
        do nn=1,nshg
           xdel=(x(nn,1)-xinlet_all)/(xoutlet_all-xinlet_all)
           xdel=1-tanh(7*xdel)
           y(nn,1)= xvel_ini + (inletVelX-xvel_ini)*xdel
           y(nn,4)= outPres1 + (pres_ini-outPres1)*xdel
        enddo
           y(:,2)=yvel_ini
           y(:,3)=zvel_ini
           y(:,5)=temp_ini
           if(nsclr.eq.1)then
              y(:,6)=evis_ini
           endif
        elseif (isetInitial.eq.2) then
            !----------------------------------
            ! set initial condition for blower
            !----------------------------------

!            blthickness=1.0e-4
!            yfloor=-0.0333375 - 1e-7
!            xthroat=-0.3048 + 1e-5
!            xSlit = 0.025            !x-coordinate of the blower slit
!
!            p0 = 121800       !total pressure in blower
!            T0 = 305          !total temperature
!
!            do nn=1,nshg
!                xcoor = x(nn,1)
!
!                if((xcoor .gt. xthroat) .and. (x(nn,2).lt.yfloor)) then
!                    y(nn,3:4)=0  !velocity
!                    y(nn,6)=1e-3 !eddy viscosity
!
!                    if(xcoor .le. xSlit) then  !test whether the point is infront of or behind the slit. If infront, use M = 0.8
!                        !4th order polynomial best fit for an exit Mach number of 0.8
!                        MEuler=0.0806961 + xcoor*(2.96302 + xcoor*(72.5101 + xcoor*(13914.6 + xcoor*1.09778e6)))
!                    else
!                        MEuler = 0.8
!                    endif
!
!                    pEuler = p0/((1 + 0.2*M*M)**3.5)
!                    TEuler = T0/(1 + 0.2*M*M)
!                    uEuler = MEuler*sqrt(1.4*287.06*TEuler)
!
!                    bldamper=min(one, d2wall(nn)/blthickness)
!                    y(nn,2)=uEuler*bldamper
!                    y(nn,1)=pEuler
!                    y(nn,5)=TEuler
!                endif
!             enddo
        endif

        return
        end

c=========================================================================
c Set the initial condition based on coordiantes
c=========================================================================

        subroutine setInitial_Duct(x)

        use readarrays
        include "common.h"
        include "mpif.h"
        include "auxmpi.h"

        real*8 x(nshg,nsd)
        real*8 xcoor,ycoor,zcoor
        integer nn
        real*8 xVel, pres, Temp
!p u v w T s(calar_1) in array qold
!x y z in array x



        a1 =       329.7
        b1 =    -0.07732
        c1 =        9.82
        a2 =       6.431
        b2 =       1.621
        c2 =      0.8998
        a3 =      -4.637
        b3 =       1.728
        c3 =      0.2168


        if (iDuctgeometryType .eq. 6) then
           xContraStart = -74*0.0254;
        elseif (iDuctgeometryType .eq. 8) then
           xContraStart = -84*0.0254;
        endif
        xContraEnd = xContraStart + 72*0.0254

        do nn=1,nshg

           xcoor=x(nn,1)
           ycoor=x(nn,2)
           zcoor=x(nn,3)
c.... [xVel,pres,Temp] are functions of [xcoor,ycoor,zcoor]

c          if(ycoor<(-0.587375+1.0e-4))then
c             ry=(ycoor+0.587375)/1.0e-4
c          elseif(ycoor>(0.587375+1.0e-4))then
c             ry=(0.587375-ycoor)/1.0e-4
c          else
c             ry=1.0
c          endif

c          if(zcoor<(-0.587375+1.0e-4))then
c             rz=(zcoor+0.587375)/1.0e-4
c          elseif(zcoor>(0.587375+1.0e-4))then
c             rz=(0.587375-zcoor)/1.0e-4
c          else
c             rz=1.0
c          endif

c          ry=max(0.0,ry)
c          rz=max(0.0,rz)

c          if(xcoor.le.0.0.and.xcoor.gt.-12*0.0254)then
c             xVel=120
c          elseif(xcoor.lt.-12*0.0254)then
c             xVel=(1.013158+120)/2+(120-1.013158)/2*
c     &       tanh((xcoor-(-84*0.0254-12*0.0254)/2)*4)
c          else
c             xVel=max(0.0,120*(1-(xcoor/(1.6*4.5*0.0254+0.85))**2))
c          endif

c          Temp=317+(-tanh(4*(xcoor+0.35+0.254))*
c     &       (13.578+5)*0.5+(13.578-5)*0.5)*ry*rz
c          pres=-tanh(4*(xcoor+0.35+0.254))*
c     &    (116600-97000)*0.5+(116600+97000)*0.5
c... above, Onkar's method

c         if(xcoor.lt.0.0 .and.
c     &    xcoor.gt.12*0.0254 .and.
c     &    ycoor .lt. -0.04446)then
c            xVel=0
c            pres=97800
c            Temp=317
c         else
c            if(xcoor .ge. -0.0254)then
c             xVel=158
c             pres=97800
c             Temp=318
c            elseif(xcoor .le. -1.0163)then
c             Temp=330
c             pres=111896
c             xVel=1.0132
c            else
c             pres=97800+0.5*(111896-97800)*
c     &       (tanh(-3.4*((xcoor+0.52085)/0.49545))+1)
c             Temp=318+0.5*(330-318)*
c     &       (tanh(-3.4*((xcoor+0.52085)/0.49545))+1)
c             xVel=158+0.5*(1.0132-158)*
c     &       (tanh(-3.4*((xcoor+0.52085)/0.49545))+1)
c            endif
c         endif
c....................


c             pres= 97000
c             xvel= 0.1

c             if (xcoor .gt. xContraEnd)then
c                 Temp = 320
c             elseif (xcoor . le. xContraEnd)then
c                 xSim = xcoor - xContraStart
c                 Temp = a1*exp(-((xSim-b1)/c1)**2) +
c     &                  a2*exp(-((xSim-b2)/c2)**2) +
c     &                  a3*exp(-((xSim-b3)/c3)**2)
c             endif
c             varSA = 0.0

c             qold(nn,1)=pres
c             qold(nn,2)=xvel
c             qold(nn,3)=0
c             qold(nn,4)=0
c             qold(nn,5)=Temp
c             qold(nn,6)=varSA

c... above Yi Chen's Method

          if(xcoor .lt. 0.0 .and.
     &       xcoor .gt. (-12*0.0254) .and.
     &       ycoor .lt. -0.044449)then ! this is the jet device
             qold(nn,1)=97000
             qold(nn,2)=0
             qold(nn,3)=0
             qold(nn,4)=0
             qold(nn,5)=320
             qold(nn,6)=1.825e-5
          endif

c... above, used for restart from M2M

        enddo ! end of loop over all nshg points
c...........................

        return
        end