/[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 3515 - (hide 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 gross 3494 ######################################################
2 gross 3458 #
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 gross 3484 from esys.escript.linearPDEs import LinearPDE
24 gross 3458 from esys.escript import unitsSI as U
25 gross 3494 from esys.escript.pdetools import Locator
26 gross 3515 from esys.escript import *
27 gross 3494 from esys.weipa import saveVTK
28 gross 3484 import math
29 gross 3458
30 gross 3510 USE_NODAL_WELL = False or True
31    
32 gross 3478 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 gross 3458
57 gross 3478 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 gross 3484
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 gross 3478
88 gross 3484 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 gross 3510 rho = rho_surf (1 + X + X*X/2) with X= C * ( p - p_ref )
108 gross 3484
109 gross 3510 with rho_surf =rho_s/B_ref * gravity
110 gross 3484 """
111 gross 3510 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 gross 3484 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 gross 3510 return self.__rho_0 * (1+ X * (1+ X/2) )
123 gross 3484
124     def getValueDifferential(self, p):
125     """
126     """
127     X= self.C * ( p - self.p_ref )
128 gross 3510 return self.__rho_0 * self.C * (1+ X)
129 gross 3497
130     def getFormationFactor(self, p):
131 gross 3510 return self.rho_surf/self.getValue(p)
132 gross 3484
133 gross 3497
134    
135 gross 3484 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 gross 3510 self.rho_surf =rho_air_s * gravity
168 gross 3484 self.tab=InterpolationTable(x=p, y=B)
169    
170     def getValue(self, p):
171     """
172     returns the density for given pressure p
173     """
174 gross 3510 return self.rho_surf/self.getFormationFactor(p)
175 gross 3484
176     def getValueDifferential(self, p):
177     """
178     """
179 gross 3497 B = self.getFormationFactor(p)
180     dBdp = self.getFormationFactorDifferential(p)
181 gross 3510 return -self.rho_surf * dBdp/(B * B)
182 gross 3497
183     def getFormationFactor(self, p):
184     return self.tab.getValue(p)
185    
186     def getFormationFactorDifferential(self, p):
187     return self.tab.getValueDifferential(p)
188 gross 3484
189 gross 3497
190 gross 3484 class InterpolationTable(MaterialPropertyWithDifferential):
191 gross 3478 """
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 gross 3484 MaterialPropertyWithDifferential.__init__(self)
201 gross 3478 if len(x) != len(y):
202     raise ValueError,"length of interpolation nodes and value lists need to be identical."
203 gross 3484 if len(x) < 1 :
204 gross 3478 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 gross 3484
215 gross 3478 def getValue(self, x):
216     """
217     returns the interpolated values of x
218     """
219 gross 3486 X=self.__x
220     Y=self.__y
221    
222     x0=X[0]
223 gross 3478 m0=whereNegative(x-x0)
224     if self.__obeyBounds:
225     if sup(m0) > 0:
226 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
227 gross 3478 out=self.__x[0]
228 gross 3486 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 gross 3478 out = out * m0 + z * (1-m0)
231 gross 3494 m0=whereNonPositive(x-X[i])
232    
233 gross 3478 if self.__obeyBounds:
234 gross 3494 if inf(m0) < 1 :
235 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
236 gross 3478 else:
237 gross 3502 out = out * m0 + X[-1] * (1-m0)
238 gross 3478 return out
239 gross 3458
240 gross 3484 def getValueDifferential(self, x):
241 gross 3486 X=self.__x
242     Y=self.__y
243    
244     x0=X[0]
245 gross 3484 m0=whereNegative(x-x0)
246     if self.__obeyBounds:
247     if sup(m0) > 0:
248 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
249 gross 3484 out=0.
250 gross 3486 for i in range(1,len(X)):
251     z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
252 gross 3484 out = out * m0 + z * (1-m0)
253 gross 3494 m0=whereNonPositive(x-X[i])
254 gross 3484
255     if self.__obeyBounds:
256     if inf(m0) < 1:
257 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
258 gross 3484 else:
259     out = out * m0
260     return out
261    
262    
263 gross 3478 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 gross 3510 def __init__(self, name, domain, Q=[0.], schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
273 gross 3478 """
274     set up well
275     """
276 gross 3510 if not len(schedule) == len(Q):
277     raise ValueError,"length of schedule and Q must match."
278 gross 3484 self.__schedule=schedule
279 gross 3510 self.__Q = Q
280 gross 3484 self.__phase=phase
281 gross 3494 self.__BHP_limit=BHP_limit
282     self.name=name
283     self.domain=domain
284 gross 3484
285 gross 3478 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 gross 3458
294 gross 3478 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 gross 3510 def getFlowRate(self,t):
304 gross 3478 """
305     returns the flow rate
306     """
307 gross 3515 out=0.
308     for i in range(len(self.__schedule)):
309     if t <= self.__schedule[i]:
310     out=self.__Q[i]
311     break
312 gross 3510 return out
313 gross 3494
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 gross 3478 def getBHP(self):
323     """
324     return bottom-hole pressure
325    
326     :note: needs to be overwritten
327     """
328 gross 3494 return self.__BHP
329    
330     def setBHP(self, BHP):
331     """
332     sets bottom-hole pressure
333     """
334     self.__BHP= BHP
335 gross 3478
336     def getPhase(self):
337     """
338     returns the pahse the well works on
339     """
340 gross 3484 return self.__phase
341    
342     class VerticalPeacemanWell(Well):
343 gross 3478 """
344     defines a well using Peaceman's formula
345     """
346 gross 3494 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 gross 3510 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5):
348 gross 3502 # reset_factor=0.1):
349 gross 3484 """
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 gross 3494 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 gross 3484 / ( (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 gross 3458
366 gross 3484 self.__D=D
367     self.r_el=r_el
368 gross 3515 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 gross 3494
376     def getLocation(self):
377     return self.X0
378    
379     def getProductivityIndex(self):
380 gross 3484 """
381     returns the productivity index of the well.
382     """
383     return self.__PI
384 gross 3494
385 gross 3484
386 gross 3478 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 gross 3510 wells=[], g=9.81*U.m/U.sec**2):
397 gross 3478 """
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 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
409 gross 3478 :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 gross 3458
440 gross 3478 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 gross 3487 self.t =0
456 gross 3510 self.g=g
457 gross 3487
458     self.__iter_max=1
459     self.__rtol=1.e-4
460 gross 3494 self.verbose=False
461     self.XI=0.5
462     def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
463 gross 3487 """
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 gross 3494 if verbose !=None: self.verbose=verbose
469     if xi !=None: self.XI=xi
470 gross 3478
471 gross 3487 def update(self, dt):
472 gross 3515 self.u, self.u_old =tuple([ v.copy() for v in self.u ] ), self.u
473 gross 3487 n=0
474 gross 3515 converged=False
475     while n < self.__iter_max and not converged:
476 gross 3487 u=self.solvePDE(dt)
477 gross 3515 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 gross 3494
487 gross 3515 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 gross 3487 n+=1
494     self.u=u
495     print "iteration completed."
496     self.t+=dt
497 gross 3458
498    
499 gross 3478 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 gross 3458
506 gross 3484 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
507 gross 3478 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 gross 3494 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
510 gross 3510 wells=[], g=9.81*U.m/U.sec**2):
511 gross 3478 """
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 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
523 gross 3478 :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 gross 3494 :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 gross 3478 """
548    
549     DualPorosity.__init__(self, domain,
550     phi_f=phi_f, phi_m=None, phi=phi,
551 gross 3484 S_fg=None, S_mg=None,
552 gross 3478 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 gross 3510 wells=wells, g=g)
557 gross 3478 self.L_g=L_g
558 gross 3494 self.sigma = sigma
559     self.A_mg = A_mg
560     self.f_rg = f_rg
561 gross 3478 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
562    
563 gross 3458
564 gross 3478 def getPDEOptions(self):
565     """
566     returns the `SolverOptions` of the underlying PDE
567     """
568     return self.__pde.getSolverOptions()
569    
570 gross 3486 def getState(self):
571 gross 3515 return self.u
572 gross 3486
573     def getOldState(self):
574 gross 3515 return self.u_old
575 gross 3486
576 gross 3510 def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0, c_mg=None):
577 gross 3484 """
578     sets the initial state
579    
580     :param p: pressure
581     :param S_fg: gas saturation in fractured rock
582 gross 3510 :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
583 gross 3484 using the gas adsorption curve.
584 gross 3485 """
585 gross 3515 if self.domain.getDim() == 2:
586     p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
587 gross 3510 else:
588     z=self.u.getDomain().getX()[0]
589     z_top=sup(z)
590     z_bottom=inf(z)
591 gross 3515 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 gross 3510 if c_mg == None:
596 gross 3515 c_mg_init=self.L_g(p_init)
597 gross 3484 else:
598 gross 3515 c_mg_init=Scalar(c_mg, Solution(self.domain))
599 gross 3484
600 gross 3515 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 gross 3487 def solvePDE(self, dt):
615 gross 3478
616 gross 3515 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 gross 3458
619 gross 3487 S_fw=1-S_fg
620 gross 3458
621 gross 3494 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 gross 3510 print "c_mg range = ",inf(c_mg),sup(c_mg)
626 gross 3515 print "wells state =",wells_state
627 gross 3494
628 gross 3487 k_fw=self.k_w(S_fw)
629 gross 3494 if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
630    
631 gross 3515
632 gross 3487 k_fg=self.k_g(S_fg)
633 gross 3494 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 gross 3487
641    
642 gross 3494 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 gross 3478
646 gross 3484 rho_fw = self.rho_w.getValue(p_f)
647 gross 3494 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 gross 3510 rho_g_surf = self.rho_g.rho_surf
652 gross 3484 drho_fgdp = self.rho_g.getValueDifferential(p_f)
653 gross 3510 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 gross 3502
657 gross 3510 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 gross 3478 A_fw = rho_fw * k_fw/mu_fw
662     A_fg = rho_fg * k_fg/mu_fg
663    
664 gross 3510 m = whereNegative(L_g - c_mg)
665 gross 3501 H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
666 gross 3478
667 gross 3484 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
668     E_fps= - phi_f * rho_fw
669 gross 3510 E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
670 gross 3484 E_fss= phi_f * rho_fg
671 gross 3478
672 gross 3494 F_fw=0.
673     F_fg=0.
674     D_fw=0.
675     D_fg=0.
676 gross 3497
677 gross 3515 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 gross 3494
684 gross 3515 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 gross 3494
689 gross 3515 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 gross 3478
695     if self.domain.getDim() == 3:
696 gross 3494 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 gross 3478 else:
699 gross 3484 F_fp_X = 0 * kronecker(self.domain)[1]
700     F_fs_X = 0 * kronecker(self.domain)[1]
701 gross 3478
702 gross 3510 F_mg_Y = H * L_g
703 gross 3458
704 gross 3494
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 gross 3510
710 gross 3494
711 gross 3510 D[0,0]=E_fpp + dt * D_fpp
712 gross 3484 D[0,1]=E_fps
713 gross 3510 D[1,0]=E_fsp + dt * D_fsp
714 gross 3494 D[1,1]=E_fss
715 gross 3510 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 gross 3478
719    
720     for i in range(self.domain.getDim()):
721 gross 3510 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 gross 3494
724 gross 3510 X[0,:]= dt * F_fp_X
725     X[1,:]= dt * F_fs_X
726 gross 3494
727 gross 3515 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 gross 3478
731 gross 3494 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 gross 3510 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 gross 3494
754 gross 3478 self.__pde.setValue(A=A, D=D, X=X, Y=Y)
755    
756     u2 = self.__pde.getSolution()
757 gross 3494 # 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 gross 3515 p_f=u2[0]
760     S_fg=clip(u2[1],minval=0, maxval=1)
761     c_mg=clip(u2[2],minval=0)
762    
763 gross 3510
764 gross 3515 # update well production and BHP
765     wells_state_new={}
766 gross 3497 for I in self.wells:
767 gross 3515 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 gross 3510
773 gross 3515 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