1c----------------------------------------------------------------------- 2c 3c This module conveys temporal BC data. Below functions read in the data 4c and interpolate it to the current time level. 5c 6c----------------------------------------------------------------------- 7 module specialBC 8 use pointer_data 9 10 real*8, allocatable :: BCt(:,:), acs(:,:), spamp(:) 11 real*8, allocatable :: ytarget(:,:) 12 real*8, allocatable :: BCperiod(:) 13 integer, allocatable :: nBCt(:), numBCt(:) 14 15 type(r2d),dimension(:),allocatable :: BCtptr 16 ! BCtptr%p(:,:) is to replace BCt(ic,:,:), in which 17 ! the second index is dynamically decided in 18 ! the setup stage. 19 20 integer ntv,nptsmax 21c$$$ integer itvn 22 end module 23 24c----------------------------------------------------------------------- 25c 26c This module conveys flow rate history for the different impedance outlets 27c over one period. Below functions read in the data and store it for the 28c current time level. 29c 30c----------------------------------------------------------------------- 31 module convolImpFlow 32 33 real*8, allocatable :: QHistImp(:,:), ValueImpt(:,:,:) 34 real*8, allocatable :: ValueListImp(:,:), ConvCoef(:,:) 35 real*8, allocatable :: ImpConvCoef(:,:), poldImp(:) 36 integer ntimeptpT,numDataImp 37 integer, allocatable :: nImpt(:), numImpt(:) 38 integer nptsImpmax 39 real*8, allocatable :: QHistTry(:,:), QHistTryF(:,:) !filter 40 integer cutfreq !filter 41 end module 42c----------------------------------------------------------------------- 43c 44c This module conveys the parameters for the different RCR outlets. 45c Below functions read in the inputs (proximal resistance, capacitance, 46c distal resistance and distal pressure) and store it for the 47c current time level. 48c 49c----------------------------------------------------------------------- 50 module convolRCRFlow 51 52 real*8, allocatable :: ValueListRCR(:,:), ValuePdist(:,:,:) !inputs 53 real*8, allocatable :: QHistRCR(:,:), HopRCR(:) !calc 54 real*8, allocatable :: RCRConvCoef(:,:), poldRCR(:) !calc 55 real*8, allocatable :: dtRCR(:) !scaled timestep: deltat/RdC 56 real*8, allocatable :: RCRic(:) !(P(0)-RQ(0)-Pd(0)) 57 integer nptsRCRmax,numDataRCR !to read inputs 58 integer, allocatable :: numRCRt(:) !to read inputs 59 end module 60 61c----------------------------------------------------------------------- 62c 63c Initialize: 64c 65c----------------------------------------------------------------------- 66 subroutine initSponge( y,x) 67 68 use specialBC 69 include "common.h" 70 71 real*8 y(nshg,nflow), x(numnp,3) 72 allocate (ytarget(nshg,nflow)) 73 74 if(matflg(5,1).eq.5) then 75 write(*,*) 'calculating IC sponge' 76 ytarget = y 77 else 78 write(*,*) 'calculating Analytic sponge' 79 80c 81c OLD style sponge pushed onto target. You need to be sure that your 82c solver.inp entries for start and stop of sponge match as well as the 83c growth rates 84c 85 vcl=datmat(1,5,1) ! velocity on centerline 86 rslc=datmat(2,5,1) ! shear layer center radius 87 bfz=datmat(3,5,1) 88 we=3.0*29./682. 89 rsteep=3.0 90 zstart=30.0 91 radst=10.0 92 radsts=radst*radst 93 do id=1,numnp 94 radsqr=x(id,2)**2+x(id,1)**2 95c if((x(id,3).gt. zstart) .or. (radsqr.gt.radsts)) then 96 rad=sqrt(radsqr) 97 radc=max(rad,radst) 98 zval=max(x(id,3),zstart) 99 utarget=(tanh(rsteep*(rslc-rad))+one)/two* 100 & (vcl-we) + we 101 Ttarget = press/(ro*Rgas) 102 ptarget= press 103 ytarget(id,1) = zero 104 ytarget(id,2) = zero 105 ytarget(id,3) = utarget 106 ytarget(id,4) = ptarget 107 ytarget(id,5) = Ttarget 108c endif 109 enddo 110 endif 111 return 112 end 113 114 115c----------------------------------------------------------------------- 116c 117c Initialize:time varying boundary condition 118c 119c----------------------------------------------------------------------- 120 subroutine initBCt( x, iBC, BC ) 121 122 use specialBC 123 include "common.h" 124 125 real*8 x(numnp,nsd), BC(nshg,ndofBC), rj1,rj2,rj3,rj4,distd,epsd 126 integer iBC(nshg) 127 character*80 card 128 real*8 distds 129 real*8 dd 130 131 integer irstart 132 real(kind=8),allocatable,dimension(:) :: bcttest 133 134 real*8 t0,tlen,t1,tstart,tend 135 integer i0,ilen,i1,nper,istart,iend 136c 137c This one should be used for boundary layer meshes where bct.dat must 138c be given to greater precision than is currently being generated. 139c 140c epsd=1.0d-12 ! this is distance SQUARED to save square root 141 142 epsd=1.0d-8 ! this is distance SQUARED to save square root 143 144 ic=0 !count the number on this processor 145 146! ************ Chun Sun: limit the memory use if required time steps 147! ************ is smaller than total bct.dat length. 148 149 ! readin the start timestep 150 151 open(unit=72,file='numstart.dat',status='old') 152 read(72,*) irstart 153 close(72) 154 155 ! use irstart-1 and nstep+1 to avoid machine tolerance issues 156 t0 = max(zero,(irstart-1)*Delt(1)) 157 tlen = (nstep(1)+1)*Delt(1) 158 ! Assumption: constant time step in one run. If you use variable 159 ! time step in one run, this should be modified. 160 t1 = t0 + tlen 161 162 if(myrank.eq.master) 163 & write(*,*) 'necessary bct timing: from ',t0,' to ',t1 164 165! ************************************************************** 166 167 168 if(any(ibits(iBC,3,3).eq.7)) then 169 if(myrank.eq.master) write(*,*) 'opening bct.dat' 170c open(unit=567, file='bct.dat',status='old') 171 open(unit=567, file='bct.dat',ACTION='READ',STATUS='old') 172 read (567,'(a80)') card 173 read (card,*) ntv, nptsmax 174c read(567,*) ntv,nptsmax 175 allocate (nBCt(numnp)) 176 allocate (numBCt(ntv)) 177 allocate (BCt(nptsmax,4)) 178 allocate (BCperiod(ntv)) 179 allocate (BCtptr(ntv)) ! dynamic bct allocation 180 do k=1,ntv 181 read(567,*) x1,x2,x3,ntpts 182c 183c Find the point on the boundary (if it is on this processor) 184c that matches this point 185c 186 do i=1,numnp 187 if(ibits(ibc(i),3,3) .eq.7) then 188 dd= distds(x1,x2,x3,x(i,1),x(i,2),x(i,3)) 189 if(dd.lt.epsd) then 190 ic=ic+1 191 nBCt(ic)=i ! the pointer to this point 192! numBCt(ic)=ntpts ! the number of time series 193 do j=1,ntpts 194c read(567,*) BCt(ic,j,4),(BCt(ic,j,n),n=1,3) 195 read(567,*) (BCt(j,n),n=1,4) 196 enddo 197 ! at this point we have all bc data of spacial point 198 ! ic/i in all time domain. now we figure out which subset 199 ! is necessary. 200 if (tlen.ge.BCt(ntpts,4)) then 201 ! whole run is larger than whole period 202 ! should take all data 203 ilen = ntpts 204 allocate(BCtptr(ic)%p(ilen,4)) 205 BCtptr(ic)%p = BCt(1:ilen,:) 206 else 207 nper = t0/BCt(ntpts,4) 208 tstart = t0-nper*BCt(ntpts,4) 209 nper = t1/BCt(ntpts,4) 210 tend = t1-nper*BCt(ntpts,4) 211 istart = iBfind(BCt(:,4),ntpts,tstart) 212 iend = iBfind(BCt(:,4),ntpts,tend) 213 if (istart>1) istart=istart-1 214 if (iend<ntpts) iend=iend+1 215 !write(*,*)'bcstart:',BCt(istart,4),tstart 216 !write(*,*)'bcend:',BCt(iend,4),tend 217 if (tstart.lt.tend) then ! does not loop 218 ilen = iend-istart+1 219 allocate(BCtptr(ic)%p(ilen,4)) 220 BCtptr(ic)%p = BCt(istart:iend,:) 221 else ! loop within two consequetive time period 222 i0 = (ntpts-istart+1) ! first segment 223 ilen = iend + i0 224 allocate(BCtptr(ic)%p(ilen,4)) 225 BCtptr(ic)%p(1:i0,:) = BCt(istart:ntpts,:) 226 BCtptr(ic)%p(i0+1:ilen,:) = BCt(1:iend,:) 227 BCtptr(ic)%p(i0+1:ilen,4) = 228 & BCtptr(ic)%p(i0+1:ilen,4) + BCt(ntpts,4) 229 endif 230 endif 231 numBCt(ic) = ilen 232 BCtptr(ic)%p(:,4) = BCtptr(ic)%p(:,4)*bcttimescale 233 BCperiod(ic) = BCt(ntpts,4)*bcttimescale 234 exit 235 endif 236 endif 237 enddo 238 if(i.eq.numnp+1) then 239c 240c if we get here the point was not found. It must be on another 241c processor so we read past this record and move on 242c 243 do j=1,ntpts 244 read(567,*) rj1,rj2,rj3,rj4 245 enddo 246 endif 247 enddo ! end of the loop over ntv 248 !BCt(:,:,4)=BCt(:,:,4)*bcttimescale 249 deallocate(BCt) 250 endif ! any 3 component nodes 251 itvn=ic 252 close(567) 253 if (ic.gt.0) then 254 write(*,*)'myrank=',myrank,' and I found ',ic,' nodes.' 255c else 256c deallocate(nBCt) 257c deallocate(numBCt) 258c deallocate(BCt) 259 endif 260 261 return 262 end 263 264!*************************************************************** 265! A Binary search return an index of a real array. 266! This array should be an ascending sorted array. 267!*************************************************************** 268 function iBfind(bArray,bLen,bVal) 269 integer,intent(in) :: bLen 270 real(kind=8),intent(in),dimension(bLen) :: bArray 271 real(kind=8),intent(in) :: bVal 272 integer mlen,indx,bstart,bend 273 bstart = 1 274 bend = bLen 275 indx = (bLen+1)/2 276 do while ((bstart+1).lt.bend) 277 if (bVal.gt.bArray(indx)) then 278 bstart = indx 279 else 280 bend = indx 281 endif 282 indx = (bstart+bend)/2 283 enddo 284 iBfind = indx 285 return 286 end function 287 288 289 subroutine BCint(timel,shp,shgl,shpb,shglb,x,BC,iBC) 290 291 use specialBC ! brings in itvn,nbct, bct, numbct, nptsmax 292 293 include "common.h" 294 295 real*8 BC(nshg,ndofBC), timel,t 296 real*8 x(numnp,nsd), 297 & shp(MAXTOP,maxsh,MAXQPT), 298 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 299 & shpb(MAXTOP,maxsh,MAXQPT), 300 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 301 302 integer iBC(numnp),nlast,i,j,nper 303 304 do i =1,itvn ! itvn is the number of varying nodes on this proc 305 306 nlast=numBCt(i) ! number of time series to interpolate from 307 nper=timel/BCperiod(i) 308 ! number of periods completed to shift off 309 310 311 t=timel-nper*BCperiod(i) ! now time in periodic domain 312 313 do j=2,nlast !loop to find the interval that we are in 314 315 if(BCtptr(i)%p(j,4).gt.t) then 316 ! this is upper bound, j-1 is lower 317 318 wr=(t-BCtptr(i)%p(j-1,4))/(BCtptr(i)%p(j,4)-BCtptr(i)%p(j-1,4)) 319 BC(nbct(i),3:5)= BCtptr(i)%p(j-1,1:3)*(one-wr) 320 & + BCtptr(i)%p(j,1:3)*wr 321 exit 322 323 endif 324 enddo 325 enddo 326 return 327 end 328 329 function distds(x1,y1,z1,x2,y2,z2) 330 real*8 distds 331 real*8 x1,y1,z1,x2,y2,z2,x,y,z 332 x=x1-x2 333 y=y1-y2 334 z=z1-z2 335 distds=x*x+y*y+z*z 336 return 337 end 338c----------------------------------------------------------------------- 339c initialize the impedance boundary condition: 340c read the data in initImpt 341c interpolate the data to match the process time step in Impint 342c----------------------------------------------------------------------- 343 subroutine initImpt() 344 345 use convolImpFlow 346 include "common.h" 347 348 open(unit=817, file='impt.dat',status='old') 349 read (817,*) nptsImpmax 350 allocate (numImpt(numImpSrfs)) 351 allocate (ValueImpt(nptsImpmax,2,numImpSrfs)) 352 ValueImpt=0 353 do k=1,numImpSrfs 354 read (817,*) numDataImp 355 numImpt(k) = numDataImp 356 do j=1,numDataImp 357 read(817,*) (ValueImpt(j,n,k),n=1,2) ! n=1 time, 2 value 358 enddo 359 enddo 360 close(817) 361 362 allocate (ValueListImp(ntimeptpT+1,numImpSrfs)) 363 ValueListImp = zero 364 ValueListImp(ntimeptpT+1,:) = ValueImpt(1,2,:) !Z(time=0), last entry 365 ValueListImp(1,:) = ValueImpt(1,2,:) !Z(time=0)=Z(time=T) 366 return 367 end 368 369 370 371 subroutine Impint(ctime,jstep) 372 373 use convolImpFlow 374 include "common.h" 375 376 real*8 ctime, ptime 377 integer nlast, nper, k, j , jstep 378 379 380 do k =1,numImpSrfs 381 nlast=numImpt(k) ! number of time series to interpolate from 382 nper=ctime/ValueImpt(nlast,1,k)!number of periods completed to shift off 383 ptime = ctime-nper*ValueImpt(nlast,1,k) ! now time in periodic domain 384 385 do j=2,nlast !loop to find the interval that we are in 386 387 if(ValueImpt(j,1,k).gt.ptime) then ! this is upper bound, j-1 is lower 388 wr=(ptime-ValueImpt(j-1,1,k)) 389 & / ( ValueImpt(j,1,k)-ValueImpt(j-1,1,k) ) 390 ValueListImp(jstep,k)= ValueImpt(j-1,2,k)*(one-wr) 391 & + ValueImpt(j,2,k)*wr 392 exit 393 endif 394 395 enddo 396 enddo 397 return 398 end 399 400c----------------------------------------------------------------------------- 401c time filter for a periodic function (sin cardinal + window function) 402c is used for the impedance and the flow rate history 403c----------------------------------------------------------------------------- 404 subroutine Filter(Filtered,DataHist,nptf,timestep,cutfreq) 405 406 include "common.h" 407 408 integer nptf, cutfreq, j, k, m, s, Filtime(nptf) 409 real*8 DataHist(nptf,numImpSrfs), Window(nptf) 410 real*8 Sinc(nptf), FilterSW(nptf), Filtered(nptf,numImpSrfs) 411 real*8 windK, timestep 412 413 windK = cutfreq*2 + 1 414 do j=1,nptf 415 Filtime(j) = j-1 416 Window(j) = 0.42+0.5*cos(2*pi*Filtime(j)/nptf) 417 & +0.08*cos(4*pi*Filtime(j)/nptf) 418 Sinc(j) = sin(pi*Filtime(j)*windK/nptf)/sin(pi*Filtime(j)/nptf) 419 enddo 420 Sinc(1) = windK 421 422 do j=1,nptf 423 FilterSW(j) = Window(nptf+1-j)*Sinc(nptf+1-j) !filter for convolution 424 enddo 425 426 Filtered = zero 427 do m=1,nptf 428 do j=1,nptf 429 s=modulo(m-nptf+j,nptf) 430 if(s.eq.zero) then 431 s=nptf 432 endif 433 Filtered(m,:) = Filtered(m,:) 434 & +FilterSW(j)*DataHist(s,:)/nptf !filter convolution 435 enddo 436 enddo 437 438 return 439 end 440c----------------------------------------------------------------------- 441c initialize the RCR boundary condition: 442c read the data in initRCRt 443c interpolate the data to match the process time step in RCRint 444c----------------------------------------------------------------------- 445 subroutine initRCRt() 446 447 use convolRCRFlow 448 include "common.h" 449 450 open(unit=818, file='rcrt.dat',status='old') 451 read (818,*) nptsRCRmax 452 allocate (numRCRt(numRCRSrfs)) 453 allocate (ValuePdist(nptsRCRmax,2,numRCRSrfs)) 454 allocate (ValueListRCR(3,numRCRSrfs)) 455 ValuePdist=0 456 ValueListRCR=0 457 do k=1,numRCRSrfs 458 read (818,*) numDataRCR 459 numRCRt(k) = numDataRCR 460 do j=1,3 461 read(818,*) ValueListRCR(j,k) ! reads Rp,C,Rd 462 enddo 463 do j=1,numDataRCR 464 read(818,*) (ValuePdist(j,n,k),n=1,2) ! n=1 time, 2 value 465 enddo 466 enddo 467 close(818) 468 469 allocate (dtRCR(numRCRSrfs)) 470 471 return 472 end 473 474 475 subroutine RCRint(ctime,Pdist) 476 477 use convolRCRFlow ! brings numRCRSrfs, ValuePdist 478 include "common.h" 479 480 real*8 ctime, ptime 481 integer nlast, nper, k, j 482 real*8 Pdist(0:MAXSURF) 483 484 do k =1,numRCRSrfs 485 nlast=numRCRt(k) ! number of time series to interpolate from 486 nper=ctime/ValuePdist(nlast,1,k)!number of periods completed to shift off 487 ptime = ctime-nper*ValuePdist(nlast,1,k) ! now time in periodic domain 488 489 do j=2,nlast !loop to find the interval that we are in 490 491 if(ValuePdist(j,1,k).gt.ptime) then ! this is upper bound, j-1 is lower 492 wr=(ptime-ValuePdist(j-1,1,k)) 493 & / ( ValuePdist(j,1,k)-ValuePdist(j-1,1,k) ) 494 Pdist(k)= ValuePdist(j-1,2,k)*(one-wr) 495 & + ValuePdist(j,2,k)*wr 496 exit 497 endif 498 499 enddo 500 enddo 501 return 502 end 503 504c----------------------------------------------------------------------- 505c returns in pold the history dependent part of the pressure in the 506c impedance/flow rate convolution for the impedance and RCR BC 507c----------------------------------------------------------------------- 508 subroutine pHist(pressHist,QHist,betas,nTimePoint,nSrfs) 509 510 include "common.h" 511 512 integer nTimePoint,nSrfs 513 real*8 pressHist(0:MAXSURF) 514 real*8 QHist(nTimePoint+1,nSrfs),betas(nTimePoint+2,nSrfs) 515 !don't need here betas(ntimePoint+2) 516 !but pb of array passing if cut at nTimePoint+1 517 pressHist=zero 518 do k=1,nSrfs 519 do j=1,nTimePoint+1 520 pressHist(k) = pressHist(k) + QHist(j,k)*betas(j,k) 521 enddo 522 enddo 523 return 524 end 525