Skip to content

Commit 82e7b13

Browse files
aleealee
committed
add irrigation flow matlab module
- on behalf of Andrew Bell <andrew.reid.bell@nyu.edu>
1 parent 90be8bd commit 82e7b13

File tree

11 files changed

+459
-1
lines changed

11 files changed

+459
-1
lines changed

‎ChannelLink.m

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
classdef ChannelLink < handle
2+
3+
properties
4+
id
5+
channelID
6+
length
7+
designFlow
8+
maintenance
9+
depreciationRate
10+
level
11+
inletNode
12+
outletNode
13+
inletWater
14+
inletTime
15+
end
16+
17+
events
18+
%none at the moment
19+
end
20+
21+
methods
22+
function A = ChannelLink(id, channelID, length, designFlow, level, inletNode, outletNode)
23+
A.id = id;
24+
A.channelID = channelID;
25+
A.length = length;
26+
A.designFlow = designFlow;
27+
A.level = level;
28+
A.inletNode = inletNode;
29+
A.outletNode = outletNode;
30+
end % Channel link
31+
32+
33+
end % methods
34+
end % classdef

‎ChannelNode.m

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
classdef ChannelNode < handle
2+
3+
properties
4+
id
5+
inletLinks
6+
outletLinks
7+
withdrawalFarms
8+
x
9+
y
10+
end
11+
12+
events
13+
%none at the moment
14+
end
15+
16+
methods
17+
function A = ChannelNode(id, x, y)
18+
%(just the constructor)
19+
A.id = id;
20+
A.x = x;
21+
A.y = y;
22+
end % Channel Node
23+
24+
25+
26+
27+
end % methods
28+
end % classdef

‎Farm.m

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
classdef Farm < handle
2+
3+
properties
4+
id
5+
withdrawal
6+
receipt
7+
wealth
8+
node
9+
channelMember
10+
x
11+
y
12+
end
13+
14+
events
15+
%none at the moment
16+
end
17+
18+
methods
19+
function A = Farm(id, withdrawal, node, x, y, channelMember)
20+
A.id = id;
21+
A.withdrawal = withdrawal;
22+
A.node = node;
23+
A.x = x;
24+
A.y = y;
25+
A.wealth = 0;
26+
A.receipt = withdrawal;
27+
A.channelMember = channelMember;
28+
29+
end % Farm
30+
31+
32+
end % methods
33+
end % classdef

‎README.md

Lines changed: 40 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,40 @@
1-
# matlab.irrigation-system
1+
PACKAGE NAME: Irrigation Flow
2+
3+
PLATFORM: MATLAB
4+
5+
AUTHOR: Andrew Bell
6+
7+
VERSION HISTORY: 1.0
8+
9+
### Description
10+
A set of objects and routines to describe an irrigation system as a set of nodes, links, and farms; and solve water flow
11+
and withdrawal through the system. Water *flow* is solved by specifying water demand at nodes and withdrawing it, with
12+
residual water allocated proportionally to all other outflow links based on their specified design flow. Thus, this
13+
simple model does not rely on physical solution (i.e., the St. Venant equations) to calculate flow rate or depth.
14+
15+
### Inputs
16+
1. An irrigation system, described in terms of channel links and nodes
17+
2. A volume of water appropriate to the time period
18+
19+
### Outputs
20+
Water receipt at each node for the time period
21+
22+
### PACKAGE CONTENTS
23+
24+
1. runModel.m: An outer script to set up the model space and call the main time loop
25+
2. setUpIrrigationSystem.m: A sample specification of an irrigation system (two watercourses off of a larger distributary) for testing purposes
26+
3. main.m: The main time loop within which solveWater.m is called
27+
4. visualizeSystem.m: A routine for visualizing the system state
28+
5. ChannelNode.m: An object class for the nodes in the system
29+
6. ChannelLink.m: An object class for the links between nodes in the system
30+
7. Farm.m: An object class for the farms withdrawing water from nodes in the system
31+
8. nodeDistance.m: A utility to calculate distance between points
32+
9. plotCircle.m: A utility called by visualizeSystem.m to plot circles
33+
34+
### INSTRUCTIONS
35+
36+
To test this package, copy all files to a folder and run runModel.m
37+
38+
### BENCHMARKS
39+
40+
Vary the level of input water supply, as well as the water demand of the farms, and observe that when supply is much greater than demand, all farms receive adequate water (appear as blue in the visualization) but as water supply decreases many downstream farms are water deficient (appear as black in the visualization).

‎main.m

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
function main(timeSteps, nodeList, farmList, linkList, channelList)
2+
%pointed to by runModel.m
3+
4+
5+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6+
%START TIME LOOP
7+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8+
9+
for indexT = 1:timeSteps.total
10+
11+
12+
inletConditions = 10 * rand();
13+
solveWater(inletConditions, indexT, nodeList);
14+
15+
visualizeSystem(timeSteps.cycle+1, farmList, nodeList, linkList, timeSteps);
16+
17+
end

‎nodeDistance.m

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
function distance = nodeDistance(node1, node2)
2+
3+
%calculate the euclidean distance between two nodes
4+
distance = ((node1.x - node2.x)^2 + (node1.y - node2.y)^2)^0.5;

‎plotCircle.m

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
function plotCircle(x,y,r,c)
2+
3+
%pretty straightforward. plot a circle at (x,y) of radius r, in color c
4+
5+
numPoints = 20;
6+
angleList = 0:2*pi/numPoints:2*pi;
7+
8+
xVec = x + r*cos(angleList);
9+
yVec = y + r*sin(angleList);
10+
11+
temp = patch(xVec,yVec,c);
12+
set(temp, 'EdgeAlpha',0.1);
13+
14+
end %plotCircle

‎runModel.m

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
%Outer script to set up and run irrigation system model
2+
3+
clear all;
4+
close all;
5+
6+
7+
%Set seed for random number generator
8+
randSeed = 14897;
9+
rand('twister',randSeed);
10+
randn('state', randSeed);
11+
12+
%Parameters
13+
timeSteps.total = 1560; %length of simulation
14+
timeSteps.cycle = 52; %e.g., year
15+
16+
%initialize system using setUpIrrigationSystem function (right now just
17+
%contains hard coding of a trivial test case
18+
[farmList, nodeList, linkList, channelList] = setUpIrrigationSystem();
19+
20+
21+
%make a first plot of the system to work from
22+
figure;
23+
visualizeSystem(timeSteps.cycle+1, farmList, nodeList, linkList, timeSteps);
24+
25+
%run model
26+
main(timeSteps, nodeList, farmList, linkList, channelList);

‎setUpIrrigationSystem.m

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
function [farmList, nodeList, linkList, channelList] = setUpIrrigationSystem()
2+
3+
%this script initializes a simple irrigation system for testing purposes,
4+
%with two watercourses running off a single distributary, each with 10
5+
%equally spaced farms along them
6+
7+
%the farm objects created here are specific to a particular application
8+
%of an irrigation system, in which farmers are making choices of what
9+
%to plant based on reliability of water receipt
10+
11+
%For now, manually build a small network of two watercourses off a single
12+
%distributary
13+
14+
nodeID = 0;
15+
farmID = 0;
16+
linkID = 0;
17+
18+
channelHierarchy(1).id = 1; channelHierarchy(1).parent = 0; channelHierarchy(1).order = 2;
19+
channelHierarchy(2).id = 2; channelHierarchy(2).parent = 1; channelHierarchy(2).order = 3;
20+
channelHierarchy(3).id = 3; channelHierarchy(3).parent = 1; channelHierarchy(3).order = 3;
21+
22+
%make the source and distributary nodes
23+
nodeID = nodeID + 1;
24+
nodeList(1) = ChannelNode(nodeID, 0, 0);
25+
26+
nodeID = nodeID + 1;
27+
nodeList(end+1) = ChannelNode(nodeID, 2, 0);
28+
29+
linkID = linkID + 1;
30+
linkList(1) = ChannelLink(linkID, channelHierarchy(1).id, nodeDistance(nodeList(1),nodeList(2)), .001, channelHierarchy(1).order, nodeList(1), nodeList(2));
31+
nodeList(1).outletLinks = linkList(1);
32+
nodeList(2).inletLinks = linkList(1);
33+
34+
channelList{1}(1) = linkList(1);
35+
36+
nodeID = nodeID + 1;
37+
nodeList(end+1) = ChannelNode(nodeID, 4, 0);
38+
39+
linkID = linkID + 1;
40+
linkList(end+1) = ChannelLink(linkID, channelHierarchy(1).id, nodeDistance(nodeList(2),nodeList(3)), .001, channelHierarchy(1).order, nodeList(2), nodeList(3));
41+
nodeList(2).outletLinks = linkList(2);
42+
nodeList(3).inletLinks = linkList(2);
43+
44+
channelList{1}(2) = linkList(2);
45+
46+
%now make the watercourse and farm nodes
47+
for indexA = 1:10
48+
49+
nodeID = nodeID + 1;
50+
farmID = farmID + 1;
51+
linkID = linkID + 1;
52+
xFarm = 1;
53+
xNode = 2;
54+
y = indexA;
55+
withdrawal = rand();
56+
57+
currentNode = ChannelNode(nodeID, xNode, y);
58+
59+
currentFarm = Farm(farmID, withdrawal, currentNode, xFarm, y, channelHierarchy([1 2]));
60+
61+
if(indexA == 1)
62+
currentLink = ChannelLink(linkID, channelHierarchy(2).id, nodeDistance(nodeList(2),currentNode), .001, channelHierarchy(2).order, nodeList(2), currentNode);
63+
nodeList(end+1) = currentNode;
64+
farmList(indexA) = currentFarm;
65+
linkList(end+1) = currentLink;
66+
nodeList(end).inletLinks = linkList(end);
67+
nodeList(end).withdrawalFarms = currentFarm;
68+
nodeList(2).outletLinks(end+1) = linkList(end);
69+
channelList{2}(1) = currentLink;
70+
71+
else
72+
currentLink = ChannelLink(linkID, channelHierarchy(2).id, nodeDistance(nodeList(end),currentNode), .001, channelHierarchy(2).order, nodeList(end), currentNode);
73+
nodeList(end+1) = currentNode;
74+
farmList(indexA) = currentFarm;
75+
linkList(end+1) = currentLink;
76+
nodeList(end).inletLinks = linkList(end);
77+
nodeList(end).withdrawalFarms = currentFarm;
78+
nodeList(end-1).outletLinks = linkList(end);
79+
channelList{2}(end+1) = currentLink;
80+
end
81+
82+
end
83+
84+
%drain for this watercourse
85+
nodeID = nodeID + 1;
86+
linkID = linkID + 1;
87+
xNode = 2;
88+
y = indexA+1;
89+
90+
currentNode = ChannelNode(nodeID, xNode, y);
91+
currentLink = ChannelLink(linkID, channelHierarchy(2).id, nodeDistance(nodeList(end),currentNode), .001, channelHierarchy(2).order, nodeList(end), currentNode);
92+
nodeList(end+1) = currentNode;
93+
linkList(end+1) = currentLink;
94+
95+
for indexA = 1:10
96+
97+
nodeID = nodeID + 1;
98+
farmID = farmID + 1;
99+
linkID = linkID + 1;
100+
xFarm = 3;
101+
xNode = 4;
102+
y = indexA;
103+
withdrawal = rand();
104+
105+
106+
currentNode = ChannelNode(nodeID, xNode, y);
107+
108+
currentFarm = Farm(farmID, withdrawal, currentNode, xFarm, y, channelHierarchy([1 3]));
109+
110+
if(indexA == 1)
111+
currentLink = ChannelLink(linkID, channelHierarchy(3).id, nodeDistance(nodeList(3),currentNode), .001, channelHierarchy(3).order, nodeList(3), currentNode);
112+
nodeList(end+1) = currentNode;
113+
farmList(end+1) = currentFarm;
114+
linkList(end+1) = currentLink;
115+
nodeList(end).inletLinks = linkList(end);
116+
nodeList(end).withdrawalFarms = currentFarm;
117+
nodeList(3).outletLinks = linkList(end);
118+
channelList{3}(1) = currentLink;
119+
else
120+
currentLink = ChannelLink(linkID, channelHierarchy(3).id, nodeDistance(nodeList(end),currentNode), .001, channelHierarchy(3).order, nodeList(end), currentNode);
121+
nodeList(end+1) = currentNode;
122+
farmList(end+1) = currentFarm;
123+
linkList(end+1) = currentLink;
124+
nodeList(end).inletLinks = linkList(end);
125+
nodeList(end).withdrawalFarms = currentFarm;
126+
nodeList(end-1).outletLinks = linkList(end);
127+
channelList{2}(end+1) = currentLink;
128+
end
129+
130+
end
131+
132+
%drain for this watercourse
133+
nodeID = nodeID + 1;
134+
linkID = linkID + 1;
135+
xNode = 4;
136+
y = indexA+1;
137+
138+
currentNode = ChannelNode(nodeID, xNode, y);
139+
currentLink = ChannelLink(linkID, channelHierarchy(3).id, nodeDistance(nodeList(end),currentNode), .001, channelHierarchy(3).order, nodeList(end), currentNode);
140+
nodeList(end+1) = currentNode;
141+
linkList(end+1) = currentLink;
142+

0 commit comments

Comments
 (0)