<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 BuildingPit Element
Step1:  Parameters
Define some parameters
Step2:  Model
Set up a model
Step3:  Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour.
Step4:  Plot an overview of the model
Step5:  Visualizations
Step6:  Model 2
Step7:  Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
# import sys
# sys.path.insert(1, "..")
import timml
Explanation: BuildingPit Element
End of explanation
kh = 2.  # m/day
f_ani = 0.05  # anisotropy factor
kv = f_ani*kh
ctop = 800.  # resistance top leaky layer in days
ztop = 0. # surface elevation
z_well = -13.  # end depth of the wellscreen
z_dw = -15.  # bottom elevation of sheetpile wall
z_extra = z_dw - 15.  # extra layer
zbot = -60.  # bottom elevation of the model
l = 40.  # length building pit in m
b = 30.   # width building pit in m
h_bem = -6.21  # m
offset = 5.  # distance groundwater extraction element from sheetpiles in m
xy = [(-l/2, -b/2), (l/2, -b/2), (l/2, b/2), (-l/2, b/2), (-l/2, -b/2)]
for (x, y) in xy:
    p2, = plt.plot(x, y, "o", label="building pit pts")
plt.axis("equal");
plt.show()
Explanation: Parameters
Define some parameters
End of explanation
z = np.array([ztop+1, ztop, z_dw, z_dw, z_extra, z_extra, zbot])
dz = z[1::2] - z[2::2]
dz
kh_arr = kh * np.ones(dz.shape)
c = np.r_[np.array([ctop]), dz[:-1]/(2*kv) + dz[1:]/(2*kv)]
c
Explanation: Model
Set up a model
End of explanation
ml = timml.ModelMaq(kaq=kh_arr, z=z, c=c, topboundary="semi", hstar=0.0, f2py=False)
layers = np.arange(np.sum(z_dw <= ml.aq.zaqbot))
last_lay_dw = layers[-1]
inhom = timml.BuildingPit(ml, xy, kaq=kh_arr, z=z[1:], topboundary="conf", 
                          c=c[1:], order=4, ndeg=3, layers=layers)
timml.HeadLineSink(ml, x1=-l/2+offset, y1=b/2-offset, x2=l/2-offset, y2=b/2-offset, hls=h_bem, 
                   layers=np.arange(last_lay_dw+1))
timml.HeadLineSink(ml, x1=-l/2+offset, y1=0, x2=l/2-offset, y2=0, hls=h_bem, 
                   layers=np.arange(last_lay_dw+1))
timml.HeadLineSink(ml, x1=-l/2+offset, y1=-b/2+offset, x2=l/2-offset, y2=-b/2+offset, hls=h_bem, 
                   layers=np.arange(last_lay_dw+1))
# ml.solve_mp(nproc=2)
ml.solve()
Qtot = 0.
for e in ml.elementlist:
    if e.name == "HeadLineSink":
        Qtot += e.discharge()
print("Debiet =", np.round(Qtot.sum(), 2), "m3/dag")
y = np.linspace(-b/2-25, b/2+1100, 201)
hl = ml.headalongline(np.zeros(201), y, layers=[0])
y_5cm = np.interp(-0.05, ml.headalongline(np.zeros(201), y, layers=0).squeeze(), y, right=np.nan)
print("Distance to 5 cm drawdown contour =", np.round(y_5cm, 2), "m")
# Q_arr[mi] = Qtot.sum()
# y5cm_arr[mi] = y_5cm
Explanation: Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour.
End of explanation
ml.plot()
Explanation: Plot an overview of the model
End of explanation
x = np.linspace(-l/2-25, l/2+1100, 201)
hl = ml.headalongline(x, np.zeros(201), layers=[last_lay_dw, last_lay_dw+1])
fig, ax = plt.subplots(1, 1, figsize=(16, 5))
ax.plot(x, hl[0].squeeze(), label="head layer {}".format(last_lay_dw))
ax.plot(x, hl[1].squeeze(), label="head layer {}".format(last_lay_dw+1))
ax.axhline(-0.05, color="r", linestyle="dashed", lw=0.75, label="-0.05 m")
ax.axhline(-0.5, color="k", linestyle="dashed", lw=0.75, label="-0.5 m")
ax.set_xlabel("x (m)")
ax.set_ylabel("head (m)")
ax.legend(loc="best")
ax.grid(b=True)
xoffset = 15
zoffset = 15
x1, x2, y1, y2 = [-l/2-xoffset, -l/2+xoffset, 0, 0]
nudge = 1e-6
n = 301
# plot head contour cross-sections
h = ml.headalongline(np.linspace(x1 + nudge, x2 - nudge, n),
                     np.linspace(y1 + nudge, y2 - nudge, n))
L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
xg = np.linspace(0, L, n)+x1
zg = 0.5 * (ml.aq.zaqbot + ml.aq.zaqtop)
zg = np.hstack((ml.aq.zaqtop[0], zg, ml.aq.zaqbot[-1]))
h = np.vstack((h[0], h, h[-1]))
levels = np.linspace(h_bem, -0.0, 51)
fig, ax = plt.subplots(1, 1, figsize=(16, 8))
ml.plot(win=[x1, x2, y1, y2], orientation="ver", newfig=False)
# ml.vcontoursf1D(x1, x2, n, levels=101, newfig=False, ax=ax, color="r")
cf = ax.contourf(xg, zg, h, levels)
cs = ax.contour(xg, zg, h, levels, colors="k", linewidths=0.5)
# cs2 = ax.contour(xg, zg2, h2, levels, colors="r", linewidths=0.5, linestyles="dashed")
# plt.clabel(cs, fmt="%.2f")
# plt.clabel(cs2, fmt="%.2f")
ax.set_ylim(z_dw-zoffset, z_dw+zoffset)
ax.set_ylabel("diepte (m NAP)");
ax.set_xlabel("m")
ax.set_aspect("equal")
plt.colorbar(cf, ax=ax)
plt.show()
Explanation: Visualizations
End of explanation
n = 11  # aantal laagjes boven en onder damwand
dz_i_top = (z_well-z_dw)/np.sum(np.arange(n+1))
dz_i_bot = (z_dw-z_extra)/np.sum(np.arange(2*n+1))
z_layers_top = np.cumsum(np.arange(n)*dz_i_top)
z_layers_bot = np.cumsum(np.arange(2*n)*dz_i_bot)
zgr = np.r_[z_dw + z_layers_top[::-1], (z_dw-z_layers_bot)[1:]]
z4 = np.r_[np.array([ztop+1, ztop, z_well, z_well]), np.repeat(zgr, 2, 0), z_extra, z_extra, zbot]
# z4 = np.r_[np.array([ztop+1, ztop, z_well]), np.repeat(zgr, 2, 0)[1:-1]]
dz4 = z4[1:-1:2] - z4[2::2]
kh_arr = kh * np.ones(dz4.shape)
c = np.r_[np.array([ctop]), dz4[:-1]/(2*kv) + dz4[1:]/(2*kv)]
kh_arr2 = kh_arr.copy()
kh_arr2[0] = 1e-5
Explanation: Model 2: Add more layers
Add more layers to the model to get a more accurate solution of the flow towards the building pit.
End of explanation
ml = timml.ModelMaq(kaq=kh_arr, z=z4, c=c, topboundary="semi", hstar=0.0)
layers = np.arange(np.sum(z_dw <= ml.aq.zaqbot))
last_lay_dw = layers[-1]
inhom = timml.BuildingPit(ml, xy, kaq=kh_arr2, z=z4[1:], topboundary="conf", 
                          c=c[1:], order=4, ndeg=3, layers=layers)
# wlayers = np.arange(np.sum(z_well <= ml.aq.zaqbot))
wlayers = np.arange(np.sum(-14 <= ml.aq.zaqbot))
wlayers=wlayers[1:]
timml.HeadLineSink(ml, x1=-l/2+offset, y1=b/2-offset, x2=l/2-offset, y2=b/2-offset, hls=h_bem, 
                   layers=wlayers)
