Skip to content

Commit

Permalink
Simplify solver - Create street canyon only if the cavity zone touche…
Browse files Browse the repository at this point in the history 
… an other block than the one giving birth to it
  • Loading branch information
@j3r3m1
j3r3m1 committed Nov 22, 2023
1 parent a90290e commit 40e25ff
Show file tree
Hide file tree
Showing 5 changed files with 139 additions and 95 deletions.
8 changes: 5 additions & 3 deletions functions/URock/CalculatesIndicators.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,8 @@ def zoneProperties(cursor, obstaclePropertiesTable, prefix = PREFIX_NAME):
1.*3*1.8*{9}/(POWER({8}/{3},0.3)*(1+0.24*{9}/{3})) AS {6},
0.22*(0.67*LEAST({3},{9})+0.33*GREATEST({3},{9})) AS {11},
0.9*(0.67*LEAST({3},{9})+0.33*GREATEST({3},{9})) AS {12},
1+0.05*{9}/{3} AS {13}
1+0.05*{9}/{3} AS {13},
{14}
FROM {7}""".format(tempoStackedLengthTab,
ID_FIELD_STACKED_BLOCK,
GEOM_FIELD,
Expand All @@ -157,7 +158,8 @@ def zoneProperties(cursor, obstaclePropertiesTable, prefix = PREFIX_NAME):
DISPLACEMENT_LENGTH_VORTEX_FIELD,
ROOFTOP_PERP_HEIGHT,
ROOFTOP_PERP_LENGTH,
ROOFTOP_WIND_FACTOR)
ROOFTOP_WIND_FACTOR,
ID_FIELD_BLOCK)
cursor.execute(query)

# Calculates the table containing the points corresponding to all polygons,
Expand Down Expand Up @@ -252,7 +254,7 @@ def zoneProperties(cursor, obstaclePropertiesTable, prefix = PREFIX_NAME):
COS_BLOCK_LEFT_AZIMUTH , SIN_BLOCK_LEFT_AZIMUTH,
COS_BLOCK_RIGHT_AZIMUTH , SIN_BLOCK_RIGHT_AZIMUTH,
tempoStackedLengthTab , stackedBlockAzimuths,
ID_FIELD_STACKED_BLOCK))
ID_FIELD_STACKED_BLOCK , ID_FIELD_BLOCK))

if not DEBUG:
# Drop intermediate tables
Expand Down
18 changes: 9 additions & 9 deletions functions/URock/InitWindField.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -381,13 +381,12 @@ def affectsPointToBuildZone(cursor, gridTable, dicOfBuildRockleZoneTable,
STREET_CANYON_NAME : f"""b.{idZone[STREET_CANYON_NAME]},
b.{UPSTREAM_HEIGHT_FIELD},
LEAST(b.{DOWNSTREAM_HEIGHT_FIELD}, b.{UPSTREAM_HEIGHT_FIELD}) AS {MAX_CANYON_HEIGHT_FIELD},
b.{DOWNSTREAM_HEIGHT_FIELD}-b.{UPSTREAM_HEIGHT_FIELD} AS CANYON_DELTAH_FIELD,
b.{DOWNSTREAM_HEIGHT_FIELD}-b.{UPSTREAM_HEIGHT_FIELD} AS {CANYON_DELTAH_FIELD},
b.{UPWIND_FACADE_ANGLE_FIELD},
a.{ID_POINT_X},
b.{BASE_HEIGHT_FIELD},
ST_YMAX(ST_INTERSECTION(a.{GEOM_FIELD},
b.{GEOM_FIELD}
) AS {Y_WALL},
b.{GEOM_FIELD})) AS {Y_WALL},
ST_LENGTH(ST_MAKELINE(ST_TOMULTIPOINT(ST_INTERSECTION(a.{GEOM_FIELD},
b.{GEOM_FIELD})
)
Expand Down Expand Up @@ -618,17 +617,17 @@ def affectsPointToBuildZone(cursor, gridTable, dicOfBuildRockleZoneTable,
DOWNSTREAM_X_RELATIVE_POSITION,
POINT_RELATIVE_POSITION_FIELD+WAKE_NAME[0]),
STREET_CANYON_NAME : """b.{0},
SIN(2*(a.{1}-PI()/2))*(0.5+(a.{9}-b.{13})*
(SIN(2*(a.{1}-PI()/2))*(0.5+(a.{9}-b.{13})*
(a.{2}-(a.{9}-b.{13}))/
(0.5*POWER(a.{2},2)))*
(0.5*POWER(a.{2},2))))*
(1+(0.6*a.{19})/(a.{12}+0.6*a.{19}))*
(POWER(a.{2}/a.{12},2)/(1+POWER(a.{2}/a.{12},2))) AS {3},
1-POWER(COS(a.{1}-PI()/2),2)*(1+(a.{9}-b.{13})*
(a.{2}-(a.{9}-b.{13}))/(POWER(0.5*a.{2},2)))*
(1-POWER(COS(a.{1}-PI()/2),2)*(1+(a.{9}-b.{13})*
(a.{2}-(a.{9}-b.{13}))/(POWER(0.5*a.{2},2))))*
(1+(0.6*a.{19})/(a.{12}+0.6*a.{19}))*
(POWER(a.{2}/a.{12},2)/(1+POWER(a.{2}/a.{12},2))) AS {4},
-ABS(0.5*(1-(a.{9}-b.{13})/(0.5*a.{2})))*
(1-(a.{2}-(a.{9}-b.{13}))/(0.5*a.{2}))*
(-ABS(0.5*(1-(a.{9}-b.{13})/(0.5*a.{2})))*
(1-(a.{2}-(a.{9}-b.{13}))/(0.5*a.{2})))*
(1+(0.6*a.{19})/(a.{12}+0.6*a.{19}))*
(POWER(a.{2}/a.{12},2)/(1+POWER(a.{2}/a.{12},2))) AS {5},
a.{6},
Expand Down Expand Up @@ -3030,6 +3029,7 @@ def setInitialWindField(cursor, initializedWindFactorTable, gridPoint,
'SELECT b.{14} - 1 AS {14}_MINUS_1,
b.{15} - 1 AS {15}_MINUS_1,
a.{2},
b.{5},
a.{8} * WIND_SPEED AS {8},
a.{6} * WIND_SPEED AS {6},
a.{9} * WIND_SPEED AS {9}
Expand Down
9 changes: 5 additions & 4 deletions functions/URock/Obstacles.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -502,7 +502,7 @@ def updateUpwindFacadeBase(cursor, upwindTable, prefix = PREFIX_NAME):
{10};
DROP TABLE IF EXISTS {3};
CREATE TABLE {3}
AS SELECT a.{1}, a.{4}, a.{5}, a.{7}, a.{8},
AS SELECT a.{1}, a.{4}, a.{5}, a.{7}, a.{8}, a.{11},
COALESCE(b.{6}, a.{6}) AS {6}
FROM {0} AS a LEFT JOIN {2} AS b ON a.{5} = b.{5}
""".format( upwindTable , ID_FIELD_STACKED_BLOCK,
Expand All @@ -514,7 +514,8 @@ def updateUpwindFacadeBase(cursor, upwindTable, prefix = PREFIX_NAME):
isSpatial=False),
DataUtil.createIndex(tableName=tempoUpwind,
fieldName=UPWIND_FACADE_FIELD,
isSpatial=False)))
isSpatial=False),
ID_FIELD_BLOCK))

