Skip to content

Commit 4957b5f

Browse files
jh6016jh6016
authored
added fac_pre_fluctuation
1 parent b5b9c49 commit 4957b5f

File tree

1 file changed

+86
-26
lines changed

1 file changed

+86
-26
lines changed

‎src/modibm.f90‎

Lines changed: 86 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -106,7 +106,7 @@ end function interp_temperature
106106

107107
type(bound_info_type) :: bound_info_u, bound_info_v, bound_info_w, bound_info_c
108108

109-
integer :: nstatfac=6, ncidfac, nrecfac=0
109+
integer :: nstatfac=7, ncidfac, nrecfac=0
110110
character(80), allocatable :: ncstatfac(:,:)
111111
character(80) :: facname = 'fac.xxx.nc'
112112
character(80),dimension(1,4) :: tncstatfac
@@ -115,12 +115,14 @@ end function interp_temperature
115115
real, allocatable :: fac_tau_y(:)
116116
real, allocatable :: fac_tau_z(:)
117117
real, allocatable :: fac_pres(:)
118+
real, allocatable :: fac_pres2(:)
118119
real, allocatable :: fac_htc(:)
119120
real, allocatable :: fac_cth(:)
120121
real, allocatable :: fac_tau_x_av(:)
121122
real, allocatable :: fac_tau_y_av(:)
122123
real, allocatable :: fac_tau_z_av(:)
123124
real, allocatable :: fac_pres_av(:)
125+
real, allocatable :: fac_pres2_av(:)
124126
real, allocatable :: fac_htc_av(:)
125127
real, allocatable :: fac_cth_av(:)
126128

@@ -172,10 +174,22 @@ subroutine initibm
172174
bound_info_v%nfctsecs = nfctsecs_v
173175
bound_info_w%nfctsecs = nfctsecs_w
174176
call initibmwallfun('fluid_boundary_u.txt', 'facet_sections_u.txt', xhat, bound_info_u)
177+
if (myid==0) then
178+
write(*,*) "u done"
179+
end if
180+
175181
call initibmwallfun('fluid_boundary_v.txt', 'facet_sections_v.txt', yhat, bound_info_v)
182+
if (myid==0) then
183+
write(*,*) "v done"
184+
end if
176185
call initibmwallfun('fluid_boundary_w.txt', 'facet_sections_w.txt', zhat, bound_info_w)
186+
if (myid==0) then
187+
write(*,*) "w done"
188+
end if
177189
end if
178190

191+
! stop 1
192+
179193
if (ltempeq .or. lmoist .or. nsv>0 .or. lwritefac) then
180194
solid_info_c%nsolpts = nsolpts_c
181195
call initibmnorm('solid_c.txt', solid_info_c)
@@ -184,6 +198,9 @@ subroutine initibm
184198
bound_info_c%nfctsecs = nfctsecs_c
185199
call initibmwallfun('fluid_boundary_c.txt', 'facet_sections_c.txt', vec0, bound_info_c)
186200

201+
if (myid==0) then
202+
write(*,*) "c done"
203+
end if
187204
allocate(mask_c(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh)); mask_c = 1.
188205
mask_c(:,:,kb-kh) = 0.
189206
call solid(solid_info_c, mask_c, rhs, 0., ih, jh, kh)
@@ -198,24 +215,28 @@ subroutine initibm
198215
allocate(fac_tau_y(1:nfcts))
199216
allocate(fac_tau_z(1:nfcts))
200217
allocate(fac_pres(1:nfcts))
218+
allocate(fac_pres2(1:nfcts))
201219
allocate(fac_htc(1:nfcts))
202220
allocate(fac_cth(1:nfcts))
203221
fac_tau_x = 0.
204222
fac_tau_y = 0.
205223
fac_tau_z = 0.
206224
fac_pres = 0.
225+
fac_pres2=0.
207226
fac_htc = 0.
208227
fac_cth = 0.
209228
allocate(fac_tau_x_av(1:nfcts))
210229
allocate(fac_tau_y_av(1:nfcts))
211230
allocate(fac_tau_z_av(1:nfcts))
212231
allocate(fac_pres_av(1:nfcts))
232+
allocate(fac_pres2_av(1:nfcts))
213233
allocate(fac_htc_av(1:nfcts))
214234
allocate(fac_cth_av(1:nfcts))
215235
fac_tau_x_av = 0.
216236
fac_tau_y_av = 0.
217237
fac_tau_z_av = 0.
218238
fac_pres_av = 0.
239+
fac_pres2_av=0.
219240
fac_htc_av = 0.
220241
fac_cth_av = 0.
221242

