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