if not DEBUG:
# Drop intermediate tables
Expand Down Expand Up @@ -591,7 +592,7 @@ def initDownwindFacades(cursor, obstaclesTable, prefix = PREFIX_NAME):
AS SELECT a.{1}, a.{2}, a.{4}, b.{8}, b.{9}, b.{10},
(ST_XMAX(b.{2}) + ST_XMIN(b.{2})) / 2 AS {12},
(ST_XMAX(b.{2}) - ST_XMIN(b.{2})) AS {13},
b.{14}, b.{15}, b.{16}, b.{17}, b.{18}
b.{14}, b.{15}, b.{16}, b.{17}, b.{18}, b.{19}
FROM {0} AS a LEFT JOIN {11} AS b
ON a.{4} = b.{4}
""".format(tempoDownwindLines , DOWNWIND_FACADE_FIELD,
Expand All @@ -608,7 +609,7 @@ def initDownwindFacades(cursor, obstaclesTable, prefix = PREFIX_NAME):
STACKED_BLOCK_X_MED , STACKED_BLOCK_WIDTH,
STACKED_BLOCK_UPSTREAMEST_X , SIN_BLOCK_LEFT_AZIMUTH,
COS_BLOCK_LEFT_AZIMUTH , COS_BLOCK_RIGHT_AZIMUTH,
SIN_BLOCK_RIGHT_AZIMUTH))
SIN_BLOCK_RIGHT_AZIMUTH , ID_FIELD_BLOCK))

if not DEBUG:
# Drop intermediate tables
Expand Down
134 changes: 86 additions & 48 deletions functions/URock/WindSolver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,15 @@ def solver(x, y, z, dx, dy, dz, u0, v0, w0, buildingCoordinates, cells4Solver, c
A = dx ** 2 / dy ** 2
B = eta ** 2 * dx ** 2 / dz ** 2

########################################################################
# # Used for debug only
# import matplotlib.pylab as plt
# fig, ax = plt.subplots()
# for i in range(buildingCoordinates[0].size):
# if buildingCoordinates[2][i] == 1:
# ax.plot(buildingCoordinates[0][i], buildingCoordinates[1][i], marker = "o")
########################################################################

# Set coefficients according to table 1 (Pardyjak et Brown, 2003)
# to modify the Equation near obstacles
e = np.ones([nx, ny, nz])
Expand All @@ -129,66 +138,95 @@ def solver(x, y, z, dx, dy, dz, u0, v0, w0, buildingCoordinates, cells4Solver, c
p = np.ones([nx, ny, nz])
q = np.ones([nx, ny, nz])

# Identify index having wall below AND front, left, right or behind
# Identify index having wall below AND (front, left, right or behind)
indBelow = pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0],
buildingCoordinates[1],
buildingCoordinates[2] + 1])))
indBelowFront = indBelow.intersection(pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0],
buildingCoordinates[1] - 1,
buildingCoordinates[2]]))))
indBelowBehind = indBelow.intersection(pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0],
buildingCoordinates[1] + 1,
buildingCoordinates[2]]))))
indBelowLeft = indBelow.intersection(pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0] + 1,
buildingCoordinates[1],
buildingCoordinates[2]]))))
indBelowRight = indBelow.intersection(pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0] - 1,
buildingCoordinates[1],
buildingCoordinates[2]]))))
indBelowAnyAround = indBelowFront.union(indBelowBehind).union(indBelowLeft).union(indBelowRight)
indBelowLeftRight = indBelowLeft.union(indBelowRight)
indBelowFrontBehind = indBelowFront.union(indBelowBehind)
indAbove = pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0],
buildingCoordinates[1],
buildingCoordinates[2] - 1])))
indFront = pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0],
buildingCoordinates[1] - 1,
buildingCoordinates[2]])))
indBehind = pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0],
buildingCoordinates[1] + 1,
buildingCoordinates[2]])))
indLeft = pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0] + 1,
buildingCoordinates[1],
buildingCoordinates[2]])))
indRight = pd.MultiIndex.from_tuples(list(zip(*[buildingCoordinates[0] - 1,
buildingCoordinates[1],
buildingCoordinates[2]])))
indBelowRight = indBelow.intersection(indRight)
indBelowLeft = indBelow.intersection(indLeft)
indBelowFront = indBelow.intersection(indFront)
indBelowBehind = indBelow.intersection(indBehind)

indE = indBelowRight.union(indRight)
indF = indBelowLeft.union(indLeft)
indG = indBelowFront.union(indFront)
indH = indBelowBehind.union(indBehind)
indM = indAbove
indN = indBelow.union(indBelowFront).union(indBelowLeft)\
.union(indBelowRight).union(indBelowBehind)
indO = indRight.union(indLeft).union(indBelowLeft).union(indBelowRight)
indP = indBehind.union(indFront).union(indBelowBehind).union(indBelowFront)
indQ = indBelow.union(indAbove).union(indBelowFront).union(indBelowLeft)\
.union(indBelowRight).union(indBelowBehind)

# Go descending order along y
if DESCENDING_Y:
e[buildingCoordinates[0] + 1, buildingCoordinates[1], buildingCoordinates[2]] = 0.
e[indBelowLeft.get_level_values(0), indBelowLeft.get_level_values(1), indBelowLeft.get_level_values(2)] = 0.
f[buildingCoordinates[0] - 1, buildingCoordinates[1], buildingCoordinates[2]] = 0.
f[indBelowRight.get_level_values(0), indBelowRight.get_level_values(1), indBelowRight.get_level_values(2)] = 0.
g[buildingCoordinates[0], buildingCoordinates[1] + 1, buildingCoordinates[2]] = 0.
g[indBelowBehind.get_level_values(0), indBelowBehind.get_level_values(1), indBelowBehind.get_level_values(2)] = 0.
h[buildingCoordinates[0], buildingCoordinates[1] - 1, buildingCoordinates[2]] = 0.
h[indBelowFront.get_level_values(0), indBelowFront.get_level_values(1), indBelowFront.get_level_values(2)] = 0.
m[buildingCoordinates[0], buildingCoordinates[1], buildingCoordinates[2] + 1] = 0.
n[buildingCoordinates[0], buildingCoordinates[1], buildingCoordinates[2] - 1] = 0.
e[indF.get_level_values(0), indF.get_level_values(1), indF.get_level_values(2)] = 0.
f[indE.get_level_values(0), indE.get_level_values(1), indE.get_level_values(2)] = 0.
g[indH.get_level_values(0), indH.get_level_values(1), indH.get_level_values(2)] = 0.
h[indG.get_level_values(0), indG.get_level_values(1), indG.get_level_values(2)] = 0.
m[indM.get_level_values(0), indM.get_level_values(1), indM.get_level_values(2)] = 0.
n[indN.get_level_values(0), indN.get_level_values(1), indN.get_level_values(2)] = 0.
else:
e[buildingCoordinates[0] - 1, buildingCoordinates[1], buildingCoordinates[2]] = 0.
e[indBelowRight.get_level_values(0), indBelowRight.get_level_values(1), indBelowRight.get_level_values(2)] = 0.
f[buildingCoordinates[0] + 1, buildingCoordinates[1], buildingCoordinates[2]] = 0.
f[indBelowLeft.get_level_values(0), indBelowLeft.get_level_values(1), indBelowLeft.get_level_values(2)] = 0.
g[buildingCoordinates[0], buildingCoordinates[1] - 1, buildingCoordinates[2]] = 0.
g[indBelowFront.get_level_values(0), indBelowFront.get_level_values(1), indBelowFront.get_level_values(2)] = 0.
h[buildingCoordinates[0], buildingCoordinates[1] + 1, buildingCoordinates[2]] = 0.
h[indBelowBehind.get_level_values(0), indBelowBehind.get_level_values(1), indBelowBehind.get_level_values(2)] = 0.
m[buildingCoordinates[0], buildingCoordinates[1], buildingCoordinates[2] - 1] = 0.
n[buildingCoordinates[0], buildingCoordinates[1], buildingCoordinates[2] + 1] = 0.
e[indE.get_level_values(0), indE.get_level_values(1), indE.get_level_values(2)] = 0.
f[indF.get_level_values(0), indF.get_level_values(1), indF.get_level_values(2)] = 0.
g[indG.get_level_values(0), indG.get_level_values(1), indG.get_level_values(2)] = 0.
h[indH.get_level_values(0), indH.get_level_values(1), indH.get_level_values(2)] = 0.
m[indM.get_level_values(0), indM.get_level_values(1), indM.get_level_values(2)] = 0.
n[indN.get_level_values(0), indN.get_level_values(1), indN.get_level_values(2)] = 0.

o[buildingCoordinates[0] - 1, buildingCoordinates[1], buildingCoordinates[2]] = 0.5
o[buildingCoordinates[0] + 1, buildingCoordinates[1], buildingCoordinates[2]] = 0.5
p[buildingCoordinates[0], buildingCoordinates[1] - 1, buildingCoordinates[2]] = 0.5
p[buildingCoordinates[0], buildingCoordinates[1] + 1, buildingCoordinates[2]] = 0.5
q[buildingCoordinates[0], buildingCoordinates[1], buildingCoordinates[2] + 1] = 0.5
q[buildingCoordinates[0], buildingCoordinates[1], buildingCoordinates[2] - 1] = 0.5
n[indBelowAnyAround.get_level_values(0), indBelowAnyAround.get_level_values(1), indBelowAnyAround.get_level_values(2)] = 0.
o[indBelowLeftRight.get_level_values(0), indBelowLeftRight.get_level_values(1), indBelowLeftRight.get_level_values(2)] = 0.5
p[indBelowFrontBehind.get_level_values(0), indBelowFrontBehind.get_level_values(1), indBelowFrontBehind.get_level_values(2)] = 0.5
q[indBelowAnyAround.get_level_values(0), indBelowAnyAround.get_level_values(1), indBelowAnyAround.get_level_values(2)] = 0.5

o[indO.get_level_values(0), indO.get_level_values(1), indO.get_level_values(2)] = 0.5
p[indP.get_level_values(0), indP.get_level_values(1), indP.get_level_values(2)] = 0.5
q[indQ.get_level_values(0), indQ.get_level_values(1), indQ.get_level_values(2)] = 0.5

for N in range(maxIterations):
print("Iteration {0} (max {1})".format( N + 1,
maxIterations))
lambdaN = np.copy(lambdaN1)


# ########################################################################
# # Only used for debug
# if DESCENDING_Y:
# for k, j, i in np.flip(cells4Solver):
# lambdaN1[i, j, k] = omega * (
# ((-1.) * (dx ** 2 * (-2. * alpha1 ** 2) * (((u0[i, j, k] - u0[i + 1, j, k]) / (dx) + (
# v0[i, j, k] - v0[i, j + 1, k]) / (dy) +
# (w0[i, j, k] - w0[i, j, k + 1]) / (dz)))) + (
# e[i, j, k] * lambdaN[i - 1, j, k] + f[i, j, k] * lambdaN1[i + 1, j, k] + A * (
# g[i, j, k] * lambdaN[i, j - 1, k] + h[i, j, k] * lambdaN1[i, j + 1, k]) + B * (
# m[i, j, k] * lambdaN[i, j, k - 1] + n[i, j, k] * lambdaN1[i, j, k + 1]))) / (
# 2. * (o[i, j, k] + A * p[i, j, k] + B * q[i, j, k]))) + (1 - omega) * lambdaN1[i, j, k]

# else:
# for i, j, k in cells4Solver:
# if i == 115 and j == 214:
# print("ok")
# lambdaN1[i, j, k] = omega * (
# ((-1.) * (dx ** 2 * (-2. * alpha1 ** 2) * (((u0[i + 1, j, k] - u0[i, j, k]) / (dx) + (
# v0[i, j + 1, k] - v0[i, j, k]) / (dy) +
# (w0[i, j, k + 1] - w0[i, j, k]) / (dz)))) + (
# e[i, j, k] * lambdaN[i + 1, j, k] + f[i, j, k] * lambdaN1[i - 1, j, k] + A * (
# g[i, j, k] * lambdaN[i, j + 1, k] + h[i, j, k] * lambdaN1[i, j - 1, k]) + B * (
# m[i, j, k] * lambdaN[i, j, k + 1] + n[i, j, k] * lambdaN1[i, j, k - 1]))) / (
# 2. * (o[i, j, k] + A * p[i, j, k] + B * q[i, j, k]))) + (1 - omega) * lambdaN1[i, j, k]
# # end of debug
# ########################################################################

lambdaN1 = calcLambda(cells4Solver, lambdaN, lambdaN1, omega, alpha1,
u0, v0, w0, dx, dy, dz, e, f, g, h, m, n, o, p, q,
DESCENDING_Y, A, B)
Expand Down
Loading

0 comments on commit 40e25ff

Please sign in to comment.