Skip to content

Commit

Permalink
#29: Added bracket operator to Waterstorage to access solute storage
Browse files Browse the repository at this point in the history
  • Loading branch information
philippkraft committed May 8, 2018
1 parent 84d6dbf commit a0f9b57
Show file tree
Hide file tree
Showing 10 changed files with 219 additions and 142 deletions.
2 changes: 1 addition & 1 deletion cmf/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
from .stopwatch import StopWatch


__version__ = '1.3.1a0.i28_numperf_reaches_tracer'
__version__ = '1.3.1a0.i29_bracket_solute'

from .cmf_core import connect_cells_with_flux as __ccwf

Expand Down
35 changes: 9 additions & 26 deletions cmf/cmf_core.py
Original file line number Diff line number Diff line change
Expand Up @@ -3328,36 +3328,13 @@ def set_adsorption(self, *args, **kwargs):
decay = _swig_property(_cmf_core.SoluteStorage_decay_get, _cmf_core.SoluteStorage_decay_set)
source = _swig_property(_cmf_core.SoluteStorage_source_get, _cmf_core.SoluteStorage_source_set)
Solute = _swig_property(_cmf_core.SoluteStorage_Solute_get)
conc = _swig_property(_cmf_core.SoluteStorage_conc_get, _cmf_core.SoluteStorage_conc_set)

def conc(self, *args, **kwargs):
"""
conc(SoluteStorage self) -> real
real conc()
const
Returns the concentration of the solute.
"""
return _cmf_core.SoluteStorage_conc(self, *args, **kwargs)


def set_conc(self, *args, **kwargs):
"""
set_conc(SoluteStorage self, real NewConcentration)
void
set_conc(real NewConcentration)
set a new concentration of dissolved tracers.
In case of adsorption functions, the isotherm is used
"""
return _cmf_core.SoluteStorage_set_conc(self, *args, **kwargs)
def __repr__(self):
return self.to_string()

__swig_destroy__ = _cmf_core.delete_SoluteStorage
SoluteStorage.set_adsorption = new_instancemethod(_cmf_core.SoluteStorage_set_adsorption, None, SoluteStorage)
SoluteStorage.conc = new_instancemethod(_cmf_core.SoluteStorage_conc, None, SoluteStorage)
SoluteStorage.set_conc = new_instancemethod(_cmf_core.SoluteStorage_set_conc, None, SoluteStorage)
_cmf_core.SoluteStorage_swigregister(SoluteStorage)
# SoluteStorage end

Expand Down Expand Up @@ -4405,9 +4382,15 @@ def create(*args, **kwargs):
def __repr__(self):
return self.to_string()


def __getitem__(self, *args, **kwargs):
"""__getitem__(WaterStorage self, solute X) -> SoluteStorage"""
return _cmf_core.WaterStorage___getitem__(self, *args, **kwargs)

__swig_destroy__ = _cmf_core.delete_WaterStorage
WaterStorage.Solute = new_instancemethod(_cmf_core.WaterStorage_Solute, None, WaterStorage)
WaterStorage.conc = new_instancemethod(_cmf_core.WaterStorage_conc, None, WaterStorage)
WaterStorage.__getitem__ = new_instancemethod(_cmf_core.WaterStorage___getitem__, None, WaterStorage)
_cmf_core.WaterStorage_swigregister(WaterStorage)
# WaterStorage end

Expand Down
204 changes: 141 additions & 63 deletions cmf/cmf_core_src/cmf_wrap.cpp

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion cmf/cmf_core_src/water/SoluteStorage.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ real SoluteStorage::dxdt( const cmf::math::Time& time )
return inflow + outflow + source_term - decay_term;
}

