Skip to content

Commit

Permalink
cleaned diffusive_utils.90 and related
Browse files Browse the repository at this point in the history
  • Loading branch information
Dong Kim committed Apr 20, 2022
1 parent 5602b03 commit 8d61d36
Show file tree
Hide file tree
Showing 7 changed files with 49 additions and 128 deletions.
78 changes: 20 additions & 58 deletions src/kernel/diffusive/diffusive.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,10 +150,6 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
double precision, dimension(mxnbathy_g, mxncomp_g, nrch_g),intent(in ) :: mann_bathy_g
double precision, dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g
double precision, dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: elv_ev_g
!TODO for DA
!integer, intent(in) :: nts_da_g
!double precision, dimension(nts_da_g, nrch_g), intent(in) :: qda_g


! Local variables
integer :: ncomp
Expand Down Expand Up @@ -206,22 +202,6 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
double precision, dimension(:,:), allocatable :: skMain
double precision, dimension(:,:), allocatable :: skRight

! test
integer :: ts
! DA
!doubleprecision :: dt_da, dmy, dmy2
!integer :: idmy
open(unit=1, file="./input/sine_tide_matagorda_city_31d.txt")
!open(unit=1, file="t-route/src/kernel/diffusive/input/tide_matagorda_city_31d.txt")
!open(unit=2, file="./input/usgsBastrop_06011200_06302355_2020.txt")
!open(unit=91, file="./output/iniq.txt")
!open(unit=101, file="./output/cn-mod simulated discharge depth elev_lowercolorado.txt")
!open(unit=102, file="./output/cn-mod q elev depth at each sim_time.txt")
allocate(dbcd(nts_db_g))
do n = 1, nts_db_g
read(1,*) dmyi, dbcd(n), dmyi
dbcd(n) = dbcd(n) + 5.0
end do
!-----------------------------------------------------------------------------
! Time domain parameters
dtini = timestep_ar_g(1) ! initial timestep duration [sec]
Expand Down Expand Up @@ -337,6 +317,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
allocate(size_bathy(mxncomp, nlinks))
allocate(usgs_da_reach(nlinks))
allocate(usgs_da(nts_da, nlinks))
allocate(dbcd(nts_db_g))

!--------------------------------------------------------------------------------------------
frnw_g = frnw_ar_g ! network mapping matrix
Expand All @@ -362,7 +343,7 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
dimensionless_Di = -999
dimensionless_Fc = -999
dimensionless_D = -999

dbcd = 0.0 !TODO: pass downstream boundary values from Python

!-----------------------------------------------------------------------------
! Identify mainstem reaches and list their ids in an array
Expand All @@ -383,13 +364,6 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt
end do
deallocate(dmy_frj)

! test
!do jm = 1, nmstem_rch !* mainstem reach only
! j = mstem_frj(jm)
! do i = 1, frnw_g(j, 1)-1
! write(91,*) i, j, iniq(i,j)
! end do
!end do
!-----------------------------------------------------------------------------
! create dx array from dx_ar_g and determine minimum dx.

Expand Down Expand Up @@ -843,17 +817,6 @@ subroutine diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qt

end do ! end of time loop

!* test
!do ts=1, ntss_ev_g
!do jm = 1, nmstem_rch !* mainstem reach only
! j = mstem_frj(jm)
! do i=1, frnw_g(j,1)
! write(101,"(f10.1, 2I10, 4F20.4)") saveInterval*real(ts-1)/60.0, i, j, q_ev_g(ts, i, j),&
! z(i,j), elv_ev_g(ts, i, j) - z(i,j), elv_ev_g(ts, i, j)
! end do
!end do
!end do

deallocate(frnw_g)
deallocate(area, bo, pere, areap, qp, z, depth, sk, co, dx)
deallocate(volRemain, froud, courant, oldQ, newQ, oldArea, newArea, oldY, newY)
Expand Down Expand Up @@ -1067,7 +1030,7 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)

