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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26