@@ -228,6 +249,7 @@ subroutine initibm
228249
call ncinfo(ncstatfac( 4,:),'pres', 'pressure', 'Pa','ft')
229250
call ncinfo(ncstatfac( 5,:),'htc', 'heat transfer coefficient', '','ft')
230251
call ncinfo(ncstatfac( 6,:),'cth', 'heat transfer coefficient (Ivo)', '','ft')
252+
call ncinfo(ncstatfac( 7,:),'pres_flc', 'pressure fluctuation', '','ft')
231253

232254
if (myid==0) then
233255
call open_nc(facname, ncidfac, nrecfac, nfcts=nfcts)
@@ -523,52 +545,73 @@ subroutine initibmwallfun(fname_bnd, fname_sec, dir, bound_info)
523545
! u
524546
if ((bound_info%recpts(n,1) < xh(bound_info%recids_u(n,1))) .or. &
525547
(bound_info%recpts(n,1) > xh(bound_info%recids_u(n,1)+1))) then
526-
write(*,*) "ERROR: x out of bounds"
527-
stop 1
548+
write(*,*) "ERROR: x out of bounds u"
549+
! write(*,*) n
550+
! stop 1
528551
end if
529552
if ((bound_info%recpts(n,2) < yf(bound_info%recids_u(n,2))) .or. &
530553
(bound_info%recpts(n,2) > yf(bound_info%recids_u(n,2)+1))) then
531-
write(*,*) "ERROR: y out of bounds"
532-
stop 1
554+
write(*,*) "ERROR: y out of bounds u"
555+
write(*,*) n
556+
write(*,*) bound_info%secfacids(n)
557+
write(*,*) bound_info%secareas(n)
558+
write(*,*) bound_info%secbndptids(n)
559+
write(*,*) bound_info%bnddst(n)
560+
!stop 1
533561
end if
534562
if ((bound_info%recpts(n,3) < zf(bound_info%recids_u(n,3))) .or. &
535563
(bound_info%recpts(n,3) > zf(bound_info%recids_u(n,3)+1))) then
536-
write(*,*) "ERROR: z out of bounds"
537-
stop 1
564+
write(*,*) "ERROR: z out of bounds u"
565+
!write(*,*) n
566+
! stop 1
538567
end if
539568

540569
! v
541570
if ((bound_info%recpts(n,1) < xf(bound_info%recids_v(n,1))) .or. &
542571
(bound_info%recpts(n,1) > xf(bound_info%recids_v(n,1)+1))) then
543-
write(*,*) "ERROR: x out of bounds"
544-
stop 1
572+
write(*,*) "ERROR: x out of bounds v"
573+
!write(*,*) n
574+
! stop 1
545575
end if
546576
if ((bound_info%recpts(n,2) < yh(bound_info%recids_v(n,2))) .or. &
547577
(bound_info%recpts(n,2) > yh(bound_info%recids_v(n,2)+1))) then
548-
write(*,*) "ERROR: y out of bounds"
549-
stop 1
578+
write(*,*) "ERROR: y out of bounds v"
579+
!write(*,*) n
580+
! stop 1
550581
end if
551582
if ((bound_info%recpts(n,3) < zf(bound_info%recids_v(n,3))) .or. &
552583
(bound_info%recpts(n,3) > zf(bound_info%recids_v(n,3)+1))) then
553-
write(*,*) "ERROR: z out of bounds"
554-
stop 1
584+
write(*,*) "ERROR: z out of bounds v"
585+
! write(*,*) n
586+
! stop 1
555587
end if
556588

557589
! w
558590
if ((bound_info%recpts(n,1) < xf(bound_info%recids_w(n,1))) .or. &
559591
(bound_info%recpts(n,1) > xf(bound_info%recids_w(n,1)+1))) then
560-
write(*,*) "ERROR: x out of bounds"
561-
stop 1
592+
write(*,*) "ERROR: x out of bounds w"
593+
write(*,*) bound_info%secfacids(n)
594+
write(*,*) bound_info%secareas(n)
595+
write(*,*) bound_info%secbndptids(n)
596+
write(*,*) bound_info%bnddst(n)
597+
!write(*,*) n
598+
! stop 1
562599
end if
563600
if ((bound_info%recpts(n,2) < yf(bound_info%recids_w(n,2))) .or. &
564601
(bound_info%recpts(n,2) > yf(bound_info%recids_w(n,2)+1))) then
565-
write(*,*) "ERROR: y out of bounds"
566-
stop 1
602+
write(*,*) "ERROR: y out of bounds w"
603+
write(*,*) bound_info%secfacids(n)
604+
write(*,*) bound_info%secareas(n)
605+
write(*,*) bound_info%secbndptids(n)
606+
write(*,*) bound_info%bnddst(n)
607+
!write(*,*) n
608+
! stop 1
567609
end if
568610
if ((bound_info%recpts(n,3) < zh(bound_info%recids_w(n,3))) .or. &
569611
(bound_info%recpts(n,3) > zh(bound_info%recids_w(n,3)+1))) then
570-
write(*,*) "ERROR: z out of bounds"
571-
stop 1
612+
!write(*,*) "ERROR: z out of bounds w"
613+
write(*,*) n
614+
! stop 1
572615
end if
573616
end if
574617
end do
@@ -1178,7 +1221,7 @@ end subroutine diffc_corr
11781221

11791222

11801223
subroutine ibmwallfun
1181-
use modglobal, only : libm, iwallmom, iwalltemp, xhat, yhat, zhat, ltempeq, lmoist, &
1224+
use modglobal, only : libm, iwallmom, iwalltemp, xhat, yhat, zhat, ltempeq, lmoist, skip_time, &
11821225
ib, ie, ih, ihc, jb, je, jh, jhc, kb, ke, kh, khc, nsv, totheatflux, totqflux, nfcts, rk3step, timee, nfcts, lwritefac, dt, dtfac, tfac, tnextfac
11831226
use modfields, only : u0, v0, w0, thl0, qt0, sv0, up, vp, wp, thlp, qtp, svp, &
11841227
tau_x, tau_y, tau_z, thl_flux
@@ -1259,17 +1302,21 @@ subroutine ibmwallfun
12591302

12601303
deallocate(rhs)
12611304

1262-
if (lwritefac .and. rk3step==3) then
1263-
1305+
! if (lwritefac .and. rk3step==3 .and. timee>=skip_time) then
1306+
if (lwritefac .and. rk3step==3 ) then
12641307
if (myid == 0) then
12651308
fac_tau_x_av = fac_tau_x_av + dt*fac_tau_x
12661309
fac_tau_y_av = fac_tau_y_av + dt*fac_tau_y
12671310
fac_tau_z_av = fac_tau_z_av + dt*fac_tau_z
12681311
fac_pres_av = fac_pres_av + dt*fac_pres
1312+
fac_pres2_av = fac_pres2_av + dt*fac_pres2
1313+
12691314
fac_htc_av = fac_htc_av + dt*fac_htc
12701315
fac_cth_av = fac_cth_av + dt*fac_cth
1271-
1272-
if (timee >= tnextfac) then
1316+
1317+
! if (timee >= tnextfac + skip_time) then
1318+
if (timee >= tnextfac) then
1319+
! tfac = timee - max(skip_time, tfac)
12731320
tfac = timee - tfac
12741321
allocate(varsfac(nfcts,nstatfac))
12751322
varsfac(:,1) = fac_tau_x_av(1:nfcts)/tfac
@@ -1278,19 +1325,25 @@ subroutine ibmwallfun
12781325
varsfac(:,4) = fac_pres_av(1:nfcts)/tfac
12791326
varsfac(:,5) = fac_htc_av(1:nfcts)/tfac
12801327
varsfac(:,6) = fac_cth_av(1:nfcts)/tfac
1328+
1329+
varsfac(:,7) = fac_pres2_av(1:nfcts)/tfac - (fac_pres_av(1:nfcts)/dtfac * fac_pres_av(1:nfcts)/tfac)
1330+
12811331
call writestat_nc(ncidfac,1,tncstatfac,(/timee/),nrecfac,.true.)
12821332
call writestat_1D_nc(ncidfac,nstatfac,ncstatfac,varsfac,nrecfac,nfcts)
12831333
deallocate(varsfac)
12841334

12851335
tfac = timee
1286-
tnextfac = NINT((timee + dtfac))*1.0
1336+
! tnextfac = NINT((timee + dtfac - skip_time))*1.0
12871337

1338+
tnextfac = NINT((timee + dtfac))*1.0
12881339
fac_tau_x_av = 0.
12891340
fac_tau_y_av = 0.
12901341
fac_tau_z_av = 0.
12911342
fac_pres_av = 0.
1343+
fac_pres2_av= 0.
12921344
fac_htc_av = 0.
12931345
fac_cth_av = 0.
1346+
12941347
end if
12951348

12961349
end if !myid
@@ -1463,12 +1516,14 @@ subroutine wallfunheat
14631516
integer i, j, k, n, m, sec, fac
14641517
real :: dist, flux, area, vol, tempvol, Tair, Tsurf, utan, cth, htc, cveg, hurel, qtair, qwall, resa, resc, ress, xrec, yrec, zrec
14651518
real, dimension(3) :: uvec, norm, span, strm
1466-
real, dimension(1:nfcts) :: fac_htc_loc, fac_cth_loc, fac_pres_loc
1519+
real, dimension(1:nfcts) :: fac_htc_loc, fac_cth_loc, fac_pres_loc, fac_pres2_loc
14671520
logical :: valid
14681521

14691522
fac_htc_loc = 0.
14701523
fac_cth_loc = 0.
14711524
fac_pres_loc = 0.
1525+
fac_pres2_loc = 0.
1526+
14721527
do sec = 1,bound_info_c%nfctsecsrank
14731528
! sec = bound_info_c%fctsecsrank(m) ! index of section
14741529
!n = bound_info_c%secbndptids(sec) ! index of boundary point
@@ -1489,6 +1544,8 @@ subroutine wallfunheat
14891544
end if
14901545

14911546
fac_pres_loc(fac) = fac_pres_loc(fac) + pres0(i,j,k) * area ! output pressure on facets
1547+
fac_pres2_loc(fac) = fac_pres2_loc(fac) + pres0(i,j,k)* pres0(i,j,k) * area
1548+
14921549

14931550
if (bound_info_c%lskipsec_loc(sec)) cycle
14941551
!if (facz0(fac) < eps1) cycle
@@ -1610,9 +1667,12 @@ subroutine wallfunheat
16101667
fac_cth_loc(1:nfcts) = fac_cth_loc(1:nfcts) / faca(1:nfcts)
16111668
fac_htc_loc(1:nfcts) = fac_htc_loc(1:nfcts) / faca(1:nfcts)
16121669
fac_pres_loc(1:nfcts) = fac_pres_loc(1:nfcts) / faca(1:nfcts)
1670+
fac_pres2_loc(1:nfcts) = fac_pres2_loc(1:nfcts) / faca(1:nfcts)
1671+
16131672
call MPI_ALLREDUCE(fac_cth_loc(1:nfcts), fac_cth(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr)
16141673
call MPI_ALLREDUCE(fac_htc_loc(1:nfcts), fac_htc(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr)
16151674
call MPI_ALLREDUCE(fac_pres_loc(1:nfcts), fac_pres(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr)
1675+
call MPI_ALLREDUCE(fac_pres2_loc(1:nfcts), fac_pres2(1:nfcts), nfcts, MY_REAL, MPI_SUM, comm3d, mpierr)
16161676
end if
16171677

16181678
end subroutine wallfunheat

0 commit comments

Comments
 (0)