/[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 3515 - (show annotations)
Thu May 19 08:20:57 2011 UTC (8 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 25635 byte(s)
some first work for DiracFunctions
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.pdetools import Locator
26 from esys.escript import *
27 from esys.weipa import saveVTK
28 import math
29
30 USE_NODAL_WELL = False or True
31
32 class MaterialProperty(object):
33 """
34 generic class for material properties depending on one (or several) status
35 variable(s)
36 """
37 def __init__(self, *args, **kwargs):
38 """
39 set up
40 """
41 pass
42 def __call__(self, *args, **kwargs):
43 """
44 return value of material property for given status
45 """
46 return self.getValue( *args, **kwargs)
47
48
49 def getValue(self, *args, **kwargs):
50 """
51 return value of material property for given status
52
53 :remark: needs to be overwritten
54 """
55 raise NotImplementedError
56
57 class MaterialPropertyWithDifferential(MaterialProperty):
58 """
59 generic class for material properties depending on one (or several) status
60 variable(s) where the derivative of the property with respect to the
61 status variables is available
62 """
63 def getValueDifferential(self, *args, **kwargs):
64 """
65 returns the value of the derivative of material property for given status variable
66
67 :remark: needs to be overwritten
68 """
69 raise NotImplementedError
70
71 class Porosity(MaterialPropertyWithDifferential):
72 """
73 defines porosity phi as function of pressure
74
75 phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
76
77 phi_p_ref is claculted from the initial porosity phi_0 at pressure p_0
78 """
79 def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
80 """
81 """
82 self.phi_p_ref=1 # will be overwritten later
83 self.p_ref=p_ref
84 self.C=C
85 # reset phi_p_ref to get phi(p_0)=phi_0
86 self.phi_p_ref=phi_0/self.getValue(p_0)
87
88 def getValue(self, p):
89 """
90 returns the porosity for given pressure p
91 """
92 X= self.C * ( p - self.p_ref )
93 return self.phi_p_ref * (1. + X * (1. + X/2 ) )
94
95 def getValueDifferential(self, p):
96 """
97 returns the porosity for given pressure p
98 """
99 X= self.C * ( p - self.p_ref )
100 return self.phi_p_ref * self.C * (1. + X )
101
102
103 class WaterDensity(MaterialPropertyWithDifferential):
104 """
105 set water density as
106
107 rho = rho_surf (1 + X + X*X/2) with X= C * ( p - p_ref )
108
109 with rho_surf =rho_s/B_ref * gravity
110 """
111 def __init__(self, B_ref=1., p_ref = 1.*U.atm, gravity=1., C = 0./U.bar, rho_s= 999.014 * U.kg/U.m**3):
112 self.rho_surf = rho_s * gravity
113 self.__rho_0 = self.rho_surf/B_ref
114 self.C=C
115 self.p_ref=p_ref
116
117 def getValue(self, p):
118 """
119 returns the density for given pressure p
120 """
121 X= self.C * ( p - self.p_ref )
122 return self.__rho_0 * (1+ X * (1+ X/2) )
123
124 def getValueDifferential(self, p):
125 """
126 """
127 X= self.C * ( p - self.p_ref )
128 return self.__rho_0 * self.C * (1+ X)
129
130 def getFormationFactor(self, p):
131 return self.rho_surf/self.getValue(p)
132
133
134
135 class WaterViscosity(MaterialProperty):
136 """
137 defines viscosity of water
138
139 mu=mu_ref /(1 + X + X*X/2)
140
141 with X=-C*(p-p_ref)
142 """
143
144 def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
145 """
146 set up
147 """
148 self.mu_ref = mu_ref
149 self.p_ref = p_ref
150 self.C=C
151 def getValue(self, p):
152 """
153 returns the viscosity for given pressure p
154 """
155 X= -self.C * ( p - self.p_ref )
156 return self.mu_ref/(1+ X * (1+ X/2) )
157
158 class GasDensity(MaterialPropertyWithDifferential):
159 """
160 set water density as
161
162 rho = gravity * rho_air_s /B(p)
163
164 where B is given by an interpolation table
165 """
166 def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
167 self.rho_surf =rho_air_s * gravity
168 self.tab=InterpolationTable(x=p, y=B)
169
170 def getValue(self, p):
171 """
172 returns the density for given pressure p
173 """
174 return self.rho_surf/self.getFormationFactor(p)
175
176 def getValueDifferential(self, p):
177 """
178 """
179 B = self.getFormationFactor(p)
180 dBdp = self.getFormationFactorDifferential(p)
181 return -self.rho_surf * dBdp/(B * B)
182
183 def getFormationFactor(self, p):
184 return self.tab.getValue(p)
185
186 def getFormationFactorDifferential(self, p):
187 return self.tab.getValueDifferential(p)
188
189
190 class InterpolationTable(MaterialPropertyWithDifferential):
191 """
192 a simple 1D interpolation table for escript Data object with non-equidistant nodes
193 """
194 def __init__(self,x=[], y=[], obeyBounds=True):
195 """
196 set up interpolation table. obeyBounds is set an exception is thrown if
197 the interpolation argument is below min(x) or larger than max(x). Otherwise
198 the value for x is set to y[0] for
199 """
200 MaterialPropertyWithDifferential.__init__(self)
201 if len(x) != len(y):
202 raise ValueError,"length of interpolation nodes and value lists need to be identical."
203 if len(x) < 1 :
204 raise ValueError,"length of interpolation nodes a list needs to at least one."
205
206 x_ref=x[0]
207 for i in range(1,len(x)):
208 if x_ref >= x[i]:
209 raise ValueError,"interpolation nodes need to be given in increasing order."
210 x_ref=x[i]
211 self.__x = x
212 self.__y = y
213 self.__obeyBounds = obeyBounds
214
215 def getValue(self, x):
216 """
217 returns the interpolated values of x
218 """
219 X=self.__x
220 Y=self.__y
221
222 x0=X[0]
223 m0=whereNegative(x-x0)
224 if self.__obeyBounds:
225 if sup(m0) > 0:
226 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
227 out=self.__x[0]
228 for i in range(1,len(X)):
229 z=(Y[i]-Y[i-1])/(X[i]-X[i-1]) * (x-X[i-1]) + Y[i-1]
230 out = out * m0 + z * (1-m0)
231 m0=whereNonPositive(x-X[i])
232
233 if self.__obeyBounds:
234 if inf(m0) < 1 :
235 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
236 else:
237 out = out * m0 + X[-1] * (1-m0)
238 return out
239
240 def getValueDifferential(self, x):
241 X=self.__x
242 Y=self.__y
243
244 x0=X[0]
245 m0=whereNegative(x-x0)
246 if self.__obeyBounds:
247 if sup(m0) > 0:
248 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
249 out=0.
250 for i in range(1,len(X)):
251 z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
252 out = out * m0 + z * (1-m0)
253 m0=whereNonPositive(x-X[i])
254
255 if self.__obeyBounds:
256 if inf(m0) < 1:
257 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
258 else:
259 out = out * m0
260 return out
261
262
263 class Well(object):
264 """
265 generic well
266
267 :var WATER: phase identifier for water well
268 :var GAS: phase identifier for gas well
269 """
270 WATER="Water"
271 GAS="Gas"
272 def __init__(self, name, domain, Q=[0.], schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
273 """
274 set up well
275 """
276 if not len(schedule) == len(Q):
277 raise ValueError,"length of schedule and Q must match."
278 self.__schedule=schedule
279 self.__Q = Q
280 self.__phase=phase
281 self.__BHP_limit=BHP_limit
282 self.name=name
283 self.domain=domain
284
285 def getGeometry(self):
286 """
287 returns a function representing the geometry of the well.
288 typically a Gaussian profile around the location of the well.
289
290 :note: needs to be overwritten
291 """
292 raise NotImplementedError
293
294 def getProductivityIndex(self):
295 """
296 returns the productivity index of the well.
297 typically a Gaussian profile around the location of the well.
298
299 :note: needs to be overwritten
300 """
301 raise NotImplementedError
302
303 def getFlowRate(self,t):
304 """
305 returns the flow rate
306 """
307 out=0.
308 for i in range(len(self.__schedule)):
309 if t <= self.__schedule[i]:
310 out=self.__Q[i]
311 break
312 return out
313
314 def getBHPLimit(self):
315 """
316 return bottom-hole pressure limit
317
318 :note: needs to be overwritten
319 """
320 return self.__BHP_limit
321
322 def getBHP(self):
323 """
324 return bottom-hole pressure
325
326 :note: needs to be overwritten
327 """
328 return self.__BHP
329
330 def setBHP(self, BHP):
331 """
332 sets bottom-hole pressure
333 """
334 self.__BHP= BHP
335
336 def getPhase(self):
337 """
338 returns the pahse the well works on
339 """
340 return self.__phase
341
342 class VerticalPeacemanWell(Well):
343 """
344 defines a well using Peaceman's formula
345 """
346 def __init__(self,name, domain, schedule = [0.], BHP_limit=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
347 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5):
348 # reset_factor=0.1):
349 """
350 set up well
351
352 :param BHP: ottom-hole pressure
353 :param Q: flow rate
354 :param r: well radius
355 :param X: location
356 :param D: dimension of well block
357 :param perm: permeabilities
358 :param s: skin factor
359 """
360 Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)
361 r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 + sqrt(perm[0]/perm[1]) * D[1]**2 )\
362 / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
363 self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
364 self.__r = r
365
366 self.__D=D
367 self.r_el=r_el
368 self.locator=Locator(ContinuousFunction(self.domain),X0[:self.domain.getDim()])
369 self.X0=self.locator.getX()
370 # this needs to be done differently
371 self.chi=whereZero(length(self.domain.getX()-self.X0))
372 self.chi*=1./integrate(self.chi)
373 if self.domain.getDim() <3:
374 self.chi*=1/self.__D[2]
375
376 def getLocation(self):
377 return self.X0
378
379 def getProductivityIndex(self):
380 """
381 returns the productivity index of the well.
382 """
383 return self.__PI
384
385
386 class DualPorosity(object):
387 """
388 generic dual porosity model
389 """
390 def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
391 S_fg=None, S_mg=None,
392 perm_f_0=None, perm_f_1=None, perm_f_2=None,
393 perm_m_0=None, perm_m_1=None, perm_m_2=None,
394 k_w =None, k_g=None, mu_w =None, mu_g =None,
395 rho_w =None, rho_g=None,
396 wells=[], g=9.81*U.m/U.sec**2):
397 """
398 set up
399
400 :param domain: domain
401 :type domain: `esys.escript.Domain`
402 :param phi_f: porosity of the fractured rock as function of the gas pressure
403 :type phi_f: `MaterialPropertyWithDifferential`
404 :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
405 :type phi_m: `MaterialPropertyWithDifferential` or None
406 :param phi: total porosity if None phi=1.
407 :type phi: scalar or None
408 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
409 :type S_fg: `MaterialPropertyWithDifferential`
410 :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
411 :type S_mg: `MaterialPropertyWithDifferential`
412 :param perm_f_0: permeability in the x_0 direction in the fractured media
413 :type perm_f_0: scalar
414 :param perm_f_1: permeability in the x_1 direction in the fractured media
415 :type perm_f_1: scalar
416 :param perm_f_2: permeability in the x_2 direction in the fractured media
417 :type perm_f_2: scalar
418 :param perm_m_0: permeability in the x_0 direction in the coal matrix
419 :type perm_m_0: scalar
420 :param perm_m_1: permeability in the x_1 direction in the coal matrix
421 :type perm_m_1: scalar
422 :param perm_m_2: permeability in the x_2 direction in the coal matrix
423 :type perm_m_2: scalar
424 :param k_w: relative permeability of water as function of water saturation
425 :type k_w: `MaterialProperty`
426 :param k_g: relative permeability of gas as function of gas saturation
427 :type k_g: `MaterialProperty`
428 :param mu_w: viscosity of water as function of water pressure
429 :type mu_w: `MaterialProperty`
430 :param mu_g: viscosity of gas as function of gas pressure
431 :type mu_g: `MaterialProperty`
432 :param rho_w: density of water as function of water pressure
433 :type rho_w: `MaterialPropertyWithDifferential`
434 :param rho_g: density of gas as function of gas pressure
435 :type rho_g: `MaterialPropertyWithDifferential`
436 :param wells : list of wells
437 :type wells: list of `Well`
438 """
439
440 self.domain = domain
441 self.phi_f = phi_f
442 self.phi_m = phi_m
443 self.phi = phi
444 self.S_fg = S_fg
445 self.S_mg = S_mg
446 self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
447 self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
448 self.k_w = k_w
449 self.k_g = k_g
450 self.mu_w = mu_w
451 self.mu_g = mu_g
452 self.rho_w = rho_w
453 self.rho_g = rho_g
454 self.wells=wells
455 self.t =0
456 self.g=g
457
458 self.__iter_max=1
459 self.__rtol=1.e-4
460 self.verbose=False
461 self.XI=0.5
462 def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
463 """
464 sets parameters to control iteration process
465 """
466 if iter_max !=None: self.__iter_max=iter_max
467 if rtol !=None: self.__rtol = rtol
468 if verbose !=None: self.verbose=verbose
469 if xi !=None: self.XI=xi
470
471 def update(self, dt):
472 self.u, self.u_old =tuple([ v.copy() for v in self.u ] ), self.u
473 n=0
474 converged=False
475 while n < self.__iter_max and not converged:
476 u=self.solvePDE(dt)
477 if self.verbose: print "iteration step %d:"
478 converged=True
479 for i in range(len(u)):
480 if isinstance(u[i], Data):
481 norm_u=Lsup(u[i])
482 norm_e=Lsup(u[i]-self.u[i])
483 else:
484 norm_e=0.
485 norm_u=1.
486
487 if norm_u > 0:
488 rerr=norm_e/norm_u
489 else:
490 rerr=norm_e
491 if rerr>self.__rtol: converged=False
492 if self.verbose: print " comp %i: relative change = %e"%(i, rerr)
493 n+=1
494 self.u=u
495 print "iteration completed."
496 self.t+=dt
497
498
499 class PorosityOneHalfModel(DualPorosity):
500 """
501 Model for gas and water in double prosity model tracking water and gas
502 pressure in the fractured rock and gas concentration in the coal matrix.
503 This corresponds to the coal bed methan model in the ECLIPSE code.
504 """
505
506 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
507 perm_f_0=None, perm_f_1=None, perm_f_2=None,
508 k_w =None, k_g=None, mu_w =None, mu_g =None,
509 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
510 wells=[], g=9.81*U.m/U.sec**2):
511 """
512 set up
513
514 :param domain: domain
515 :type domain: `esys.escript.Domain`
516 :param phi_f: porosity of the fractured rock as function of the gas pressure
517 :type phi_f: `MaterialPropertyWithDifferential`
518 :param phi: total porosity if None phi=1.
519 :type phi: scalar or None
520 :param L_g: gas adsorption as function of gas pressure
521 :type L_g: `MaterialPropertyWithDifferential`
522 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
523 :type S_fg: `MaterialPropertyWithDifferential`
524 :param perm_f_0: permeability in the x_0 direction in the fractured media
525 :type perm_f_0: scalar
526 :param perm_f_1: permeability in the x_1 direction in the fractured media
527 :type perm_f_1: scalar
528 :param perm_f_2: permeability in the x_2 direction in the fractured media
529 :type perm_f_2: scalar
530 :param k_w: relative permeability of water as function of water saturation
531 :type k_w: `MaterialProperty`
532 :param k_g: relative permeability of gas as function of gas saturation
533 :type k_g: `MaterialProperty`
534 :param mu_w: viscosity of water as function of water pressure
535 :type mu_w: `MaterialProperty`
536 :param mu_g: viscosity of gas as function of gas pressure
537 :type mu_g: `MaterialProperty`
538 :param rho_w: density of water as function of water pressure
539 :type rho_w: `MaterialPropertyWithDifferential`
540 :param rho_g: density of gas as function of gas pressure
541 :type rho_g: `MaterialPropertyWithDifferential`
542 :param wells : list of wells
543 :type wells: list of `Well`
544 :param sigma: shape factor for gas matrix diffusion
545 :param A_mg: diffusion constant for gas adsorption
546 :param f_rg: gas re-adsorption factor
547 """
548
549 DualPorosity.__init__(self, domain,
550 phi_f=phi_f, phi_m=None, phi=phi,
551 S_fg=None, S_mg=None,
552 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
553 perm_m_0=None , perm_m_1=None, perm_m_2=None,
554 k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
555 rho_w =rho_w, rho_g=rho_g,
556 wells=wells, g=g)
557 self.L_g=L_g
558 self.sigma = sigma
559 self.A_mg = A_mg
560 self.f_rg = f_rg
561 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
562
563
564 def getPDEOptions(self):
565 """
566 returns the `SolverOptions` of the underlying PDE
567 """
568 return self.__pde.getSolverOptions()
569
570 def getState(self):
571 return self.u
572
573 def getOldState(self):
574 return self.u_old
575
576 def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0, c_mg=None):
577 """
578 sets the initial state
579
580 :param p: pressure
581 :param S_fg: gas saturation in fractured rock
582 :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
583 using the gas adsorption curve.
584 """
585 if self.domain.getDim() == 2:
586 p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
587 else:
588 z=self.u.getDomain().getX()[0]
589 z_top=sup(z)
590 z_bottom=inf(z)
591 p_init=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
592
593 S_fg_init=Scalar(S_fg, Solution(self.domain))
594
595 if c_mg == None:
596 c_mg_init=self.L_g(p_init)
597 else:
598 c_mg_init=Scalar(c_mg, Solution(self.domain))
599
600 wells_init={}
601 for I in self.wells:
602 BHP_limit_I=I.getBHPLimit()
603 p_f_I=I.locator(p_init)
604 q_I = I.getFlowRate(0)
605 if abs(q_I)>0:
606 raise ValueError,"initial flux must be zero."
607
608 if p_f_I < BHP_limit_I: # fix me!
609 raise ValueError,"initial pressure must be greater than BHP limit."
610 wells_init[I.name]={ "bhp" : p_f_I, "q_gas": 0., "q_water": 0.}
611
612 self.u=(p_init,S_fg_init, c_mg_init, wells_init)
613
614 def solvePDE(self, dt):
615
616 p_f, S_fg, c_mg, wells_state =self.getState()
617 p_f_old, S_fg_old, c_mg_old, wells_sate_old=self.getOldState()
618
619 S_fw=1-S_fg
620
621 if self.verbose:
622 print "p_f range = ",inf(p_f),sup(p_f)
623 print "S_fg range = ",inf(S_fg),sup(S_fg)
624 print "S_fw range = ",inf(S_fw),sup(S_fw)
625 print "c_mg range = ",inf(c_mg),sup(c_mg)
626 print "wells state =",wells_state
627
628 k_fw=self.k_w(S_fw)
629 if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
630
631
632 k_fg=self.k_g(S_fg)
633 if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
634
635 mu_fw=self.mu_w(p_f)
636 if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
637
638 mu_fg=self.mu_g(p_f)
639 if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
640
641
642 phi_f =self.phi_f.getValue(p_f)
643 dphi_fdp=self.phi_f.getValueDifferential(p_f)
644 if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
645
646 rho_fw = self.rho_w.getValue(p_f)
647 drho_fwdp = self.rho_w.getValueDifferential(p_f)
648 if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
649
650 rho_fg = self.rho_g.getValue(p_f)
651 rho_g_surf = self.rho_g.rho_surf
652 drho_fgdp = self.rho_g.getValueDifferential(p_f)
653 if self.verbose:
654 print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
655 print "rho_fg surf = ",rho_g_surf
656
657 L_g = self.L_g(p_f)
658 dL_gdp = self.L_g.getValueDifferential(p_f)
659 if self.verbose: print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
660
661 A_fw = rho_fw * k_fw/mu_fw
662 A_fg = rho_fg * k_fg/mu_fg
663
664 m = whereNegative(L_g - c_mg)
665 H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
666
667 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
668 E_fps= - phi_f * rho_fw
669 E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
670 E_fss= phi_f * rho_fg
671
672 F_fw=0.
673 F_fg=0.
674 D_fw=0.
675 D_fg=0.
676
677 for I in self.wells:
678 Pi_I = I.getProductivityIndex()
679 p_f_I=I.locator(p_f)
680 S_fg_I=I.locator(S_fg)
681 BHP_I=wells_state[I.name]["bhp"]
682 if self.verbose: print "well %s: BHP = %e"%(I.name, BHP_I)
683
684 D_fw = D_fw + Pi_I * I.chi
685 F_fw = F_fw + Pi_I * BHP_I *I.chi
686 D_fg = D_fg + Pi_I * I.chi
687 F_fg = F_fg + Pi_I * BHP_I *I.chi
688
689 F_fp_Y = A_fw * F_fw
690 F_fs_Y = A_fg * F_fg
691 D_fpp = A_fw * D_fw
692 D_fsp = A_fg * D_fg
693
694
695 if self.domain.getDim() == 3:
696 F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
697 F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
698 else:
699 F_fp_X = 0 * kronecker(self.domain)[1]
700 F_fs_X = 0 * kronecker(self.domain)[1]
701
702 F_mg_Y = H * L_g
703
704
705 D=self.__pde.createCoefficient("D")
706 A=self.__pde.createCoefficient("A")
707 Y=self.__pde.createCoefficient("Y")
708 X=self.__pde.createCoefficient("X")
709
710
711 D[0,0]=E_fpp + dt * D_fpp
712 D[0,1]=E_fps
713 D[1,0]=E_fsp + dt * D_fsp
714 D[1,1]=E_fss
715 D[1,2]= rho_g_surf # rho_g_surf
716 D[0,2]= -dt * H * dL_gdp * 0
717 D[2,2]= 1 + dt * H
718
719
720 for i in range(self.domain.getDim()):
721 A[0,i,0,i] = dt * A_fw * self.perm_f[i]
722 A[1,i,1,i] = dt * A_fg * self.perm_f[i]
723
724 X[0,:]= dt * F_fp_X
725 X[1,:]= dt * F_fs_X
726
727 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old + dt * F_fp_Y
728 Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + rho_g_surf * c_mg_old + dt * F_fs_Y
729 Y[2] = c_mg_old + dt * ( F_mg_Y - H * dL_gdp * p_f *0 )
730
731 print "HHH D[0,0] = ",D[0,0]
732 print "HHH D[0,1] = ",D[0,1]
733 print "HHH D[0,2] = ",D[0,2]
734 print "HHH D[1,0] = ",D[1,0]
735 print "HHH D[1,1] = ",D[1,1]
736 print "HHH D[1,2] = ",D[1,2]
737 print "HHH D[2,0] = ",D[2,0]
738 print "HHH D[2,1] = ",D[2,1]
739 print "HHH D[2,2] = ",D[2,2]
740 print "HHH A_fw = ",A_fw
741 print "HHH A_fg = ",A_fg
742 print "HHH A[0,0,0,0] = ",A[0,0,0,0]
743 print "HHH A[0,1,0,1] = ",A[0,1,0,1]
744 print "HHH A[1,0,1,0] = ",A[1,0,1,0]
745 print "HHH A[1,1,1,1] = ",A[1,1,1,1]
746 print "HHH Y[0] ",Y[0]
747 print "HHH Y[1] ",Y[1]
748 print "HHH Y[2] ",Y[2]
749 print "HHH F_fp_Y ",F_fp_Y
750 print "HHH F_fs_Y ",F_fs_Y
751 print "HHH F_mg_Y ",F_mg_Y
752 print "HHH H = ",H
753
754 self.__pde.setValue(A=A, D=D, X=X, Y=Y)
755
756 u2 = self.__pde.getSolution()
757 # this is deals with round-off errors to maintain physical meaningful values
758 # we should do this in a better way to detect values that are totally wrong
759 p_f=u2[0]
760 S_fg=clip(u2[1],minval=0, maxval=1)
761 c_mg=clip(u2[2],minval=0)
762
763
764 # update well production and BHP
765 wells_state_new={}
766 for I in self.wells:
767 q_I = I.getFlowRate(self.t+dt)
768 Pi_I = I.getProductivityIndex()
769 p_f_I=I.locator(p_f)
770 S_fg_I=I.locator(S_fg)
771 BHP_limit_I=I.getBHPLimit()
772
773 A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)*self.rho_w(p_f_I)
774 A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)*self.rho_g(p_f_I)
775 if I.getPhase() == Well.WATER:
776 BHP= p_f_I - self.rho_w.rho_surf*q_I/(A_fw_I * Pi_I)
777 else:
778 BHP= p_f_I - self.rho_g.rho_surf*q_I/(A_fg_I * Pi_I)
779 BHP=max(BHP, BHP_limit_I)
780 if self.verbose: print "well %s: BHP = %e"%(I.name, BHP)
781
782 wells_state_new[I.name]={ "bhp" : BHP, "q_gas": A_fg_I * Pi_I*(p_f_I-BHP)/self.rho_g.rho_surf, "q_water": A_fw_I * Pi_I*(p_f_I-BHP)/self.rho_w.rho_surf }
783
784 print wells_state_new
785 return (p_f,S_fg, c_mg, wells_state_new)

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26