xref: /phasta/phSolver/common/dtn.f (revision 6b966dd8220c388f9c18e3b5a7cb164734f542c6)

      module dtnmod
      integer, allocatable :: ifeature(:)
      end module


      subroutine initDtN
      use dtnmod
      include "common.h"
      allocate (ifeature(nshg))
      end

      subroutine finalizeDtN
      use dtnmod
      include "common.h"
      deallocate(ifeature)
      end

      subroutine DtN(iBC,BC,y)
      use dtnmod
      include "common.h"
      real*8 BC(nshg,ndofBC),y(nshg,ndof),tmp(nsclr)
      integer iBC(nshg)
      do i=1,nshg
         itype=ifeature(i)
         if(btest(iBC(i),13)) then
            do j=1,nsclr
               tmp(j)=y(i,5+j)
            end do
            call Dirichlet2Neumann(nsclr, itype, tmp)
c
c  put the value in the position a Dirichlet value would be in BC.
c  later we will localize this value to the BCB array.
c  this is not dangerous because we should NEVER need to set Dirichlet
c  on the same node as a DtN condition
c
            do j=1,nsclr
               BC(i,6+j)=-tmp(j)
            end do
         endif
      end do
      return
      end

      subroutine Dirichlet2Neumann_faux(nscalar, itype, tmp)
c
c This is a fake routine, designed to do nothing but feed back
c fluxes as if there were a fast reaction at the surface and the
c thickness of the stagnant BL were given by "distance".
c It can handle up to 24 different scalars.
c
c If itype is zero, the flux is arbitrarily set to half what it would
c be for any other itype.
c
c The assumption of "units" is that the initial concentrations are in
c moles/cubic-meter and the fluxes are in moles/(sec square-meter)
c
c The listed diffusivities are characteristic of metal-ions in room
c temperature water, in units of square-meter/sec.
c
c
      integer itype, nscalar
      real*8 tmp(nscalar)

      integer i
      real*8 distance
      real*8 units

c  Completely fake diffusivities
      real*8 D(24)
      data D/
     &       5.0d-05, 1.0d-5, 8.0d-10, 7.0d-10, 6.0d-10, 5.0d-10,
     &       4.0d-10, 3.0d-10, 2.0d-10, 1.0d-10, 0.9d-10, 0.8d-10,
     &       1.0d-10, 9.0d-10, 8.0d-10, 7.0d-10, 6.0d-10, 5.0d-10,
     &       4.0d-10, 3.0d-10, 2.0d-10, 1.0d-10, 0.9d-10, 0.8d-10/

      do i=1,nscalar
         tmp(i) = 0.0d0
      enddo
      return
      distance = 10.0d-03
      units = 1.0d-3
      units = 1.0
      if(nscalar.gt.24) then
         write(*,*) 'Problem in Dir2Neu: nscalar larger than 24!'
         stop
      endif

      do i=1,nscalar
         tmp(i) = D(i) * ( tmp(i) - 0.0 ) / distance  * units
         if(itype.eq.2) tmp(i) = tmp(i) / 2.0d+00
      enddo
c      tmp(1)=1.0d-5
      return
      end









      subroutine dtnl(iBC,BC,ienb,iBCB,BCB)
      include "common.h"
      integer ienb(npro,nshl), iBC(nshg),iBCB(npro,ndiBCB)
      real*8  BCB(npro,nshlb,ndBCB), tmpBCB(npro,nshlb,nsclr),
     &        BC(nshg,ndofBC),        tmpBC(nshg,nsclr)

      nstart=ndofBC-nsclr
c      tmpBC=zero
c      do i=1,nshg
c         if(btest(iBC(i),13)) then
            do j=1,nsclr
c               tmpBC(i,j)=BC(i,nstart+j)
               tmpBC(:,j)=BC(:,nstart+j)
            enddo
c         endif
c      enddo

      call localb(tmpBC,tmpBCB,ienb,nsclr,'gather  ')

      do i=1,npro
         do j=1,nsclr
            if(iBCB(i,2).lt.0) then  !this is a face with dtn
               do k=1,nshlb
                  BCB(i,k,6+j)=tmpBCB(i,k,j)
               enddo
            endif
         enddo
      enddo
      return
      end

c
c This routine just calls the appropriate version of D2N for the number
c of scalars used
c
      subroutine Dirichlet2Neumann(nscalar, itype, tmp)
      integer nscalar, itype
      real*8 tmp(nscalar),foo

c Just short circuit the routine for a little bit.
c      tmp(1)=0.0d0
c      return
      if(nscalar .eq. 1) then
c         write(*,*) 'Entering D2N1'
c          foo= rand(0)
         call Dirichlet2Neumann_1(nscalar,itype,tmp)
c         write(*,*) 'Returning from D2N after DTN1'
c         return
      elseif(nscalar.eq.2) then
         call Dirichlet2Neumann_2(nscalar,itype,tmp)
      else
         write(*,*) 'FATAL ERROR: cannont handle ',nscalar,' scalars'
         stop
      endif

      return
      end

      subroutine Dirichlet2Neumann_2(nscalar, itype, tmp)
c
c This is an interface routine, designed to call return a value for
c the flux to a point on the wafer due to electrochemical deposition
c to Ken Jansen's PHASTA given a boundary conditions and an index for
c a particular feature.
c
c There is an inherent assumption that we are going to be doing
c electroplating. This routine sets up the filenames and the
c top-of-the-domain boundary conditions.
c
      implicit none

      integer maxdata,maxtypes
      parameter(maxdata=100,maxtypes=5)

      integer itype, nscalar
      real*8 tmp(nscalar)
c For each table up to maxtypes, we have 4 pieces of data--two independent,
c two dependent--for each point, up to maxdata+1.
      real*8 table(4,0:maxdata,0:maxdata,maxtypes)
      save table

      integer i,j,n
      logical readfile(maxtypes)
      save readfile
      data (readfile(i),i=1,maxtypes) / maxtypes*.false./

      real*8  dx(2,maxtypes)
      integer numdata(2,maxtypes)
      save dx
      save numdata

      real*8  x,y, z(3,2)
c We can only deal with two parameter models for now.
      if(nscalar .ne. 2) then
         write(*,*) 'Sorry, Dirichlet2Neumann handles 2 scalars!'
         write(*,*) 'You asked for ', nscalar
         write(*,*) 'STOPPING...'
         stop
      endif

c If we haven't read in our parameters for this featuretype yet...

      if( .not. readfile(itype)) then
         readfile(itype) = .true.
         call readtable_2(itype,table,numdata,dx,
     &        maxdata,maxtypes)
      endif

      x = tmp(1)
      y = tmp(2)

      if(.false.) then
         if( x .gt. table(1,0,0,itype) .or.
     &        x .lt. table(1,numdata(1,itype)-1,0,itype) ) then
            write(*,*) 'Sorry, concentration 1 asked for: ', x
            write(*,*) '  is out of the table bounds.'
            write(*,*)  '#1  [ ',table(1,0,0,itype), ' , ',
     &           table(1,numdata(1,itype)-1,0,itype), ' ] ',
     &           numdata(1,itype)-1

            write(*,*) '  STOPPING...'
            stop
         endif
         if( y .gt. table(2,0,0,itype) .or.
     &        y .lt. table(2,0,numdata(2,itype)-1,itype) ) then
            write(*,*) 'Sorry, concentration 2 asked for: ', y
            write(*,*) '  is out of the table bounds.'
            write(*,*)  '#2   [ ',table(2,0,0,itype), ' , ',
     &           table(2,0,numdata(2,itype)-1,itype), ' ] ',
     &           numdata(2,itype)-1
            write(*,*) '  STOPPING...'
            stop
         endif
      endif

      i = int ( (x - table(1,0,0,itype) ) / dx(1,itype))
      j = int ( (y - table(2,0,0,itype) ) / dx(2,itype))
c      write(*,*) 'i,j,x,y: ',i,j,x,y
      if(i .lt. 0) then
         i = 0
c         x = table(1,0,0,itype)
c         write(*,*) 'Reseting i low: ',i,j,x,y
      endif
      if(j .lt. 0) then
         j = 0
         y = table(2,0,0,itype)
c         write(*,*) 'Reseting j low: ',i,j,x,y
      endif
      if(i .ge. numdata(1,itype)) then
         i = numdata(1,itype)-2
