xref: /phasta/phSolver/common/aveprep.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      module aveall

      real*8, allocatable :: ylist(:)
      integer, allocatable :: ifathe(:,:)

      end module

c------------------------------------------------------------------------------

      subroutine setave

      use aveall

      include "common.h"

      idim = numel*intg(1,1)
      allocate ( ylist(idim) )
      allocate ( ifathe(numel,maxsh) )

c.... return

      return
      end

c------------------------------------------------------------------------------

      subroutine aveprep(shp,x)

      use pointer_data
      use aveall

      include "common.h"

      dimension shp(MAXTOP,maxsh,MAXQPT)
      dimension x(numnp,3)
      integer   nlist

      nlist  = 0
      ylist  = zero
      lfathe = 0
      do iblk = 1,nelblk
        lcsyst = lcblk(3,iblk)
        iel  = lcblk(1,iblk) !Element number where this block begins
        npro = lcblk(1,iblk+1) - iel
        lelCat = lcblk(2,iblk)
        nenl = lcblk(5,iblk)
        nshl = lcblk(10,iblk)
        inum  = iel + npro - 1

        ngauss = nint(lcsyst)

        call getylist( ylist, ifathe(iel:inum,:), shp(lcsyst,1:nshl,:),
     &       x, mien(iblk)%p)


      enddo

      return
      end















      subroutine getnu (ien, strl, xmudmi, cdelsq, lfathe)

      include "common.h"

      dimension  ien(npro,nshl),       strl(npro,ngauss),
     &           xmudmi(npro,ngauss),   cdelsq(nlist),
     &           lfathe(npro,ngauss)

      do intp = 1, ngauss
      do iel  = 1, npro

         xmudmi(iel,intp) = cdelsq( lfathe(iel,intp) )*
     &        strl(iel,intp)

      enddo
      enddo

c
c  local clipping
c
      rmu=datmat(1,2,1)
cc      xmudmi=min(xmudmi,1000.0*rmu) !don't let it get larger than 1000 mu
cc      xmudmi=max(xmudmi, -rmu) ! don't let (xmudmi + mu) < 0

      do int = 1, ngauss
      xmudmi(:,int) =
     &     min(xmudmi(:,int),1000.0*rmu) !don't let it get larger than 1000 mu
      xmudmi(:,int) =
     &     max(xmudmi(:,int), -rmu)    ! don't let (xmudmi + mu) < 0
      enddo

      return
      end
      subroutine getxnut (fres, xden, xnum, ien, shp)

      include "common.h"

      dimension ien(npro,nshl),          fres(nshg,22),
     &          shp(nshl,maxsh),         fresli(npro,22),
     &          fresl(npro,nshl,22),     strnrmi(npro),
     &          xnutl(npro),             fwr(npro),
     &          xmij(npro,6),            xlij(npro,6),
     &          xdenli(npro),            xnumli(npro),
     &          xdenl(npro),             xnuml(npro),
     &          xden(1),                 xnum(1)


      call local (fres, fresl, ien, 22, 'gather  ')

      xnuml  = zero
      xdenl  = zero

      do intp = 1, ngauss

         fresli = zero
         do i = 1, nenl
            do j = 1, 22
               fresli(:,j) = fresli(:,j) + shp(i,intp)*fresl(:,i,j)
            enddo
         enddo

         strnrmi(:) =  sqrt(
     &   two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2)
     &   + four * ( fresli(:,13)**2 + fresli(:,14)**2 +
     &   fresli(:,15)**2 ) )

         fwr = fwr1 * fresli(:,22) * strnrmi(:)

         xmij(:,1) = -fwr
     &        * fresli(:,10) + fresli(:,16)

         xmij(:,2) = -fwr
     &        * fresli(:,11) + fresli(:,17)

         xmij(:,3) = -fwr
     &        * fresli(:,12) + fresli(:,18)

         xmij(:,4) = -fwr * fresli(:,13) + fresli(:,19)
         xmij(:,5) = -fwr * fresli(:,14) + fresli(:,20)
         xmij(:,6) = -fwr * fresli(:,15) + fresli(:,21)

         fresli(:,22) = one/fresli(:,22)

         xlij(:,1) = fresli(:,4) - fresli(:,1) *
     &        fresli(:,1) * fresli(:,22)
         xlij(:,2) = fresli(:,5) - fresli(:,2) *
     &        fresli(:,2) * fresli(:,22)
         xlij(:,3) = fresli(:,6) - fresli(:,3) *
     &        fresli(:,3) * fresli(:,22)
         xlij(:,4) = fresli(:,7) - fresli(:,1) *
     &        fresli(:,2) * fresli(:,22)
         xlij(:,5) = fresli(:,8) - fresli(:,1) *
     &        fresli(:,3) * fresli(:,22)
         xlij(:,6) = fresli(:,9) - fresli(:,2) *
     &        fresli(:,3) * fresli(:,22)

         xnumli(:) =    xlij(:,1) * xmij(:,1)
     &        + xlij(:,2) * xmij(:,2)
     &        + xlij(:,3) * xmij(:,3)
     &        + two * (xlij(:,4) * xmij(:,4)
     &        + xlij(:,5) * xmij(:,5)
     &        + xlij(:,6) * xmij(:,6))
         xdenli(:) =        xmij(:,1) * xmij(:,1)
     &        + xmij(:,2) * xmij(:,2)
     &        + xmij(:,3) * xmij(:,3)
     &        + two * (xmij(:,4) * xmij(:,4)
     &        + xmij(:,5) * xmij(:,5)
     &        + xmij(:,6) * xmij(:,6))

c         xdenli = xdenli

         xdenl(:) = xdenl(:) + xdenli(:)
         xnuml(:) = xnuml(:) + xnumli(:)

      enddo

      do nel = 1, npro
         xden(1) = xden(1) + xdenl(nel)
         xnum(1) = xnum(1) + xnuml(nel)
      enddo

      return

      end









      subroutine getylist (ylist, lfathe, shp, x, ien)

      include "common.h"

      dimension lfathe(npro,maxsh),    shp(nshl,maxsh),
     &          x(numnp,3),            ien(npro,nshl),
     &          xl(npro,nenl,nsd),     ylist(idim)

      integer intp

      call localx (x,      xl,     ien,    3,  'gather  ')


      do nel  = 1, npro
      do intp = 1, ngauss

c... Compute the y-coordinate of the current quad. pt. and call it yint.

         yint = zero
         xint = zero
         zint = zero
         do n = 1, nenl
            yint = yint + xl(nel,n,2)*shp(n,intp)
            xint = xint + xl(nel,n,1)*shp(n,intp)
            zint = zint + xl(nel,n,3)*shp(n,intp)
         enddo

c... Check ylist to see if yint is already included in ylist.

         if(nlist .eq. 0)then ! In this case yint is definitely not in ylist.
            nlist            = nlist + 1
            ylist(nlist)     = yint
            lfathe(nel,intp) = nlist

         else

            do ilist = 1, nlist
               imatch = ilist
               if( abs(yint-ylist(ilist)) .lt. 0.00001 ) then
                  lfathe(nel,intp) = imatch
                  goto 5
               endif
            enddo

            nlist            = nlist + 1
            ylist(nlist)     = yint
            lfathe(nel,intp) = nlist

         endif

 5       yint = zero

c         if( nlist .eq. 5)then
c            write(*,*)'ylist yint=',ylist(4),ylist(5)
c            stop
c         endif

      enddo ! End loop over quad. pts.
      enddo ! End loop over elements in current block.

      return
      end
      subroutine projdmc (y,      shgl,      shp,
     &                   iper,   ilwork,       x)

      use pointer_data

      use aveall   ! This module helps us average cdelsq computed at quad pts

      use quadfilt   ! This module gives us shglf(maxtp,nsd,maxsh,ngaussf),
c                    shpf(maxtp,maxsh,ngaussf), and Qwtf(maxtp,ngaussf).
c                    Shpf and shglf are the shape funciotns and their
c                    gradient evaluated using the quadrature rule desired
c                    for computing the dmod. Qwt contains the weights of the
c                    quad. points.

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

c
      dimension fres(nshg,24),         fwr(nshg),
     &          strnrm(nshg),          cdelsq(nshg),
     &          strl(numel,ngauss),
     &          y(nshg,5),            yold(nshg,5),
     &          iper(nshg),
     &          ilwork(nlwork),
     &          x(numnp,3),
     &          shgl(MAXTOP,nsd,maxsh,MAXQPT), shp(MAXTOP,maxsh,MAXQPT),
     &          xden(idim),                    xnum(idim),
     &          xdent(idim),                   xnumt(idim)

      fres = zero
      yold(:,1)=y(:,4)
      yold(:,2:4)=y(:,1:3)