real SoluteStorage::conc() const
real SoluteStorage::get_conc() const
{
real V = m_water->get_volume();
real c = this->adsorption->freesolute(get_state(),V) /V;
Expand Down
2 changes: 1 addition & 1 deletion cmf/cmf_core_src/water/SoluteStorage.h
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ namespace cmf {
/// @brief The solute, which is stored in this
const cmf::water::solute& Solute;
/// @brief Returns the concentration of the solute
real conc() const;
real get_conc() const;
/// @brief set a new concentration of dissolved tracers.
///
/// In case of adsorption functions, the isotherm is used
Expand Down
4 changes: 2 additions & 2 deletions cmf/cmf_core_src/water/WaterStorage.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,15 +40,15 @@ WaterStorage::WaterStorage(cmf::project& _project, const std::string& _Name,
}


SoluteStorage& WaterStorage::Solute( const solute& _Solute )
SoluteStorage& WaterStorage::Solute( const solute _Solute )
{
return *m_Concentrations[_Solute.Id];
}


real WaterStorage::conc(const solute& _Solute) const
{
return Solute(_Solute).conc();
return Solute(_Solute).get_conc();
}

std::shared_ptr<WaterStorage> WaterStorage::from_node( flux_node::ptr node )
Expand Down
15 changes: 13 additions & 2 deletions cmf/cmf_core_src/water/WaterStorage.h
Original file line number Diff line number Diff line change
Expand Up @@ -85,8 +85,19 @@ namespace cmf {

static std::shared_ptr<WaterStorage> from_node(cmf::water::flux_node::ptr node);
/// @brief Returns the water quality of the water storage.
SoluteStorage& Solute(const cmf::water::solute& _Solute);
const SoluteStorage& Solute(const cmf::water::solute& _Solute) const {return *m_Concentrations[_Solute.Id];}
SoluteStorage& Solute(const cmf::water::solute _Solute);
const SoluteStorage& Solute(const cmf::water::solute _Solute) const {return *m_Concentrations[_Solute.Id];}
#ifndef SWIG
SoluteStorage& operator[](cmf::water::solute _Solute) {
return *m_Concentrations[_Solute.Id];
}
const SoluteStorage& operator[](cmf::water::solute _Solute) const {
return *m_Concentrations[_Solute.Id];
}

#endif // !SWIG


/// @brief Returns the concentration of the given solute
real conc(const cmf::water::solute& _Solute) const;
/// @brief Returns the current WaterQuality (concentration of all solutes)
Expand Down
9 changes: 9 additions & 0 deletions cmf/cmf_core_src/water/water.i
Original file line number Diff line number Diff line change
Expand Up @@ -77,8 +77,11 @@
}
}

%attribute(cmf::water::SoluteStorage, real, conc, get_conc, set_conc);

%include "water/SoluteStorage.h"
%extend__repr__(cmf::water::SoluteStorage);


// flux_connection and Node
%feature("ref") cmf::water::flux_connection ""
Expand Down Expand Up @@ -180,6 +183,12 @@ namespace cmf{namespace water {class flux_connection;}}
}}
%extend__repr__(cmf::water::WaterStorage)

%extend cmf::water::WaterStorage {
cmf::water::SoluteStorage& __getitem__(cmf::water::solute X) {
return (*$self)[X];
}
}

%include "water/simple_connections.h"

%include "water/collections.i"
Expand Down
30 changes: 15 additions & 15 deletions demo/adsorption.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
import cmf
# project with three tracers
p = cmf.project('No Linear Freundlich Langmuir')
N,X,Y,Z = p.solutes
N, X, Y, Z = p.solutes

# A water storage without specific properties
ws = p.NewStorage('storage',0,0,0)
Expand All @@ -17,41 +17,41 @@

# The storage should contain 1 unit of every tracer for the beginning
for s in p.solutes:
ws.conc(s,1.0)
ws[s].conc = 1.0

# Tracer N does not adsorb, nothing to change
# Tracer X has a linear isotherm xa/m=Kc, with K = 1 and sorbent mass m = 1
ws.Solute(X).set_adsorption(cmf.LinearAdsorption(1,1))
ws[X].set_adsorption(cmf.LinearAdsorption(1, 1))
# Tracer Y has a Freundlich isotherm xa/m=Kc^n,
# with K = 1 and n=0.5 and sorbent mass m = 1
ws.Solute(Y).set_adsorption(cmf.FreundlichAdsorbtion(1.,1.,1.0))
ws[Y].set_adsorption(cmf.FreundlichAdsorbtion(1., 1., 1.0))
# Tracer Y has a Langmuir isotherm xa/m=Kc/(1+Kc),
# with K = 1 and sorbent mass m = 1
ws.Solute(Z).set_adsorption(cmf.LangmuirAdsorption(1.,1.))
ws[Z].set_adsorption(cmf.LangmuirAdsorption(1., 1.))

# Now we put a constant clean water flux into the storage
inflow=cmf.NeumannBoundary.create(ws)
inflow = cmf.NeumannBoundary.create(ws)
inflow.flux = 1.0 # 1 m3/day
for s in p.solutes:
inflow.concentration[s] = 0.0
# And an outlet, a linear storage term with a retention time of 1 day
outlet = p.NewOutlet('out',0,0,0)
cmf.kinematic_wave(ws,outlet,1.)
outlet = p.NewOutlet('out', 0, 0, 0)
cmf.LinearStorageConnection(ws, outlet, 1.)

# Create a solver
solver = cmf.ExplicitEuler_fixed(p)

# Rinse the storage for 1 week and get the outlet concentration at every hour
# and the remaining tracer in the storage
result = [ [outlet.conc(t,s) for s in p.solutes]
+ [ws.Solute(s).state for s in p.solutes]
for t in solver.run(solver.t,2*cmf.week,cmf.h)
]
result = [[outlet.conc(t, s) for s in p.solutes]
+ [ws[s].state for s in p.solutes]
for t in solver.run(solver.t, 2 * cmf.week, cmf.h)
]
# Plot result ]
from pylab import *
result= array(result)
conc=result[:,:len(p.solutes)]
load = result[:,len(p.solutes):]
result = array(result)
conc = result[:, :len(p.solutes)]
load = result[:, len(p.solutes):]
subplot(211)
plot(conc)
grid()
Expand Down
58 changes: 27 additions & 31 deletions demo/adsorption_richards.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,44 +8,38 @@
import cmf
# project with three tracers
p = cmf.project('No Linear Freundlich Langmuir')
N,X,Y,Z = p.solutes
N, X, Y, Z = p.solutes