timml.HeadLineSink(ml, x1=-l/2+offset, y1=0, x2=l/2-offset, y2=0, hls=h_bem, 
                   layers=wlayers, order=5)
timml.HeadLineSink(ml, x1=-l/2+offset, y1=-b/2+offset, x2=l/2-offset, y2=-b/2+offset, hls=h_bem, 
                   layers=wlayers);
# ml.solve_mp(nproc=2)
ml.solve()
Qtot = 0.
for e in ml.elementlist:
    if e.name == "HeadLineSink":
        Qtot += e.discharge()
print("Debiet =", np.round(Qtot.sum(), 2), "m3/dag")
y = np.linspace(-b/2-25, b/2+1100, 201)
hl = ml.headalongline(np.zeros(201), y, layers=[0])
y_5cm = np.interp(-0.05, ml.headalongline(np.zeros(201), y, layers=0).squeeze(), y, right=np.nan)
print("Distance to 5 cm drawdown contour =", np.round(y_5cm, 2), "m")
last_lay_dw = layers[-1]
x = np.linspace(-l/2-25, l/2+1, 201)
# x = np.linspace(-l/2-25, l/2+1100, 201)
hl = ml.headalongline(x, np.zeros(201), layers=[0, last_lay_dw, last_lay_dw+1])
fig, ax = plt.subplots(1, 1, figsize=(16, 5))
ax.plot(x, hl[0].squeeze(), label="head layer 0")
ax.plot(x, hl[1].squeeze(), label="head layer {}".format(last_lay_dw))
ax.plot(x, hl[2].squeeze(), label="head layer {}".format(last_lay_dw+1))
ax.axhline(-0.05, color="r", linestyle="dashed", lw=0.75, label="-0.05 m")
ax.axhline(-0.5, color="k", linestyle="dashed", lw=0.75, label="-0.5 m")
ax.set_xlabel("x (m)")
ax.set_ylabel("head (m)")
ax.legend(loc="best")
ax.grid(b=True)
xoffset = 50
zoffset = 15
x1, x2, y1, y2 = [-l/2-xoffset, -l/2+xoffset, 0, 0]
nudge = 1e-6
n = 301
# plot head contour cross-sections
h = ml.headalongline(np.linspace(x1 + nudge, x2 - nudge, n),
                     np.linspace(y1 + nudge, y2 - nudge, n))
L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
xg = np.linspace(0, L, n) + x1
zg = 0.5 * (ml.aq.zaqbot + ml.aq.zaqtop)
zg = np.hstack((ml.aq.zaqtop[0], zg, ml.aq.zaqbot[-1]))
h = np.vstack((h[0], h, h[-1]))
levels = np.linspace(h_bem-.1, -0.0, 51)
fig, ax = plt.subplots(1, 1, figsize=(16, 8))
ml.plot(win=[x1, x2, y1, y2], orientation="ver", newfig=False)
# ml.vcontoursf1D(x1, x2, n, levels=101, newfig=False, ax=ax, color="r")
cf = ax.contourf(xg, zg, h, levels)
cs = ax.contour(xg, zg, h, levels, colors="k", linewidths=0.5)
# cs2 = ax.contour(xg, zg2, h2, levels, colors="r", linewidths=0.5, linestyles="dashed")
# plt.clabel(cs, fmt="%.2f")
# plt.clabel(cs2, fmt="%.2f")
ax.set_ylim(z_dw-zoffset, z_dw+zoffset)
ax.set_ylabel("depth (m NAP)");
ax.set_xlabel("m");
# ax.set_aspect("equal")
cb = plt.colorbar(cf, ax=ax)
cb.ax.set_ylabel("stijghoogte (m)")
# ax.set_ylim(-20, -5)
# ax.set_xlim(-70, -57)
Explanation: Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour.
End of explanation
<END_TASK>