c
c  hack in an interesting velocity field (uncomment to test dmod)
c
c      yold(:,5) = 1.0
c      yold(:,2) = 2.0*x(:,1) - 3*x(:,2)
c      yold(:,3) = 3.0*x(:,1) + 4.0*x(:,2)
c      yold(:,4) = 4.0*x(:,1) + x(:,2) + x(:,3)
c      yold(:,1) = Rgas * yold(:,5) ! Necessary to make model suitable
c                               suitable for the


      intrul=intg(1,itseq)
      intind=intpt(intrul)

      do iblk = 1,nelblk
        lcsyst = lcblk(3,iblk)
        iel  = lcblk(1,iblk) !Element number where this block begins
        npro = lcblk(1,iblk+1) - iel
        lelCat = lcblk(2,iblk)
        nenl = lcblk(5,iblk)
        nshl = lcblk(10,iblk)
        inum  = iel + npro - 1

        ngauss  = nint(lcsyst)
        ngaussf = nintf(lcsyst)

        call asithf (yold, x, strl(iel:inum,:), mien(iblk)%p, fres,
     &               shglf(lcsyst,:,1:nshl,:),
     &               shpf(lcsyst,1:nshl,:),Qwtf(lcsyst,1:ngaussf))

      enddo
c

      if(numpe>1) call commu (fres, ilwork, 24, 'in ')
c
c account for periodicity in filtered variables
c
      do j = 1,nshg
        i = iper(j)
        if (i .ne. j) then
           fres(i,:) = fres(i,:) + fres(j,:)
        endif
      enddo

      do j = 1,nshg
        i = iper(j)
        if (i .ne. j) then
           fres(j,:) = fres(i,:)
        endif
      enddo

      if(numpe>1)then
         call commu (fres, ilwork, 24, 'out')
      endif

      fres(:,23) = one / fres(:,23)
      do j = 1,22
        fres(:,j) = fres(:,j) * fres(:,23)
      enddo

      xden   = zero
      xdent  = zero
      xnum   = zero
      xnumt  = zero

      do iblk = 1,nelblk
        lcsyst = lcblk(3,iblk)
        iel  = lcblk(1,iblk) !Element number where this block begins
        npro = lcblk(1,iblk+1) - iel
        lelCat = lcblk(2,iblk)
        nenl = lcblk(5,iblk)
        inum  = iel + npro - 1
        ngauss = nint(lcsyst)
        call smagcoeff (fres(:,1:22), xden, xnum, mien(iblk)%p,
     &                  shp(lcsyst,1:nshl,:), ifathe(iel:inum,:))

      enddo

      if(numpe.gt.1)then
         call drvAllreduce( xden, xdent, idim )
         call drvAllreduce( xnum, xnumt, idim )
      else
         xdent  = xden
         xnumt  = xnum
      endif

      do i = 1, nlist
         cdelsq(i) = pt5*xnumt(i)/xdent(i)
      enddo

c...  Uncomment for averaging over all directions

c      sumc1 = 0.0
c      sumc2 = 0.0
c      do i = 1, nlist
c         sumc1 = sumc1 + xnumt(i)
c         sumc2 = sumc2 + xdent(i)
c      enddo
c      xfact = pt5*sumc1/sumc2
c      cdelsq(:) = xfact
c      if(myrank .eq. master)then
c         write(22,*)cdelsq(1)
c      endif

