Nested for loop and 3D Plot within Class Object

Question

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)

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)

Unfortunately, my table doesnt work and the desired 3D Plot looks strange. I am still a total beginner in Python and need some advices.

Can anyone help?

Answer 1

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.

---

• 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
  

  If you change xsteps and ysteps to 10, and call:

  b1 = Bohrloch(2,3,0*24*3600,6.0/1000)
  print b1.getPressure(t)  
  

  Your output will be:

  [ 5.44152501 4.40905986        nan        nan  2.87481753  2.64950827
    2.46756653 2.31503845 2.18379093 2.06866598]
  

  Why would you want to replace the 3rd and 4th elements with nan?

---

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:

---

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

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: