/[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 3527 - (show annotations)
Tue May 24 06:59:20 2011 UTC (8 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 25239 byte(s)
modifications to the coal seam gas 
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, X0=[0.,0.,0.], *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 self.locator=Locator(DiracDeltaFunctions(self.domain),X0[:self.domain.getDim()])
285 self.X0=self.locator.getX()
286
287 def getLocation(self):
288 return self.X0
289
290 def getProductivityIndex(self):
291 """
292 returns the productivity index of the well.
293 typically a Gaussian profile around the location of the well.
294
295 :note: needs to be overwritten
296 """
297 raise NotImplementedError
298
299 def getFlowRate(self,t):
300 """
301 returns the flow rate
302 """
303 out=0.
304 for i in range(len(self.__schedule)):
305 if t <= self.__schedule[i]:
306 out=self.__Q[i]
307 break
308 return out
309
310 def getBHPLimit(self):
311 """
312 return bottom-hole pressure limit
313
314 :note: needs to be overwritten
315 """
316 return self.__BHP_limit
317
318 def getPhase(self):
319 """
320 returns the pahse the well works on
321 """
322 return self.__phase
323
324 class VerticalPeacemanWell(Well):
325 """
326 defines a well using Peaceman's formula
327 """
328 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],
329 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5):
330 # reset_factor=0.1):
331 """
332 set up well
333
334 :param BHP: ottom-hole pressure
335 :param Q: flow rate
336 :param r: well radius
337 :param X: location
338 :param D: dimension of well block
339 :param perm: permeabilities
340 :param s: skin factor
341 """
342 Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit, X0=X0)
343 r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 + sqrt(perm[0]/perm[1]) * D[1]**2 )\
344 / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
345
346 r_el=0.11271331774384821*D[0]
347 self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
348 self.__r = r
349
350 self.__D=D
351 self.r_el=r_el
352
353 if self.domain.getDim() == 3:
354 self.geo_fac=1
355 else:
356 self.geo_fac=1./D[2]
357
358
359
360 def getProductivityIndex(self):
361 """
362 returns the productivity index of the well.
363 """
364 return self.__PI
365
366
367 class DualPorosity(object):
368 """
369 generic dual porosity model
370 """
371 def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
372 S_fg=None, S_mg=None,
373 perm_f_0=None, perm_f_1=None, perm_f_2=None,
374 perm_m_0=None, perm_m_1=None, perm_m_2=None,
375 k_w =None, k_g=None, mu_w =None, mu_g =None,
376 rho_w =None, rho_g=None,
377 wells=[], g=9.81*U.m/U.sec**2):
378 """
379 set up
380
381 :param domain: domain
382 :type domain: `esys.escript.Domain`
383 :param phi_f: porosity of the fractured rock as function of the gas pressure
384 :type phi_f: `MaterialPropertyWithDifferential`
385 :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
386 :type phi_m: `MaterialPropertyWithDifferential` or None
387 :param phi: total porosity if None phi=1.
388 :type phi: scalar or None
389 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
390 :type S_fg: `MaterialPropertyWithDifferential`
391 :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
392 :type S_mg: `MaterialPropertyWithDifferential`
393 :param perm_f_0: permeability in the x_0 direction in the fractured media
394 :type perm_f_0: scalar
395 :param perm_f_1: permeability in the x_1 direction in the fractured media
396 :type perm_f_1: scalar
397 :param perm_f_2: permeability in the x_2 direction in the fractured media
398 :type perm_f_2: scalar
399 :param perm_m_0: permeability in the x_0 direction in the coal matrix
400 :type perm_m_0: scalar
401 :param perm_m_1: permeability in the x_1 direction in the coal matrix
402 :type perm_m_1: scalar
403 :param perm_m_2: permeability in the x_2 direction in the coal matrix
404 :type perm_m_2: scalar
405 :param k_w: relative permeability of water as function of water saturation
406 :type k_w: `MaterialProperty`
407 :param k_g: relative permeability of gas as function of gas saturation
408 :type k_g: `MaterialProperty`
409 :param mu_w: viscosity of water as function of water pressure
410 :type mu_w: `MaterialProperty`
411 :param mu_g: viscosity of gas as function of gas pressure
412 :type mu_g: `MaterialProperty`
413 :param rho_w: density of water as function of water pressure
414 :type rho_w: `MaterialPropertyWithDifferential`
415 :param rho_g: density of gas as function of gas pressure
416 :type rho_g: `MaterialPropertyWithDifferential`
417 :param wells : list of wells
418 :type wells: list of `Well`
419 """
420
421 self.domain = domain
422 self.phi_f = phi_f
423 self.phi_m = phi_m
424 self.phi = phi
425 self.S_fg = S_fg
426 self.S_mg = S_mg
427 self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
428 self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
429 self.k_w = k_w
430 self.k_g = k_g
431 self.mu_w = mu_w
432 self.mu_g = mu_g
433 self.rho_w = rho_w
434 self.rho_g = rho_g
435 self.wells=wells
436 self.t =0
437 self.g=g
438
439 self.__iter_max=1
440 self.__rtol=1.e-4
441 self.verbose=False
442 self.XI=0.5
443 def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
444 """
445 sets parameters to control iteration process
446 """
447 if iter_max !=None: self.__iter_max=iter_max
448 if rtol !=None: self.__rtol = rtol
449 if verbose !=None: self.verbose=verbose
450 if xi !=None: self.XI=xi
451
452 def update(self, dt):
453 self.u, self.u_old =tuple([ v.copy() for v in self.u ] ), self.u
454 n=0
455 converged=False
456 while n < self.__iter_max and not converged:
457 u=self.solvePDE(dt)
458 if self.verbose: print "iteration step %d:"%n
459 converged=True
460 for i in range(len(u)):
461 if isinstance(u[i], Data):
462 norm_u=Lsup(u[i])
463 norm_e=Lsup(u[i]-self.u[i])
464 else:
465 norm_e=0.
466 norm_u=1.
467
468 if norm_u > 0:
469 rerr=norm_e/norm_u
470 else:
471 rerr=norm_e
472 if norm_e>self.__rtol * norm_u + 1.e-10: converged=False
473 if self.verbose: print " comp %i: change = %e (value = %e)"%(i, norm_e,norm_u)
474 n+=1
475 self.u=u
476 print "iteration completed."
477 self.t+=dt
478
479
480 class PorosityOneHalfModel(DualPorosity):
481 """
482 Model for gas and water in double prosity model tracking water and gas
483 pressure in the fractured rock and gas concentration in the coal matrix.
484 This corresponds to the coal bed methan model in the ECLIPSE code.
485 """
486
487 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
488 perm_f_0=None, perm_f_1=None, perm_f_2=None,
489 k_w =None, k_g=None, mu_w =None, mu_g =None,
490 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
491 wells=[], g=9.81*U.m/U.sec**2):
492 """
493 set up
494
495 :param domain: domain
496 :type domain: `esys.escript.Domain`
497 :param phi_f: porosity of the fractured rock as function of the gas pressure
498 :type phi_f: `MaterialPropertyWithDifferential`
499 :param phi: total porosity if None phi=1.
500 :type phi: scalar or None
501 :param L_g: gas adsorption as function of gas pressure
502 :type L_g: `MaterialPropertyWithDifferential`
503 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
504 :type S_fg: `MaterialPropertyWithDifferential`
505 :param perm_f_0: permeability in the x_0 direction in the fractured media
506 :type perm_f_0: scalar
507 :param perm_f_1: permeability in the x_1 direction in the fractured media
508 :type perm_f_1: scalar
509 :param perm_f_2: permeability in the x_2 direction in the fractured media
510 :type perm_f_2: scalar
511 :param k_w: relative permeability of water as function of water saturation
512 :type k_w: `MaterialProperty`
513 :param k_g: relative permeability of gas as function of gas saturation
514 :type k_g: `MaterialProperty`
515 :param mu_w: viscosity of water as function of water pressure
516 :type mu_w: `MaterialProperty`
517 :param mu_g: viscosity of gas as function of gas pressure
518 :type mu_g: `MaterialProperty`
519 :param rho_w: density of water as function of water pressure
520 :type rho_w: `MaterialPropertyWithDifferential`
521 :param rho_g: density of gas as function of gas pressure
522 :type rho_g: `MaterialPropertyWithDifferential`
523 :param wells : list of wells
524 :type wells: list of `Well`
525 :param sigma: shape factor for gas matrix diffusion
526 :param A_mg: diffusion constant for gas adsorption
527 :param f_rg: gas re-adsorption factor
528 """
529
530 DualPorosity.__init__(self, domain,
531 phi_f=phi_f, phi_m=None, phi=phi,
532 S_fg=None, S_mg=None,
533 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
534 perm_m_0=None , perm_m_1=None, perm_m_2=None,
535 k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
536 rho_w =rho_w, rho_g=rho_g,
537 wells=wells, g=g)
538 self.L_g=L_g
539 self.sigma = sigma
540 self.A_mg = A_mg
541 self.f_rg = f_rg
542 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
543
544
545 def getPDEOptions(self):
546 """
547 returns the `SolverOptions` of the underlying PDE
548 """
549 return self.__pde.getSolverOptions()
550
551 def getState(self):
552 return self.u
553
554 def getOldState(self):
555 return self.u_old
556
557 def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0, c_mg=None):
558 """
559 sets the initial state
560
561 :param p: pressure
562 :param S_fg: gas saturation in fractured rock
563 :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
564 using the gas adsorption curve.
565 """
566 if self.domain.getDim() == 2:
567 p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
568 else:
569 z=self.u.getDomain().getX()[0]
570 z_top=sup(z)
571 z_bottom=inf(z)
572 p_init=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
573
574 S_fg_init=Scalar(S_fg, Solution(self.domain))
575
576 if c_mg == None:
577 c_mg_init=self.L_g(p_init)
578 else:
579 c_mg_init=Scalar(c_mg, Solution(self.domain))
580
581 q_gas=Scalar(0., DiracDeltaFunctions(self.domain))
582 q_water=Scalar(0., DiracDeltaFunctions(self.domain))
583 BHP=interpolate(p_init, DiracDeltaFunctions(self.domain))
584
585 self.u=(p_init,S_fg_init, c_mg_init, BHP, q_gas, q_water)
586
587 def solvePDE(self, dt):
588
589 p_f, S_fg, c_mg, BHP_old, q_gas_old, q_water_old =self.getState()
590 p_f_old, S_fg_old, c_mg_old, BHP, q_gas, q_water =self.getOldState()
591
592 S_fw=1-S_fg
593
594 if self.verbose:
595 print "p_f range = ",inf(p_f),sup(p_f)
596 print "S_fg range = ",inf(S_fg),sup(S_fg)
597 print "S_fw range = ",inf(S_fw),sup(S_fw)
598 print "c_mg range = ",inf(c_mg),sup(c_mg)
599 print "BHP =",BHP
600 print "q_gas =",q_gas
601 print "q_water =",q_water
602
603 k_fw=self.k_w(S_fw)
604 if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
605
606
607 k_fg=self.k_g(S_fg)
608 if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
609
610 mu_fw=self.mu_w(p_f)
611 if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
612
613 mu_fg=self.mu_g(p_f)
614 if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
615
616
617 phi_f =self.phi_f.getValue(p_f)
618 dphi_fdp=self.phi_f.getValueDifferential(p_f)
619 if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
620
621 rho_fw = self.rho_w.getValue(p_f)
622 drho_fwdp = self.rho_w.getValueDifferential(p_f)
623 if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
624
625 rho_fg = self.rho_g.getValue(p_f)
626 rho_g_surf = self.rho_g.rho_surf
627 drho_fgdp = self.rho_g.getValueDifferential(p_f)
628 if self.verbose:
629 print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
630 print "rho_fg surf = ",rho_g_surf
631
632 L_g = self.L_g(p_f)
633 dL_gdp = self.L_g.getValueDifferential(p_f)
634 if self.verbose: print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
635
636 A_fw = rho_fw * k_fw/mu_fw
637 A_fg = rho_fg * k_fg/mu_fg
638
639 m = whereNegative(L_g - c_mg)
640 H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
641
642 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
643 E_fps= - phi_f * rho_fw
644 E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
645 E_fss= phi_f * rho_fg
646
647
648
649 F_fw=0.
650 F_fg=0.
651 D_fw=0.
652 D_fg=0.
653
654 prod_index=Scalar(0., DiracDeltaFunctions(self.domain))
655 geo_fac=Scalar(0., DiracDeltaFunctions(self.domain))
656 for I in self.wells:
657 prod_index.setTaggedValue(I.name, I.getProductivityIndex() )
658 geo_fac.setTaggedValue(I.name, I.geo_fac)
659
660 F_fp_Y = A_fw * prod_index * BHP * geo_fac
661 F_fs_Y = A_fg * prod_index * BHP * geo_fac
662 D_fpp = A_fw * prod_index * geo_fac
663 D_fsp = A_fg * prod_index * geo_fac
664
665
666 if self.domain.getDim() == 3:
667 F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
668 F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
669 else:
670 F_fp_X = 0 * kronecker(self.domain)[1]
671 F_fs_X = 0 * kronecker(self.domain)[1]
672
673 F_mg_Y = H * L_g
674
675
676 D=self.__pde.createCoefficient("D")
677 A=self.__pde.createCoefficient("A")
678 Y=self.__pde.createCoefficient("Y")
679 X=self.__pde.createCoefficient("X")
680
681 d_dirac=self.__pde.createCoefficient("d_dirac")
682 y_dirac=self.__pde.createCoefficient("y_dirac")
683
684
685
686 D[0,0]=E_fpp
687 D[0,1]=E_fps
688 d_dirac[0,0]=dt * D_fpp
689
690 D[1,0]=E_fsp
691 D[1,1]=E_fss
692 D[1,2]= rho_g_surf
693 d_dirac[1,0]=dt * D_fsp
694
695 D[0,2]= -dt * H * dL_gdp * 0
696 D[2,2]= 1 + dt * H
697
698
699 for i in range(self.domain.getDim()):
700 A[0,i,0,i] = dt * A_fw * self.perm_f[i]
701 A[1,i,1,i] = dt * A_fg * self.perm_f[i]
702
703 X[0,:]= dt * F_fp_X
704 X[1,:]= dt * F_fs_X
705
706 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old
707 Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + rho_g_surf * c_mg_old
708 Y[2] = c_mg_old + dt * ( F_mg_Y - H * dL_gdp * p_f *0 )
709
710 y_dirac[0] =dt * F_fp_Y
711 y_dirac[1] =dt * F_fs_Y
712
713 print "HHH D[0,0] = ",D[0,0]
714 print "HHH D[0,1] = ",D[0,1]
715 print "HHH D[0,2] = ",D[0,2]
716 print "HHH D[1,0] = ",D[1,0]
717 print "HHH D[1,1] = ",D[1,1]
718 print "HHH D[1,2] = ",D[1,2]
719 print "HHH D[2,0] = ",D[2,0]
720 print "HHH D[2,1] = ",D[2,1]
721 print "HHH D[2,2] = ",D[2,2]
722 print "HHH A_fw = ",A_fw
723 print "HHH A_fg = ",A_fg
724 print "HHH A[0,0,0,0] = ",A[0,0,0,0]
725 print "HHH A[0,1,0,1] = ",A[0,1,0,1]
726 print "HHH A[1,0,1,0] = ",A[1,0,1,0]
727 print "HHH A[1,1,1,1] = ",A[1,1,1,1]
728 print "HHH Y[0] ",Y[0]
729 print "HHH Y[1] ",Y[1]
730 print "HHH Y[2] ",Y[2]
731 print "HHH F_fp_Y ",F_fp_Y
732 print "HHH F_fs_Y ",F_fs_Y
733 print "HHH F_mg_Y ",F_mg_Y
734 print "HHH H = ",H
735
736 self.__pde.setValue(A=A, D=D, X=X, Y=Y, d_dirac=d_dirac , y_dirac=y_dirac)
737
738 u2 = self.__pde.getSolution()
739 # this is deals with round-off errors to maintain physical meaningful values
740 # we should do this in a better way to detect values that are totally wrong
741 p_f=u2[0]
742 S_fg=clip(u2[1],minval=0, maxval=1)
743 c_mg=clip(u2[2],minval=0)
744
745
746
747 q=Scalar(0., DiracDeltaFunctions(self.domain))
748 BHP_limit=Scalar(0., DiracDeltaFunctions(self.domain))
749 water_well_mask=Scalar(0., DiracDeltaFunctions(self.domain))
750 for I in self.wells:
751 q.setTaggedValue(I.name, I.getFlowRate(self.t+dt))
752 BHP_limit.setTaggedValue(I.name, I.getBHPLimit())
753 if I.getPhase() == Well.WATER:
754 water_well_mask.setTaggedValue(I.name, 1)
755
756 p_f_wells=interpolate(p_f, DiracDeltaFunctions(self.domain))
757 S_fg_wells=interpolate(S_fg, DiracDeltaFunctions(self.domain))
758 A_fw_wells= self.k_w(1-S_fg_wells)/self.mu_w(p_f_wells)*self.rho_w(p_f_wells)
759 A_fg_wells= self.k_g(S_fg_wells)/self.mu_g(p_f_wells)*self.rho_g(p_f_wells)
760
761 BHP=clip( p_f_wells - q/prod_index * ( m * self.rho_w.rho_surf/A_fw_wells + (1-m)* self.rho_g.rho_surf/A_fg_wells ), minval = BHP_limit)
762 BHP=clip( p_f_wells - q/prod_index * self.rho_w.rho_surf/A_fw_wells, minval = BHP_limit)
763 q_gas = prod_index * A_fg_wells * (p_f_wells - BHP)/self.rho_g.rho_surf
764 q_water = prod_index * A_fw_wells * (p_f_wells - BHP)/self.rho_w.rho_surf
765
766 return (p_f,S_fg, c_mg, BHP, q_gas, q_water )

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26