1 module dtnmod 2 integer, allocatable :: ifeature(:) 3 end module 4 5 6 subroutine initDtN 7 use dtnmod 8 include "common.h" 9 allocate (ifeature(nshg)) 10 end 11 12 13 subroutine DtN(iBC,BC,y) 14 use dtnmod 15 include "common.h" 16 real*8 BC(nshg,ndofBC),y(nshg,ndof),tmp(nsclr) 17 integer iBC(nshg) 18 do i=1,nshg 19 itype=ifeature(i) 20 if(btest(iBC(i),13)) then 21 do j=1,nsclr 22 tmp(j)=y(i,5+j) 23 end do 24 call Dirichlet2Neumann(nsclr, itype, tmp) 25c 26c put the value in the position a Dirichlet value would be in BC. 27c later we will localize this value to the BCB array. 28c this is not dangerous because we should NEVER need to set Dirichlet 29c on the same node as a DtN condition 30c 31 do j=1,nsclr 32 BC(i,6+j)=-tmp(j) 33 end do 34 endif 35 end do 36 return 37 end 38 39 subroutine Dirichlet2Neumann_faux(nscalar, itype, tmp) 40c 41c This is a fake routine, designed to do nothing but feed back 42c fluxes as if there were a fast reaction at the surface and the 43c thickness of the stagnant BL were given by "distance". 44c It can handle up to 24 different scalars. 45c 46c If itype is zero, the flux is arbitrarily set to half what it would 47c be for any other itype. 48c 49c The assumption of "units" is that the initial concentrations are in 50c moles/cubic-meter and the fluxes are in moles/(sec square-meter) 51c 52c The listed diffusivities are characteristic of metal-ions in room 53c temperature water, in units of square-meter/sec. 54c 55c 56 integer itype, nscalar 57 real*8 tmp(nscalar) 58 59 integer i 60 real*8 distance 61 real*8 units 62 63c Completely fake diffusivities 64 real*8 D(24) 65 data D/ 66 & 5.0d-05, 1.0d-5, 8.0d-10, 7.0d-10, 6.0d-10, 5.0d-10, 67 & 4.0d-10, 3.0d-10, 2.0d-10, 1.0d-10, 0.9d-10, 0.8d-10, 68 & 1.0d-10, 9.0d-10, 8.0d-10, 7.0d-10, 6.0d-10, 5.0d-10, 69 & 4.0d-10, 3.0d-10, 2.0d-10, 1.0d-10, 0.9d-10, 0.8d-10/ 70 71 do i=1,nscalar 72 tmp(i) = 0.0d0 73 enddo 74 return 75 distance = 10.0d-03 76 units = 1.0d-3 77 units = 1.0 78 if(nscalar.gt.24) then 79 write(*,*) 'Problem in Dir2Neu: nscalar larger than 24!' 80 stop 81 endif 82 83 do i=1,nscalar 84 tmp(i) = D(i) * ( tmp(i) - 0.0 ) / distance * units 85 if(itype.eq.2) tmp(i) = tmp(i) / 2.0d+00 86 enddo 87c tmp(1)=1.0d-5 88 return 89 end 90 91 92 93 94 95 96 97 98 99 subroutine dtnl(iBC,BC,ienb,iBCB,BCB) 100 include "common.h" 101 integer ienb(npro,nshl), iBC(nshg),iBCB(npro,ndiBCB) 102 real*8 BCB(npro,nshlb,ndBCB), tmpBCB(npro,nshlb,nsclr), 103 & BC(nshg,ndofBC), tmpBC(nshg,nsclr) 104 105 nstart=ndofBC-nsclr 106c tmpBC=zero 107c do i=1,nshg 108c if(btest(iBC(i),13)) then 109 do j=1,nsclr 110c tmpBC(i,j)=BC(i,nstart+j) 111 tmpBC(:,j)=BC(:,nstart+j) 112 enddo 113c endif 114c enddo 115 116 call localb(tmpBC,tmpBCB,ienb,nsclr,'gather ') 117 118 do i=1,npro 119 do j=1,nsclr 120 if(iBCB(i,2).lt.0) then !this is a face with dtn 121 do k=1,nshlb 122 BCB(i,k,6+j)=tmpBCB(i,k,j) 123 enddo 124 endif 125 enddo 126 enddo 127 return 128 end 129 130c 131c This routine just calls the appropriate version of D2N for the number 132c of scalars used 133c 134 subroutine Dirichlet2Neumann(nscalar, itype, tmp) 135 integer nscalar, itype 136 real*8 tmp(nscalar),foo 137 138c Just short circuit the routine for a little bit. 139c tmp(1)=0.0d0 140c return 141 if(nscalar .eq. 1) then 142c write(*,*) 'Entering D2N1' 143c foo= rand(0) 144 call Dirichlet2Neumann_1(nscalar,itype,tmp) 145c write(*,*) 'Returning from D2N after DTN1' 146c return 147 elseif(nscalar.eq.2) then 148 call Dirichlet2Neumann_2(nscalar,itype,tmp) 149 else 150 write(*,*) 'FATAL ERROR: cannont handle ',nscalar,' scalars' 151 stop 152 endif 153 154 return 155 end 156 157 subroutine Dirichlet2Neumann_2(nscalar, itype, tmp) 158c 159c This is an interface routine, designed to call return a value for 160c the flux to a point on the wafer due to electrochemical deposition 161c to Ken Jansen's PHASTA given a boundary conditions and an index for 162c a particular feature. 163c 164c There is an inherent assumption that we are going to be doing 165c electroplating. This routine sets up the filenames and the 166c top-of-the-domain boundary conditions. 167c 168 implicit none 169 170 integer maxdata,maxtypes 171 parameter(maxdata=100,maxtypes=5) 172 173 integer itype, nscalar 174 real*8 tmp(nscalar) 175c For each table up to maxtypes, we have 4 pieces of data--two independent, 176c two dependent--for each point, up to maxdata+1. 177 real*8 table(4,0:maxdata,0:maxdata,maxtypes) 178 save table 179 180 integer i,j,n 181 logical readfile(maxtypes) 182 save readfile 183 data (readfile(i),i=1,maxtypes) / maxtypes*.false./ 184 185 real*8 dx(2,maxtypes) 186 integer numdata(2,maxtypes) 187 save dx 188 save numdata 189 190 real*8 x,y, z(3,2) 191c We can only deal with two parameter models for now. 192 if(nscalar .ne. 2) then 193 write(*,*) 'Sorry, Dirichlet2Neumann handles 2 scalars!' 194 write(*,*) 'You asked for ', nscalar 195 write(*,*) 'STOPPING...' 196 stop 197 endif 198 199c If we haven't read in our parameters for this featuretype yet... 200 201 if( .not. readfile(itype)) then 202 readfile(itype) = .true. 203 call readtable_2(itype,table,numdata,dx, 204 & maxdata,maxtypes) 205 endif 206 207 x = tmp(1) 208 y = tmp(2) 209 210 if(.false.) then 211 if( x .gt. table(1,0,0,itype) .or. 212 & x .lt. table(1,numdata(1,itype)-1,0,itype) ) then 213 write(*,*) 'Sorry, concentration 1 asked for: ', x 214 write(*,*) ' is out of the table bounds.' 215 write(*,*) '#1 [ ',table(1,0,0,itype), ' , ', 216 & table(1,numdata(1,itype)-1,0,itype), ' ] ', 217 & numdata(1,itype)-1 218 219 write(*,*) ' STOPPING...' 220 stop 221 endif 222 if( y .gt. table(2,0,0,itype) .or. 223 & y .lt. table(2,0,numdata(2,itype)-1,itype) ) then 224 write(*,*) 'Sorry, concentration 2 asked for: ', y 225 write(*,*) ' is out of the table bounds.' 226 write(*,*) '#2 [ ',table(2,0,0,itype), ' , ', 227 & table(2,0,numdata(2,itype)-1,itype), ' ] ', 228 & numdata(2,itype)-1 229 write(*,*) ' STOPPING...' 230 stop 231 endif 232 endif 233 234 i = int ( (x - table(1,0,0,itype) ) / dx(1,itype)) 235 j = int ( (y - table(2,0,0,itype) ) / dx(2,itype)) 236c write(*,*) 'i,j,x,y: ',i,j,x,y 237 if(i .lt. 0) then 238 i = 0 239c x = table(1,0,0,itype) 240c write(*,*) 'Reseting i low: ',i,j,x,y 241 endif 242 if(j .lt. 0) then 243 j = 0 244 y = table(2,0,0,itype) 245c write(*,*) 'Reseting j low: ',i,j,x,y 246 endif 247 if(i .ge. numdata(1,itype)) then 248 i = numdata(1,itype)-2 249c x = table(1,i+1,0,itype) 250c write(*,*) 'Reseting i high: ',i,j,x,y 251 endif 252 if(j .ge. numdata(2,itype)) then 253 j = numdata(2,itype)-2 254 y = table(1,0,j+1,itype) 255c write(*,*) 'Reseting j high: ',i,j,x,y 256 endif 257 258 do n=3,4 259 260 z(1,1) = table(n,i,j,itype) 261 z(3,1) = table(n,i+1,j,itype) 262 z(1,2) = table(n,i,j+1,itype) 263 z(3,2) = table(n,i+1,j+1,itype) 264 265 z(2,1) = (z(3,1) - z(1,1)) / dx(1,itype) 266 & *(x-table(1,i,j,itype)) + z(1,1) 267 z(2,2) = (z(3,2) - z(1,2)) / dx(1,itype) 268 & *(x-table(1,i,j,itype)) + z(1,2) 269 tmp(n-2) = (z(2,2) - z(2,1))/dx(2,itype) 270 & *(y-table(2,i,j,itype)) + z(2,1) 271 272 enddo 273c write(*,*) 'Interpolation from ',x,y,' to:', tmp(1),tmp(2) 274 return 275 276 end 277c-------------------------------------------------------------------- 278 279 subroutine readtable_2(islot,table,numdata,dx,maxdata,maxslots) 280c 281c Read a table of ordered quadruplets and place them into the slot in 282c TABLE that is assosciated with ISLOT. Store the number of 283c data in NUMDATA and the spacing in DX. The file to be read 284c is 'TABLE.[islot]' The data but be in a rectangular regular grid. 285c 286c AUTHOR: Max Bloomfield, May 2000 287c 288 implicit none 289c 290 integer islot 291 integer maxslots,numdata(2,maxslots), maxdata 292c 293 real*8 table(4,0:maxdata,0:maxdata,maxslots), dx(2,maxslots) 294 real*8 x1,x2,y1,y2,x2old 295c 296 character*250 linein,filename 297c 298 integer i,j,k 299 logical verbose 300 301 verbose = .true. 302 303 i=0 304 j=0 305 write(filename,1066) islot 306 1066 format('TABLE.',i1) 307 308 open(file=filename,unit=1066) 309 310 write(*,*) 'Opening ', filename 311 312 1 continue 313 read (unit=1066,fmt='(a)',end=999) linein 314 315 if (linein(1:1).eq.'#') then 316 write (*,'(a)') linein 317 goto 1 318 endif 319c 320 if (i.gt.maxdata**2+maxdata+1) then 321 write(*,*) 322 & 'reached the maximum number of data points allowed' 323 write(*,*) 'FATAL ERROR: stopping' 324 stop 325 endif 326c 327 read (linein,*) x1,x2,y1,y2 328 if(i .gt. 0 .and. x2 .ne. x2old) then 329c increment the outer index in this nested loop. That is, go on 330c to the next "row" (not in fortran speak, but in table speak.) 331 j = j + 1 332 i=0 333 endif 334 335 table(1,i,j,islot) = x1*1.d-0 336 table(2,i,j,islot) = x2*1.d-0 337 table(3,i,j,islot) = y1*1.d-0 338 table(4,i,j,islot) = y2*1.d-0 339c 340 i=i+1 341 x2old = x2 342 343 goto 1 344c 345 999 continue 346c 347 numdata(1,islot) = I 348 numdata(2,islot) = j+1 349c 350 dx(1,islot) = table(1,2,1,islot) - table(1,1,1,islot) 351 dx(2,islot) = table(2,1,2,islot) - table(2,1,1,islot) 352 353 if(verbose) then 354 write(*,*) 'Table is ',i,' by ',j+1 355 356 write(*,*) 'there are ',i*(j+1),' flux data points' 357 write(*,*) 'closing unit 1066' 358 close(1066) 359c 360 write(*,*) 'The abscissa are ', 361 & dx(1,islot),' and ',dx(2,islot),' apart.' 362 363 write(*,*) 'the flux data are ' 364 do i=0,numdata(1,islot)-1 365 do j=0,numdata(2,islot)-1 366 write(*,*) i,j,(table(k,i,j,islot), k=1,4) 367 end do 368 end do 369 370 endif 371 return 372 end 373 374 375 376 subroutine Dirichlet2Neumann_1(nscalar, itype, tmp) 377c 378c This is an interface routine, designed to call return a value for 379c the flux to a point on the wafer due to electrochemical deposition 380c to Ken Jansen's PHASTA given a boundary conditions and an index for 381c a particular feature. 382c 383c There is an inherent assumption that we are going to be doing 384c electroplating. This routine sets up the filenames and the 385c top-of-the-domain boundary conditions. 386c 387 implicit none 388 389 integer maxdata,maxtypes 390 parameter(maxdata=200,maxtypes=2) 391 392 integer itype, nscalar 393 real*8 tmp(nscalar) 394 real*8 table(0:maxdata,2,maxtypes) 395 save table 396 397 integer i 398 logical readfile(maxtypes) 399 save readfile 400 data (readfile(i),i=1,maxtypes) / maxtypes*.false./ 401 402 real*8 dx(maxtypes) 403 save dx 404 integer numdata(maxtypes) 405 save numdata 406 407 real*8 dt, conc_BC, flux_BC 408c We can only deal with one parameter models for now. 409 410 if(nscalar .ne. 1) then 411 write(*,*) 'Sorry, Dirichlet2Neumann can only handle 1 scalar!' 412 write(*,*) 'You asked for ', nscalar 413 write(*,*) 'STOPPING...' 414 stop 415 endif 416 417c If we haven't read in our parameters for this featuretype yet... 418 419 if( .not. readfile(itype)) then 420 readfile(itype) = .true. 421 call readtable_1(itype,table,numdata(itype),dx(itype), 422 & maxdata,maxtypes) 423c write(*,*) 'back from readtable' 424 if(dx(itype) .eq. 0.0d0) then 425 write(*,*) 'DX for table ',itype,' is zero (I think!)' 426 stop 427 endif 428 endif 429c write(*,*) 'returning from D2N' 430 431 conc_BC = tmp(1) 432 433c if( conc_BC .lt. table(0,1,itype) .or. 434c & conc_BC .gt. table(numdata(itype),1,itype) ) then 435c write(*,*) 'Sorry, concentration asked for: ', conc_BC 436c write(*,*) ' is out of the table bounds.' 437c write(*,*) '[',table(0,1,itype),', 438c & ',table(numdata(itype),1,itype),']' 439c write(*,*) ' STOPPING...' 440c stop 441c endif 442 443 i = int ( (conc_BC - table(0,1,itype) ) / dx(itype)) 444 445 if( conc_BC .lt. table(0,1,itype))then 446 i = 0 447 conc_BC = table(i,1,itype) 448 elseif( conc_BC .gt. table(numdata(itype),1,itype)) then 449 i = numdata(itype) 450 conc_BC = table(i,1,itype) 451 endif 452 453 454 dt = conc_BC - table(i,1,itype) 455 flux_BC = dt * (table(i+1,2,itype) - table(i,2,itype)) + 456 & table(i,2,itype) 457 458 459 tmp(1) = flux_BC 460 461 462 end 463c-------------------------------------------------------------------- 464 465 subroutine readtable_1(islot,table,numdata,dx,maxdata,maxslots) 466c 467c Read a table of ordered pairs and place them into the slot in 468c TABLE that is assosciated with ISLOT. Store the number of 469c data in NUMDATA and the spacing in DX. The file to be read 470c is 'TABLE.[islot]' 471c 472c AUTHOR: Max Bloomfield, May 2000 473c 474 implicit none 475c 476 integer islot 477 integer numdata, maxdata, maxslots 478c 479 real*8 table(0:maxdata,2,maxslots),dx 480c 481 character*80 linein,filename 482c 483 integer i,j 484 logical verbose 485 verbose = .true. 486 487 i=-1 488 489 write(filename,1066) islot 490 1066 format('TABLE.',i1) 491 open(file=filename,unit=1066) 492 if(verbose) write(*,*) 'Opening ', filename 493 494 1 continue 495 read (unit=1066,fmt='(a)',end=999) linein 496 497 if (linein(1:1).eq.'#') then 498 write (*,'(a)') linein 499 goto 1 500 endif 501c 502 i=i+1 503 if (i.ge.maxdata) then 504 write(*,*) 505 & 'reached the maximum number of data points allowed' 506 write(*,*) 'FATAL ERROR: stopping' 507 stop 508 endif 509c 510 read (linein,*) table(i,1,islot), table(i,2,islot) 511 table(i,1,islot)= table(i,1,islot)*1.0d-0 512 table(i,2,islot)= table(i,2,islot)*1.0d-0 513c 514 goto 1 515c 516 999 continue 517c 518 numdata = i 519 dx = table(1,1,islot)-table(0,1,islot) 520c 521 if(verbose) then 522 write(*,*) 'there are ',numdata,' flux data points' 523 write(*,*) 'closing unit 1066' 524 close(1066) 525c 526 write(*,*) 'the flux data are ' 527 do 101 j=0,i 528 write(*,*) j,table(j,1,islot), table(j,2,islot) 529 101 continue 530 endif 531 return 532 end 533