c         x = table(1,i+1,0,itype)
c         write(*,*) 'Reseting i high: ',i,j,x,y
      endif
      if(j .ge. numdata(2,itype)) then
         j = numdata(2,itype)-2
         y = table(1,0,j+1,itype)
c         write(*,*) 'Reseting j high: ',i,j,x,y
      endif

      do n=3,4

         z(1,1) = table(n,i,j,itype)
         z(3,1) = table(n,i+1,j,itype)
         z(1,2) = table(n,i,j+1,itype)
         z(3,2) = table(n,i+1,j+1,itype)

         z(2,1) = (z(3,1) - z(1,1)) / dx(1,itype)
     &        *(x-table(1,i,j,itype)) + z(1,1)
         z(2,2) = (z(3,2) - z(1,2)) / dx(1,itype)
     &        *(x-table(1,i,j,itype)) + z(1,2)
         tmp(n-2) = (z(2,2) - z(2,1))/dx(2,itype)
     &        *(y-table(2,i,j,itype)) + z(2,1)

      enddo
c      write(*,*) 'Interpolation from ',x,y,' to:', tmp(1),tmp(2)
      return

      end
c--------------------------------------------------------------------

      subroutine readtable_2(islot,table,numdata,dx,maxdata,maxslots)
c
c    Read a table of ordered quadruplets and place them into the slot in
c    TABLE that is assosciated with ISLOT. Store the number of
c    data in NUMDATA and the spacing in DX.  The file to be read
c    is 'TABLE.[islot]' The data but be in a rectangular regular grid.
c
c     AUTHOR: Max Bloomfield, May 2000
c
      implicit none
c
      integer islot
      integer maxslots,numdata(2,maxslots), maxdata
c
      real*8 table(4,0:maxdata,0:maxdata,maxslots), dx(2,maxslots)
      real*8 x1,x2,y1,y2,x2old
c
      character*250 linein,filename
c
      integer i,j,k
      logical verbose

      verbose = .true.

      i=0
      j=0
      write(filename,1066) islot
 1066 format('TABLE.',i1)

      open(file=filename,unit=1066)

      write(*,*) 'Opening ', filename

 1    continue
      read (unit=1066,fmt='(a)',end=999) linein

      if (linein(1:1).eq.'#') then
         write (*,'(a)') linein
         goto 1
      endif
c
      if (i.gt.maxdata**2+maxdata+1) then
         write(*,*)
     &        'reached the maximum number of data points allowed'
         write(*,*) 'FATAL ERROR: stopping'
         stop
      endif
c
      read (linein,*) x1,x2,y1,y2
      if(i .gt. 0 .and. x2 .ne. x2old) then
c        increment the outer index in this nested loop. That is, go on
c        to the next "row" (not in fortran speak, but in table speak.)
         j = j + 1
         i=0
      endif

      table(1,i,j,islot) = x1*1.d-0
      table(2,i,j,islot) = x2*1.d-0
      table(3,i,j,islot) = y1*1.d-0
      table(4,i,j,islot) = y2*1.d-0
c
      i=i+1
      x2old = x2

      goto 1
c
 999  continue
c
      numdata(1,islot) = I
      numdata(2,islot) = j+1
c
      dx(1,islot) = table(1,2,1,islot) - table(1,1,1,islot)
      dx(2,islot) = table(2,1,2,islot) - table(2,1,1,islot)

      if(verbose) then
         write(*,*) 'Table is ',i,' by ',j+1

         write(*,*) 'there are ',i*(j+1),' flux data points'
         write(*,*) 'closing unit 1066'
         close(1066)
c
         write(*,*) 'The abscissa are ',
     &        dx(1,islot),' and ',dx(2,islot),' apart.'

         write(*,*) 'the flux data are '
         do i=0,numdata(1,islot)-1
            do j=0,numdata(2,islot)-1
               write(*,*) i,j,(table(k,i,j,islot), k=1,4)
            end do
         end do

      endif
      return
      end



      subroutine Dirichlet2Neumann_1(nscalar, itype, tmp)