! Local variables
integer :: ncomp
integer :: i
integer :: i, irow, flag_da
double precision :: a1, a2, a3, a4
double precision :: b1, b2, b3, b4
double precision :: dd1, dd2, dd3, dd4
Expand Down Expand Up @@ -1220,13 +1183,21 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)
varr_da(n) = usgs_da(n, j)
end do
qp(ncomp,j) = intp_y(nts_da, tarr_da, varr_da, t)
if (qp(ncomp,j) <= 0.0) then
flag_da = 1
! check usgs_da value is in good quality
irow = locate(tarr_da, t)
if (irow == nts_da) then
irow = irow-1
endif
if ((varr_da(irow)<= -4443.999).or.(varr_da(irow+1)<= -4443.999)) then
! when usgs data is missing or in poor quality
qp(ncomp,j) = eei(ncomp) * qp_ghost + ffi(ncomp)
flag_da = 0
endif
deallocate(varr_da)
else
qp(ncomp,j) = eei(ncomp) * qp_ghost + ffi(ncomp)
flag_da = 0
endif
qpx(ncomp,j) = exi(ncomp) *qpx_ghost + fxi(ncomp)

Expand All @@ -1236,7 +1207,7 @@ subroutine mesh_diffusive_forward(dtini_given, t0, t, tfin, saveInterval,j)
end do

! when a reach hasn't been applied to DA
if ((usgs_da_reach(j) == 0).or.(qp(ncomp,j) <= 0.0)) then
if ((usgs_da_reach(j) == 0).or.(flag_da == 0)) then
qp(1, j) = newQ(1, j)
qp(1, j) = qp(1, j) + allqlat
endif
Expand Down Expand Up @@ -1494,30 +1465,24 @@ function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds)
integer :: xcolID, ycolID
double precision :: x1, x2, df, dxx, dxold, f, fh, fl, temp, xh, xl
double precision :: y_norm, y_ulm_multi, y_llm_multi, elv_norm, y_old
double precision :: rtsafe

!open(unit=201, file="./output/depth computing.txt")

double precision :: rtsafe

y_ulm_multi = 2.0
y_llm_multi = 0.1

xcolID = 10
ycolID = 1
elv_norm = intp_xsec_tab(i, j, nel, xcolID, ycolID, abs(Q_cur)) ! normal elevation not depth
y_norm = elv_norm - z(i, j)
y_old = oldY(i, j) - z(i, j)

!write(201,*) "very 0: ", i, j, z(i,j), Q_cur, Q_ds, elv_norm, oldY(i,j), y_norm, y_old

y_old = oldY(i, j) - z(i, j)

! option 1 for initial point
!x1 = y_norm * y_llm_multi
!x2 = y_norm * y_ulm_multi
! option 2 for initial point
x1 = 0.5 * (y_norm + y_old) * y_llm_multi
x2 = 0.5 * (y_norm + y_old) * y_ulm_multi

!write(201,*) "initial x:", i, j, x1, x2, y_norm


call funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, x1, y_ds, fl, df)

call funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, x2, y_ds, fh, df)
Expand All @@ -1529,11 +1494,9 @@ function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds)

if (fl == 0.0) then
rtsafe = x1
!write(201,*) "return 1: ", i, j, x1, x2, rtsafe
return
elseif (fh == 0.0) then
rtsafe = x2
!write(201,*) "return 2: ", i, j, x1, x2, rtsafe
return
elseif (fl < 0.0) then ! orient the search so that f(xl) < 0.
xl = x1
Expand Down Expand Up @@ -1564,7 +1527,7 @@ function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds)
rtsafe = rtsafe - dxx
if (temp == rtsafe) return
end if
!write(201,*) "iteration:", i, j, xl, xh, rtsafe

if (abs(dxx) < xacc) return ! convergence criterion.

! one new function evaluation per iteration.
Expand All @@ -1579,8 +1542,7 @@ function rtsafe(i, j, Q_cur, Q_ds, z_cur, z_ds, y_ds)

! when root is not converged:
rtsafe = y_norm
!write(201,*) "no covg: ", i, j, rtsafe


end function rtsafe

subroutine funcd_diffdepth(i, j, Q_cur, Q_ds, z_cur, z_ds, y_cur, y_ds, f, df)
Expand Down
16 changes: 2 additions & 14 deletions src/kernel/diffusive/pydiffusive.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,6 @@ module diffusive_interface

