In [1]:
import numpy as np
import pandas as pd
In [2]:
import os,sys
#connect echoRD Tools
pathdir='../' #path to echoRD
lib_path = os.path.abspath(pathdir)
sys.path.append(lib_path)
import vG_conv as vG
from hydro_tools import plotparticles_t,hydroprofile,plotparticles_column
In [3]:
Q63=pd.read_csv('./ref_Q63/Q63_data.csv')
#index set to time
#Q63.index=pd.to_datetime(Q63['Datum Uhrzeit'])
#index set to experiment
Q63.index=Q63['minute 0=start']
In [4]:
#clean erroneous values
Q63hc=Q63.iloc[:,11:-2]
Q63hc[Q63>1000.]=np.nan
#tension is given in cm
In [5]:
#define tens -> theta conversion
def phi_theta(tens):
#calculate theta from tension in cm
th_s=0.336
th_r=0.071
alpha=0.015
n=15
m=0.93
return th_r+((th_s-th_r)/(1+(alpha*tens)**n)**m)
def phi_rtheta(tens):
#calculate relative saturation from tension in cm
th_s=0.336
th_r=0.071
alpha=0.015
n=15
m=0.93
return (1/(1+(alpha*tens)**n)**m)
In [6]:
def ref_plot_t(tens,ti):
from scipy.interpolate import griddata
import matplotlib.pylab as plt
#positions of tensiometers
z_tens=(110.-np.array([97.5,82.5,67.5,52.5,37.5,22.5]))/100.
x_tens=np.array([2.,6.,14.,30.])/100.
#grid setup
x_mx=x_tens.repeat(len(z_tens)).reshape((len(x_tens),len(z_tens))).T
x_mx=np.append(x_mx,-x_mx)
z_mx=-z_tens[::-1].repeat(len(x_tens)).reshape((len(z_tens),len(x_tens)))
z_mx=np.append(z_mx,z_mx)
N = 100j
extent = (-0.5,0.5,-1.1,0)
xs,ys = np.mgrid[extent[0]:extent[1]:N, extent[2]:extent[3]:N]
#interpolation
resampled = griddata((x_mx, z_mx), np.tile(phi_rtheta(tens),(2, 1)).ravel(), (xs, ys),method='cubic')
#plot
fig, ax = plt.subplots()
#im = ax.imshow(resampled.T, extent=extent, origin='lower',vmin=0., vmax=1.,cmap='Blues')
con = ax.contour(resampled.T, extent=extent, levels=[0.2,0.5,0.8],origin='lower',vmin=0., vmax=1.,cmap='Blues')#colors='k')
ax.clabel(con, fontsize=9, inline=1)
ax.scatter(x_mx.ravel(), z_mx.ravel(), c=np.tile(Q63hc.iloc[3,:].values,(2, 1)).ravel(), s=10,
vmin=0., vmax=1.)
ax.set_title(''.join(['Theta* Observed at min=',str(ti)]))
ax.set_aspect('equal')
#fig.colorbar(con)
plt.show()
# plt.savefig(''.join(['Q63_c',str(int(round(ti))),'.png']))
# plt.close(fig)
In [7]:
tt=200. #time of experiment in minutes
idx=np.where(Q63hc.index>tt)[0]
ref_plot_t(Q63hc.iloc[idx[0],:],tt)
In [8]:
#read break through data
g63_flow=pd.read_csv('./ref_Q63/g63_flow.dat', delimiter='\t', na_values='NA')
g63_flow.t=g63_flow.t_day.values*(24.*60.*60.)
g63_flow.index=g63_flow.t.astype(int)
#translate into recovery profile
P_breakthrough=g63_flow.discharge_l/np.amax(g63_flow.discharge_l)
#corresponding velocity of fast hump
ref=np.append(0.001,1./g63_flow.index.values[1:]) #m/s
In [9]:
#slow hump not drained
#take observed moisture profile
z=np.array([0.025,0.175,0.325,0.475,0.625,0.775])
theta=np.array([0.335999585,0.34,0.31,0.21,0.119014148,0.126440703])
#plot(theta,-z)
#thus the velocity distribution inside the column (referring to theta=0.07 at start)
dtheta=theta-0.07
ptheta=np.cumsum(dtheta/(sum(dtheta)))
ref_slow=z/(24.*3600.)
In [10]:
#combine the two distributions through the ratio of water drained and retained
ratio=g63_flow.discharge_l.iloc[-1]/g63_flow.inflow_l.iloc[-1]
z_r=z*(1-ratio)
adv=(1-ratio)+P_breakthrough*ratio
basis=np.append(z_r,adv.values[4:]) #dropped zero values to avoid inconsistency
basis_v=np.append(ref_slow,ref[4:])
In [11]:
n=1000
part=np.random.random(n)
p_v=np.arange(n)*0.
for i in np.arange(n):
idx=np.where(basis>=part[i])[0][0]
weight=(part[i]-basis[idx-1])/np.diff(basis[[idx-1,idx]])
p_v[i]=basis_v[idx-1]*weight + basis_v[idx]*(1-weight)
p_bin=np.append(np.arange(20),(p_v*(6*3600)*20).astype(int))
fill=np.bincount(p_bin)
plot((fill-1)/1000.,-np.arange(20)/20.)
Out[11]:
In [12]:
#now write out the concentrations to be used as caos_py input
(fill-1)/1000.
Out[12]:
In [13]:
import run_echoRD as rE
[dr,mc,mcp,pdyn,cinf,vG]=rE.loadconnect(pathdir='../',mcinif='mcini_g63')
#mc=rE.preproc_echoRD(mc, dr, mcp,'XI_new.pickle')
mcp.mcpick_out(mc,'mc_dump_g63ax.pickle')
mc.advectref='Shipitalo'
[mc,particles,npart]=dr.particle_setup(mc)
precTS=pd.read_csv(mc.precf, sep=',',skiprows=3)
In [14]:
runname='echoRD_column'
maccoat=1.
saveDT=30
#select if column geometry is referred to or not:
#mc.prects='column'
mc.prects=False
mc.colref=False
precTS.tstart=60 #modify start time to get results more quickly
#precTS.intense=0.063*60./1000.# intensity in m3/s
#precTS.intense*=2.*mc.particleD/(0.5*np.pi*0.5**2)
#precTS.intense*=0.001 #conversion BUT WHY
precTS.intense*=5
precTS
Out[14]:
In [14]:
In [15]:
#define macropore as coarse sand
mc.soilgrid[:,97:100]=13
imshow(mc.soilgrid)
Out[15]:
In [16]:
#update matrix definitions for "macropore"
mc.soilmatrix.loc[13]=[14,-mc.a_velocity_real[-1]*10.,0.3,0.01,0.01,27.,1500,-2.3,'g63_mac',1.-1./27.]
mc.soilmatrix
Out[16]:
In [17]:
import os, sys
#connect echoRD Tools
pathdir='../' #path to echoRD
lib_path = os.path.abspath(pathdir)
sys.path.append(lib_path)
import vG_conv as vG
from hydro_tools import plotparticles_t
#[mc,particles,npart]=rE.particle_setup_obs((mc.soilgrid/mc.soilgrid)[:,0]*0.1,mc,vG,dr,pdyn)
[thS,npart]=pdyn.gridupdate_thS(particles.lat,particles.z,mc)
#[A,B]=plotparticles_t(particles,thS/100.,mc,vG,store=True)
In [18]:
#1: MDA
#2: MED
infiltmeth='MDA'
#infiltmeth='MED'
#3: RWdiff
#4: Ediss
#exfiltmeth='RWdiff'
exfiltmeth='Ediss'
#5: film_uconst
#6: dynamic u
film=True
#film=False
#7: maccoat1
#8: maccoat10
#9: maccoat100
macscale=1. #scale the macropore coating
This will start the model. Be aware that the setup runs with a large number of particles on a very fine grid with great moisture contrasts. Depending on your machine it may run for more days. An alternative brief dummy test is given as comment below. This runs about 30 s.
Alternatively, you may use the prepared runfile run_g63b.py on a powerful machine. This is prepared to require only the basic packages numpy, pandas, scipy and cPickle. However intended for HPC it is not parallelized (yet) and will also run for quite a while. The resulting TSstore is pickled and can be used as given below.
In [19]:
clogswitch=False
t_end=3600.
drained=pd.DataFrame(np.array([]))
leftover=0
output=60. #mind to set also in TXstore.index definition
dummy=np.floor(t_end/output)
#prepare store arrays
TSstore=np.zeros((int(dummy),np.shape(thS)[0],np.shape(thS)[1]))
NPstore=np.zeros((int(dummy),int(np.abs(mc.soildepth*100)+1)))
colnames=['part1','part2','part3','part4','part5']
TXstore=pd.DataFrame(np.zeros((int(dummy),len(colnames))),columns=colnames)
t=0.
#loop through plot cycles
for i in np.arange(dummy.astype(int)):
#plot and store states
plotparticles_column(particles,mc,pdyn,vG,runname,t,i,saving=True,relative=False)
[TXstore.iloc[i],NPstore[i,:]]=plotparticles_t(particles,thS/100.,mc,vG,runname,t,i,saving=True,store=True)
#store thS
TSstore[i,:,:]=thS
#[particles,npart,thS,leftover,drained,t]=rE.CAOSpy_rundx(i*output,(i+1)*output,mc,pdyn,cinf,precTS,particles,leftover,drained,6.,splitfac=4,prec_2D=True,maccoat=macscale,saveDT=saveDT,clogswitch=clogswitch,infilt_method=infiltmeth,exfilt_method=exfiltmeth,film=film)
[particles,npart,thS,leftover,drained,t]=rE.CAOSpy_rund_diffonly(i*output,(i+1)*output,mc,pdyn,cinf,precTS,particles,leftover,drained,6.,splitfac=4,prec_2D=True,saveDT=saveDT)
In [ ]:
#!ffmpeg -r 2 -i ./results/echoRD_g63a_diff_TESTagain_coat_t_%03d.png -c:v libx264 -r 5 -pix_fmt yuv420p ./results/echoRD_g63a_diff_TESTagain_coat.mp4
In [19]:
#imshow(TSstore[10,:,:])
#colorbar()
In [23]:
#pickle TSstore array - use only if a new model run was created (will overwrite file)
import cPickle as pickle
f = open(''.join(['./results',runname,'TS.pick']),"wb")
pickle.dump(pickle.dumps([pickle.dumps(TSstore),pickle.dumps(TXstore),pickle.dumps(NPstore)]), f, protocol=2)
f.close()
In [24]:
#load prepared data set with plane geometry
runname='NechoRD_2D_Column_AdQ5x_'
# --> result withradial symmetric geometry below
In [25]:
#unpickle TSstore array
import cPickle as pickle
f = open(''.join(['./results/',runname,'TS.pick']),"rb")
#f = open('./resultsechoRD_g63a_diff_TS.pick',"rb")
[TS_store,TX_store,NP_store] = pickle.loads(pickle.load(f))
f.close()
TSstore=pickle.loads(TS_store)
TXstore=pickle.loads(TX_store)
NPstore=pickle.loads(NP_store)
In [41]:
#plot definition to qualitatively compare experiment with simulation
def comp_plot_t(tens,sim,ti,imfile=False,fsize=(10,11),saving=False,extent2=-0.5):
from scipy.interpolate import griddata
import matplotlib.pylab as plt
import scipy.ndimage as spn
#positions of tensiometers
z_tens=(110.-np.array([97.5,82.5,67.5,52.5,37.5,22.5]))/100.
x_tens=np.array([2.,6.,14.,30.])/100.
#get reference data
idx=np.where(tens.index>ti)[0]
tens=tens.iloc[idx[0],:]
#grid setup
x_mx=x_tens.repeat(len(z_tens)).reshape((len(x_tens),len(z_tens))).T
x_mx=np.append(x_mx,-x_mx)
z_mx=-z_tens[::-1].repeat(len(x_tens)).reshape((len(z_tens),len(x_tens)))
z_mx=np.append(z_mx,z_mx)
N = 100j
extent = (-0.5,0.5,-1.,0)
#extent = (-0.5,0.5,-.5,0)
xs,ys = np.mgrid[extent[0]:extent[1]:N, extent[2]:extent[3]:N]
#interpolation
resampled = griddata((x_mx, z_mx), np.tile(phi_rtheta(tens),(2, 1)).ravel(), (xs, ys),method='cubic')
#plot
plt.figure(figsize=fsize)
fig, ax = plt.subplots()
TSsmooth=spn.gaussian_filter(sim[int(ti),:,:]/50.,2.)
im = ax.imshow(TSsmooth, extent=extent, origin='upper',vmin=0.1, vmax=0.8,cmap='Blues')
#con1 = ax.contour(TSsmooth, extent=extent,levels=[0.1,0.2,0.3,0.5,0.8],origin='upper',vmin=0.1, vmax=0.9,cmap='YlGn')#colors='k')
if imfile:
ima=plt.imread(imfile)
#im1 = ax.imshow(ima[:,500:],extent=(0., 0.5, -1.09, 0.01), origin='upper')
im1 = ax.imshow(ima[:int(extent2*-1000),500:],extent=(0., 0.5, extent2, 0.01), origin='upper')
#ax.text(-0.3,-1.05,'Model')
#ax.text(0.1,-1.05,'Obs Photo')
ax.text(-0.3,extent2+0.03,'Model',fontsize=14)
ax.text(0.1,extent2+0.03,'Obs Photo',fontsize=14)
con = ax.contour(resampled.T, extent=extent, levels=[0.2,0.5,0.8],origin='lower',vmin=0., vmax=1.,cmap='Blues')#colors='k')
ax.clabel(con, fontsize=14, inline=1)
ax.scatter(x_mx.ravel(), z_mx.ravel(), c=np.tile(Q63hc.iloc[3,:].values,(2, 1)).ravel(), s=10,
vmin=0., vmax=1.,alpha=0.5, linewidths=0.1)
ax.set_title(''.join(['Theta* Model vs. Observed\ntime: ',str(ti-1),'min']))
ax.set_aspect('equal')
ax.set_ylim(extent2,0.)
ax.set_xlabel('width [m]')
ax.set_ylabel('depth [m]')
plt.rcParams.update({'font.size': 14})
plt.rcParams['figure.figsize'] = 5, 5
#fig.colorbar(im)
#plt.show()
if saving:
plt.savefig(saving)
#plt.close(fig)
In [27]:
#the plot function assumes that TSstore is recorded every minute.
imfile=False
#imfile='./ref_Q63/Q63_5.png'
#comp_plot_t(Q63hc,TSstore,5.,imfile)
imfile='./ref_Q63/Q63_10.png'
comp_plot_t(Q63hc,TSstore,11.,imfile,saving='./results/col10_plane.pdf')
imfile='./ref_Q63/Q63_20.png'
comp_plot_t(Q63hc,TSstore,21.,imfile,saving='./results/col20_plane.pdf')
imfile='./ref_Q63/Q63_50.png'
comp_plot_t(Q63hc,TSstore,51.,imfile,saving='./results/col50_plane.pdf')
In [28]:
import scipy.ndimage as spn
plt.figure(figsize=(18,6))
allplt=181
for i in np.arange(8):
k=allplt+i
plt.subplot(k)
plt.imshow(spn.filters.median_filter(TSstore[1+i**2,:,:],size=mc.smooth),vmin=0., vmax=50., cmap='Blues',origin='upper')
plt.title('t='+str(i**2)+'min',fontsize=16)
#plt.subplot(k+1)
#plt.colorbar()
plt.tight_layout()
In [29]:
runname='NechoRD_2D_Column_AdP4x_'
#load prepared data set with plane geometry
In [30]:
#unpickle TSstore array
import cPickle as pickle
f = open(''.join(['./results/',runname,'TS.pick']),"rb")
#f = open('./resultsechoRD_g63a_diff_TS.pick',"rb")
[TS_store,TX_store,NP_store] = pickle.loads(pickle.load(f))
f.close()
TSstore2=pickle.loads(TS_store)
TXstore2=pickle.loads(TX_store)
NPstore2=pickle.loads(NP_store)
In [31]:
shape(TSstore2)
Out[31]:
In [43]:
#the plot function assumes that TSstore is recorded every minute.
imfile=False
#imfile='./ref_Q63/Q63_5.png'
#comp_plot_t(Q63hc,TSstore,5.,imfile)
imfile='./ref_Q63/Q63_10.png'
comp_plot_t(Q63hc,TSstore2,11.,imfile,saving='./results/col10_col.pdf')
imfile='./ref_Q63/Q63_20.png'
comp_plot_t(Q63hc,TSstore2,21.,imfile,saving='./results/col20_col.pdf')
imfile='./ref_Q63/Q63_50.png'
comp_plot_t(Q63hc,TSstore2,51.,imfile,saving='./results/col50_col.pdf')
imfile='./ref_Q63/Q63_200.png'
comp_plot_t(Q63hc,TSstore2,201.,imfile,saving='./results/col200_col.pdf',extent2=-0.8)
In [68]:
plt.figure(figsize=(18,6))
allplt=191
for i in np.arange(9):
k=allplt+i
plt.subplot(k)
plt.imshow(spn.filters.median_filter(TSstore2[1+i**2.5,:,:],size=mc.smooth),vmin=0., vmax=50., cmap='Blues',origin='upper')
plt.title('t='+str(int(np.round(i**2.5)))+'min',fontsize=16)
#plt.subplot(k+1)
#plt.colorbar()
plt.tight_layout()
In [ ]: