612 lines
24 KiB
Python
612 lines
24 KiB
Python
# SPDX-License-Identifier: GPL-2.0-or-later
|
|
|
|
from time import time
|
|
import unittest
|
|
import sys
|
|
import os
|
|
from random import random as rand, shuffle
|
|
import numpy as np
|
|
|
|
numexpr_available = False
|
|
|
|
|
|
def getmemsize():
|
|
return 0.0
|
|
|
|
|
|
def getptime():
|
|
return time()
|
|
|
|
|
|
class Grid:
|
|
|
|
def __init__(self, size=10, dtype=np.single):
|
|
self.center = np.zeros([size, size], dtype)
|
|
self.water = None
|
|
self.sediment = None
|
|
self.scour = None
|
|
self.flowrate = None
|
|
self.sedimentpct = None
|
|
self.sedimentpct = None
|
|
self.capacity = None
|
|
self.avalanced = None
|
|
self.minx = None
|
|
self.miny = None
|
|
self.maxx = None
|
|
self.maxy = None
|
|
self.zscale = 1
|
|
self.maxrss = 0.0
|
|
self.sequence = [0, 1, 2, 3]
|
|
self.watermax = 1.0
|
|
self.flowratemax = 1.0
|
|
self.scourmax = 1.0
|
|
self.sedmax = 1.0
|
|
self.scourmin = 1.0
|
|
|
|
def init_water_and_sediment(self):
|
|
if self.water is None:
|
|
self.water = np.zeros(self.center.shape, dtype=np.single)
|
|
if self.sediment is None:
|
|
self.sediment = np.zeros(self.center.shape, dtype=np.single)
|
|
if self.scour is None:
|
|
self.scour = np.zeros(self.center.shape, dtype=np.single)
|
|
if self.flowrate is None:
|
|
self.flowrate = np.zeros(self.center.shape, dtype=np.single)
|
|
if self.sedimentpct is None:
|
|
self.sedimentpct = np.zeros(self.center.shape, dtype=np.single)
|
|
if self.capacity is None:
|
|
self.capacity = np.zeros(self.center.shape, dtype=np.single)
|
|
if self.avalanced is None:
|
|
self.avalanced = np.zeros(self.center.shape, dtype=np.single)
|
|
|
|
def __str__(self):
|
|
return ''.join(self.__str_iter__(fmt="%.3f"))
|
|
|
|
def __str_iter__(self, fmt):
|
|
for row in self.center[::]:
|
|
values = []
|
|
for v in row:
|
|
values.append(fmt % v)
|
|
yield ' '.join(values) + '\n'
|
|
|
|
@staticmethod
|
|
def fromFile(filename):
|
|
if filename == '-':
|
|
filename = sys.stdin
|
|
g = Grid()
|
|
g.center = np.loadtxt(filename, np.single)
|
|
return g
|
|
|
|
def toFile(self, filename, fmt="%.3f"):
|
|
if filename == '-':
|
|
filename = sys.stdout.fileno()
|
|
with open(filename, "w") as f:
|
|
for line in self.__str_iter__(fmt):
|
|
f.write(line)
|
|
|
|
def raw(self, format="%.3f"):
|
|
fstr = format + " " + format + " " + format + " "
|
|
a = self.center / self.zscale
|
|
minx = 0.0 if self.minx is None else self.minx
|
|
miny = 0.0 if self.miny is None else self.miny
|
|
maxx = 1.0 if self.maxx is None else self.maxx
|
|
maxy = 1.0 if self.maxy is None else self.maxy
|
|
dx = (maxx - minx) / (a.shape[0] - 1)
|
|
dy = (maxy - miny) / (a.shape[1] - 1)
|
|
for row in range(a.shape[0] - 1):
|
|
row0 = miny + row * dy
|
|
row1 = row0 + dy
|
|
for col in range(a.shape[1] - 1):
|
|
col0 = minx + col * dx
|
|
col1 = col0 + dx
|
|
yield (fstr % (row0, col0, a[row][col]) +
|
|
fstr % (row0, col1, a[row][col + 1]) +
|
|
fstr % (row1, col0, a[row + 1][col]) + "\n")
|
|
yield (fstr % (row0, col1, a[row][col + 1]) +
|
|
fstr % (row1, col0, a[row + 1][col]) +
|
|
fstr % (row1, col1, a[row + 1][col + 1]) + "\n")
|
|
|
|
def toRaw(self, filename, infomap=None):
|
|
with open(filename if type(filename) == str else sys.stdout.fileno(), "w") as f:
|
|
f.writelines(self.raw())
|
|
if infomap:
|
|
with open(os.path.splitext(filename)[0] + ".inf" if type(filename) == str else sys.stdout.fileno(), "w") as f:
|
|
f.writelines("\n".join("%-15s: %s" % t for t in sorted(infomap.items())))
|
|
|
|
@staticmethod
|
|
def fromRaw(filename):
|
|
"""initialize a grid from a Blender .raw file.
|
|
currently supports just rectangular grids of all triangles
|
|
"""
|
|
g = Grid.fromFile(filename)
|
|
# we assume tris and an axis aligned grid
|
|
g.center = np.reshape(g.center, (-1, 3))
|
|
g._sort()
|
|
return g
|
|
|
|
def _sort(self, expfact):
|
|
# keep unique vertices only by creating a set and sort first on x then on y coordinate
|
|
# using rather slow python sort but couldn't wrap my head around np.lexsort
|
|
verts = sorted(list({tuple(t) for t in self.center[::]}))
|
|
x = set(c[0] for c in verts)
|
|
y = set(c[1] for c in verts)
|
|
nx = len(x)
|
|
ny = len(y)
|
|
self.minx = min(x)
|
|
self.maxx = max(x)
|
|
self.miny = min(y)
|
|
self.maxy = max(y)
|
|
xscale = (self.maxx - self.minx) / (nx - 1)
|
|
yscale = (self.maxy - self.miny) / (ny - 1)
|
|
# note: a purely flat plane cannot be scaled
|
|
if (yscale != 0.0) and (abs(xscale / yscale) - 1.0 > 1e-3):
|
|
raise ValueError("Mesh spacing not square %d x %d %.4f x %4.f" % (nx, ny, xscale, yscale))
|
|
self.zscale = 1.0
|
|
if abs(yscale) > 1e-6:
|
|
self.zscale = 1.0 / yscale
|
|
|
|
# keep just the z-values and null any offset
|
|
# we might catch a reshape error that will occur if nx*ny != # of vertices
|
|
# (if we are not dealing with a heightfield but with a mesh with duplicate
|
|
# x,y coords, like an axis aligned cube
|
|
self.center = np.array([c[2] for c in verts], dtype=np.single).reshape(nx, ny)
|
|
self.center = (self.center - np.amin(self.center)) * self.zscale
|
|
if self.rainmap is not None:
|
|
rmscale = np.max(self.center)
|
|
self.rainmap = expfact + (1 - expfact) * (self.center / rmscale)
|
|
|
|
@staticmethod
|
|
def fromBlenderMesh(me, vg, expfact):
|
|
g = Grid()
|
|
g.center = np.asarray(list(tuple(v.co) for v in me.vertices), dtype=np.single)
|
|
g.rainmap = None
|
|
if vg is not None:
|
|
for v in me.vertices:
|
|
vg.add([v.index], 0.0, 'ADD')
|
|
g.rainmap = np.asarray(list((v.co[0], v.co[1], vg.weight(v.index)) for v in me.vertices), dtype=np.single)
|
|
g._sort(expfact)
|
|
return g
|
|
|
|
def setrainmap(self, rainmap):
|
|
self.rainmap = rainmap
|
|
|
|
def _verts(self, surface):
|
|
a = surface / self.zscale
|
|
minx = 0.0 if self.minx is None else self.minx
|
|
miny = 0.0 if self.miny is None else self.miny
|
|
maxx = 1.0 if self.maxx is None else self.maxx
|
|
maxy = 1.0 if self.maxy is None else self.maxy
|
|
dx = (maxx - minx) / (a.shape[0] - 1)
|
|
dy = (maxy - miny) / (a.shape[1] - 1)
|
|
for row in range(a.shape[0]):
|
|
row0 = miny + row * dy
|
|
for col in range(a.shape[1]):
|
|
col0 = minx + col * dx
|
|
yield (row0, col0, a[row][col])
|
|
|
|
def _faces(self):
|
|
nrow, ncol = self.center.shape
|
|
for row in range(nrow - 1):
|
|
for col in range(ncol - 1):
|
|
vi = row * ncol + col
|
|
yield (vi, vi + ncol, vi + 1)
|
|
yield (vi + 1, vi + ncol, vi + ncol + 1)
|
|
|
|
def toBlenderMesh(self, me):
|
|
# pass me as argument so that we don't need to import bpy and create a dependency
|
|
# the docs state that from_pydata takes iterators as arguments but it will
|
|
# fail with generators because it does len(arg)
|
|
me.from_pydata(list(self._verts(self.center)), [], list(self._faces()))
|
|
|
|
def toWaterMesh(self, me):
|
|
# pass me as argument so that we don't need to import bpy and create a dependency
|
|
# the docs state that from_pydata takes iterators as arguments but it will
|
|
# fail with generators because it does len(arg)
|
|
me.from_pydata(list(self._verts(self.water)), [], list(self._faces()))
|
|
|
|
def peak(self, value=1):
|
|
nx, ny = self.center.shape
|
|
self.center[int(nx / 2), int(ny / 2)] += value
|
|
|
|
def shelf(self, value=1):
|
|
nx, ny = self.center.shape
|
|
self.center[:nx / 2] += value
|
|
|
|
def mesa(self, value=1):
|
|
nx, ny = self.center.shape
|
|
self.center[nx / 4:3 * nx / 4, ny / 4:3 * ny / 4] += value
|
|
|
|
def random(self, value=1):
|
|
self.center += np.random.random_sample(self.center.shape) * value
|
|
|
|
def neighborgrid(self):
|
|
self.up = np.roll(self.center, -1, 0)
|
|
self.down = np.roll(self.center, 1, 0)
|
|
self.left = np.roll(self.center, -1, 1)
|
|
self.right = np.roll(self.center, 1, 1)
|
|
|
|
def zeroedge(self, quantity=None):
|
|
c = self.center if quantity is None else quantity
|
|
c[0, :] = 0
|
|
c[-1, :] = 0
|
|
c[:, 0] = 0
|
|
c[:, -1] = 0
|
|
|
|
def diffuse(self, Kd, IterDiffuse, numexpr):
|
|
self.zeroedge()
|
|
c = self.center[1:-1, 1:-1]
|
|
up = self.center[:-2, 1:-1]
|
|
down = self.center[2:, 1:-1]
|
|
left = self.center[1:-1, :-2]
|
|
right = self.center[1:-1, 2:]
|
|
if(numexpr and numexpr_available):
|
|
self.center[1:-1, 1:-1] = ne.evaluate('c + Kd * (up + down + left + right - 4.0 * c)')
|
|
else:
|
|
self.center[1:-1, 1:-1] = c + (Kd / IterDiffuse) * (up + down + left + right - 4.0 * c)
|
|
self.maxrss = max(getmemsize(), self.maxrss)
|
|
return self.center
|
|
|
|
def avalanche(self, delta, iterava, prob, numexpr):
|
|
self.zeroedge()
|
|
c = self.center[1:-1, 1:-1]
|
|
up = self.center[:-2, 1:-1]
|
|
down = self.center[2:, 1:-1]
|
|
left = self.center[1:-1, :-2]
|
|
right = self.center[1:-1, 2:]
|
|
where = np.where
|
|
|
|
if(numexpr and numexpr_available):
|
|
self.center[1:-1, 1:-1] = ne.evaluate('c + where((up -c) > delta ,(up -c -delta)/2, 0) \
|
|
+ where((down -c) > delta ,(down -c -delta)/2, 0) \
|
|
+ where((left -c) > delta ,(left -c -delta)/2, 0) \
|
|
+ where((right-c) > delta ,(right-c -delta)/2, 0) \
|
|
+ where((up -c) < -delta,(up -c +delta)/2, 0) \
|
|
+ where((down -c) < -delta,(down -c +delta)/2, 0) \
|
|
+ where((left -c) < -delta,(left -c +delta)/2, 0) \
|
|
+ where((right-c) < -delta,(right-c +delta)/2, 0)')
|
|
else:
|
|
sa = (
|
|
# incoming
|
|
where((up - c) > delta, (up - c - delta) / 2, 0)
|
|
+ where((down - c) > delta, (down - c - delta) / 2, 0)
|
|
+ where((left - c) > delta, (left - c - delta) / 2, 0)
|
|
+ where((right - c) > delta, (right - c - delta) / 2, 0)
|
|
# outgoing
|
|
+ where((up - c) < -delta, (up - c + delta) / 2, 0)
|
|
+ where((down - c) < -delta, (down - c + delta) / 2, 0)
|
|
+ where((left - c) < -delta, (left - c + delta) / 2, 0)
|
|
+ where((right - c) < -delta, (right - c + delta) / 2, 0)
|
|
)
|
|
randarray = np.random.randint(0, 100, sa.shape) * 0.01
|
|
sa = where(randarray < prob, sa, 0)
|
|
self.avalanced[1:-1, 1:-1] = self.avalanced[1:-1, 1:-1] + sa / iterava
|
|
self.center[1:-1, 1:-1] = c + sa / iterava
|
|
|
|
self.maxrss = max(getmemsize(), self.maxrss)
|
|
return self.center
|
|
|
|
def rain(self, amount=1, variance=0, userainmap=False):
|
|
self.water += (1.0 - np.random.random(self.water.shape) * variance) * \
|
|
(amount if ((self.rainmap is None) or (not userainmap)) else self.rainmap * amount)
|
|
|
|
def spring(self, amount, px, py, radius):
|
|
# px, py and radius are all fractions
|
|
nx, ny = self.center.shape
|
|
rx = max(int(nx * radius), 1)
|
|
ry = max(int(ny * radius), 1)
|
|
px = int(nx * px)
|
|
py = int(ny * py)
|
|
self.water[px - rx:px + rx + 1, py - ry:py + ry + 1] += amount
|
|
|
|
def river(self, Kc, Ks, Kdep, Ka, Kev, numexpr):
|
|
zeros = np.zeros
|
|
where = np.where
|
|
min = np.minimum
|
|
max = np.maximum
|
|
abs = np.absolute
|
|
arctan = np.arctan
|
|
sin = np.sin
|
|
|
|
center = (slice(1, -1, None), slice(1, -1, None))
|
|
up = (slice(None, -2, None), slice(1, -1, None))
|
|
down = (slice(2, None, None), slice(1, -1, None))
|
|
left = (slice(1, -1, None), slice(None, -2, None))
|
|
right = (slice(1, -1, None), slice(2, None, None))
|
|
|
|
water = self.water
|
|
rock = self.center
|
|
sediment = self.sediment
|
|
height = rock + water
|
|
|
|
# !! this gives a runtime warning for division by zero
|
|
verysmallnumber = 0.0000000001
|
|
water += verysmallnumber
|
|
sc = where(water > verysmallnumber, sediment / water, 0)
|
|
|
|
sdw = zeros(water[center].shape)
|
|
svdw = zeros(water[center].shape)
|
|
sds = zeros(water[center].shape)
|
|
angle = zeros(water[center].shape)
|
|
for d in (up, down, left, right):
|
|
if(numexpr and numexpr_available):
|
|
hdd = height[d]
|
|
hcc = height[center]
|
|
dw = ne.evaluate('hdd-hcc')
|
|
inflow = ne.evaluate('dw > 0')
|
|
wdd = water[d]
|
|
wcc = water[center]
|
|
# nested where() represent min() and max()
|
|
dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0')
|
|
sdw = ne.evaluate('sdw + dw')
|
|
scd = sc[d]
|
|
scc = sc[center]
|
|
rockd = rock[d]
|
|
rockc = rock[center]
|
|
sds = ne.evaluate('sds + dw * where(inflow, scd, scc)')
|
|
svdw = ne.evaluate('svdw + abs(dw)')
|
|
angle = ne.evaluate('angle + arctan(abs(rockd-rockc))')
|
|
else:
|
|
dw = (height[d] - height[center])
|
|
inflow = dw > 0
|
|
dw = where(inflow, min(water[d], dw), max(-water[center], dw)) / 4.0
|
|
sdw = sdw + dw
|
|
sds = sds + dw * where(inflow, sc[d], sc[center])
|
|
svdw = svdw + abs(dw)
|
|
angle = angle + np.arctan(abs(rock[d] - rock[center]))
|
|
|
|
if numexpr and numexpr_available:
|
|
wcc = water[center]
|
|
scc = sediment[center]
|
|
rcc = rock[center]
|
|
water[center] = ne.evaluate('wcc + sdw')
|
|
sediment[center] = ne.evaluate('scc + sds')
|
|
sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)')
|
|
fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw')
|
|
ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)')
|
|
rock[center] = ne.evaluate('rcc - ds')
|
|
# there isn't really a bottom to the rock but negative values look ugly
|
|
rock[center] = ne.evaluate('where(rcc<0,0,rcc)')
|
|
sediment[center] = ne.evaluate('scc + ds')
|
|
else:
|
|
wcc = water[center]
|
|
scc = sediment[center]
|
|
rcc = rock[center]
|
|
water[center] = wcc * (1 - Kev) + sdw
|
|
sediment[center] = scc + sds
|
|
sc = where(wcc > 0, scc / wcc, 2 * Kc)
|
|
fKc = Kc * svdw
|
|
ds = where(fKc > sc, (fKc - sc) * Ks, (fKc - sc) * Kdep) * wcc
|
|
self.flowrate[center] = svdw
|
|
self.scour[center] = ds
|
|
self.sedimentpct[center] = sc
|
|
self.capacity[center] = fKc
|
|
sediment[center] = scc + ds + sds
|
|
|
|
def flow(self, Kc, Ks, Kz, Ka, numexpr):
|
|
zeros = np.zeros
|
|
where = np.where
|
|
min = np.minimum
|
|
max = np.maximum
|
|
abs = np.absolute
|
|
arctan = np.arctan
|
|
sin = np.sin
|
|
|
|
center = (slice(1, -1, None), slice(1, -1, None))
|
|
rock = self.center
|
|
ds = self.scour[center]
|
|
rcc = rock[center]
|
|
rock[center] = rcc - ds * Kz
|
|
# there isn't really a bottom to the rock but negative values look ugly
|
|
rock[center] = where(rcc < 0, 0, rcc)
|
|
|
|
def rivergeneration(
|
|
self,
|
|
rainamount,
|
|
rainvariance,
|
|
userainmap,
|
|
Kc,
|
|
Ks,
|
|
Kdep,
|
|
Ka,
|
|
Kev,
|
|
Kspring,
|
|
Kspringx,
|
|
Kspringy,
|
|
Kspringr,
|
|
numexpr,
|
|
):
|
|
self.init_water_and_sediment()
|
|
self.rain(rainamount, rainvariance, userainmap)
|
|
self.zeroedge(self.water)
|
|
self.zeroedge(self.sediment)
|
|
self.river(Kc, Ks, Kdep, Ka, Kev, numexpr)
|
|
self.watermax = np.max(self.water)
|
|
|
|
def fluvial_erosion(
|
|
self,
|
|
rainamount,
|
|
rainvariance,
|
|
userainmap,
|
|
Kc,
|
|
Ks,
|
|
Kdep,
|
|
Ka,
|
|
Kspring,
|
|
Kspringx,
|
|
Kspringy,
|
|
Kspringr,
|
|
numexpr,
|
|
):
|
|
self.flow(Kc, Ks, Kdep, Ka, numexpr)
|
|
self.flowratemax = np.max(self.flowrate)
|
|
self.scourmax = np.max(self.scour)
|
|
self.scourmin = np.min(self.scour)
|
|
self.sedmax = np.max(self.sediment)
|
|
|
|
def analyze(self):
|
|
self.neighborgrid()
|
|
# just looking at up and left to avoid needless double calculations
|
|
slopes = np.concatenate((np.abs(self.left - self.center), np.abs(self.up - self.center)))
|
|
return '\n'.join(["%-15s: %.3f" % t for t in [
|
|
('height average', np.average(self.center)),
|
|
('height median', np.median(self.center)),
|
|
('height max', np.max(self.center)),
|
|
('height min', np.min(self.center)),
|
|
('height std', np.std(self.center)),
|
|
('slope average', np.average(slopes)),
|
|
('slope median', np.median(slopes)),
|
|
('slope max', np.max(slopes)),
|
|
('slope min', np.min(slopes)),
|
|
('slope std', np.std(slopes))
|
|
]])
|
|
|
|
|
|
class TestGrid(unittest.TestCase):
|
|
|
|
def test_diffuse(self):
|
|
g = Grid(5)
|
|
g.peak(1)
|
|
self.assertEqual(g.center[2, 2], 1.0)
|
|
g.diffuse(0.1, numexpr=False)
|
|
for n in [(2, 1), (2, 3), (1, 2), (3, 2)]:
|
|
self.assertAlmostEqual(g.center[n], 0.1)
|
|
self.assertAlmostEqual(g.center[2, 2], 0.6)
|
|
|
|
def test_diffuse_numexpr(self):
|
|
g = Grid(5)
|
|
g.peak(1)
|
|
g.diffuse(0.1, numexpr=False)
|
|
h = Grid(5)
|
|
h.peak(1)
|
|
h.diffuse(0.1, numexpr=True)
|
|
self.assertEqual(list(g.center.flat), list(h.center.flat))
|
|
|
|
def test_avalanche_numexpr(self):
|
|
g = Grid(5)
|
|
g.peak(1)
|
|
g.avalanche(0.1, numexpr=False)
|
|
h = Grid(5)
|
|
h.peak(1)
|
|
h.avalanche(0.1, numexpr=True)
|
|
print(g)
|
|
print(h)
|
|
np.testing.assert_almost_equal(g.center, h.center)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
|
|
import argparse
|
|
|
|
parser = argparse.ArgumentParser(description='Erode a terrain while assuming zero boundary conditions.')
|
|
parser.add_argument('-I', dest='iterations', type=int, default=1, help='the number of iterations')
|
|
parser.add_argument('-Kd', dest='Kd', type=float, default=0.01, help='Diffusion constant')
|
|
parser.add_argument('-Kh', dest='Kh', type=float, default=6, help='Maximum stable cliff height')
|
|
parser.add_argument('-Kp', dest='Kp', type=float, default=0.1, help='Avalanche probability for unstable cliffs')
|
|
parser.add_argument('-Kr', dest='Kr', type=float, default=0.1, help='Average amount of rain per iteration')
|
|
parser.add_argument('-Kspring', dest='Kspring', type=float, default=0.0,
|
|
help='Average amount of wellwater per iteration')
|
|
parser.add_argument('-Kspringx', dest='Kspringx', type=float, default=0.5, help='relative x position of spring')
|
|
parser.add_argument('-Kspringy', dest='Kspringy', type=float, default=0.5, help='relative y position of spring')
|
|
parser.add_argument('-Kspringr', dest='Kspringr', type=float, default=0.02, help='radius of spring')
|
|
parser.add_argument('-Kdep', dest='Kdep', type=float, default=0.1, help='Sediment deposition constant')
|
|
parser.add_argument('-Ks', dest='Ks', type=float, default=0.1, help='Soil softness constant')
|
|
parser.add_argument('-Kc', dest='Kc', type=float, default=1.0, help='Sediment capacity')
|
|
parser.add_argument('-Ka', dest='Ka', type=float, default=2.0, help='Slope dependency of erosion')
|
|
parser.add_argument(
|
|
'-ri',
|
|
action='store_true',
|
|
dest='rawin',
|
|
default=False,
|
|
help='use Blender raw format for input')
|
|
parser.add_argument(
|
|
'-ro',
|
|
action='store_true',
|
|
dest='rawout',
|
|
default=False,
|
|
help='use Blender raw format for output')
|
|
parser.add_argument('-i', action='store_true', dest='useinputfile', default=False,
|
|
help='use an inputfile (instead of just a synthesized grid)')
|
|
parser.add_argument(
|
|
'-t',
|
|
action='store_true',
|
|
dest='timingonly',
|
|
default=False,
|
|
help='do not write anything to an output file')
|
|
parser.add_argument('-infile', type=str, default="-", help='input filename')
|
|
parser.add_argument('-outfile', type=str, default="-", help='output filename')
|
|
parser.add_argument('-Gn', dest='gridsize', type=int, default=20, help='Gridsize (always square)')
|
|
parser.add_argument('-Gp', dest='gridpeak', type=float, default=0, help='Add peak with given height')
|
|
parser.add_argument('-Gs', dest='gridshelf', type=float, default=0, help='Add shelve with given height')
|
|
parser.add_argument('-Gm', dest='gridmesa', type=float, default=0, help='Add mesa with given height')
|
|
parser.add_argument('-Gr', dest='gridrandom', type=float, default=0,
|
|
help='Add random values between 0 and given value')
|
|
parser.add_argument('-m', dest='threads', type=int, default=1, help='number of threads to use')
|
|
parser.add_argument('-u', action='store_true', dest='unittest', default=False, help='perform unittests')
|
|
parser.add_argument('-a', action='store_true', dest='analyze', default=False,
|
|
help='show some statistics of input and output meshes')
|
|
parser.add_argument('-d', action='store_true', dest='dump', default=False,
|
|
help='show sediment and water meshes at end of run')
|
|
parser.add_argument('-n', action='store_true', dest='usenumexpr', default=False, help='use numexpr optimizations')
|
|
|
|
args = parser.parse_args()
|
|
print("\nInput arguments:")
|
|
print("\n".join("%-15s: %s" % t for t in sorted(vars(args).items())), file=sys.stderr)
|
|
|
|
if args.unittest:
|
|
unittest.main(argv=[sys.argv[0]])
|
|
sys.exit(0)
|
|
|
|
if args.useinputfile:
|
|
if args.rawin:
|
|
grid = Grid.fromRaw(args.infile)
|
|
else:
|
|
grid = Grid.fromFile(args.infile)
|
|
else:
|
|
grid = Grid(args.gridsize)
|
|
|
|
if args.gridpeak > 0:
|
|
grid.peak(args.gridpeak)
|
|
if args.gridmesa > 0:
|
|
grid.mesa(args.gridmesa)
|
|
if args.gridshelf > 0:
|
|
grid.shelf(args.gridshelf)
|
|
if args.gridrandom > 0:
|
|
grid.random(args.gridrandom)
|
|
|
|
if args.analyze:
|
|
print('\nstatistics of the input grid:\n\n', grid.analyze(), file=sys.stderr, sep='')
|
|
t = getptime()
|
|
for g in range(args.iterations):
|
|
if args.Kd > 0:
|
|
grid.diffuse(args.Kd, args.usenumexpr)
|
|
if args.Kh > 0 and args.Kp > rand():
|
|
grid.avalanche(args.Kh, args.usenumexpr)
|
|
if args.Kr > 0 or args.Kspring > 0:
|
|
grid.fluvial_erosion(
|
|
args.Kr,
|
|
args.Kc,
|
|
args.Ks,
|
|
args.Kdep,
|
|
args.Ka,
|
|
args.Kspring,
|
|
args.Kspringx,
|
|
args.Kspringy,
|
|
args.Kspringr,
|
|
args.usenumexpr,
|
|
)
|
|
t = getptime() - t
|
|
print("\nElapsed time: %.1f seconds, max memory %.1f Mb.\n" % (t, grid.maxrss), file=sys.stderr)
|
|
if args.analyze:
|
|
print('\nstatistics of the output grid:\n\n', grid.analyze(), file=sys.stderr, sep='')
|
|
|
|
if not args.timingonly:
|
|
if args.rawout:
|
|
grid.toRaw(args.outfile, vars(args))
|
|
else:
|
|
grid.toFile(args.outfile)
|
|
|
|
if args.dump:
|
|
print("sediment\n", np.array_str(grid.sediment, precision=3), file=sys.stderr)
|
|
print("water\n", np.array_str(grid.water, precision=3), file=sys.stderr)
|
|
print("sediment concentration\n", np.array_str(grid.sediment / grid.water, precision=3), file=sys.stderr)
|