implicit none
contains
!subroutine c_diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, &
! nts_qtrib_g, mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, &
! tw_ar_g, twcc_ar_g, mann_ar_g, manncc_ar_g, so_ar_g, dx_ar_g, &
! iniq, frnw_col, frnw_g, qlat_g, ubcd_g, dbcd_g, qtrib_g, &
! paradim, para_ar_g, mxnbathy_g, x_bathy_g, z_bathy_g, &
! mann_bathy_g, size_bathy_g, q_ev_g, elv_ev_g) bind(c)
subroutine c_diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qtrib_g, nts_da_g, &
mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g, mann_ar_g, &
manncc_ar_g, so_ar_g, dx_ar_g, &
Expand Down Expand Up @@ -43,14 +37,8 @@ subroutine c_diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_
real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: x_bathy_g
real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: z_bathy_g
real(c_double), dimension(mxnbathy_g, mxncomp_g, nrch_g), intent(in ) :: mann_bathy_g
real(c_double), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g, elv_ev_g

!call diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, &
! nts_qtrib_g, mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, &
! tw_ar_g, twcc_ar_g, mann_ar_g, manncc_ar_g, so_ar_g, dx_ar_g, &
! iniq, frnw_col, frnw_g, qlat_g, ubcd_g, dbcd_g, qtrib_g, &
! paradim, para_ar_g, mxnbathy_g, x_bathy_g, z_bathy_g, &
! mann_bathy_g, size_bathy_g, q_ev_g, elv_ev_g)
real(c_double), dimension(ntss_ev_g, mxncomp_g, nrch_g), intent(out) :: q_ev_g, elv_ev_g