c
c This is an interface routine, designed to call return a value for
c the flux to a point on the wafer due to electrochemical deposition
c to Ken Jansen's PHASTA given a boundary conditions and an index for
c a particular feature.
c
c There is an inherent assumption that we are going to be doing
c electroplating. This routine sets up the filenames and the
c top-of-the-domain boundary conditions.
c
      implicit none

      integer maxdata,maxtypes
      parameter(maxdata=200,maxtypes=2)

      integer itype, nscalar
      real*8 tmp(nscalar)
      real*8 table(0:maxdata,2,maxtypes)
      save table

      integer i
      logical readfile(maxtypes)
      save readfile
      data (readfile(i),i=1,maxtypes) / maxtypes*.false./

      real*8  dx(maxtypes)
      save dx
      integer numdata(maxtypes)
      save numdata

      real*8 dt, conc_BC, flux_BC
c We can only deal with one parameter models for now.

      if(nscalar .ne. 1) then
         write(*,*) 'Sorry, Dirichlet2Neumann can only handle 1 scalar!'
         write(*,*) 'You asked for ', nscalar
         write(*,*) 'STOPPING...'
         stop
      endif

c If we haven't read in our parameters for this featuretype yet...

      if( .not. readfile(itype)) then
         readfile(itype) = .true.
         call readtable_1(itype,table,numdata(itype),dx(itype),
     &        maxdata,maxtypes)
c         write(*,*) 'back from readtable'
         if(dx(itype) .eq. 0.0d0) then
            write(*,*) 'DX for table ',itype,' is zero (I think!)'
            stop
         endif
      endif
c      write(*,*) 'returning from D2N'

      conc_BC = tmp(1)

c      if( conc_BC .lt. table(0,1,itype) .or.
c     &    conc_BC .gt. table(numdata(itype),1,itype) ) then
c         write(*,*) 'Sorry, concentration asked for: ', conc_BC
c         write(*,*) '  is out of the table bounds.'
c         write(*,*) '[',table(0,1,itype),',
c     &        ',table(numdata(itype),1,itype),']'
c         write(*,*) '  STOPPING...'
c         stop
c      endif

      i = int ( (conc_BC - table(0,1,itype) ) / dx(itype))

      if( conc_BC .lt. table(0,1,itype))then
         i = 0
         conc_BC =  table(i,1,itype)
      elseif( conc_BC .gt. table(numdata(itype),1,itype)) then
         i = numdata(itype)
         conc_BC =  table(i,1,itype)
      endif


      dt = conc_BC - table(i,1,itype)
      flux_BC = dt * (table(i+1,2,itype) - table(i,2,itype)) +
     &     table(i,2,itype)


      tmp(1) = flux_BC


      end
c--------------------------------------------------------------------

      subroutine readtable_1(islot,table,numdata,dx,maxdata,maxslots)
c
c     Read a table of ordered pairs and place them into the slot in
c     TABLE that is assosciated with ISLOT. Store the number of
c     data in NUMDATA and the spacing in DX.  The file to be read
c     is 'TABLE.[islot]'
c
c     AUTHOR: Max Bloomfield, May 2000
c
      implicit none
c
      integer islot
      integer numdata, maxdata, maxslots
c
      real*8 table(0:maxdata,2,maxslots),dx
c
      character*80 linein,filename
c
      integer i,j
      logical verbose
      verbose = .true.

      i=-1

      write(filename,1066) islot
 1066 format('TABLE.',i1)
      open(file=filename,unit=1066)
      if(verbose) write(*,*) 'Opening ', filename

 1    continue
      read (unit=1066,fmt='(a)',end=999) linein

      if (linein(1:1).eq.'#') then
         write (*,'(a)') linein
         goto 1
      endif
c
      i=i+1
      if (i.ge.maxdata) then
         write(*,*)
     &        'reached the maximum number of data points allowed'
         write(*,*) 'FATAL ERROR: stopping'
         stop
      endif
c
      read (linein,*) table(i,1,islot), table(i,2,islot)
      table(i,1,islot)= table(i,1,islot)*1.0d-0
      table(i,2,islot)= table(i,2,islot)*1.0d-0
c
      goto 1
c
 999  continue
c
      numdata = i
      dx = table(1,1,islot)-table(0,1,islot)
c
      if(verbose) then
         write(*,*) 'there are ',numdata,' flux data points'
         write(*,*) 'closing unit 1066'
         close(1066)
c
         write(*,*) 'the flux data are '
         do 101 j=0,i
            write(*,*) j,table(j,1,islot), table(j,2,islot)
 101     continue
      endif
      return
      end