Merge from upstream and repaired cmf_wrap.cpp

.gitignore

@@ -15,8 +15,7 @@ __pycache__/
 *.VC.opendb
 
 # Distribution / packaging
-.Python
-build/
+.Pythonbuild/
 develop-eggs/
 dist/
 downloads/
@@ -100,3 +99,5 @@ cmake-build-*/
 *.sbn
 *.sbx
 nohup.out
+/demo/results/
+/build/

cmf/cmf_core_src/upslope/surfacewater.cpp

@@ -67,6 +67,30 @@ void DiffusiveSurfaceRunoff::connect_cells( cmf::upslope::Cell& cell1,cmf::upslo
 		new DiffusiveSurfaceRunoff(sw_upper,lower.get_surfacewater(),w);
 	}
 
+}
+/// Returns the square root of the gradient with singularity prevention
+///
+/// \param grad Gradient in m/m
+/// \return sign(grad) * sqrt(grad) but with a flatted singularity at grad = 0
+real sqrt_slope(real grad) {
+
+    const real s_sqrt = sign(grad) * sqrt(std::abs(grad));
+    // linear slope width is a value in which slope range the slope should be altered
+    // to prevent a singularity in dq/ds
+    if (cmf::options.diffusive_slope_singularity_protection > 0.0) {
+        // Only a shortcut for faster writing
+        const real &s0 = cmf::options.diffusive_slope_singularity_protection;
+        // Weight of linear part
+        real w_lin = exp(-square((grad / s0)));
+        // linear part using the slope at s0/4
+        real s_lin = grad / (2. * sqrt(s0 / 4));
+        // Weighted sum of sqrt(s) and a*s
+        return w_lin * s_lin + (1 - w_lin) * s_sqrt;
+    } else {
+        return s_sqrt;
+
+    }
+
 }
 
 real DiffusiveSurfaceRunoff::calc_q( cmf::math::Time t )
@@ -101,21 +125,9 @@ real DiffusiveSurfaceRunoff::calc_q( cmf::math::Time t )
 	// Get slope
 	real grad = dPsi/m_distance;
 
-	// Get signed square root for 
-	real s_sqrt = sign(grad) * sqrt(std::abs(grad));
+	// Get signed square root for gradient
+	real s_sqrt = sqrt_slope(grad);
 
-	// linear slope width is a value in which slope range the slope should be altered
-	// to prevent a singularity in dq/ds
-	if (cmf::options.diffusive_slope_singularity_protection > 0.0) {
-        // Only a shortcut for faster writing
-		const real & s0 = cmf::options.diffusive_slope_singularity_protection;
-        // Weight of linear part
-		real w_lin = exp(-square((grad/s0)));
-        // linear part using the slope at s0/4
-        real s_lin =grad/(2.*sqrt(s0/4));
-        // Weighted sum of sqrt(s) and a*s
-		s_sqrt = w_lin * s_lin + (1-w_lin) * s_sqrt;
-	}
 
 	return prevent_negative_volume(m_flowwidth * pow(d,5/3.) * s_sqrt/left->get_nManning() * 86400.);
 }
@@ -126,3 +138,4 @@ void DiffusiveSurfaceRunoff::NewNodes()
 	wright = SurfaceWater::cast(right_node());
 	owright=OpenWaterStorage::cast(right_node());
 }
+

cmf/cmf_core_src/upslope/surfacewater.h

