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