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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3487 - (show annotations)
Mon Mar 28 03:44:22 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 20645 byte(s)
some compile fixes
1 #######################################################
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.escript.linearPDEs import LinearPDE
24 from esys.escript import unitsSI as U
25 from esys.escript import sqrt, log, whereNegative, sup, inf
26 import math
27
28 class MaterialProperty(object):
29 """
30 generic class for material properties depending on one (or several) status
31 variable(s)
32 """
33 def __init__(self, *args, **kwargs):
34 """
35 set up
36 """
37 pass
38 def __call__(self, *args, **kwargs):
39 """
40 return value of material property for given status
41 """
42 return self.getValue( *args, **kwargs)
43
44
45 def getValue(self, *args, **kwargs):
46 """
47 return value of material property for given status
48
49 :remark: needs to be overwritten
50 """
51 raise NotImplementedError
52
53 class MaterialPropertyWithDifferential(MaterialProperty):
54 """
55 generic class for material properties depending on one (or several) status
56 variable(s) where the derivative of the property with respect to the
57 status variables is available
58 """
59 def getValueDifferential(self, *args, **kwargs):
60 """
61 returns the value of the derivative of material property for given status variable
62
63 :remark: needs to be overwritten
64 """
65 raise NotImplementedError
66
67 class Porosity(MaterialPropertyWithDifferential):
68 """
69 defines porosity phi as function of pressure
70
71 phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
72
73 phi_p_ref is claculted from the initial porosity phi_0 at pressure p_0
74 """
75 def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
76 """
77 """
78 self.phi_p_ref=1 # will be overwritten later
79 self.p_ref=p_ref
80 self.C=C
81 # reset phi_p_ref to get phi(p_0)=phi_0
82 self.phi_p_ref=phi_0/self.getValue(p_0)
83
84 def getValue(self, p):
85 """
86 returns the porosity for given pressure p
87 """
88 X= self.C * ( p - self.p_ref )
89 return self.phi_p_ref * (1. + X * (1. + X/2 ) )
90
91 def getValueDifferential(self, p):
92 """
93 returns the porosity for given pressure p
94 """
95 X= self.C * ( p - self.p_ref )
96 return self.phi_p_ref * self.C * (1. + X )
97
98
99 class WaterDensity(MaterialPropertyWithDifferential):
100 """
101 set water density as
102
103 rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )
104
105 with rho_ref =rho_s/B_ref * gravity
106 """
107 def __init__(self, B_ref=1., p_ref = 1.*U.atm, gravity=1., C = 0./U.bar, rho_s= 998.2 * U.kg/U.m**3):
108 self.rho_ref = rho_s/B_ref * gravity
109 self.C=C
110 self.p_ref=p_ref
111
112 def getValue(self, p):
113 """
114 returns the density for given pressure p
115 """
116 X= self.C * ( p - self.p_ref )
117 return self.rho_ref * (1+ X * (1+ X/2) )
118
119 def getValueDifferential(self, p):
120 """
121 """
122 X= self.C * ( p - self.p_ref )
123 return self.rho_ref * self.C * (1+ X)
124
125 class WaterViscosity(MaterialProperty):
126 """
127 defines viscosity of water
128
129 mu=mu_ref /(1 + X + X*X/2)
130
131 with X=-C*(p-p_ref)
132 """
133
134 def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
135 """
136 set up
137 """
138 self.mu_ref = mu_ref
139 self.p_ref = p_ref
140 self.C=C
141 def getValue(self, p):
142 """
143 returns the viscosity for given pressure p
144 """
145 X= -self.C * ( p - self.p_ref )
146 return self.mu_ref/(1+ X * (1+ X/2) )
147
148 class GasDensity(MaterialPropertyWithDifferential):
149 """
150 set water density as
151
152 rho = gravity * rho_air_s /B(p)
153
154 where B is given by an interpolation table
155 """
156 def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
157 self.rho_ref =rho_air_s * gravity
158 self.tab=InterpolationTable(x=p, y=B)
159
160 def getValue(self, p):
161 """
162 returns the density for given pressure p
163 """
164 return self.rho_ref/self.tab.getValue(p)
165
166 def getValueDifferential(self, p):
167 """
168 """
169 B = self.tab.getValue(p)
170 dBdp = self.tab.getValueDifferential(p)
171 return -self.rho_ref * dBdp/(B * B)
172
173 class InterpolationTable(MaterialPropertyWithDifferential):
174 """
175 a simple 1D interpolation table for escript Data object with non-equidistant nodes
176 """
177 def __init__(self,x=[], y=[], obeyBounds=True):
178 """
179 set up interpolation table. obeyBounds is set an exception is thrown if
180 the interpolation argument is below min(x) or larger than max(x). Otherwise
181 the value for x is set to y[0] for
182 """
183 MaterialPropertyWithDifferential.__init__(self)
184 if len(x) != len(y):
185 raise ValueError,"length of interpolation nodes and value lists need to be identical."
186 if len(x) < 1 :
187 raise ValueError,"length of interpolation nodes a list needs to at least one."
188
189 x_ref=x[0]
190 for i in range(1,len(x)):
191 if x_ref >= x[i]:
192 raise ValueError,"interpolation nodes need to be given in increasing order."
193 x_ref=x[i]
194 self.__x = x
195 self.__y = y
196 self.__obeyBounds = obeyBounds
197
198 def getValue(self, x):
199 """
200 returns the interpolated values of x
201 """
202 X=self.__x
203 Y=self.__y
204
205 x0=X[0]
206 m0=whereNegative(x-x0)
207 if self.__obeyBounds:
208 if sup(m0) > 0:
209 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
210 out=self.__x[0]
211 for i in range(1,len(X)):
212 z=(Y[i]-Y[i-1])/(X[i]-X[i-1]) * (x-X[i-1]) + Y[i-1]
213 out = out * m0 + z * (1-m0)
214 m0=whereNegative(x-X[i])
215
216 if self.__obeyBounds:
217 if inf(m0) < 1:
218 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
219 else:
220 out = out * m0 + V[-1] * (1-m0)
221 return out
222
223 def getValueDifferential(self, x):
224 X=self.__x
225 Y=self.__y
226
227 x0=X[0]
228 m0=whereNegative(x-x0)
229 if self.__obeyBounds:
230 if sup(m0) > 0:
231 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
232 out=0.
233 for i in range(1,len(X)):
234 z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
235 out = out * m0 + z * (1-m0)
236 m0=whereNegative(x-X[i])
237
238 if self.__obeyBounds:
239 if inf(m0) < 1:
240 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
241 else:
242 out = out * m0
243 return out
244
245
246 class Well(object):
247 """
248 generic well
249
250 :var WATER: phase identifier for water well
251 :var GAS: phase identifier for gas well
252 """
253 WATER="Water"
254 GAS="Gas"
255 def __init__(self, Q=0., schedule=[0.], phase="Water", *args, **kwargs):
256 """
257 set up well
258 """
259 self.__schedule=schedule
260 self.__phase=phase
261 self.__Q = Q
262
263 def getGeometry(self):
264 """
265 returns a function representing the geometry of the well.
266 typically a Gaussian profile around the location of the well.
267
268 :note: needs to be overwritten
269 """
270 raise NotImplementedError
271
272 def getProductivityIndex(self):
273 """
274 returns the productivity index of the well.
275 typically a Gaussian profile around the location of the well.
276
277 :note: needs to be overwritten
278 """
279 raise NotImplementedError
280
281 def getFlowRate(self, t):
282 """
283 returns the flow rate
284 """
285 if self.isOpen(t):
286 return self.__Q
287 else:
288 return 0.
289
290 def getBHP(self):
291 """
292 return bottom-hole pressure
293
294 :note: needs to be overwritten
295 """
296 raise NotImplementedError
297
298 def getPhase(self):
299 """
300 returns the pahse the well works on
301 """
302 return self.__phase
303
304 def isOpen(self, t):
305 """
306 returns True is the well is open at time t
307 """
308 out = False
309 t0=min(t, self.__schedule[0])
310 i=0
311 while t < t0:
312 if out:
313 out=False
314 else:
315 out=True
316 i+=1
317 if i < len(self.__schedule[i]):
318 t0=self.__schedule[i]
319 else:
320 t0=t
321 return out
322
323 class VerticalPeacemanWell(Well):
324 """
325 defines a well using Peaceman's formula
326 """
327 def __init__(self,name, schedule = [0.], BHP=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
328 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
329 """
330 set up well
331
332 :param BHP: ottom-hole pressure
333 :param Q: flow rate
334 :param r: well radius
335 :param X: location
336 :param D: dimension of well block
337 :param perm: permeabilities
338 :param s: skin factor
339 """
340 Well.__init__(self, Q=Q, schedule=schedule, phase=phase,)
341 r_el=0.28 * ( (perm[1]/perm[0])**0.5 * D[0]**2 + (perm[0]/perm[1])**0.5 * D[1]**2 ) \
342 / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
343 self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
344 self.__r = r
345
346 self.__BHP = BHP
347 self.__D=D
348 self.r_el=r_el
349 self.X0=X0
350
351 def getGeometry(self):
352 """
353 returns a function representing the geometry of the well.
354 typically a Gaussian profile around the location of the well.
355
356 :note: needs to be overwritten
357 """
358 raise NotImplementedError
359
360 def getProductivityIndex(self,t):
361 """
362 returns the productivity index of the well.
363 typically a Gaussian profile around the location of the well.
364 """
365 return self.__PI
366
367 def getBHP(self):
368 """
369 return bottom-hole pressure
370 """
371 return self.__BHP
372
373 class DualPorosity(object):
374 """
375 generic dual porosity model
376 """
377 def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
378 S_fg=None, S_mg=None,
379 perm_f_0=None, perm_f_1=None, perm_f_2=None,
380 perm_m_0=None, perm_m_1=None, perm_m_2=None,
381 k_w =None, k_g=None, mu_w =None, mu_g =None,
382 rho_w =None, rho_g=None,
383 wells=[]):
384 """
385 set up
386
387 :param domain: domain
388 :type domain: `esys.escript.Domain`
389 :param phi_f: porosity of the fractured rock as function of the gas pressure
390 :type phi_f: `MaterialPropertyWithDifferential`
391 :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
392 :type phi_m: `MaterialPropertyWithDifferential` or None
393 :param phi: total porosity if None phi=1.
394 :type phi: scalar or None
395 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
396 :type S_fg: `MaterialPropertyWithDifferential`
397 :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
398 :type S_mg: `MaterialPropertyWithDifferential`
399 :param perm_f_0: permeability in the x_0 direction in the fractured media
400 :type perm_f_0: scalar
401 :param perm_f_1: permeability in the x_1 direction in the fractured media
402 :type perm_f_1: scalar
403 :param perm_f_2: permeability in the x_2 direction in the fractured media
404 :type perm_f_2: scalar
405 :param perm_m_0: permeability in the x_0 direction in the coal matrix
406 :type perm_m_0: scalar
407 :param perm_m_1: permeability in the x_1 direction in the coal matrix
408 :type perm_m_1: scalar
409 :param perm_m_2: permeability in the x_2 direction in the coal matrix
410 :type perm_m_2: scalar
411 :param k_w: relative permeability of water as function of water saturation
412 :type k_w: `MaterialProperty`
413 :param k_g: relative permeability of gas as function of gas saturation
414 :type k_g: `MaterialProperty`
415 :param mu_w: viscosity of water as function of water pressure
416 :type mu_w: `MaterialProperty`
417 :param mu_g: viscosity of gas as function of gas pressure
418 :type mu_g: `MaterialProperty`
419 :param rho_w: density of water as function of water pressure
420 :type rho_w: `MaterialPropertyWithDifferential`
421 :param rho_g: density of gas as function of gas pressure
422 :type rho_g: `MaterialPropertyWithDifferential`
423 :param wells : list of wells
424 :type wells: list of `Well`
425 """
426
427 self.domain = domain
428 self.phi_f = phi_f
429 self.phi_m = phi_m
430 self.phi = phi
431 self.S_fg = S_fg
432 self.S_mg = S_mg
433 self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
434 self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
435 self.k_w = k_w
436 self.k_g = k_g
437 self.mu_w = mu_w
438 self.mu_g = mu_g
439 self.rho_w = rho_w
440 self.rho_g = rho_g
441 self.wells=wells
442 self.t =0
443
444 self.__iter_max=1
445 self.__rtol=1.e-4
446 self.__verbose=False
447 def setIterationControl(self, iter_max=None, rtol=None, verbose=None):
448 """
449 sets parameters to control iteration process
450 """
451 if iter_max !=None: self.__iter_max=iter_max
452 if rtol !=None: self.__rtol = rtol
453 if verbose !=None: self.__verbose=verbose
454
455 def update(self, dt):
456 self.u, self.u_old = self.u.copy(), self.u
457 n=0
458 rerr=1.
459 while n < self.__iter_max and rerr > self.__rtol:
460 u=self.solvePDE(dt)
461 norm_u=Lsup(u)
462 norm_e=Lsup(u-self.u)
463 if norm_u > 0:
464 rerr=norm_e/norm_u
465 else:
466 rerr=norm_e
467 if self.__verbose: print "iteration step %d: relative change = %e"%(n, rerr)
468 n+=1
469 self.u=u
470 print "iteration completed."
471 self.t+=dt
472
473
474 class PorosityOneHalfModel(DualPorosity):
475 """
476 Model for gas and water in double prosity model tracking water and gas
477 pressure in the fractured rock and gas concentration in the coal matrix.
478 This corresponds to the coal bed methan model in the ECLIPSE code.
479 """
480
481 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
482 perm_f_0=None, perm_f_1=None, perm_f_2=None,
483 k_w =None, k_g=None, mu_w =None, mu_g =None,
484 rho_w =None, rho_g=None,
485 wells=[]):
486 """
487 set up
488
489 :param domain: domain
490 :type domain: `esys.escript.Domain`
491 :param phi_f: porosity of the fractured rock as function of the gas pressure
492 :type phi_f: `MaterialPropertyWithDifferential`
493 :param phi: total porosity if None phi=1.
494 :type phi: scalar or None
495 :param L_g: gas adsorption as function of gas pressure
496 :type L_g: `MaterialPropertyWithDifferential`
497 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
498 :type S_fg: `MaterialPropertyWithDifferential`
499 :param perm_f_0: permeability in the x_0 direction in the fractured media
500 :type perm_f_0: scalar
501 :param perm_f_1: permeability in the x_1 direction in the fractured media
502 :type perm_f_1: scalar
503 :param perm_f_2: permeability in the x_2 direction in the fractured media
504 :type perm_f_2: scalar
505 :param k_w: relative permeability of water as function of water saturation
506 :type k_w: `MaterialProperty`
507 :param k_g: relative permeability of gas as function of gas saturation
508 :type k_g: `MaterialProperty`
509 :param mu_w: viscosity of water as function of water pressure
510 :type mu_w: `MaterialProperty`
511 :param mu_g: viscosity of gas as function of gas pressure
512 :type mu_g: `MaterialProperty`
513 :param rho_w: density of water as function of water pressure
514 :type rho_w: `MaterialPropertyWithDifferential`
515 :param rho_g: density of gas as function of gas pressure
516 :type rho_g: `MaterialPropertyWithDifferential`
517 :param wells : list of wells
518 :type wells: list of `Well`
519 """
520
521 DualPorosity.__init__(self, domain,
522 phi_f=phi_f, phi_m=None, phi=phi,
523 S_fg=None, S_mg=None,
524 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
525 perm_m_0=None , perm_m_1=None, perm_m_2=None,
526 k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
527 rho_w =rho_w, rho_g=rho_g,
528 wells=wells)
529 self.L_g=L_g
530 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
531
532
533 def getPDEOptions(self):
534 """
535 returns the `SolverOptions` of the underlying PDE
536 """
537 return self.__pde.getSolverOptions()
538
539 def getState(self):
540 return self.u[0], self.u[1], self.u[2]
541
542 def getOldState(self):
543 return self.u_old[0], self.u_old[1], self.u_old[2]
544
545 def setInitialState(self, p=1.*U.atm, S_fg=0, C_mg=None):
546 """
547 sets the initial state
548
549 :param p: pressure
550 :param S_fg: gas saturation in fractured rock
551 :param C_mg: gas concentration in coal matrix. if not given it is calculated
552 using the gas adsorption curve.
553 """
554 self.u=self.__pde.createSolution()
555 self.u[0]=p
556 self.u[1]=S_fg
557 if C_mg == None:
558 self.u[2]= self.L_g(p)
559 else:
560 self.u[2]=C_mg
561
562 def solvePDE(self, dt):
563
564 p_f, S_fg, C_mg=self.getState()
565 p_f_old, S_fg_old, C_mg_old=self.getOldState()
566
567 S_fw=1-S_fg
568
569 k_fw=self.k_w(S_fw)
570 k_fg=self.k_g(S_fg)
571 mu_fw=self.mu_w(p_f)
572 mu_fg=self.mu_g(p_f)
573
574
575 phi_f =self.phi_f.getValue(S_fg)
576 dphi_fdp=self.phi_f.getValueDifferential(S_fg)
577
578 S_fw= 1-S_fg
579
580 rho_fw = self.rho_w.getValue(p_f)
581 drho_fwdp = self.rho_w.getValueDifferential( p_f)
582 rho_fg = self.rho_g.getValue(p_f)
583 drho_fgdp = self.rho_g.getValueDifferential(p_f)
584
585 L_g = self.L_g.getValue(p_f)
586 dL_gdp = self.rho_w.getValueDifferential(p_f)
587
588
589 A_fw = rho_fw * k_fw/mu_fw
590 A_fg = rho_fg * k_fg/mu_fg
591
592 m = whereNegative(L_g - C_mg)
593 B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )
594
595
596 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
597 E_fps= - phi_f * rho_fw
598 E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f * drho_fgdp )
599 E_fss= phi_f * rho_fg
600
601
602
603 Q_w=0.
604 Q_g=0.
605 R_w=0.
606 R_g=0.
607
608 for I in self.wells:
609 chi_I= I.getGeometry()
610 Pi_I = I.getProductivityIndex()
611 P_I = I.getBHP()
612 Q_I = I.getFlowRate()
613
614 if self.getPhase == Well.WATER:
615 Q_w = Q_w + Q_I * chi_I
616 Q_g = Q_g + Pi_I * P_I * chi_I
617 R_g = R_g + Pi_I * chi_I
618 else:
619 Q_g = Q_g + Q_I * chi_I
620 Q_w = Q_w + Pi_I * P_I * chi_I
621 R_w = R_w + Pi_I * chi_I
622
623 F_fp_Y = A_fw * Q_w
624 F_fs_Y = A_fg * Q_g
625 D_fss = - A_fg * R_g
626 D_fpp = - A_fw * R_w
627
628 if self.domain.getDim() == 3:
629 F_fp_X = A_fw * self.g * rho_fw * self.perm_f[2] * kronecker(self.domain)[2]
630 F_fs_X = A_fg * self.g * rho_fg * self.perm_f[2] * kronecker(self.domain)[2]
631 else:
632 F_fp_X = 0 * kronecker(self.domain)[1]
633 F_fs_X = 0 * kronecker(self.domain)[1]
634
635 F_mg = B * ( L_g - dL_gdp * S_fg )
636
637
638 D=self.__pde.getNewCoefficient("D")
639 A=self.__pde.getNewCoefficient("A")
640 Y=self.__pde.getNewCoefficient("Y")
641 X=self.__pde.getNewCoefficient("X")
642
643 dtXI = dt*self.XI
644 dtcXI = dt*(1-self.XI)
645
646 D[0,0]=E_fpp + dtXI * D_fpp
647 D[0,1]=E_fps
648 D[1,0]=E_fsp
649 D[1,1]=E_fss + dtXI * D_fss
650 D[1,2]=rho_fg
651 D[2,1]= dtXI * B * dL_gdp
652 D[2,2]= 1 + dtXI * B
653
654 H_fw = dtcXI * A_fw * grad( p_f_old)
655 H_fg = dtcXI * A_fg * grad(S_fg_old)
656
657 for i in range(self.domain.getDim()):
658 A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
659 A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
660
661 X[0,i]= dt * F_fp_X - H_fw * self.perm_f[i]
662 X[1,i]= dt * F_fs_X - H_fg * self.perm_f[i]
663
664 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old + dt * F_fp_Y - dtcXI * D_fpp * p_f_old
665 Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + rho_fg * C_mg_old + dt * F_fs_Y - dtcXI * D_fss * S_fg_old
666 Y[2] = C_mg_old + dt * F_mg_Y - dtcXI * ( dL_gdp * S_fg_old + C_mg_old)
667
668 self.__pde.setValue(A=A, D=D, X=X, Y=Y)
669
670 u2 = self.__pde.getSolution()
671 return u2

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26