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