cmf.set_parallel_threads(1)
# Create a single cell c with a surfacewater storage, which references 3 solute storages
c = p.NewCell(0,0,0,1000,with_surfacewater = True)
c = p.NewCell(0, 0, 0, 1000, with_surfacewater=True)
# Create 50 layers with 2cm thickness
for i in range(50):
# Add a layer. Each layer will reference 3 solute storages
c.add_layer((i+1)*0.02, cmf.VanGenuchtenMualem())
l = c.layers[-1]
l = c.add_layer((i+1)*0.02, cmf.VanGenuchtenMualem())
# Tracer N does not adsorb, nothing to change
l.Solute(N).state = 1.
l[N].state = 1.

# Tracer X has a linear isotherm xa/m=Kc, with K = 1 and sorbent mass m = 1
l.Solute(X).set_adsorption(cmf.LinearAdsorption(1,5))
l.Solute(X).state = 1.
l[X].set_adsorption(cmf.LinearAdsorption(1,5))
l[X].state = 1.
# Tracer Y has a Freundlich isotherm xa/m=Kc^n,
# with K = 1 and n=0.5 and sorbent mass m = 1
c.layers[-1].Solute(Y).set_adsorption(cmf.FreundlichAdsorbtion(2.,.5,1.0,1e-9))
l.Solute(Y).state = 1.
l[Y].set_adsorption(cmf.FreundlichAdsorbtion(2., .5, 1.0, 1e-9))
l[Y].state = 1.

# Tracer Y has a Langmuir isotherm xa/m=Kc/(1+Kc),
# with K = 1 and sorbent mass m = 1
l.Solute(Z).set_adsorption(cmf.LangmuirAdsorption(1.,5.))
l.Solute(Z).state = 1.
l[Z].set_adsorption(cmf.LangmuirAdsorption(1.,5.))
l[Z].state = 1.

# Use Richards equation
c.install_connection(cmf.Richards)
# Make a groundwater boundary condition
gw = p.NewOutlet('gw',0,0,-1.5)
cmf.Richards(c.layers[-1],gw)
gw = p.NewOutlet('gw', 0, 0, -1.5)
cmf.Richards(c.layers[-1], gw)

# Template for the water solver
wsolver = cmf.CVodeIntegrator(1e-9)
# Template for the solute solver
ssolver = cmf.ImplicitEuler(1e-9)
# Creating the SWI, the storage objects of the project are internally assigned to the correct solver
solver = cmf.SoluteWaterIntegrator(p.solutes, wsolver,ssolver,p)
solver = cmf.CVodeIntegrator(p, 1e-9)

c.saturated_depth = 1.5

Expand All @@ -60,19 +54,19 @@
tracer_state = []
wetness = []
# save groundwater recharge
recharge=[]
recharge = []
# save concentration of recharge
crecharge=[]
crecharge = []

# Run for one week with hourly time step
for t in solver.run(solver.t,solver.t + cmf.week,cmf.h):
for t in solver.run(solver.t, solver.t + cmf.week, cmf.h):
# Get concentration of all layers
tracer_state.append([[l.Solute(T).state for T in p.solutes] for l in c.layers])
tracer_state.append([[l[T].state for T in p.solutes] for l in c.layers])
conc.append([[l.conc(T) for T in p.solutes] for l in c.layers])
wetness.append([l.wetness for l in c.layers])
# Get water balance of groundwater
recharge.append(gw.waterbalance(t))
crecharge.append([gw.conc(t,T) for T in p.solutes])
crecharge.append([gw.conc(t, T) for T in p.solutes])
print(t)


Expand All @@ -81,17 +75,19 @@
import pylab as plt

pc = len(p.solutes) + 2
ax1 = plt.subplot(pc,1,1)
plt.plot(recharge,'k:')
plt.legend(['water'],loc=2)
ax1 = plt.subplot(pc, 1, 1)
plt.plot(recharge, 'k:')
plt.legend(['water'], loc=2)
plt.twinx()
plt.plot(crecharge)
plt.legend([str(s) for s in p.solutes],loc=1)
plt.subplot(pc,1,2,sharex=ax1)
plt.imshow(np.transpose(wetness),cmap=plt.cm.RdYlBu,aspect='auto')
plt.subplot(pc, 1, 2, sharex=ax1)
plt.imshow(np.transpose(wetness), cmap=plt.cm.RdYlBu, aspect='auto')
plt.ylabel('wetness')
for i,s in enumerate(p.solutes):
plt.subplot(pc,1,3+i,sharex=ax1)
plt.imshow(np.transpose(tracer_state)[i],cmap=plt.cm.viridis,aspect='auto',vmax=1.0, vmin=0.0)
plt.subplot(pc, 1, 3+i, sharex=ax1)
plt.imshow(np.transpose(tracer_state)[i],
cmap=plt.cm.viridis, aspect='auto',
vmax=1.0, vmin=0.0)
plt.ylabel(s)
plt.show()

0 comments on commit a0f9b57

Please sign in to comment.