As presented in section 6.2.1 in Jackisch and Zehe (WRR, 2015), this is the testcase to reproduce a plot scale sprinkler experiment with the 2D version of echoRD.
Make sure to have numpy, pandas, scipy and echoRD installed.
In [1]:
import numpy as np
import pandas as pd
In [2]:
import scipy as sp
import matplotlib.pyplot as plt
import os, sys
In [3]:
#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_specht
In [4]:
obs=pd.read_csv('brspecht.dat',delimiter='\t')
obs.index=np.arange(-0.05,-1.,-0.1)
plt.figure(figsize=(18,6))
plt.subplot(131)
plt.imshow(obs.iloc[:,:10],cmap='Blues')
plt.colorbar()
plt.title('Br- Tracer Recovery\nProfile 1')
plt.subplot(132)
plt.imshow(obs.iloc[:,10:],cmap='Blues')
plt.title('Br- Tracer Recovery\nProfile 2')
plt.colorbar()
plt.subplot(133)
plt.plot(obs.iloc[:,:10].mean(axis=1),obs.index,label='Profile 1',c='b')
plt.errorbar(obs.iloc[:,:10].mean(axis=1),obs.index,xerr=obs.iloc[:,:10].std(axis=1),c='b')
plt.plot(obs.iloc[:,10:].mean(axis=1),obs.index,label='Profile 2',c='r')
plt.errorbar(obs.iloc[:,10:].mean(axis=1),obs.index,xerr=obs.iloc[:,10:].std(axis=1),c='r')
plt.legend(loc=4)
plt.xlabel('c(Br-)')
plt.ylabel('depth [m]')
plt.title('Br- Tracer Recovery [mg]\nLateral Mean with Var over Depth')
plt.tight_layout()
In [5]:
#connect to echoRD
import run_echoRD as rE
#connect and load project
[dr,mc,mcp,pdyn,cinf,vG]=rE.loadconnect(pathdir='../',mcinif='mcini_specht4y')
mcp.mcpick_out(mc,'specht4y_x2.pickle')
mc.advectref='Shipitalo'
[mc,particles,npart]=dr.particle_setup(mc)
precTS=pd.read_csv(mc.precf, sep=',',skiprows=3)
In [6]:
mc.prects=False
#theta=mc.zgrid[:,1]*0.+0.273
#[mc,particles,npart]=rE.particle_setup_obs(theta,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 [7]:
mc.prects=False
precTS.tstart=60
precTS.tend=60+2.3*3600
precTS.total=0.02543
precTS.intense=precTS.total/(precTS.tend-precTS.tstart)
precTS
Out[7]:
In [7]:
In [8]:
t_end=6.*3600.
saveDT=20
mc.FC[-1]=60
In [9]:
#1: MDA
#2: MED
infiltmeth='MDA'
#3: RWdiff
#4: Ediss
#exfiltmeth='RWdiff'
exfiltmeth='Ediss'
#5: film_uconst
#6: dynamic u
film=True
#7: maccoat1
#8: maccoat10
#9: maccoat100
macscale=1. #scale the macropore coating
clogswitch=False
In [ ]:
runname='echoRD_2D_Weiherbach_sprinkler_1457h_'
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),len(mc.zgrid[:,1])+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_specht(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
#check against br- profile
#if (i==60) | (i==120) | (i==180) | (i==240) | (i==300):
# z1=np.append(particles.loc[((particles.age>0.) & (particles.lat<=0.1)),'z'].values,obs.index)
# z2=np.append(particles.loc[((particles.age>0.) & (particles.lat>0.1) & (particles.lat<=0.2)),'z'].values,obs.index)
# z3=np.append(particles.loc[((particles.age>0.) & (particles.lat>0.2) & (particles.lat<=0.3)),'z'].values,obs.index)
# advect_dummy1=(precTS.conc.values*mc.particleA)*(np.bincount(np.ceil((z1*-10).astype(float)).astype(np.int))[1:]-1)/mc.particleA
# advect_dummy2=(precTS.conc.values*mc.particleA)*(np.bincount(np.ceil((z2*-10).astype(float)).astype(np.int))[1:]-1)/mc.particleA
# advect_dummy3=(precTS.conc.values*mc.particleA)*(np.bincount(np.ceil((z3*-10).astype(float)).astype(np.int))[1:]-1)/mc.particleA
# sim=pd.DataFrame([advect_dummy1,advect_dummy2,advect_dummy3],columns=obs.index).T
# hydroprofile(obs,sim,fsize=(3, 6),xbound=[0.,0.03,3],ybound=[-1.,0.,5],ptitle='Bromid Recovery',saving=''.join(['./results/',runname,str(i),'BR.pdf']))
[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=False,maccoat=macscale,saveDT=saveDT,clogswitch=clogswitch,infilt_method=infiltmeth,exfilt_method=exfiltmeth,film=film)
In [16]:
#pickle TSstore array - use only if a new model run was created (will overwrite file)
import cPickle as pickle
f = open(''.join([runname,'TS.pick']),"wb")
pickle.dump(pickle.dumps([pickle.dumps(TSstore),pickle.dumps(TXstore),pickle.dumps(NPstore)]), f, protocol=2)
f.close()
In [10]:
#in case there are different runs stored:
#runname='../resultsechoRD_2D_Weiherbach_1458z_'
runname='NechoRD_2D_Weiherbach_R457A_'
#This is a 26h run of the given Weiherbach setup.
#Infiltration: MDA
#Macropore-Matrix Exchange: Energy Dissipation
#Dynamic Film Flow
#No Coating.
In [11]:
#unpickle TSstore array
import cPickle as pickle
f = open(''.join(['./results/',runname,'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 [12]:
np.shape(TSstore)
#stored every minute
Out[12]:
In [13]:
import scipy.ndimage as spn
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(TSstore[1+i**3,:,:],size=mc.smooth),vmin=0., vmax=100., cmap='Blues',origin='upper')
plt.title('t='+str(i**3)+'min',fontsize=16)
#plt.subplot(k+1)
#plt.colorbar()
plt.tight_layout()
plt.savefig(''.join(['./results/',runname,'moistdev.pdf']))
In [14]:
sim=pd.DataFrame(NPstore[3*60,:]*precTS.conc.values*mc.particlemass*1000.,np.linspace(0.,mc.soildepth,len(NPstore[0,:])))
sim.columns=['simulation 3h']
sim['simulation 24h']=pd.DataFrame(NPstore[24*60,:]*precTS.conc.values*mc.particlemass*1000.,np.linspace(0.,mc.soildepth,len(NPstore[0,:])))
#hydroprofile(obs,sim,fsize=(3, 6),xbound=[0.,0.025,3],ybound=[-1.,0.,5],ptitle='Bromid Recovery\nSimulation after 6h\nObservation after 24h\n',saving=''.join(['./results/',runname,str(360),'BR.pdf']))
hydroprofile(obs,sim,fsize=(3, 6),xbound=[0.,0.03,3],ybound=[-1.,0.,5],ptitle='Bromid Recovery\n',xlab='c(Br-) [mg/l]',saving=''.join(['./results/',runname,str(360),'BR.pdf']))
In [60]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
ts=[61,121,241,1441]
k=0
fig, axes = plt.subplots(nrows=1, ncols=4, sharex=True, sharey=True)
for ax in axes.flat:
dummy=spn.filters.median_filter(TSstore[ts[k],:,:]-TSstore[1,:,:],size=mc.smooth)
im = ax.imshow(dummy, vmin=0, vmax=50., cmap='Blues')
ax.set_ylim(100,0)
k+=1
cax,kw = mpl.colorbar.make_axes([ax for ax in axes.flat])
plt.colorbar(im, cax=cax, **kw)
plt.show()
#plt.savefig('./results/weiherbach_diffmoist2.pdf')
In [59]:
fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(151)
ax1.imshow(spn.filters.median_filter(TSstore[61,:,:]-TSstore[1,:,:],size=mc.smooth),vmin=0., vmax=50., cmap='Blues',origin='upper')
ax1.set_title('t=60min')
ax2 = fig.add_subplot(152)
ax2.imshow(spn.filters.median_filter(TSstore[121,:,:]-TSstore[1,:,:],size=mc.smooth),vmin=0., vmax=50., cmap='Blues',origin='upper')
ax2.set_title('t=120min')
ax3 = fig.add_subplot(153)
ax3.imshow(spn.filters.median_filter(TSstore[241,:,:]-TSstore[1,:,:],size=mc.smooth),vmin=0., vmax=50., cmap='Blues',origin='upper')
ax3.set_title('t=240min')
ax4 = fig.add_subplot(154)
ax4.imshow(spn.filters.median_filter(TSstore[1441,:,:]-TSstore[1,:,:],size=mc.smooth),vmin=0., vmax=50., cmap='Blues',origin='upper')
ax4.set_title('t=24h')
ylabels = ['0','0','-0.2','-0.4','-0.6','-0.8']
ax1.set_yticklabels(ylabels)
plt.setp(ax2.get_yticklabels(), visible=False)
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)
ax1.set_xticks([0,10,20,30])
ax2.set_xticks([0,10,20,30])
ax3.set_xticks([0,10,20,30])
ax4.set_xticks([0,10,20,30])
xlabels = ['0','0.1','0.2','0.3']
ax1.set_xticklabels(xlabels)
ax2.set_xticklabels(xlabels)
ax3.set_xticklabels(xlabels)
ax4.set_xticklabels(xlabels)
plt.tight_layout()
#plt.savefig('./results/weiherbach_diffmoist.pdf')
In [ ]: