在类对象内嵌套for循环和3D绘图 [英] Nested for loop and 3D Plot within Class Object
问题描述
我确定这是一个容易处理的问题，但我不能弄清楚。我创建了一个钻孔类，并计算我的每个钻孔/井周围的孔隙压力。沿着单个轴，我的代码如下所示：
i am sure this is an easy problem to deal with, but i cant figure it out. I created a Borehole Class and want to compute my pore pressure around each Borehole/Well. Along a single axis, my code looks like this:
from scipy.special import *
import matplotlib.pyplot as plt
import numpy as np
from math import *
## Globale Variablen ##
rhof = 1000 # Dichte Flüssigkeit [kg/m³]
lameu = 11.2*10**9 # Lamé-Parameter, undrained [GPa]
lame = 8.4*10**9 # Lamé-Parameter, drained [GPa]
pi # durch Pythonmodul "math" gegeben
alpha = 0.65 # Biot-Willis-Koeffizient
G = 8.4*10**9 # Schermodul [GPa]
k = 1.0e-15 # Permeabilität [m²] bzw. [Darcy]
eta = 0.001 # Viskosität des Fluids [Pa*s]
## Berechnung der Parameter ##
kappa = k/eta
c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))
## Wertebereich ##
xmin = 0
xmax = 100
xsteps = 1.0
x = np.arange(xmin, xmax, xsteps)
## Class ##
class Bohrloch(object):
loch_zaehler = 0
def __init__(self, xlage, tstart, q): # Funktion, um BL zu erzeugen
self.xlage = xlage
#self.ylage = ylage # Lage der Bohrung
self.tstart = tstart # Start der Injektion/Produktion
self.q = q # Fluidmenge
## Druck ##
def getPressure(self, t): # gibt nach Zeit t die zugehörigen Druckwerte aus
if (t-self.tstart<0): # Fehlermeldung, falls Startpunkt nach t liegt
return ()
print "Startpunkt liegt außerhalb des Förderzeitraumes!"
else:
self.r = np.sqrt((x-self.xlage)**2)
self.P = (self.q/(rhof*4*pi*kappa))*(expn(1,self.r**2/(4*c*(t-self.tstart))))
#self.P[self.xlage] = 0 # gibt Bohrlochlage wieder
self.z = self.P/1e6
return self.z # Druckwerte in [MPa]
def pressureTable (self, t, xschritt): # erstellt Wertetabelle
self.getPressure(t)
for i in range (xmin, xmax, xschritt):
print i, " ", self.z[i]
t = 1000*24*3600
b1 = Bohrloch(50,0*24*3600,6.0/1000)
b1.pressureTable(t,1)
使用这种方法，我得到我想要的压力表。现在我想有一个压力表的x和y值，包括一个3D绘图。这是我的代码到目前为止：
With this method i get my desired pressure table. Now i want to have a pressure table for x and y values, including an 3D Plot. This is my code so far:
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from scipy.special import *
import matplotlib.pyplot as plt
import numpy as np
from math import *
## Globale Variablen ##
rhof = 1000 # Dichte Flüssigkeit [kg/m³]
lameu = 11.2*10**9 # Lamé-Parameter, undrained [GPa]
lame = 8.4*10**9 # Lamé-Parameter, drained [GPa]
pi # durch Pythonmodul "math" gegeben
alpha = 0.65 # Biot-Willis-Koeffizient
G = 8.4*10**9 # Schermodul [GPa]
k = 1.0e-15 # Permeabilität [m²] bzw. [Darcy]
eta = 0.001 # Viskosität des Fluids [Pa*s]
## Berechnung der Parameter ##
kappa = k/eta
c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))
## Wertebereich ##
xmin = 0
xmax = 100
xsteps = 1.0
x = np.arange(xmin,xmax,xsteps)
ymin = 0
ymax = 100
ysteps = 1.0
y = np.arange(ymin,ymax,ysteps)
## Klassendefinition ##
class Bohrloch(object):
loch_zaehler = 0
def __init__(self, xlage, ylage, tstart, q): # Funktion, um BL zu erzeugen
self.xlage = xlage # x-Lage der Bohrung
self.ylage = ylage # y-Lage der Bohrung
self.tstart = tstart # Start der Injektion/Produktion
self.q = q # Fluidmenge
## Druck ##
def getPressure(self, t):
if (t-self.tstart<0):
return ()
print "Startpunkt liegt außerhalb des Förderzeitraumes!"
else:
self.r = np.sqrt((x-self.xlage)**2+(y-self.ylage)**2)
self.P = (self.q/(rhof*4*pi*kappa))*(expn(1,self.r**2/(4*c*(t-self.tstart))))
self.P[self.xlage] = np.nan
self.P[self.ylage] = np.nan
self.z = self.P/1e6
return self.z # Druckwerte in [MPa]
def pressureTable (self, t, xschritt, yschritt):
self.getPressure(t)
for k in range (xmin, xmax, xschritt):
for l in range (ymin, ymax, yschritt):
# my mistake should be here?
print k, " ", l, " ", self.z[k][l]
def pressurePlot3D (self, t):
self.getPressure(t)
Z = self.z
X, Y = np.meshgrid(x,y)
Z[Z == np.inf] = np.nan
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0,
antialiased=False, vmin=np.nanmin(Z), vmax=np.nanmax(Z))
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.set_xlim(xmin,xmax) # x-Achsenskala vorgeben
ax.set_ylim(ymin,ymax) # y-Achsenskala vorgeben
ax.set_title('Druckverteilung')
ax.set_xlabel('x-Richtung [m]')
ax.set_ylabel('y-Richtung Well [m]')
ax.set_zlabel('Druck in [MPa]')
plt.show()
t = 1000*24*3600
b1 = Bohrloch(50,50,0*24*3600,6.0/1000)
b1.pressureTable(t,1)
b1.pressurePlot3D(t)
不幸的是，我的表不工作，剧情看起来很奇怪。我仍然是Python的初学者，需要一些建议。
Unfortunately, my table doesnt work and the desired 3D Plot looks strange. I am still a total beginner in Python and need some advices.
任何人都可以帮助？
推荐答案
问题是 self.z 不是一个二维数组/列表。因此，尝试访问 self.z [k] [l] 会导致 IndexError：无效的标量变量索引
The problem is that self.z is not a two-dimensional array/list. Therefore, trying to access self.z[k][l] results in IndexError: invalid index to scalar variable.
-
我不太明白要如何实现第二个维度。您介绍y位置，但是，您只需使用
I do not quite understand how you want to implement the second dimension. You introduce the y-position, but then, you just calculate a 1D radius array by using both the x- and y-location in
self.r = np.sqrt((x-self.xlage)**2+(y-self.ylage)**2)
下一个问题是，您打算做什么：
The next question is, what do you intend with:
self.P[self.xlage] = np.nan
self.P[self.ylage] = np.nan
b $ b
如果将 xsteps 和 ysteps 更改为 10 ，然后致电：
If you change xsteps and ysteps to 10, and call:
b1 = Bohrloch(2,3,0*24*3600,6.0/1000)
print b1.getPressure(t)
您的输出将是：
[ 5.44152501 4.40905986 nan nan 2.87481753 2.64950827
2.46756653 2.31503845 2.18379093 2.06866598]
为什么要用 nan 替换第3和第4个元素？
Why would you want to replace the 3rd and 4th elements with nan?
这些问题也是您制图程序的基础。因为您现在在数组中有 np.nan 值，这些将不会显示在绘图中。因为 self.z 不是二维的，你可能没有得到你可能期待的表面：
These issues are also at the basis of your plotting routine. Because you now have np.nan values in your array, these won't show in the plot. Because self.z is not two-dimensional, you are probably not getting the surface you may be expecting:
这里是一个简单的2D实现的方法。我不太熟悉你想要做什么，但它得到了想法：
Here's a simple way of coming up with a 2D implementation. I am not familiar enough with what you are trying to do, but it gets the idea across:
def getPressure(self, t):
if (t-self.tstart<0):
return ()
print "Startpunkt liegt außerhalb des Förderzeitraumes!"
else:
# you need to initialize r, P and z as list of lists
# make this dependent on your x coordinates
# the second dimension will grow dynamically
self.r = [[] for ri in range(len(x))]
self.P = [[] for ri in range(len(x))]
self.z = [[] for ri in range(len(x))]
# iterate through both x and y independently
for ii in range(len(x)):
for jj in range(len(y)):
# append to the list that corresponds to the current x -value
# also, use x[ii] and y[jj] to call one x-, y-value at a time
self.r[ii].append(np.sqrt((x[ii]-self.xlage)**2+(y[jj]-self.ylage)**2))
# calling r[ii][-1] ensures you are using the value that was last added to the list:
self.P[ii].append((self.q/(rhof*4*pi*kappa))*(expn(1,self.r[ii][-1]**2/(4*c*(t-self.tstart)))))
self.z[ii].append(self.P[ii][-1]/1e6)
# now, you can use xlage and ylage to blank one value
# do this for both P and z, because z is now calculated inside the loop
self.P[self.xlage][self.ylage] = np.nan
self.z[self.xlage][self.ylage] = np.nan
return self.z
从绘制例程中删除此行： Z [Z == np.inf] = np.nan ，使用原始命令：
From your plotting routine, remove this line: Z[Z == np.inf] = np.nan, use your original command:
b1 = Bohrloch(50,50,0*24*3600,6.0/1000)
b1.pressurePlot3D(t)
，您现在将得到这个情节：
and you will now get this plot:
这篇关于在类对象内嵌套for循环和3D绘图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
