/[escript]/trunk/finley/test/python/coalgas.py
ViewVC logotype

Annotation of /trunk/finley/test/python/coalgas.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3478 - (hide annotations)
Wed Mar 23 04:02:39 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 13599 byte(s)
coal gas stuff

1 gross 3458 #######################################################
2     #
3     # Copyright (c) 2003-2010 by University of Queensland
4     # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6     #
7     # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10     #
11     ########################################################
12     """
13     Gas in Coal Seam (fully coupled version)
14     """
15     __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16     Earth Systems Science Computational Center (ESSCC)
17     http://www.uq.edu.au/esscc
18     Primary Business: Queensland, Australia"""
19     __license__="""Licensed under the Open Software License version 3.0
20     http://www.opensource.org/licenses/osl-3.0.php"""
21     __url__="https://launchpad.net/escript-finley"
22    
23     from esys.finley import Rectangle, Brick
24     from esys.escript import unitsSI as U
25    
26 gross 3478 class MaterialProperty(object):
27     """
28     generic class for material properties depending on one (or several) status
29     variable(s)
30     """
31     def __init__(self, *args, **kwargs):
32     """
33     set up
34     """
35     pass
36     def __call__(self, *args, **kwargs):
37     """
38     return value of material property for given status
39     """
40     return self.getValue( *args, **kwargs)
41    
42    
43     def getValue(self, *args, **kwargs):
44     """
45     return value of material property for given status
46    
47     :remark: needs to be overwritten
48     """
49     raise NotImplementedError
50 gross 3458
51 gross 3478 class MaterialPropertyWithDifferential(MaterialProperty):
52     """
53     generic class for material properties depending on one (or several) status
54     variable(s) where the derivative of the property with respect to the
55     status variables is available
56     """
57     def getValueDifferential(self, *args, **kwargs):
58     """
59     returns the value of the derivative of material property for given status variable
60    
61     :remark: needs to be overwritten
62     """
63     raise NotImplementedError
64    
65     class InterpolationTable(object):
66     """
67     a simple 1D interpolation table for escript Data object with non-equidistant nodes
68     """
69     def __init__(self,x=[], y=[], obeyBounds=True):
70     """
71     set up interpolation table. obeyBounds is set an exception is thrown if
72     the interpolation argument is below min(x) or larger than max(x). Otherwise
73     the value for x is set to y[0] for
74     """
75     if len(x) != len(y):
76     raise ValueError,"length of interpolation nodes and value lists need to be identical."
77     if len(x) > 0:
78     raise ValueError,"length of interpolation nodes a list needs to at least one."
79    
80     x_ref=x[0]
81     for i in range(1,len(x)):
82     if x_ref >= x[i]:
83     raise ValueError,"interpolation nodes need to be given in increasing order."
84     x_ref=x[i]
85     self.__x = x
86     self.__y = y
87     self.__obeyBounds = obeyBounds
88     def getValue(self, x):
89     """
90     returns the interpolated values of x
91     """
92     x0=self.__x[0]
93     m0=whereNegative(x-x0)
94     if self.__obeyBounds:
95     if sup(m0) > 0:
96     raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
97     out=self.__x[0]
98     for i in range(1,len(self.__x)):
99     x1=self.__x[i]
100     z=(y[i]-y[i-1])/(x[i]-x[i-1]) * (x-x[i-1]) + y[i-1]
101     out = out * m0 + z * (1-m0)
102     m0=whereNegative(x-x[i])
103    
104     if self.__obeyBounds:
105     if inf(m0) < 1:
106     raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
107     else:
108     out = out * m0 + y[-1] * (1-m0)
109     return out
110 gross 3458
111 gross 3478 class Well(object):
112     """
113     generic well
114    
115     :var WATER: phase identifier for water well
116     :var GAS: phase identifier for gas well
117     """
118     WATER="Water"
119     GAS="Gas"
120     def __init__(self, *args, **kwargs):
121     """
122     set up well
123     """
124     pass
125     def getGeometry(self):
126     """
127     returns a function representing the geometry of the well.
128     typically a Gaussian profile around the location of the well.
129    
130     :note: needs to be overwritten
131     """
132     raise NotImplementedError
133 gross 3458
134 gross 3478 def getProductivityIndex(self):
135     """
136     returns the productivity index of the well.
137     typically a Gaussian profile around the location of the well.
138    
139     :note: needs to be overwritten
140     """
141     raise NotImplementedError
142    
143     def getFlowRate(self):
144     """
145     returns the flow rate
146     """
147     return 0.
148    
149     def getBHP(self):
150     """
151     return bottom-hole pressure
152    
153     :note: needs to be overwritten
154     """
155     raise NotImplementedError
156    
157     def getPhase(self):
158     """
159     returns the pahse the well works on
160     """
161     return self.WATER
162 gross 3458
163 gross 3478 class PeacemanWell(Well):
164     """
165     defines a well using Peaceman's formula
166     """
167     pass
168 gross 3458
169 gross 3478 class DualPorosity(object):
170     """
171     generic dual porosity model
172     """
173     def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
174     S_fg=None, S_mg=None,
175     perm_f_0=None, perm_f_1=None, perm_f_2=None,
176     perm_m_0=None, perm_m_1=None, perm_m_2=None,
177     k_w =None, k_g=None, mu_w =None, mu_g =None,
178     rho_w =None, rho_g=None,
179     wells=[]):
180     """
181     set up
182    
183     :param domain: domain
184     :type domain: `esys.escript.Domain`
185     :param phi_f: porosity of the fractured rock as function of the gas pressure
186     :type phi_f: `MaterialPropertyWithDifferential`
187     :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
188     :type phi_m: `MaterialPropertyWithDifferential` or None
189     :param phi: total porosity if None phi=1.
190     :type phi: scalar or None
191     :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=p_fg-p_fw
192     :type S_fg: `MaterialPropertyWithDifferential`
193     :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
194     :type S_mg: `MaterialPropertyWithDifferential`
195     :param perm_f_0: permeability in the x_0 direction in the fractured media
196     :type perm_f_0: scalar
197     :param perm_f_1: permeability in the x_1 direction in the fractured media
198     :type perm_f_1: scalar
199     :param perm_f_2: permeability in the x_2 direction in the fractured media
200     :type perm_f_2: scalar
201     :param perm_m_0: permeability in the x_0 direction in the coal matrix
202     :type perm_m_0: scalar
203     :param perm_m_1: permeability in the x_1 direction in the coal matrix
204     :type perm_m_1: scalar
205     :param perm_m_2: permeability in the x_2 direction in the coal matrix
206     :type perm_m_2: scalar
207     :param k_w: relative permeability of water as function of water saturation
208     :type k_w: `MaterialProperty`
209     :param k_g: relative permeability of gas as function of gas saturation
210     :type k_g: `MaterialProperty`
211     :param mu_w: viscosity of water as function of water pressure
212     :type mu_w: `MaterialProperty`
213     :param mu_g: viscosity of gas as function of gas pressure
214     :type mu_g: `MaterialProperty`
215     :param rho_w: density of water as function of water pressure
216     :type rho_w: `MaterialPropertyWithDifferential`
217     :param rho_g: density of gas as function of gas pressure
218     :type rho_g: `MaterialPropertyWithDifferential`
219     :param wells : list of wells
220     :type wells: list of `Well`
221     """
222 gross 3458
223 gross 3478 self.domain = domain
224     self.phi_f = phi_f
225     self.phi_m = phi_m
226     self.phi = phi
227     self.S_fg = S_fg
228     self.S_mg = S_mg
229     self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
230     self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
231     self.k_w = k_w
232     self.k_g = k_g
233     self.mu_w = mu_w
234     self.mu_g = mu_g
235     self.rho_w = rho_w
236     self.rho_g = rho_g
237     self.wells=wells
238    
239 gross 3458
240    
241    
242    
243 gross 3478 class PorosityOneHalfModel(DualPorosity):
244     """
245     Model for gas and water in double prosity model tracking water and gas
246     pressure in the fractured rock and gas concentration in the coal matrix.
247     This corresponds to the coal bed methan model in the ECLIPSE code.
248     """
249 gross 3458
250 gross 3478 def __init__(self, domain, phi_f=None, phi=None, L_g=None, S_fg=None,
251     perm_f_0=None, perm_f_1=None, perm_f_2=None,
252     k_w =None, k_g=None, mu_w =None, mu_g =None,
253     rho_w =None, rho_g=None,
254     wells=[]):
255     """
256     set up
257    
258     :param domain: domain
259     :type domain: `esys.escript.Domain`
260     :param phi_f: porosity of the fractured rock as function of the gas pressure
261     :type phi_f: `MaterialPropertyWithDifferential`
262     :param phi: total porosity if None phi=1.
263     :type phi: scalar or None
264     :param L_g: gas adsorption as function of gas pressure
265     :type L_g: `MaterialPropertyWithDifferential`
266     :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=p_fg-p_fw
267     :type S_fg: `MaterialPropertyWithDifferential`
268     :param perm_f_0: permeability in the x_0 direction in the fractured media
269     :type perm_f_0: scalar
270     :param perm_f_1: permeability in the x_1 direction in the fractured media
271     :type perm_f_1: scalar
272     :param perm_f_2: permeability in the x_2 direction in the fractured media
273     :type perm_f_2: scalar
274     :param k_w: relative permeability of water as function of water saturation
275     :type k_w: `MaterialProperty`
276     :param k_g: relative permeability of gas as function of gas saturation
277     :type k_g: `MaterialProperty`
278     :param mu_w: viscosity of water as function of water pressure
279     :type mu_w: `MaterialProperty`
280     :param mu_g: viscosity of gas as function of gas pressure
281     :type mu_g: `MaterialProperty`
282     :param rho_w: density of water as function of water pressure
283     :type rho_w: `MaterialPropertyWithDifferential`
284     :param rho_g: density of gas as function of gas pressure
285     :type rho_g: `MaterialPropertyWithDifferential`
286     :param wells : list of wells
287     :type wells: list of `Well`
288     """
289    
290     DualPorosity.__init__(self, domain,
291     phi_f=phi_f, phi_m=None, phi=phi,
292     S_fg=S_fg, S_mg=None,
293     perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
294     perm_m_0=None , perm_m_1=None, perm_m_2=None,
295     k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
296     rho_w =rho_w, rho_g=rho_g,
297     wells=wells)
298     self.L_g=L_g
299     self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
300    
301 gross 3458
302 gross 3478 def getPDEOptions(self):
303     """
304     returns the `SolverOptions` of the underlying PDE
305     """
306     return self.__pde.getSolverOptions()
307    
308     def solvePDE(self):
309    
310     p_fw=self.u[0]
311     p_fg=self.u[1]
312     p_fc=p_fg-p_fw
313     C_mg=self.u[3]
314    
315    
316     p_fw_old=self.u_old[0]
317     p_fg_old=self.u_old[1]
318     C_mg_old=self.u_old[3]
319 gross 3458
320    
321 gross 3478 phi_f =self.phi_f.getValue(p_fg)
322     dphi_fdp=self.phi_f.getValueDifferential(p_fg)
323    
324     S_fg= self.S_fg.getValue(p_fc)
325     dS_fgdp= self.S_fg.getValueDifferential(p_fc)
326     S_fw= 1-S_fg
327    
328     rho_fw = self.rho_w.getValue(p_fw)
329     drho_fwdp = self.rho_w.getValueDifferential(p_fw)
330     rho_fg = self.rho_g.getValue(p_fg)
331     drho_fgdp = self.rho_g.getValueDifferential(p_fg)
332    
333     L_g = self.getValue(p_fg)
334     dL_gdp_fg = self.rho_w.getValueDifferential(p_fg)
335    
336     A_fw = rho_fw * k_fw/mu_fw
337     A_fg = rho_fg * k_fg/mu_fg
338    
339     m = whereNegative(L_g - C_mg)
340     B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )
341    
342    
343     E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )
344     E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )
345     E_fgw= - phi_f * rho_fg * dS_fgdp
346     E_fgg= S_fg * rho_fg * dphi_fdp + phi_f * rho_fg * dS_fgdp + phi_f * S_fg * drho_fgdp
347 gross 3458
348    
349 gross 3478
350     Q_w=0.
351     Q_g=0.
352     R_w=0.
353     R_g=0.
354    
355     for I in self.wells:
356     chi_I= I.getGeometry()
357     Pi_I = I.getProductivityIndex()
358     P_I = I.getBHP()
359     Q_I = I.getFlowRate()
360    
361     if self.getPhase == Well.WATER:
362     Q_w = Q_w + Q_I * chi_I
363     Q_g = Q_g + Pi_I * P_I * chi_I
364     R_g = R_g + Pi_I * chi_I
365     else:
366     Q_g = Q_g + Q_I * chi_I
367     Q_w = Q_w + Pi_I * P_I * chi_I
368     R_w = R_w + Pi_I * chi_I
369    
370     F_fw_Y = A_fw * Q_w
371     F_fg_Y = A_fg * Q_g
372     D_fgg = - A_fg * R_g
373     D_fww = - A_fw * R_w
374    
375     if self.domain.getDim() == 3:
376     F_fw_X = A_fw * self.g * rho_fw * self.perm_f[2] * kronecker(self.domain)[2]
377     F_fg_X = A_fg * self.g * rho_fg * self.perm_f[2] * kronecker(self.domain)[2]
378     else:
379     F_fw_X = 0 * kronecker(self.domain)[1]
380     F_fg_X = 0 * kronecker(self.domain)[1]
381    
382     F_mg = B * ( L_g - dL_gdp_fg * p_fg )
383 gross 3458
384 gross 3478
385     D=self.__pde.getNewCoefficient("D")
386     A=self.__pde.getNewCoefficient("A")
387     Y=self.__pde.getNewCoefficient("Y")
388     X=self.__pde.getNewCoefficient("X")
389    
390     dtXI = dt*self.XI
391     dtcXI = dt*(1-self.XI)
392    
393     D[0,0]=E_fww + dtXI * D_fww
394     D[0,1]=E_fwg
395     D[1,0]=E_fgw
396     D[1,1]=E_fgg + dtXI * D_fgg
397     D[1,2]=rho_fg
398     D[2,1]= dtXI * B * dL_gdp_fg
399     D[2,2]= 1 + dtXI * B
400    
401     H_fw = dtcXI * A_fw * grad(p_fw_old)
402     H_fg = dtcXI * A_fg * grad(p_fg_old)
403    
404     for i in range(self.domain.getDim()):
405     A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
406     A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
407    
408     X[0,i]= dt * F_fw_X - H_fw * self.perm_f[i]
409     X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]
410    
411     Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old + dt * F_fw_Y - dtcXI * D_fww * p_fw_old
412     Y[1] = E_fgw * p_fw_old + E_fgg * p_fg_old + rho_fg * C_mg_old + dt * F_fg_Y - dtcXI * D_fgg * p_fg_old
413     Y[2] = C_mg_old + dt * F_mg_Y - dtcXI * ( dL_gdp_fg * p_fg_old + C_mg_old)
414    
415     self.__pde.setValue(A=A, D=D, X=X, Y=Y)
416    
417     u2 = self.__pde.getSolution()
418     return u2

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26