i have coded 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 *
import csv
## Globale Variablen ##
# kommen später ins main file, MÜSSEN vor dem import von diesem Modul definiert werden!)
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 ##
# Berechnung der Diffusivität
kappa = k/eta
# Berechnung der Diffusivität
c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))
## Wertebereich ##
xmin = 0
xmax = 50
xsteps = 1
x = np.arange(xmin,xmax,xsteps)
ymin = 0
ymax = 50
ysteps = 1
y = np.arange(ymin,ymax,ysteps)
#X, Y = np.meshgrid(x,y)
## Klassendefinition ##
# hier wird eine Bohrloch-Klasse definiert
class Bohrloch(object):
# __init__ erzeugt das eigentliche Objekt mit seinen Parametern
def __init__(self, xlage, ylage, tstart, q):
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
#==============================================================================
# Drücke
#==============================================================================
# getPressure erzeugt einen Array mit allen Werten für x,y im Zeitraum t
def getPressure(self, t):
if (t-self.tstart<0): # Förderzeit muss in Förderzeitraum liegen
print "Startpunkt liegt außerhalb des Förderzeitraumes!"
return ()
else:
# erzeugen des Abstandsvektors
self.r = np.sqrt((x-self.xlage)**2+(y-self.ylage)**2)
# Druckformel nach Rudnicki (1986)
self.P = (self.q/(rhof*4*pi*kappa))*(expn(1,self.r**2/(4*c*(t-self.tstart))))
# Druck wird direkt bei der Bohrung als "not a number" definiert
self.P[self.xlage] = np.nan
self.P[self.ylage] = np.nan
# Umrechnung des Druckes in [MPa]
self.z = self.P/1e6
# return gibt die Druckwerte aus (und speichert sie zwischen)
return self.z
The issue i have is that i get back an 1 dimensional array of self.z. The idea behind the code is the following:
I have an x- and and y-axis, e.g. with the ranges above. I try to get pressure values (depending from x,y) for each possible combination of these 1d-arrays, e.g. values for x,y = 1,1 ; 1,2 ; 1,3 and so on until 50 on each axis. If i test the dimension of my array self.z with ndim, i see that i only have covered one dimension. The desired output should be a two dimensional array with 50x50 entrys (or 49x49), or am i wrong? I'm really stuck to this problem, i think it should be simple to solve but i can't get it properly. Later on i want to plot a 3D Plot with these values, i know i have to use np.meshgrid for this issue. Do i have to use it also here?
I hope you guys get my idea and can help!
According to your explanation, the problem is in your definition of 'self.r' You want a list of the distances r between (x,y) and (xlage,ylage) for all combinations of x and y, i.e. a matrix. However, you have a vector which contains r belonging to (x1,y1), (x2,y2),(x3,y3) ... , i.e. you are missing the combinations (x1,y2), ... etc.
The solution is
# put this into section 'Wertebereich'
xx,yy = np.meshgrid(x,y)
# and change the function to
def getPressure(self, t):
if (t-self.tstart<0): # Förderzeit muss in Förderzeitraum liegen
print "Startpunkt liegt außerhalb des Förderzeitraumes!"
return ()
else:
# erzeugen des Abstandsvektors
self.r = np.sqrt((xx-self.xlage)**2+(yy-self.ylage)**2)
# Druckformel nach Rudnicki (1986)
self.P = (self.q/(rhof*4*pi*kappa))*(expn(1,self.r**2/(4*c*(t-self.tstart))))
# Druck wird direkt bei der Bohrung als "not a number" definiert
self.P = np.where((xx==self.xlage) & (yy==self.ylage),self.P,NaN)
# Umrechnung des Druckes in [MPa]
self.z = self.P/1e6
# return gibt die Druckwerte aus (und speichert sie zwischen)
return self.z
I can't check the code right now but it should do the job.
You can then plot it by using xx,yy and the output z, as described here