@@ -133,9 +133,10 @@ namespace cmf {
 			private:
 				static void connect_cells(cmf::upslope::Cell& c1,cmf::upslope::Cell& c2,ptrdiff_t dummy);
 				std::weak_ptr<cmf::upslope::SurfaceWater> wleft;
+                real m_distance, m_flowwidth;
+            protected:
 				virtual real calc_q(cmf::math::Time t);
 				void NewNodes();
-				real m_distance, m_flowwidth;
 			public:
 				/// Creates a KinematicSurfaceRunoff between a SurfaceWater (left) with another (right) node.
 				///
@@ -187,10 +188,11 @@ namespace cmf {
 				std::weak_ptr<cmf::upslope::SurfaceWater> wleft;
 				std::weak_ptr<cmf::river::OpenWaterStorage> owright;
 				std::weak_ptr<cmf::upslope::SurfaceWater> wright;
-				real calc_q(cmf::math::Time t) override;
-				void NewNodes() override;
 				real m_distance, m_flowwidth;
-			public:
+            protected:
+                real calc_q(cmf::math::Time t) override;
+                void NewNodes() override;
+            public:
 				DiffusiveSurfaceRunoff(cmf::upslope::SurfaceWater::ptr left,cmf::water::flux_node::ptr right,real flowwidth,real distance=-1)
 					: flux_connection(left,right,"Diffusive surface runoff"), m_distance(distance), m_flowwidth(flowwidth) {
 

cmf/geometry/irregular_grid.py

@@ -0,0 +1,202 @@
+import numpy as np
+import scipy.spatial as scsp
+from shapely.geometry import Polygon
+
+from .. import project as Project
+from . import create_cell, mesh_project
+
+def delaunay_cells(project: Project, x, y, z):
+    """
+    Creates a mesh using delaunay triangles between the given coordinates as cells.
+    cell.x, cell.y and cell.z represent the centroid of the cell. This results in loss of information.
+    In most cases one would use the voronoi_cells function
+
+    x, y and z must have the same length
+
+    :param project:
+    :param x: sequence of x coordinates
+    :param y: sequence of y coordinates
+    :param z: sequence of heights
+    :return: list of cells
+    """
+    points = np.array([x, y]).T
+    tri = scsp.Delaunay(points)
+
+    def get_edge_length(tri: scsp.Delaunay, tri1: int, tri2: int):
+        try:
+            n1, n2 = set(tri.simplices[tri1]) & set(tri.simplices[tri2])
+        except ValueError:
+            raise ValueError(
+                f'Triangle {tri1}:{tri.simplices[tri1]} and {tri2}:{tri.simplices[tri2]} do not share an edge')
+        c1, c2 = tri.points[[n1, n2]]
+        return np.sqrt(((c1 - c2) ** 2).sum())
+
+    def get_polygons(tri: scsp.Delaunay):
+
+        return [
+            Polygon(tri.points[s])
+            for s in tri.simplices
+        ]
+
+    def get_heights(tri: scsp.Delaunay, z):
+        return z[tri.simplices].mean(1)
+
+    polys = get_polygons(tri)
+    z_cell = get_heights(tri, z)
+
+    cells = [
+        create_cell(project, polygon, height)
+        for polygon, height in zip(polys, z_cell)
+    ]
+
+    for i, c in enumerate(cells):
+        for n in tri.neighbors[i]:
+            if n >= 0:
+                l = get_edge_length(tri, i, n)
+                c.topology.AddNeighbor(cells[n], l)
+
+    return cells
+
+
+# %%
+def __voronoi_finite_polygons_2d(vor, radius=None):
+    """
+    Reconstruct infinite voronoi regions in a 2D diagram to finite
+    regions.
+
+    Parameters
+    ----------
+    vor : Voronoi
+        Input diagram
+    radius : float, optional
+        Distance to 'points at infinity'.
+
+    Returns
+    -------
+    regions : list of tuples
+        Indices of vertices in each revised Voronoi regions.
+    vertices : list of tuples
+        Coordinates for revised Voronoi vertices. Same as coordinates
+        of input vertices, with 'points at infinity' appended to the
+        end.
+
+    """
+
+    if vor.points.shape[1] != 2:
+        raise ValueError("Requires 2D input")
+
+    new_regions = []
+    new_vertices = vor.vertices.tolist()
+
+    center = vor.points.mean(axis=0)
+    if radius is None:
+        radius = vor.points.ptp(axis=0).max()
+
+    # Construct a map containing all ridges for a given point
+    all_ridges = {}
+    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
+        all_ridges.setdefault(p1, []).append((p2, v1, v2))
+        all_ridges.setdefault(p2, []).append((p1, v1, v2))
+
+    # Reconstruct infinite regions
+    for p1, region in enumerate(vor.point_region):
+        vertices = vor.regions[region]
+
+        if all(v >= 0 for v in vertices):
+            # finite region
+            new_regions.append(vertices)
+            continue
+        elif p1 not in all_ridges:
+            continue
+        # reconstruct a non-finite region
+        ridges = all_ridges[p1]
+        new_region = [v for v in vertices if v >= 0]
+
+        for p2, v1, v2 in ridges:
+            if v2 < 0:
+                v1, v2 = v2, v1
+            if v1 >= 0:
+                # finite ridge: already in the region
+                continue
+
+            # Compute the missing endpoint of an infinite ridge
+
+            t = vor.points[p2] - vor.points[p1]  # tangent
+            t /= np.linalg.norm(t)
+            n = np.array([-t[1], t[0]])  # normal
+
+            midpoint = vor.points[[p1, p2]].mean(axis=0)
+            direction = np.sign(np.dot(midpoint - center, n)) * n
+            far_point = vor.vertices[v2] + direction * radius
+
+            new_region.append(len(new_vertices))
+            new_vertices.append(far_point.tolist())
+
+        # sort region counterclockwise
+        vs = np.asarray([new_vertices[v] for v in new_region])
+        c = vs.mean(axis=0)
+        angles = np.arctan2(vs[:, 1] - c[1], vs[:, 0] - c[0])
+        new_region = np.array(new_region)[np.argsort(angles)]
+
+        # finish
+        new_regions.append(new_region.tolist())
+
+    return new_regions, np.asarray(new_vertices)
+
+def voronoi_polygons(x, y, buffer=0.0, mask=None):
+    """
+    Creates closed voronoi polygons for x and y coordinates. With voronoi_cells, cmf cells are easily created
+    from x, y and z coordinates. Use this function directly only, if you need to process the voronoi polygons further,
+    eg. with unions of landuse or soil maps.
+
+    If this throws an unclear error, make sure you do not have any duplicate coordinates
+
+    :param x: a sequence of x values
+    :param y: a sequence of y values
+    :param buffer: a float indicating the buffer around the hull. Expressed as a fraction of the sqrt of the hull's area
+    :param mask: a polygon to use as a mask. If given, the convex hull is ignored
+    :return: A list of polygons
+    """
+    points = np.array([x, y]).T
+    vor = scsp.Voronoi(points)
+    if not mask:
+        hull = Polygon(points[scsp.ConvexHull(points).vertices])
+        mask = hull.buffer(buffer * hull.area ** 0.5)
+    regions, vertices = __voronoi_finite_polygons_2d(vor)
+    return [Polygon(vertices[r]).intersection(mask) for r in regions]
+
+
+def voronoi_cells(project: Project, x, y, z, buffer=0.0, mask=None):
+    """
+    Creates for each (x, y, z) coordinate a voronoi (nearest neighbor) cell.
+    The cells are cut of at a buffer around the convex hull of the point or using a mask.
+
+    For each coordinate a cell is created.
+
+    If this throws an unclear error, make sure you do not have any duplicate coordinates
+
+    :param project: The cmf project
+    :param x: a sequence of x values
+    :param y: a sequence of x values
+    :param z: a sequence of x values
+    :param buffer: a float indicating the buffer around the hull. Expressed as a fraction of the sqrt of the hull's area
+    :param mask: a polygon to use as a mask. If given, the convex hull is ignored
+    :return: list of cells
+    """
+    def make_cell(project, x, y, z, shape, Id=None, with_surfacewater=False):
+        cell = project.NewCell(x, y, z, shape.area, with_surfacewater=with_surfacewater)
+        cell.geometry = shape
+        if Id is not None:
+            cell.Id = Id
+        return cell
+
+    vpolys = voronoi_polygons(x, y, z, buffer, mask)
+    cells = [
+        make_cell(project, x[i], y[i], z[i], vpolys[i], Id=i, with_surfacewater=True)
+        for i in range(len(x))
+    ]
+
+    mesh_project(project)
+
+    return cells
+

demo/jacobian.py

@@ -0,0 +1,38 @@
+import cmf
+from matplotlib import pyplot as plt
+
+p = cmf.project('X')
+X, = p.solutes
+
+V1 = p.NewStorage('V1', 0, 0, 0)
+V2 = p.NewStorage('V2', 1, 0, 0)
+V3 = p.NewStorage('V3', 0, 0, 1)
+V4 = p.NewStorage('V4', 1, 0, 1)
+
+NB = p.NewNeumannBoundary('NB', V1)
+NB.flux = 1
+NB.concentration[X] = 1
+
+DB = p.NewOutlet('DB', -1, 0, 0)
+
+q_nb_1 = NB.connection_to(V1)
+q_1_2 = cmf.LinearStorageConnection(V1, V2, 0.1)
+q_1_3 = cmf.LinearStorageConnection(V1, V3, 0.5)
+q_2_4 = cmf.LinearStorageConnection(V2, V4, 1)
+q_4_3 = cmf.LinearStorageConnection(V4, V3, 0.5)
+q_3_out = cmf.LinearStorageConnection(V3, DB, 1)
+
+V2[X].decay = 0.5
+V4[X].decay = 0.5
+V3[X].decay = 0.1
+
+for S in [V1, V2, V3, V4]:
+    S.volume = 0.0
+
+solver = cmf.CVodeDiag(p)
+solver.t = cmf.Time()
+solver.integrate_until(cmf.year)
+print(solver.get_info())
+t = solver.t
+plt.imshow(solver.get_jacobian(), vmin=-1, vmax=1, cmap=plt.cm.coolwarm)
+plt.colorbar()

demo/pump_trigger.py

@@ -10,9 +10,9 @@ class SchmittTriggerSourceConnection:
     It is similar to the SchmittTriggerTargetConnection,
     but that is triggered by the state of the target.
 
-    If the volume of the water source is greater then the upper regulation level,
+    If th source is greater then the upper regulation level,
     the flux is switched on, if the level is smaller then the lower
-    regulation level, it is switched off. Between those levels, the
+    regulation level, it is swie volume of the watertched off. Between those levels, the
     flux stays in the same state as before.
 
     A trigger where the flow is turned on if the source has more than

demo/cmf2d.py → demo/subsurface/cmf2d.py

@@ -9,7 +9,7 @@
 import numpy as np
 from argparse import ArgumentParser
 from datetime import datetime, timedelta
-
+from matplotlib import cm
 try:
     import pylab
     import cmf.draw
@@ -176,7 +176,7 @@ def run(p, outlet, until, dt=cmf.day):
     outflow = cmf.timeseries(solver.t, dt)
     for t in solver.run(solver.t, solver.t + until, dt):
         outflow.add(outlet(t))
-        print("%20s - %6.1f l/day" % (t, outlet(t) * 1e3))
+        print("%20s - %6.1f l/day - %6.2f Wetness" % (t, outlet(t) * 1e3, p.cells[0].layers[0].Wetness))
     return outflow, solver.info
 
 
@@ -186,12 +186,12 @@ def animate(p, outlet, until, dt=cmf.day):
 
     solver = cmf.CVodeKrylov(p, 1e-9)
     solver.t = cmf.Time(1, 1, 1980)
-    hp = HillPlot(p, solver.t)
+    hp = HillPlot(p, solver.t, cmap=cm.viridis_r)
     hp.scale = 1000
 
     def integration():
         for t in solver.run(solver.t, solver.t + until, dt):
-            print(t)
+            print("%20s - %6.1f l/day - %6.2f Wetness" % (t, outlet(t) * 1e3, p.cells[0].layers[0].wetness))
             yield t
 
     anim = hp.get_animator(integration)