# Diff of /trunk/finley/test/python/coalgas.py

revision 3485 by gross, Thu Mar 24 22:44:40 2011 UTC revision 3486 by gross, Sun Mar 27 22:27:45 2011 UTC
22
23  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
24  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
25  from esys.escript import sqrt, log  from esys.escript import sqrt, log, whereNegative, sup, inf
26  import math  import math
27
28  class MaterialProperty(object):  class MaterialProperty(object):
# Line 199  class InterpolationTable(MaterialPropert Line 199  class InterpolationTable(MaterialPropert
199        """        """
200        returns the interpolated values of x        returns the interpolated values of x
201        """        """
202        x0=self.__x[0]        X=self.__x
203          Y=self.__y
204
205          x0=X[0]
206        m0=whereNegative(x-x0)        m0=whereNegative(x-x0)
207        if self.__obeyBounds:        if self.__obeyBounds:
208       if sup(m0) > 0:       if sup(m0) > 0:
209          raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])          raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
210        out=self.__x[0]        out=self.__x[0]
211        for i in range(1,len(self.__x)):        for i in range(1,len(X)):
212        x1=self.__x[i]        z=(Y[i]-Y[i-1])/(X[i]-X[i-1]) * (x-X[i-1]) + Y[i-1]
z=(y[i]-y[i-1])/(x[i]-x[i-1]) * (x-x[i-1]) + y[i-1]
213        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
214        m0=whereNegative(x-x[i])        m0=whereNegative(x-X[i])
215
216        if self.__obeyBounds:        if self.__obeyBounds:
217          if inf(m0) < 1:          if inf(m0) < 1:
218             raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])             raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
219        else:        else:
220          out = out * m0 + y[-1] * (1-m0)          out = out * m0 + V[-1] * (1-m0)
221        return out        return out
222
223     def getValueDifferential(self, x):     def getValueDifferential(self, x):
224        x0=self.__x[0]        X=self.__x
225          Y=self.__y
226
227          x0=X[0]
228        m0=whereNegative(x-x0)        m0=whereNegative(x-x0)
229        if self.__obeyBounds:        if self.__obeyBounds:
230       if sup(m0) > 0:       if sup(m0) > 0:
231          raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])          raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
232        out=0.        out=0.
233        for i in range(1,len(self.__x)):        for i in range(1,len(X)):
234        x1=self.__x[i]        z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
z=(y[i]-y[i-1])/(x[i]-x[i-1])
235        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
236        m0=whereNegative(x-x[i])        m0=whereNegative(x-X[i])
237
238        if self.__obeyBounds:        if self.__obeyBounds:
239          if inf(m0) < 1:          if inf(m0) < 1:
240             raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])             raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
241        else:        else:
242          out = out * m0          out = out * m0
243        return out          return out
# Line 436  class DualPorosity(object): Line 440  class DualPorosity(object):
440        self.rho_g = rho_g            self.rho_g = rho_g
441        self.wells=wells        self.wells=wells
442
443
444
445
446
# Line 505  class PorosityOneHalfModel(DualPorosity) Line 509  class PorosityOneHalfModel(DualPorosity)
509       """       """
510       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
511
512        def setInitialState(p=1.*U.atm, S_fg=0,  C_mg=None):        def getState(self):
513         return self.u[0], self.u[1],  self.u[2]
514
515          def getOldState(self):
516         return self.u_old[0], self.u_old[1],  self.u_old[2]
517
518          def setInitialState(self, p=1.*U.atm, S_fg=0,  C_mg=None):
519          """          """
520          sets the initial state          sets the initial state
521
# Line 514  class PorosityOneHalfModel(DualPorosity) Line 524  class PorosityOneHalfModel(DualPorosity)
524          :param C_mg: gas concentration in coal matrix. if not given it is calculated          :param C_mg: gas concentration in coal matrix. if not given it is calculated
526          """              """
527          self.u=self.__pde.getNewCoefficient("u")          self.u=self.__pde.createSolution()
528          u[0]=p          self.u[0]=p
529          u[1]=S_fg          self.u[1]=S_fg
530          if C_mg == None:          if C_mg == None:
531            u[2]= self.L_g(p)            self.u[2]= self.L_g(p)
532          else:          else:
533            u[2]=C_mg            self.u[2]=C_mg
534
535        def solvePDE(self):        def solvePDE(self):
536
537       p_f=self.u[0]       p_f, S_fg, C_mg=self.getState()
538       S_fg=self.u[1]       p_f_old, S_fg_old, C_mg_old=self.getOldState()
C_mg=self.u[3]

p_f_old=self.u_old[0]
S_fg_old=self.u_old[1]
C_mg_old=self.u_old[3]
539
540
541       phi_f   =self.phi_f.getValue(S_fg)       phi_f   =self.phi_f.getValue(S_fg)

Legend:
 Removed from v.3485 changed lines Added in v.3486