c...  End of averaging over all directions.

      do iblk = 1,nelblk
        lcsyst = lcblk(3,iblk)
        iel  = lcblk(1,iblk) !Element number where this block begins
        npro = lcblk(1,iblk+1) - iel
        lelCat = lcblk(2,iblk)
        nenl = lcblk(5,iblk)
        nshl = lcblk(10,iblk)
        inum  = iel + npro - 1

        ngauss = nint(lcsyst)
        ngaussf = nintf(lcsyst)

        if (ngauss .ne. ngaussf) then
        call getstrl (yold, x,      mien(iblk)%p,
     &               strl(iel:inum,:), shgl(lcsyst,:,1:nshl,:),
     &               shp(lcsyst,1:nshl,:))
        endif

      enddo

      do iblk = 1,nelblk
         lcsyst = lcblk(3,iblk)
         iel  = lcblk(1,iblk)
         npro = lcblk(1,iblk+1) - iel
         lelCat = lcblk(2,iblk)
         inum  = iel + npro - 1

         ngauss = nint(lcsyst)

         call getnu (mien(iblk)%p, strl(iel:inum,:),
     &        mxmudmi(iblk)%p,cdelsq,ifathe(iel:inum,:))
      enddo


      return
      end



      subroutine smagcoeff (fres, xden, xnum, ien, shp, lfathe)

      include "common.h"

      dimension ien(npro,nshl),          fres(nshg,22),
     &          shp(nshl,maxsh),         fresli(npro,22),
     &          fresl(npro,nshl,22),     strnrmi(npro),
     &          xnutl(npro),             fwr(npro),
     &          xmij(npro,6),            xlij(npro,6),
     &          xdenli(npro,ngauss),      xnumli(npro,ngauss),
     &          xden(idim),              xnum(idim),
     &          lfathe(npro,maxsh)

      call local (fres, fresl, ien, 22, 'gather  ')

      do intp = 1, ngauss

         fresli = zero
         do i = 1, nenl
            do j = 1, 22
               fresli(:,j) = fresli(:,j) + shp(i,intp)*fresl(:,i,j)
            enddo
         enddo

         strnrmi(:) =  sqrt(
     &   two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2)
     &   + four * ( fresli(:,13)**2 + fresli(:,14)**2 +
     &   fresli(:,15)**2 ) )

         fwr = fwr1 * fresli(:,22) * strnrmi(:)

         xmij(:,1) = -fwr
     &        * fresli(:,10) + fresli(:,16)

         xmij(:,2) = -fwr
     &        * fresli(:,11) + fresli(:,17)

         xmij(:,3) = -fwr
     &        * fresli(:,12) + fresli(:,18)

         xmij(:,4) = -fwr * fresli(:,13) + fresli(:,19)
         xmij(:,5) = -fwr * fresli(:,14) + fresli(:,20)
         xmij(:,6) = -fwr * fresli(:,15) + fresli(:,21)

         fresli(:,22) = one/fresli(:,22)

         xlij(:,1) = fresli(:,4) - fresli(:,1) *
     &        fresli(:,1) * fresli(:,22)
         xlij(:,2) = fresli(:,5) - fresli(:,2) *
     &        fresli(:,2) * fresli(:,22)
         xlij(:,3) = fresli(:,6) - fresli(:,3) *
     &        fresli(:,3) * fresli(:,22)
         xlij(:,4) = fresli(:,7) - fresli(:,1) *
     &        fresli(:,2) * fresli(:,22)
         xlij(:,5) = fresli(:,8) - fresli(:,1) *
     &        fresli(:,3) * fresli(:,22)
         xlij(:,6) = fresli(:,9) - fresli(:,2) *
     &        fresli(:,3) * fresli(:,22)

         xnumli(:,intp) =    xlij(:,1) * xmij(:,1)
     &        + xlij(:,2) * xmij(:,2)
     &        + xlij(:,3) * xmij(:,3)
     &        + two * (xlij(:,4) * xmij(:,4)
     &        + xlij(:,5) * xmij(:,5)
     &        + xlij(:,6) * xmij(:,6))
         xdenli(:,intp) =        xmij(:,1) * xmij(:,1)
     &        + xmij(:,2) * xmij(:,2)
     &        + xmij(:,3) * xmij(:,3)
     &        + two * (xmij(:,4) * xmij(:,4)
     &        + xmij(:,5) * xmij(:,5)
     &        + xmij(:,6) * xmij(:,6))


      enddo ! End loop over quad pts

c... For debugging

c      do nel = 1, npro
c      do intp = 1, ngauss

c         if ( mod(lfathe(nel,intp),2) .eq. 0 )then
c            xnumli(nel,intp) = 2.0
c            xdenli(nel,intp) = 3.0
c         else
c            xnumli(nel,intp) = 5.0
c            xdenli(nel,intp) = 2.5
c         endif

c      enddo
c      enddo

c...  End of debugging stuff.


      do intp = 1, ngauss
      do nel  = 1, npro

         xden( lfathe(nel,intp) ) = xden( lfathe(nel,intp) ) +
     &        xdenli(nel,intp)

         xnum( lfathe(nel,intp) ) = xnum( lfathe(nel,intp) ) +
     &        xnumli(nel,intp)

      enddo
      enddo

      return

      end