call diffnw(timestep_ar_g, nts_ql_g, nts_ub_g, nts_db_g, ntss_ev_g, nts_qtrib_g, nts_da_g, &
mxncomp_g, nrch_g, z_ar_g, bo_ar_g, traps_ar_g, tw_ar_g, twcc_ar_g, mann_ar_g, &
manncc_ar_g, so_ar_g, dx_ar_g, &
Expand Down
4 changes: 2 additions & 2 deletions src/troute-network/troute/nhd_network_utilities_v02.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -754,8 +754,8 @@ def build_data_assimilation_lastobs(data_assimilation_parameters):
da_parameter_dict["da_decay_coefficient"] = data_assimilation_parameters.get(
"da_decay_coefficient",
120)
da_parameter_dict["diffusive_streamflow_DA"] = streamflow_da_parameters.get(
"diffusive_streamflow_DA", False)
da_parameter_dict["diffusive_streamflow_nudging"] = streamflow_da_parameters.get(
"diffusive_streamflow_nudging", False)

# TODO: Add parameters here for interpolation length (14/59), QC threshold (1.0)

Expand Down
3 changes: 2 additions & 1 deletion src/troute-nwm/src/nwm_routing/__main__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,13 +127,14 @@ def nwm_route(

LOG.debug("MC computation complete in %s seconds." % (time.time() - start_time_mc))
start_time_diff = time.time()

#import pdb; pdb.set_trace()
# call diffusive wave simulation and append results to MC results
results.extend(
compute_diffusive_routing(
results,
diffusive_network_data,
cpu_pool,
t0,
dt,
nts,
q0,
Expand Down
4 changes: 2 additions & 2 deletions src/troute-routing/troute/routing/compute.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1108,7 +1108,7 @@ def compute_diffusive_routing(
diffusive_parameters,
waterbodies_df,
topobathy_data,
):
):

results_diffusive = []
for tw in diffusive_network_data: # <------- TODO - by-network parallel loop, here.
Expand All @@ -1134,7 +1134,7 @@ def compute_diffusive_routing(
topobathy_data_bytw = pd.DataFrame()

# diffusive streamflow DA activation switch
if da_parameter_dict['diffusive_streamflow_DA']==True:
if da_parameter_dict['diffusive_streamflow_nudging']==True:
diff_usgs_df = usgs_df
else:
diff_usgs_df = pd.DataFrame()
Expand Down
50 changes: 10 additions & 40 deletions src/troute-routing/troute/routing/diffusive_utils.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
import numpy as np
from functools import partial, reduce
import troute.nhd_network as nhd_network
from datetime import datetime, timedelta
import pandas as pd


def adj_alt1(
Expand Down Expand Up @@ -531,24 +533,21 @@ def fp_da_map(
nts_da_g -- (int) number of DA timesteps
usgs_da_g -- (float) usgs oberved streamflow data [cms]
usgs_da_reach_g -- (int) indices of stream reaches where DA is applied
"""
from datetime import datetime, timedelta
import pandas as pd

"""
nts_da_g = int((tfin_g - t0_g) * 3600.0 / dt_da_g) + 1 # include initial time 0 to the final time
usgs_da_g = -444*np.ones((nts_da_g, nrch_g))
usgs_da_g = -4444.0*np.ones((nts_da_g, nrch_g))
usgs_da_reach_g = np.zeros(nrch_g, dtype='i4')

if not usgs_df.empty:
dt_timeslice = timedelta(minutes=dt_da_g/60.0)
tfin = t0 + dt_timeslice*nsteps
timestamps = pd.date_range(t0, tfin, freq=dt_timeslice)

usgs_df_complete = usgs_df.replace(np.nan, -444)
usgs_df_complete = usgs_df.replace(np.nan, -4444.0)

for i in range(len(timestamps)):
if timestamps[i] not in usgs_df.columns:
usgs_df_complete.insert(i, timestamps[i], -444*np.ones(len(usgs_df)), allow_duplicates=False)
usgs_df_complete.insert(i, timestamps[i], -4444.0*np.ones(len(usgs_df)), allow_duplicates=False)

frj = -1
for x in range(mx_jorder, -1, -1):
Expand All @@ -560,28 +559,8 @@ def fp_da_map(
segID = seg_list[seg]
if segID in usgs_df_complete.index:
usgs_da_g[:,frj] = usgs_df_complete.loc[segID].values
usgs_da_reach_g[frj] = frj + 1 # Fortran-Python index relationship, that is Python i = Fortran i+1

'''
# test
frj= -1
with open("./output/usgs_da_g.txt",'a') as usgs_da:
for x in range(mx_jorder,-1,-1):
for head_segment, reach in ordered_reaches[x]:
seg_list= reach["segments_list"]
ncomp= reach["number_segments"]
frj= frj+1
#for seg in range(0,ncomp):
# segID= seg_list[seg]
for tsi in range (0,nts_da_g):
t_min= dt_da_g/60.0*float(tsi)
dmy1= usgs_da_g[tsi, frj]
usgs_da.write("%s %s %s\n" %\
(str(frj+1), str(t_min), str(dmy1)) )
# <- +1 is added to frj for accommodating Fortran indexing.
# test
np.savetxt("./output/usgs_da_reach_g.txt", usgs_da_reach_g, fmt='%15.5f')
'''
usgs_da_reach_g[frj] = frj + 1 # Fortran-Python index relationship, that is Python i = Fortran i+1

return nts_da_g, usgs_da_g, usgs_da_reach_g

def diffusive_input_data_v02(
Expand Down Expand Up @@ -631,7 +610,7 @@ def diffusive_input_data_v02(
# upstream boundary condition timestep (sec)
dt_ub_g = dt
# downstream boundary condition timestep (sec)
dt_db_g = 360.0 # dt * qts_subdivisions
dt_db_g = dt * qts_subdivisions
# tributary inflow timestep (sec)
dt_qtrib_g = dt
# usgs_df time step used for data assimilation. The original timestep of USGS streamflow is 15 min but later interpolated at every dt in min.
Expand Down Expand Up @@ -802,16 +781,7 @@ def diffusive_input_data_v02(

# covert data type from integer to float for frnw
dfrnw_g = frnw_g.astype('float')

np.savetxt("./output/frnw_ar.txt", frnw_g, fmt='%10i')
frj= -1
with open("./output/pynw.txt",'a') as pynwtxt:
for x in range(mx_jorder,-1,-1):
for head_segment, reach in ordered_reaches[x]:
frj= frj+1
pynwtxt.write("%s %s\n" %\
(str(frj+1), str(pynw[frj])))


# ---------------------------------------------------------------------------------
# Step 0-5
# Prepare channel geometry data
Expand Down
Loading

0 comments on commit 8d61d36

Please sign in to comment.