In [1]:
import numpy as np
import pandas as pd
In [2]:
import scipy as sp
import matplotlib.pyplot as plt
import os, sys
import cPickle as pickle
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_colpach
In [4]:
obs=pd.read_csv('brprofileXI.dat',delimiter='\t')
obs.index=np.arange(-0.025,-0.74,-0.05)
plt.figure(figsize=(12,6))
plt.subplot(121)
plt.imshow(obs,cmap='Blues')
plt.xlabel('Grid a 5cm')
plt.ylabel('Grid a 5cm')
plt.colorbar()
plt.title('Br- Tracer Recovery\nProfile 1')
plt.subplot(122)
plt.plot(obs.mean(axis=1),obs.index,label='Profile 1')
plt.errorbar(obs.mean(axis=1),obs.index,xerr=obs.std(axis=1),label='Profile 1')
xlabel('c(Br-)')
ylabel('depth [m]')
plt.title('Br- Tracer Recovery\nLateral Sum over Depth')
plt.tight_layout()
In [5]:
soilsXI=pd.read_csv('./XI_soilsamples.csv',delimiter=',')
soilsXI['rDepth']=np.round(soilsXI.Depth*10.)/10.
soilrefXI=soilsXI.groupby(['rDepth']).mean()[['Ks','ths','thr','alpha','n','BD','Depth']]
soilrefXI.columns=['ks','ts','tr','alpha','n','roh','z']
soilrefXI.index=np.arange(5)+1
soilrefXI.z=-soilrefXI.z
soilrefXI.alpha=soilrefXI.alpha
soilrefXI.to_csv('./matrix_XI.dat',sep=' ',index_label='no',float_format='%.4f')
soilrefXI
Out[5]:
In [5]:
In [6]:
#connect to echoRD
import run_echoRD as rE
#connect and load project
[dr,mc,mcp,pdyn,cinf,vG]=rE.loadconnect(pathdir='../',mcinif='mcini_XI')
#mc=rE.preproc_echoRD(mc, dr, mcp,'XI_new.pickle')
mcp.mcpick_out(mc,'XI_new.pickle')
mc.advectref='Shipitalo'
[mc,particles,npart]=dr.particle_setup(mc)
precTS=pd.read_csv(mc.precf, sep=',',skiprows=3)
In [7]:
mc.prects=False
[thS,npart]=pdyn.gridupdate_thS(particles.lat,particles.z,mc)
theta=thS[:,1]/100.
theta[theta<0.22]=0.22
theta[:30]+=0.15
#update particles to be correct
[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 [8]:
np.shape(mc.md_contact)[1]-1
Out[8]:
In [9]:
precTS.tstart-=340
precTS.tend-=340
precTS.intense*=10.
precTS
Out[9]:
In [10]:
#tracer per particle
n_part_tot=precTS.intense.values*3600.*mc.mgrid.width.values/mc.particleA
tracer_part=mc.tracer_appl_Br*mc.particleD/n_part_tot
In [11]:
t_end=6.*3600.
saveDT=0.5
clogswitch=False
#1: MDA
#2: MED
infiltmeth='MED'
#3: RWdiff
#4: Ediss
exfiltmeth='Ediss'
#5: film_uconst
#6: dynamic u
film=True
#7: maccoat1
#8: maccoat10
#9: maccoat100
macscale=1000. #scale the macropore coating [1...inf] - large means high coating, less infiltration
infiltscale=0.12 #scale the particles provided to the matrix domain at top boundary [1...0] - small means more water to macropores
In [ ]:
runname='echoRD_2D_Colpach_sprinkler_infk1_'
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_colpach(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=False,maccoat=macscale,saveDT=saveDT,clogswitch=clogswitch,infilt_method=infiltmeth,exfilt_method=exfiltmeth,film=film,infiltscale=infiltscale)
if i/10.==np.round(i/10.):
f = open(''.join(['./results/N',runname,'TS.pick']),"wb")
pickle.dump(pickle.dumps([pickle.dumps(TSstore),pickle.dumps(TXstore),pickle.dumps(NPstore)]), f, protocol=2)
f.close()
f = open(''.join(['./results/N',runname,'_particles.pick']),"wb")
pickle.dump(pickle.dumps(pickle.dumps(particles)), f, protocol=2)
f.close()
In [12]:
mc.a_velocity_real
Out[12]:
In [9]:
#pickle TSstore array - use only if a new model run was created (will overwrite file)
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 [13]:
#in case there are different runs stored:
#runname='NechoRD_2D_Colpach_R457x_'
#runname='NechoRD_2D_Colpach_Vsk2x_'
runname='NechoRD_2D_Colpach_Vsk5x_'
#runname='NechoRD_2D_Colpach_VsU7x_'
#This is a 26h run of the given Colpach setup.
#Infiltration: MDA
#Macropore-Matrix Exchange: Energy Dissipation
#Dynamic Film Flow
#No Coating.
In [14]:
#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 [15]:
np.shape(TSstore)
#stored every minute
Out[15]:
In [17]:
import scipy.ndimage as spn
t_inst=[1,6,11,21,41,61,121,181,601]
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[t_inst[i],:,:],size=mc.smooth),vmin=0., vmax=100., cmap='Blues',origin='upper')
plt.title('t='+str(t_inst[i]-1)+'min',fontsize=16)
#plt.subplot(k+1)
#plt.colorbar()
plt.tight_layout()
#plt.savefig(''.join(['./results/',runname,'moistdev.pdf']))
In [20]:
sim=pd.DataFrame(NPstore[121,:]*precTS.conc.values*mc.particlemass*1000.,np.linspace(0.,mc.soildepth,len(NPstore[0,:])))
sim.columns=['simulation 2h']
sim['simulation 1h']=pd.DataFrame(NPstore[61,:]*precTS.conc.values*mc.particlemass*1000.,np.linspace(0.,mc.soildepth,len(NPstore[0,:])))
#sim['simulation 3h']=pd.DataFrame(NPstore[181,:]*precTS.conc.values*mc.particlemass*1000.,np.linspace(0.,mc.soildepth,len(NPstore[0,:])))
#sim['simulation 40min']=pd.DataFrame(NPstore[41,:]*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*50000.,fsize=(4, 10),xbound=[0.,600.,5],ybound=[-1.2,0.,5],ptitle='Bromide Recovery\n',xlab='c(Br-) [mg/l]',saving='colpach_sx5.pdf')
In [ ]:
In [ ]: