diff --git a/AquaME.lua b/AquaME.lua index 4324ab3..9870b5b 100644 --- a/AquaME.lua +++ b/AquaME.lua @@ -22,6 +22,8 @@ AquaME = Model { -- Default and constraints for actual lenght of the area modeled, in meters: xAreaSize = Choice {min = 0.001}, -- meters yAreaSize = Choice {min = 0.001}, -- meters + -- Neighborhood strategy. + strategy = Choice {"moore", "vonneumann"}, -- Default and constraints for initial relaxation parameter for all cells, representing the flow resistance. alpha = Choice {min = 0, max = 1, default = 0.988}, -- Default and constraints for alpha0, a constant for calculating alpha for rivers. @@ -34,45 +36,58 @@ AquaME = Model { infiltrationBias = Choice {min = 0}, -- meters/hours -- Load scenario terrain data. - loadTerrain = function (_,_) end, + loadTerrain = function () end, -- Load scenario precipitation data. - loadPrecipitation = function (_,_,_) return 0 end, + loadPrecipitation = function () return 0 end, - -- Return optional parms specific for instance. - otherParms = function (_) return end, + -- Return optional parms specific for scenario instance. + loadScenarioParms = function () end, -- Initialize model: init = function(model) - -- Maximum duration, time steps. - model.finalTime = math.ceil(model.duration / model.stepHours) + 1 + -- Loads into model any addtional parm needed for scenario instance: + model:loadScenarioParms() - model.xdim = math.ceil(model.xAreaSize / model.xCellSize) -- Number or cells in the x dimension. - model.ydim = math.ceil(model.yAreaSize / model.yCellSize) -- Number or cells in the y dimension. + -- Calculates the maximum duration in time steps: + model.finalTime = math.ceil(model.duration / model.stepHours) + 1 - -- Sequence of the 9 partitions for runoff calculation (Rinaldi at al, Figure 5). - model.sequence = {{2,2},{1,1},{2,1}, - {1,2},{2,3},{1,3}, - {3,2},{3,1},{3,3}} + -- Calculates the number or cells in the x and y dimensions: + model.xdim = math.ceil(model.xAreaSize / model.xCellSize) + model.ydim = math.ceil(model.yAreaSize / model.yCellSize) + + -- Define the sequence of partitions for runoff calculation according to the strategy and to the center cell x and y coordinates. + if model.strategy == "vonneumann" + then model.sequence = { {1,2,3,4,5}, + {4,5,1,2,3}, + {2,3,4,5,1}, + {5,1,2,3,4}, + {3,4,5,1,2} + } + else model.sequence = { {2,4,6}, -- (Rinaldi at al, Figure 5) + {3,1,5}, + {8,7,9} + } + end -- No outflow at the beginning of the simulation. - model.outflow = 0 + model.outflow = {} model.cell = Cell { init = function(cell) - -- Set defaults for cell using nil for lower memory consumption: - -- cell.open -- Cell has no open border. - -- cell.isRiver -- Cell is not on a river. - -- cell.infiltration -- No initial infiltration. - -- cell.infiltrationBias -- Typical infiltration bias for terrain from model. - -- cell.alpha = model.alpha -- Typical alpha for terrain from model. + -- Many defaults use nil for lower memory consumption: + -- cell.open -- Cell has no open border. + -- cell.isRiver -- Cell is not on a river. + -- cell.infiltration -- No initial infiltration. + -- cell.infiltrationBias -- Typical infiltration bias for terrain comes from model. + -- cell.alpha -- Typical alpha for terrain comes from model. -- Set defaults for cell using atibuttes because they are usually updated in each time step: cell.height = 0 -- All plain. cell.water = 0 -- No initial water. - -- Fix cells' values from default ones according to scenario: + -- Fix cells' height and alpha from default ones according to scenario: model:loadTerrain(cell) end @@ -80,45 +95,29 @@ AquaME = Model { --** print(sessionInfo().time, ": Creates cellular space.") model.cs = CellularSpace {xdim = model.xdim, ydim = model.ydim, instance = model.cell} - --** print(sessionInfo().time, ": Creates cellular space neighborhood with Moore strategy (3x3 grid).") - model.cs:createNeighborhood {} + --** print(sessionInfo().time, ": Creates cellular space neighborhood with model strategy.") + model.cs:createNeighborhood {strategy = model.strategy} - --** print(sessionInfo().time, ": Creates 9 trajectories representing the 9 partitions with their grids.") - model.partition = {}; + --** print(sessionInfo().time, ": Creates a grid around each cell.") + model.grid = {} + forEachCell(model.cs, function(cell) - for t = 1, #model.sequence do + -- Calculates the index (i.e., the order of evaluation) of the partition that the grid around the cell belongs to: + local partition = model.sequence[1 + cell.x % #model.sequence][1 + cell.y % #model.sequence] - local x = model.sequence[t][1] -- The x reference for the center cells of this partition. - local y = model.sequence[t][2] -- The y reference for the center cells of this partition. + -- Creates the partition, if it hasn't been created yet: + model.grid[partition] = model.grid[partition] or {} - model.partition[t] = Trajectory { - -- Each partition contains cells defined by (x,y) in such a way that they are separated by 2 cells above and below, - -- so there is no neighbor at the same partition. - target = model.cs, select = function(cell) return (1 + cell.x % 3) == x and (1 + cell.y % 3) == y end - } - end + -- Initiates the grid around the cell with the cell itself, on the partition the grid belongs to: + model.grid[partition][cell] = {cell} - --** print(sessionInfo().time, ": Creates grids around each cell.") - model.grid = {} - forEachCell(model.cs, function(cell) - -- Creates a grid around the cell: - model.grid[cell] = {cell} - -- Add cell's neighbors to the grid: - forEachNeighbor(cell, function(neighbor) table.insert(model.grid[cell], neighbor) end) - -- Order the cells index in the grid according to the terrain height, in ascending order (Rinaldi et al. 2007, equation 1). - table.sort(model.grid[cell], function(lower, higher) return lower.height < higher.height end) - end) + -- Adds cell's neighbors to this grid: + forEachNeighbor(cell, function(neighbor) table.insert(model.grid[partition][cell], neighbor) end) - --** print(sessionInfo().time, ": Define chart.") + -- Orders the cells in the grid according to the terrain height, in ascending order (Rinaldi et al. 2007, equation 1). + table.sort(model.grid[partition][cell], function(lower, higher) return lower.height < higher.height end) - -- [[** To turning chart off. - model.chart = Chart { - title = model.name, - target = model, --** Alternative: model.cs:get(2, model.ydim) - select = "outflow", --** Alternatives: water, infiltration, alpha, height. - color = "blue" - } - --]] + end) --** print(sessionInfo().time, ": Define timed events.") model.timer = Timer { @@ -129,12 +128,12 @@ AquaME = Model { forEachCell(model.cs, function(cell) - -- Adds precipitation to cell water, if there is any defined in model. + -- Adds precipitation to cell water. cell.water = cell.water + model:loadPrecipitation(cell, event:getTime()) --** print(sessionInfo().time, ": Calculates cell infiltration.") - -- Using nil for defaults for lower memory consumption. + -- Replace nil by default values: local infiltration = cell.infiltration or 0 local infiltrationBias = cell.infiltrationBias or model.infiltrationBias local beta = cell.beta or model.beta @@ -153,12 +152,12 @@ AquaME = Model { local alpha = 1 - model.alpha0 * cell.water ^ model.nExp -- Alpha can't be lower than zero. if alpha < 0 then alpha = 0 end - --[[** Alpha for rivers shouldn't be higher than alpha for terrain, but authors oriented this way: + --[[** Alpha for rivers shouldn't be higher than alpha for terrain, but Rinaldi et al. (2007) oriented this way: if alpha > model.alpha then alpha = model.alpha end --]] - -- Alpha is recalculated and saved for rivers only: for other cells, alpha is not saved and always assumed to be model.alpha. cell.alpha = alpha end + end) end}, @@ -166,72 +165,65 @@ AquaME = Model { --** print(sessionInfo().time, ": Calculates runoff.") - -- For all partitions ordered by model.sequence (Rinaldi et al. 2007, figure 5). - for t = 1, #model.sequence do + -- For each partition in their order (Rinaldi et al. 2007, figure 5): + for partition = 1, #model.grid do + + -- For each grid of the partition, obtains the center cell and its grid: + for cell, grid in pairs(model.grid[partition]) do + --[[ The steps in this loop were carefully crafted in order to enhance performance, + as this is the most repeated loop in the algorithm. + --]] - forEachCell(model.partition[t], function(cell) + local N = #grid -- The size of the grid. Depends on the strategy and if the cell is at a corner or border. + local h = {}; h[N] = grid[N].height -- Will store the height of each cell. + local s = 0 -- Will store the sum of h[i] from h[1] to h[N-1] - -- Local variable for improved readability and performance: - local grid = model.grid[cell] -- Obtain the gridaroung the cell previously ordered by height. - local N = #grid -- The size of the grid ordered by height. Mostly 3x3 = 9, but could be 4 at corners and 6 at borders. - local wOld = {} -- Current water level of the ith cell of the grid ordered by height. - local h = {} -- Cell height of the ith cell of the ordered grid ordered by height. - local W = 0 -- Current water volume contained in the grid divided by the unit-cell area. + local W = grid[N].water -- Will store the current water volume contained in the grid divided by the unit-cell area (Rinaldi et al. 2007, equation 3). + local alpha = cell.alpha or model.alpha -- Resistance to runoff. Using nil as default for lower memory consumption. - --For each ith cell in the ordered grid: - for i = 1, N do - wOld[i] = grid[i].water; h[i] = grid[i].height - -- W is the current water volume contained in the grid divided by the unit-cell area (Rinaldi et al. 2007, equation 3). - W = W + wOld[i] + -- Store the next water level. Adjusted a few lines below for lower cells that receive water from higher cells. + grid[N].water = alpha * grid[N].water + + for i = 1, N-1 do + h[i] = grid[i].height + s = s + h[i] + W = W + grid[i].water + grid[i].water = alpha * grid[i].water end - -- If there is no water in the grid, there is nothing to calculate. + --[[** If there is no water in the grid, there is nothing else to calculate, + but testing this would only help performance if there was a large number of cells with no water each time step. if W == 0 then return end + --]] -- Obtain k, the number of cells that remains wet after the water drains down (Rinaldi et al. 2007, equation 4): - local k = N - repeat - local s = 0; for i = 1, k do s = s + h[k] - h[i] end - k = k - 1 - until W >= s or k == 0 - k = k + 1 - - -- Calculate the equilibrium surface water height, H, as the height the water would reach (Rinaldi et al. 2007, equation 2): - local H = W/k; for i = 1, k do H = H + h[i]/k end - - -- Calculate the new water level of ith cell of the grid as a linear combination of wOld[i] and the drained water level (Rinaldi et al. 2007, equation 6): - -- The drained water level wDrain[i] is zero for cells below k and H - h[i] for cells above (Rinaldi et al. 2007, equation 5). - -- The value of alpha corresponds to the center cell of the grid. - local alpha = cell.alpha or model.alpha -- Using nil as default for lower memory consumption. - for i = 1, k do grid[i].water = alpha * wOld[i] + (1 - alpha) * (H - h[i]) end - for i = k + 1, N do grid[i].water = alpha * wOld[i] end - - --[[** Checks conservation of water at the grid: - local Wnew = 0; for i = 1, k do Wnew = Wnew + H - h[i] end - if W > 0 and (W - Wnew)/W > 10^(-8) then - print("W changed at: (",cell.x,",",cell.y,") in:", 100*(Wnew-W)/W, "%.") - end - --**]] + -- Obtain also "s", the sum of h[i] until h[k-1], a shortcut for obtaining H (see below). + local k = N; while W + s < (k-1)*h[k] and k >= 2 do k = k - 1; s = s - h[k] end + + -- Calculate H, the equilibrium surface height the water would reach if enough time was provided (Rinaldi et al. 2007, equation 2): + local H = (W + s + h[k])/k + + -- The added drained water level wDrain[i] is H - h[i] for cells below k (Rinaldi et al. 2007, equation 5). + -- The new water level of ith cell of the grid as a linear combination of current water and drained water levels (Rinaldi et al. 2007, equation 6): + for i = 1, k do grid[i].water = grid[i].water + (1 - alpha) * (H - h[i]) end - end) + end end end}, - Event {priority = 3, action = function(_) + Event {priority = 3, action = function(event) --** print(sessionInfo().time, ": Calculates outflow.") local outflow = 0 + -- If cell is at an open border, let water ouflow from cell . forEachCell(model.cs, function(cell) - - -- If cell at an open border,let water ouflow from cell . if cell.open then outflow = outflow + cell.water; cell.water = 0 end - end) - -- Calculates outflow in m3/sec: - model.outflow = outflow + -- Calculates and save outflow in m3/sec for the time step: + model.outflow[event:getTime()] = outflow * model.xCellSize -- meters per cell, x dimension * model.yCellSize -- meters per cell, y dimension / model.stepHours -- hours per time step @@ -239,14 +231,16 @@ AquaME = Model { end}, - -- Updates chart. - Event {priority = 4, action = model.chart}, - - --[[** - Event {priority = 5, period = 1000, action = function(event) + -- [[** + Event {priority = 5, period = 5000, action = function(event) print(sessionInfo().time, ": End of time step ", event:getTime(), ".") - end + end}, --]] + + Event {priority = 6, period = model.finalTime, action = function(_) + --** print(sessionInfo().time, ": Shows chart.") + Chart{target = DataFrame{outflow = model.outflow}, select = "outflow", color = "blue" } + end}, } end } @@ -254,6 +248,7 @@ AquaME = Model { -- Creates a AquaME scenario redefining model parameters and functions. simplePlain = AquaME { + strategy = "vonneumann", name = "Simple plain", duration = 40, -- Duration of the simulation in hours. stepHours = 1/1000, -- Hours per time step. @@ -264,102 +259,73 @@ simplePlain = AquaME { alpha = 0.988, -- Flow resistance for terrain. alpha0 = 0.012, -- Constant for calculating alpha for rivers. - -- Add parameters for functions specific for scenario. - otherParms = function (_) return { - plainSlope = 0.0085, -- Plain slope at the ydim only, in meters per meter. - channelSlope = 0.0025, -- Channel slope at the xdim only, in meters per meter. - channelWidth = 5, -- Channel width, in meters. Adjusted to grids (ceil). - -- If absent or zero, no channel and opens the lower border of plain for outflow. - channelDepth = 2.5, -- Additional depth at the highest border of the channel, in meters. - rainBegin = 0, -- Instant when when rain begins after the beginning of the simulatio, in hours. - rainDuration = 10, -- Rain duration, in hours. - rainRate = 0.015, -- Rate of rain, in mm per hour. - rainAll = false -- True if it rains also on the channel. - } + -- Adds specific parameters to the scenario. + loadScenarioParms = function (model) + + model.plainSlope = 0.0085 -- Plain slope at the ydim only, in meters per meter. + model.channelSlope = 0.0025 -- Channel slope at the xdim only, in meters per meter. + model.channelWidth = 5 -- Channel width, in meters. It will be adjusted to an integer number of grids (ceil). + + model.channelDepth = 2.5 -- Additional depth at the highest border of the channel, in meters. + model.rainBegin = 0 -- Instant when when rain begins after the beginning of the simulatio, in hours. + model.rainDuration = 5 -- Rain duration, in hours. + model.rainRate = 0.005 -- Rate of rain, in mm per hour. + model.rainAll = false -- True if it rains also on the channel. + + -- Assumes an integer number of grids for channel greater than zero. + model.channelGridWidth = math.ceil(model.channelWidth / model.yCellSize) + + -- Calculates rain begin and end in time steps. + model.stepRainBegin = model.rainBegin // model.stepHours + 1 + model.stepRainEnd = model.stepRainBegin - 1 + model.rainDuration // model.stepHours + + -- Calculates rain rate in mm per time step. + model.stepRainRate = model.rainRate * model.stepHours + end, - -- Loads height and defines rivers for scenario. + -- Loads height and defines rivers for scenario (Rinaldi, 2007, figure 14). loadTerrain = function(model, cell) - -- Other parms specific for scenario. Use of local parms for improved performance. - local parm = model:otherParms() - - -- Defines cell.heigh and channel width (Rinaldi, 2007, figure 14). - - -- Calculates channel width in number of cells. - local channelGridWidth = 0 - if not parm.channelWidth then - -- If no channel is defined, than lower border of the plain is all open for outflow. - cell.open = (cell.x == 0) - elseif parm.channelWidth == 0 then - -- If channel defined as zero, then assumes 1 grid. - channelGridWidth = 1 - elseif parm.channelWidth ~= 0 then - -- If channel defined and not zero, assumes an integer number of grids (ceil). - channelGridWidth = math.ceil(parm.channelWidth / model.yCellSize) - end - -- If cell is on the channel: - if cell.y >= model.ydim - channelGridWidth then - -- Channel behave as a river. - cell.isRiver = true - -- Height of cell in the channel calculated from its slope on xdim and doesn't depend on cell.y - cell.height = cell.x * parm.channelSlope * model.xCellSize - -- In this scenario, the lowest end of the channel is an open border. - cell.open = (cell.x == 0) - - -- If cell is not on the channel, cell is on the plain: - else cell.height = - -- Plain reachs its lowest level at the border of the channel and doesn't depende on cell.x: - (model.ydim - 1 - channelGridWidth - cell.y) * parm.plainSlope * model.yCellSize + - -- The lowest plain border is above the highest channel border... - (model.xdim - 1) * parm.channelSlope * model.xCellSize + - -- ... and an additional depth at the highest border of the channel may be present. - (parm.channelDepth or 0) + if cell.y >= model.ydim - model.channelGridWidth + then + -- Channel behave as a river. + cell.isRiver = true + -- Height of cell in the channel calculated from its slope on xdim and doesn't depend on cell.y + cell.height = cell.x * model.channelSlope * model.xCellSize + -- The lowest end of the channel is an open border. + cell.open =(cell.x == 0) + -- If cell is not on the channel, cell is on the plain: + else cell.height = + -- Plain reachs its lowest level at the border of the channel and doesn't depende on cell.x: + (model.ydim - 1 - model.channelGridWidth - cell.y) * model.plainSlope * model.yCellSize + + -- The lowest plain border is above the highest channel border... + (model.xdim - 1) * model.channelSlope * model.xCellSize + + -- ... and an additional depth at the highest border of the channel may be present. + (model.channelDepth or 0) end + end, -- Loads precipitation for the scenario. loadPrecipitation = function (model, cell, t) - -- Other parms specific for scenario. Use of local parms for improved performance. - local parm = model:otherParms() - - -- Calculate channel width in number of cells. - local channelGridWidth = 0 - if not parm.channelWidth then - -- If no channel is defined, than lower border of the plain is all open for outflow. - cell.open = (cell.x == 0) - elseif parm.channelWidth == 0 then - -- If channel defined as zero, then assumes 1 grid. - channelGridWidth = 1 - elseif parm.channelWidth ~= 0 then - -- If channel defined and not zero, assumes an integer number of cells (ceil). - channelGridWidth = math.ceil(parm.channelWidth / model.yCellSize) - end - - -- Calculate rain begin and end in time steps. - local rainBegin = parm.rainBegin // model.stepHours + 1 - local rainEnd = rainBegin - 1 + parm.rainDuration // model.stepHours - - -- Calculate rain rate in mm per time step. - local rainRate = parm.rainRate * model.stepHours - - if t >= rainBegin and t <= rainEnd and - -- Depending on a scenario option, it rains all over the scenario, ... - (parm.rainAll or - -- ... or only on the plain cells - cell.y <= model.ydim - 1 - channelGridWidth) + -- Makes it rain during the period set. + if t <= model.stepRainEnd and t >= model.stepRainBegin and + -- Depending on a scenario option, it rains all over the scenario or only on the plain. + (cell.y <= model.ydim - 1 - model.channelGridWidth or model.rainAll) -- Fixed rain rate during the period for all cells in the region: - then return rainRate + then return model.stepRainRate -- No rain out of the time interval: else return 0 end + end, } --** print(sessionInfo().time, ": SimplePlain instance created.") -print(sessionInfo().time, ": Run scenario.") +print(sessionInfo().time, ": Runs scenario.") simplePlain:run()