Test wells with resistance

This notebook checks the implementation for wells with specified screen resistance.

import numpy as np

import timflow.steady as tfs

# model parameters single layer
k = 10  # m/day
c = 1000.0  # resistance top leaky layer in days
z = [1.0, 0, -20]

res = 0
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res)
ml.solve()
x = 100
y = 0
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is {ml.head(x, y).item():.5f}")
print(f"head inside well: {w.headinside().item():.5f}")

Number of elements, Number of equations: 2 , 0
No unknowns. Solution complete
resistance: res=0.00000
head at x=100, y=0 is -0.13103
head inside well: -0.67812

res = 0.01
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res)
ml.solve()
hin = w.headinside()
hout = ml.head(w.xc, w.yc, w.layers)
Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res
x = 100
y = 0
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is {ml.head(x, y).item():.5f}")
print(f"head inside well: {w.headinside().item():.5f}")
print("discharge: ", w.discharge())
print("discharge from equation", Qcheck)

Number of elements, Number of equations: 2 , 0
No unknowns. Solution complete
resistance: res=0.01000
head at x=100, y=0 is -0.13103
head inside well: -0.75770
discharge:  [100.]
discharge from equation [100.]

# model parameters multiple layers
k = [10, 20, 40]  # m/day
c = [1000, 500, 800]  # resistance top leaky layer in days
z = [1.0, 0, -10, -12, -18, -20, -32]

res = 0
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res, layers=[1, 2])
ml.solve()
x = 100
y = 0
layer = 2
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is", ml.head(x, y))
print("head inside well is", w.headinside())
print("discharge:", w.discharge())

Number of elements, Number of equations: 2 , 2
.
.

solution complete
resistance: res=0.00000
head at x=100, y=0 is [-0.02343071 -0.05355064 -0.06735737]
head inside well is [-0.24736602 -0.24736602]
discharge: [ 0.         21.37172528 78.62827472]

res = 0.01
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.Well(ml, 0, 0, Qw=100, rw=0.1, res=res, layers=[1, 2])
ml.solve()
hin = w.headinside()
hout = ml.head(w.xc, w.yc, w.layers)
Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res
x = 100
y = 0
layer = 2
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is", ml.head(x, y))
print("head inside well is", w.headinside())
print("discharge:", w.discharge())
print("discharge from equation", Qcheck)

Number of elements, Number of equations: 2 , 2
.
.

solution complete
resistance: res=0.01000
head at x=100, y=0 is [-0.02451537 -0.05693459 -0.06594459]
head inside well is [-0.34014217 -0.34014217]
discharge: [ 0.         24.17931547 75.82068453]
discharge from equation [24.17931547 75.82068453]

HeadWell

# model parameters
k = 10  # m/day
c = 1000.0  # resistance top leaky layer in days
z = [1.0, 0, -20]

res = 0
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res)
ml.solve()
x = 100
y = 0
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is {ml.head(x, y).item():.5f}")
print(f"head inside well: {w.headinside().item():.5f}")
print("discharge:", w.discharge())

Number of elements, Number of equations: 2 , 1
.
.

solution complete
resistance: res=0.00000
head at x=100, y=0 is -0.38645
head inside well: -2.00000
discharge: [294.93146472]

res = 0.01
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res)
ml.solve()
hin = w.headinside()
hout = ml.head(w.xc, w.yc, w.layers)
Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res
x = 100
y = 0
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is {ml.head(x, y).item():.5f}")
print(f"head inside well: {w.headinside().item():.5f}")
print("discharge:", w.discharge())
print("discharge from equation", Qcheck)

Number of elements, Number of equations: 2 , 1
.
.

solution complete
resistance: res=0.01000
head at x=100, y=0 is -0.34587
head inside well: -2.00000
discharge: [263.95632199]
discharge from equation [263.95632199]

# model parameters multiple layers
k = [10, 20, 40]  # m/day
c = [1000, 500, 800]  # resistance top leaky layer in days
z = [1.0, 0, -10, -12, -18, -20, -32]

res = 0
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res, layers=[1, 2])
ml.solve()
x = 100
y = 0
layer = 2
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is", ml.head(x, y))
print("head inside well is", w.headinside())
print("discharge:", w.discharge())

Number of elements, Number of equations: 2 , 2
.
.

solution complete
resistance: res=0.00000
head at x=100, y=0 is [-0.18944163 -0.43296678 -0.54459675]
head inside well is [-2. -2.]
discharge: [  0.         172.79434715 635.72412705]

res = 0.01
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.HeadWell(ml, 0, 0, hw=-2, rw=0.1, res=res, layers=[1, 2])
ml.solve()
hin = w.headinside()
hout = ml.head(w.xc, w.yc, w.layers)
Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res
x = 100
y = 0
layer = 2
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is", ml.head(x, y))
print("head inside well is", w.headinside())
print("discharge:", w.discharge())
print("discharge from equation", Qcheck)

Number of elements, Number of equations: 2 , 2
.
.

solution complete
resistance: res=0.01000
head at x=100, y=0 is [-0.14414778 -0.33476934 -0.38774719]
head inside well is [-2. -2.]
discharge: [  0.         142.17181677 445.81760322]
discharge from equation [142.17181677 445.81760322]

TargetHeadWell

# model parameters single layer
k = 10  # m/day
c = 1000.0  # resistance top leaky layer in days
z = [1.0, 0, -20]

res = 0
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.TargetHeadWell(ml, 0, 0, rw=0.1, res=res, hcp=-0.4, xcp=100, ycp=200)
ml.solve()
x = 100
y = 200
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is {ml.head(x, y).item():.5f}")
print(f"head inside well: {w.headinside().item():.5f}")
print("discharge:", w.discharge())

Number of elements, Number of equations: 2 , 1
.
.

solution complete
resistance: res=0.00000
head at x=100, y=200 is -0.40000
head inside well: -3.68731
discharge: [543.7521144]

res = 0.01
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.TargetHeadWell(ml, 0, 0, rw=0.1, res=res, hcp=-0.4, xcp=100, ycp=200)
ml.solve()
hin = w.headinside()
hout = ml.head(w.xc, w.yc, w.layers)
Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res
x = 100
y = 200
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is {ml.head(x, y).item():.5f}")
print(f"head inside well: {w.headinside().item():.5f}")
print("discharge:", w.discharge())
print("discharge from equation", Qcheck)

Number of elements, Number of equations: 2 , 1
.
.

solution complete
resistance: res=0.01000
head at x=100, y=200 is -0.40000
head inside well: -4.12002
discharge: [543.7521144]
discharge from equation [543.7521144]

# model parameters multiple layers
k = [10, 20, 40]  # m/day
c = [1000, 500, 800]  # resistance top leaky layer in days
z = [1.0, 0, -10, -12, -18, -20, -32]

res = 0
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.TargetHeadWell(
    ml, 0, 0, rw=0.1, res=res, layers=[1, 2], hcp=-0.4, xcp=100, ycp=200, lcp=1
)
ml.solve()
x = 100
y = 200
layer = 2
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is", ml.head(x, y))
print("head inside well is", w.headinside())
print("discharge:", w.discharge())

Number of elements, Number of equations: 2 , 2
.
.

solution complete
resistance: res=0.00000
head at x=100, y=200 is [-0.22679647 -0.4        -0.54195309]
head inside well is [-2.86902936 -2.86902936]
discharge: [  0.         247.87602799 911.95559409]

res = 0.01
ml = tfs.ModelMaq(kaq=k, c=c, z=z, topboundary="semi", hstar=0)
w = tfs.TargetHeadWell(
    ml, 0, 0, rw=0.1, res=res, layers=[1, 2], hcp=-0.4, xcp=100, ycp=200, lcp=1
)
ml.solve()
hin = w.headinside()
hout = ml.head(w.xc, w.yc, w.layers)
Qcheck = ml.aq.Haq[w.layers] * 2 * np.pi * w.rw * (hout - hin) / res
x = 100
y = 200
layer = 2
print(f"resistance: res={res:.5f}")
print(f"head at x={x}, y={y} is", ml.head(x, y))
print("head inside well is", w.headinside())
print("discharge:", w.discharge())
print("discharge from equation", Qcheck)

Number of elements, Number of equations: 2 , 2
.
.

solution complete
resistance: res=0.01000
head at x=100, y=200 is [-0.226388   -0.4        -0.51338359]
head inside well is [-3.8014717 -3.8014717]
discharge: [  0.         270.23106888 847.38150063]
discharge from equation [270.23106888 847.38150063]