/[escript]/trunk/finley/test/python/runcoalgas.py
ViewVC logotype

Annotation of /trunk/finley/test/python/runcoalgas.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3485 - (hide annotations)
Thu Mar 24 22:44:40 2011 UTC (8 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 25386 byte(s)
some compile fixes
1 gross 3478 #######################################################
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     Coal Seam gasL ECLIPSE test case
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 import *
24 gross 3484 from esys.escript import unitsSI as U
25 gross 3478 from coalgas import *
26 gross 3484 import time
27     from esys.finley import Rectangle
28 gross 3478
29 gross 3484 L_X=168*U.km
30     L_Y=168*U.km
31     L_Z=10*U.m
32 gross 3478
33 gross 3484 N_X=21
34     N_Y=21
35 gross 3478
36    
37 gross 3484 PERM_F_X = 100 * U.mDarcy
38     PERM_F_Y = 100 * U.mDarcy
39     PERM_F_Z = 1e-4 * U.mDarcy
40     PHI_F_0 = 0.01
41     P_F_0 = 69 * U.bar
42 gross 3478
43 gross 3484 # these object correspond to the ECLIPSE input files
44     PVTW={ "p_ref" : 1000 * U.bar ,
45     "B_ref" : 0.997 ,
46     "C" : 3.084E-06 /U.bar,
47     "mu_ref" : 0.68673 * U.cPoise,
48     "C_v" : 0/U.bar
49    
50     }
51     GRAVITY = { "water" : 1.0,
52     "gas" : .553 }
53     ROCK = { "p_ref" : 1000 * U.bar ,
54     "C" : 3.3E-004 * 1./U.bar }
55    
56     LANGMUIR = [
57     [ 0 * U.bar , 0.00000000 ],
58     [ 100 * U.bar , 0.00213886 ],
59     [ 200 * U.bar , 0.00383259 ],
60     [ 300 * U.bar , 0.00520706 ],
61     [ 400 * U.bar , 0.00634474 ],
62     [ 500 * U.bar , 0.00730199 ],
63     [ 600 * U.bar , 0.00811857 ],
64     [ 700 * U.bar , 0.00882336 ],
65     [ 800 * U.bar , 0.00943786 ],
66     [ 900 * U.bar , 0.00997836 ],
67     [ 1000 * U.bar , 0.01045748 ],
68     [ 1200 * U.bar , 0.01126912 ] ]
69    
70     PVDG = [
71     [ 14.70 * U.bar ,200.3800 , 0.012025 * U.cPoise ] ,
72     [ 20.00 * U.bar ,146.0600 , 0.012030 * U.cPoise ] ,
73     [ 25.00 * U.bar ,116.1461 , 0.012034 * U.cPoise ] ,
74     [ 30.00 * U.bar ,96.3132 , 0.012038 * U.cPoise ] ,
75     [ 35.00 * U.bar ,82.2113 , 0.012043 * U.cPoise ] ,
76     [ 49.33 * U.bar ,57.7891 , 0.012055 * U.cPoise ] ,
77     [ 59.00 * U.bar ,48.0866 , 0.012064 * U.cPoise ] ,
78     [ 69.00 * U.bar ,40.9441 , 0.012073 * U.cPoise ] ,
79     [ 75.00 * U.bar ,37.5839 , 0.012078 * U.cPoise ] ,
80     [ 83.00 * U.bar ,33.8685 , 0.012085 * U.cPoise ] ,
81     [ 90.00 * U.bar ,31.1661 , 0.012092 * U.cPoise ] ,
82     [ 95.00 * U.bar ,29.4827 , 0.012097 * U.cPoise ] ,
83     [ 100.00 * U.bar ,27.9698 , 0.012101 * U.cPoise ] ,
84     [ 105.00 * U.bar ,26.6028 , 0.012106 * U.cPoise ] ,
85     [ 118.60 * U.bar ,23.4749 , 0.012119 * U.cPoise ] ,
86     [ 120.00 * U.bar ,23.1937 , 0.012120 * U.cPoise ] ,
87     [ 140.00 * U.bar ,19.7977 , 0.012140 * U.cPoise ] ,
88     [ 153.23 * U.bar ,18.0443 , 0.012153 * U.cPoise ] ,
89     [ 160.00 * U.bar ,17.2607 , 0.012159 * U.cPoise ] ,
90     [ 170.00 * U.bar ,16.2188 , 0.012169 * U.cPoise ] ,
91     [ 187.86 * U.bar ,14.6373 , 0.012188 * U.cPoise ] ,
92     [ 222.49 * U.bar ,12.3027 , 0.012224 * U.cPoise ] ,
93     [ 257.13 * U.bar ,10.6038 , 0.012262 * U.cPoise ] ,
94     [ 291.76 * U.bar ,9.3134 , 0.012301 * U.cPoise ] ,
95     [ 326.39 * U.bar ,8.3001 , 0.012341 * U.cPoise ] ,
96     [ 361.02 * U.bar ,7.4835 , 0.012383 * U.cPoise ] ,
97     [ 395.66 * U.bar ,6.8114 , 0.012425 * U.cPoise ] ,
98     [ 430.29 * U.bar ,6.2491 , 0.012470 * U.cPoise ] ,
99     [ 464.92 * U.bar ,5.7715 , 0.012515 * U.cPoise ] ,
100     [ 499.55 * U.bar ,5.3610 , 0.012562 * U.cPoise ] ,
101     [ 534.19 * U.bar ,5.0043 , 0.012610 * U.cPoise ] ,
102     [ 568.82 * U.bar ,4.6917 , 0.012659 * U.cPoise ] ,
103     [ 603.45 * U.bar ,4.4154 , 0.012710 * U.cPoise ] ,
104     [ 638.08 * U.bar ,4.1695 , 0.012762 * U.cPoise ] ,
105     [ 672.72 * U.bar ,3.9491 , 0.012815 * U.cPoise ] ,
106     [ 707.35 * U.bar ,3.7507 , 0.012869 * U.cPoise ] ,
107     [ 741.98 * U.bar ,3.5711 , 0.012925 * U.cPoise ] ,
108     [ 776.61 * U.bar ,3.4076 , 0.012982 * U.cPoise ] ,
109     [ 811.25 * U.bar ,3.2583 , 0.013041 * U.cPoise ] ,
110     [ 845.88 * U.bar ,3.1214 , 0.013100 * U.cPoise ] ,
111     [ 880.51 * U.bar ,2.9953 , 0.013161 * U.cPoise ] ,
112     [ 915.14 * U.bar ,2.8790 , 0.013223 * U.cPoise ] ,
113     [ 949.78 * U.bar ,2.7712 , 0.013287 * U.cPoise ] ,
114     [ 984.41 * U.bar ,2.6711 , 0.013352 * U.cPoise ] ,
115     [ 1019.00 * U.bar ,2.5781 , 0.013418 * U.cPoise ] ,
116     [ 1053.70 * U.bar ,2.4909 , 0.013486 * U.cPoise ] ,
117     [ 1088.30 * U.bar ,2.4096 , 0.013554 * U.cPoise ] ,
118     [ 1122.90 * U.bar ,2.3334 , 0.013624 * U.cPoise ] ,
119     [ 1157.60 * U.bar ,2.2616 , 0.013696 * U.cPoise ] ,
120     [ 1192.20 * U.bar ,2.1942 , 0.013768 * U.cPoise ] ,
121     [ 1226.80 * U.bar ,2.1307 , 0.013842 * U.cPoise ] ,
122     [ 1261.50 * U.bar ,2.0705 , 0.013917 * U.cPoise ] ,
123     [ 1296.10 * U.bar ,2.0138 , 0.013994 * U.cPoise ] ,
124     [ 1330.70 * U.bar ,1.9600 , 0.014072 * U.cPoise ] ,
125     [ 1365.40 * U.bar ,1.9089 , 0.014151 * U.cPoise ] ]
126    
127    
128     SGFN = [
129     [ 0 , 0 , 0 ],
130     [ 0.05 , 0 , 0 ],
131     [ 0.1333 , 0.00610 , 0 ],
132     [ 0.2167 , 0.02990 , 0 ],
133     [ 0.3 , 0.0759 , 0 ],
134     [ 0.3833 , 0.1471 , 0 ],
135     [ 0.46667 , 0.2458 , 0 ],
136     [ 0.55 , 0.3739 , 0 ],
137     [ 0.6333 , 0.53300 , 0 ],
138     [ 0.7167 , 0.7246 , 0 ],
139     [ 0.8 , 0.95 , 0 ] ]
140     SWFN = [
141     [ 0.20000 , 0.00000, 0 ],
142     [ 0.28330 , 0.03280, 0 ],
143     [ 0.36670 , 0.09270, 0 ],
144     [ 0.45000 , 0.17030, 0 ],
145     [ 0.53330 , 0.26220, 0 ],
146     [ 0.61670 , 0.36650, 0 ],
147     [ 0.70000 , 0.48170, 0 ],
148     [ 0.78330 , 0.60710, 0 ],
149     [ 0.86670 , 0.74170, 0 ],
150     [ 0.95000 , 0.88500, 0 ],
151     [ 1.00000 , 1.00000, 0 ] ]
152    
153    
154     wellspecs = {
155     'P1' : { "X0" : [106, 106],
156     "r" : 0.253*U.m,
157     "s" : 0,
158     "Q" : 2000*U.Barrel/U.day*GRAVITY["water"],
159     "BHP" : 75*U.psi,
160     "schedule" : [0.*U.yr, 4*U.yr]
161     }
162     }
163    
164     CELL_X=L_X/N_X
165     CELL_Y=L_Y/N_Y
166     # print input
167     print("<%s> Execution started."%time.asctime())
168    
169     print "length x-direction = %f km"%(L_X/U.km)
170     print "number of cells in x direction = %d"%N_X
171     print "cell size in x direction = %f m"%(CELL_X/U.m)
172     print "length y-direction = %f km"%(L_Y/U.km)
173     print "number of cells in y direction = %d"%N_Y
174     print "cell size in y direction = %f m"%(CELL_Y/U.m)
175     print "fracture permeability in x direction= %f mD"%(PERM_F_X/(U.mDarcy))
176     print "fracture permeability in y direction= %f mD"%(PERM_F_Y/(U.mDarcy))
177     print "fracture permeability in z direction= %f mD"%(PERM_F_Z/(U.mDarcy))
178     print "initial porosity in fractured rock= %f"%PHI_F_0
179    
180     domain=Rectangle(N_X, N_Y, l0=L_X, l1=L_Y)
181     print("<%s> Mesh set up completed."%time.asctime())
182    
183    
184     model = PorosityOneHalfModel(domain,
185     phi_f=Porosity(phi_0=PHI_F_0, p_0=P_F_0, p_ref=ROCK["p_ref"], C = ROCK["C"]),
186     L_g=InterpolationTable([ l[0] for l in LANGMUIR ], [ l[1] for l in LANGMUIR ] ),
187     perm_f_0=PERM_F_X,
188     perm_f_1=PERM_F_Y,
189     perm_f_2=PERM_F_Z,
190     k_w =InterpolationTable([ l[0] for l in SWFN ], [ l[1] for l in SWFN ] ),
191     k_g= InterpolationTable([ l[0] for l in SGFN ], [ l[1] for l in SGFN ] ),
192     mu_w = WaterViscosity(mu_ref = PVTW["mu_ref"], p_ref=PVTW["p_ref"], C=PVTW["C_v"]),
193     mu_g = InterpolationTable([ l[0] for l in PVDG ], [ l[2] for l in PVDG ] ),
194     rho_w = WaterDensity(B_ref=PVTW["B_ref"], p_ref = PVTW["p_ref"], C=PVTW["C"], gravity=GRAVITY["water"]),
195     rho_g=GasDensity( p = [ l[0] for l in PVDG ], B = [ l[1] for l in PVDG ], gravity=GRAVITY["gas"]),
196     wells=[ VerticalPeacemanWell(i,BHP=wellspecs[i]["BHP"],
197     Q=wellspecs[i]["Q"],
198     r=wellspecs[i]["r"],
199     X0=[ (wellspecs[i]["X0"][0]+0.5)*CELL_X, (wellspecs[i]["X0"][1]+0.5)*CELL_Y],
200     D=[CELL_X, CELL_Y, L_Z],
201     perm=[PERM_F_X, PERM_F_Y, PERM_F_Z],
202     schedule=wellspecs[i]["schedule"],
203     s=wellspecs[i]["s"]) for i in wellspecs]
204     )
205    
206     model.setInitialState(p=P_F_0, S_fg=0, C_mg=None)
207    
208 gross 3485
209 gross 3484 1/0
210    
211    
212 gross 3478 #=======================================================
213    
214     from esys.escript.linearPDEs import LinearPDE
215     from esys.finley import Rectangle, Brick
216     from esys.escript import unitsSI as U
217     from esys.weipa import saveVTK
218     from math import pi, ceil
219     import sys
220     import time
221    
222     # ======================= Default Values ==================================================
223     H=2.*U.m # height
224     L=H # length
225     NE_H=5 # number of elements in H-direction.
226     NE_L=int(ceil(L*NE_H/H))
227     DIM=2
228     T_END=40000.0*U.sec # end time at max strain = 0.5
229     DT=1.*U.sec
230    
231     CHI=1.0 # =0 = explicit Euler, 0.5 = Crank-Nicholson, 1. = implicit Euler
232    
233     G=-9.81*U.m/U.sec**2
234     S_W_MIN = 0.01
235     S_W_MAX = 0.75
236    
237     M=3./2. *0+1
238     GAMMA=1.e-4*(1./U.Pa)
239    
240     V_L=0.015*U.m**3/U.kg *0
241     P_L=6.109*U.Mega*U.Pa
242     P_A=101.325*U.Kilo*U.Pa
243     RHO_AG=0.717*U.kg/U.m**3
244     P_0G=P_A
245     RHO_0G=RHO_AG
246     RHO_W = 1000.*U.kg/U.m**3
247     RHO_COAL=1250*U.kg/U.m**3
248     ETA_W=1.e-3*U.kg/U.m/U.sec
249     ETA_G=1.1e-5*U.kg/U.m/U.sec
250     LAM=1./3.
251     PHI_0=0.00804
252     EPS_L=0.02295
253     K_S=8425*U.Mega*U.Pa
254     ALPHA=0.005*1./U.Pa *0
255     EPS_V=0
256     K_0=3.7996e-17*U.m**3
257     D_G=20e-10*0
258    
259    
260     TOL=1.e-4
261     OUT_DIR="ReservoirResults"
262    
263    
264     class Porosity(object):
265     def __init__(self, phi0=1., alpha=0., Ks=1., eps_v=0, p0=0., eps_L=0, p_L=0):
266     self.phi0=phi0
267     self.alpha=alpha
268     self.eps_v=eps_v
269     self.p_0=p0
270     self.p_L=p_L
271     self.eps_L = eps_L
272     self.Ks = Ks
273    
274     def getValue(self,p=0):
275     return self.phi0+ self.alpha*(self.eps_v + (p-self.p_0)/self.Ks + self.eps_L * self.p_L * (self.p_0-p)/((self.p_L + self.p_0) * (p + self.p_L)))
276    
277     def __call__(self, p=0):
278     return self.getValue(p)
279    
280     def getDp(self, p=0):
281     return self.alpha * (1./self.Ks - self.eps_L * self.p_L/(p-self.p_L)**2 )
282    
283     class Density(object):
284     def __init__(self, p_a =0, rho_a=0, p0 =0, rho0=0):
285     self.p_a=p_a
286     self.rho_a=rho_a
287     self.p0=p0
288     self.rho0=rho0
289     def getValue(self,p=0):
290     return (self.rho_a - self.rho0)/(self.p_a - self.p0)*(p-self.p0) + self.rho0
291     def __call__(self, p=0):
292     return self.getValue(p)
293    
294     def getDp(self, p=0):
295     return (self.rho_a - self.rho0)/(self.p_a - self.p0)
296    
297    
298     class GasStorageCapacity(object):
299     def __init__(self):
300     pass
301    
302     def __call__(self, p=0):
303     return self.getValue(p)
304    
305     def getValue(self,p=0):
306     return 0.*p
307     def getDp(self, p=0):
308     return 0.*p
309    
310     class GasStorageCapacity(object):
311     def __init__(self):
312     raise NotImplementedError
313    
314     def getValue(self,p=0):
315     raise NotImplementedError
316     def getDp(self, p=0):
317     raise NotImplementedError
318     def __call__(self, p=0):
319     return self.getValue(p)
320    
321    
322    
323     class LangmuirIsotherm(GasStorageCapacity):
324     def __init__(self, V_L=0, p_L=0., rho_a=1. ):
325     self.V_L=V_L
326     self.p_L=p_L
327     self.rho_a=rho_a
328    
329     def getValue(self,p=0):
330     return (self.V_L*self.rho_a)*p/(self.p_L +p)
331    
332     def getDp(self, p=0):
333     return (self.V_L*self.rho_a)*p/(self.p_L +p)**2
334    
335    
336     class WaterVolumeFraction(object):
337     """
338     class defining water volume fraction S_w as function of suction pressure p
339     """
340     def __init__(self):
341     raise NotImplementedError
342     def __call__(self, p=0):
343     return self.getValue(p)
344     def getValue(self,p=0):
345     """
346     returns water volume fraction for given suction pressure p
347     """
348     raise NotImplementedError
349     def getDp(self, p=0):
350     """
351     returns derivative of water volume fraction with respect to suction pressure for given suction pressure p
352     """
353     raise NotImplementedError
354    
355     class VanGenuchten(WaterVolumeFraction):
356     def __init__(self, S_min=0., S_max=1., m =1, gamma=0):
357     self.S_min = S_min
358     self.S_max = S_max
359     self.m = m
360     self.gamma = gamma
361     def getEffFraction(self,s):
362     return (s-self.S_min)/(self.S_max-self.S_min)
363     def getFraction(self,s_eff):
364     return self.S_min+s_eff*(self.S_max-self.S_min)
365    
366     def getValue(self,p=0):
367     s_hat = (1+(abs(p) * self.gamma)**self.m)**(1./self.m-1)
368     return self.getFraction(s_hat)
369     def getDp(self, p=0):
370     return (self.S_max-self.S_min) * (self.m-1)/self.m*(1+(self.gamma * abs(p))**self.m)**(1./self.m) * self.gamma**self.m * abs(p)**(self.m-1)*sign(p)
371    
372     if DIM==2:
373     dom=Rectangle(NE_L,NE_H,l0=L,l1=H,order=1,optimize=True)
374     else:
375     dom=Brick(NE_L,NE_L,NE_H,l0=L,l1=L,l2=H,order=1,optimize=True)
376    
377    
378    
379     S_w=VanGenuchten(S_min=S_W_MIN, S_max=S_W_MAX, m =M, gamma=GAMMA)
380     C_c=LangmuirIsotherm(V_L=V_L, p_L=P_L, rho_a=RHO_AG )
381    
382     Rho_g=Density(p_a = P_A, rho_a= RHO_AG, p0 = 2*P_A, rho0=RHO_0G)
383     Phi=Porosity(phi0=PHI_0, alpha=ALPHA, Ks=K_S, eps_v=EPS_V, p0=P_A, eps_L=EPS_L, p_L=P_L)
384    
385     pde=LinearPDE(dom,numEquations=3, numSolutions=3)
386     pde.getSolverOptions().setVerbosityOn()
387     pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
388    
389     t=0
390     n=0
391     dt=DT
392     mkDir(OUT_DIR)
393    
394     z=dom.getX()[DIM-1]
395     z_top=sup(z)
396     z_bot=inf(z)
397    
398    
399     suction=(RHO_W-RHO_AG)*G*(z-z_top)
400     p_g=RHO_AG*G*(z-z_top)+P_A
401     m_g=( 1-S_w(suction)) * ( Rho_g(p_g) * Phi(p_g) - RHO_COAL * C_c(p_g) )
402    
403     q=pde.createCoefficient("q")
404     r=pde.createCoefficient("r")
405     # fix pressures and movable mass on the surface
406     q[0]=whereZero(z-z_top)+whereZero(z-z_bot)
407     q[1]=whereZero(z-z_top)+whereZero(z-z_bot)
408     q[2]=whereZero(z-z_top)+whereZero(z-z_bot)
409     r[0]=suction
410     r[1]=m_g
411     r[2]=p_g
412     pde.setValue(q=q, r=r)
413    
414     norm_p_g=Lsup(p_g)
415     norm_suction=Lsup(suction)
416    
417     print " m_g %e"%(Lsup(m_g))
418     print " p_g %e"%(norm_p_g)
419     print " suction %e"%(norm_suction)
420    
421     while t< T_END:
422     suction_old=suction
423     p_g_old=p_g
424     m_g_old=m_g
425     norm_dp_g = norm_p_g
426     norm_dsuction = norm_suction
427    
428     print "Time step %d at %e sec"%(n+1,t+dt)
429     niter=0
430    
431     while norm_dp_g > norm_p_g * TOL or norm_dsuction > norm_suction * TOL:
432     print " Iteration step %d:"%(niter,)
433    
434     s_w=S_w(suction)
435     DS_wDp=S_w.getDp(suction)
436    
437     rho_g=Rho_g(p_g)
438     DRho_gDp=Rho_g.getDp(p_g)
439    
440     c = C_c(p_g)
441     DC_cDp=C_c.getDp(p_g)
442    
443     phi=Phi(p_g)
444     DPhiDp=Phi.getDp(p_g)
445    
446     print " porosity range phi:",inf(phi), sup(phi)
447     print " water volume fraction :",inf(s_w), sup(s_w)
448    
449     D=pde.createCoefficient("D")
450     A=pde.createCoefficient("A")
451     X=pde.createCoefficient("X")
452     Y=pde.createCoefficient("Y")
453    
454    
455     # p_g equation:
456     Hp =phi*DRho_gDp + rho_g*DPhiDp - RHO_COAL*DC_cDp
457     Homega=DS_wDp*(rho_g*phi - RHO_COAL*c)
458     D[2,0]=-Homega
459     D[2,1]=-1.
460     D[2,2]= Hp
461     Y[2]=Hp * p_g_old - m_g_old - Homega*suction_old
462    
463     # revise :::
464     k=K_0*(phi/PHI_0)**3
465     # s_w_eff=S_w.getEffFraction(s_w)
466     s_w_eff=s_w
467     k_hat=sqrt(s_w_eff) * (1-(1-s_w_eff**(1./LAM))**LAM)**2
468     k_w=k*k_hat
469     k_g=k*(1-k_hat)
470    
471     K00=CHI*dt*k_w/ETA_W
472     K11=CHI*dt*rho_g*D_G
473     K12=CHI*dt*rho_g*k_g/ETA_G
474    
475     # suction equation:
476     Y[0]=phi*DS_wDp*p_g_old+s_w*DPhiDp*suction_old
477     X[0,:] = - dt * k_w/ETA_W * ( (1-CHI) * grad(suction+p_g) - G * RHO_W * kronecker(DIM)[DIM-1] )
478     print Y[0]
479     print X[0,:]
480     D[0,0]=phi*DS_wDp
481     D[0,2]=s_w*DPhiDp
482     for i in range(DIM):
483     A[0,i,0,i]=K00
484     A[0,i,2,i]=K00
485    
486    
487     # p_g equation:
488     Y[1]=m_g_old
489     X[1,:] = - dt * rho_g * ((1-CHI) * ( D_G * grad(m_g_old) + rho_g * k_g/ETA_G * ( grad(p_g_old) - G * rho_g * kronecker(DIM)[DIM-1]) ))
490     print X[1,:]
491     D[1,1]=1.
492     for i in range(DIM):
493     A[1,i,1,i]=K11
494     A[1,i,2,i]=K12
495    
496    
497     pde.setValue(A=A, D=D, X=X, Y=Y)
498     u=pde.getSolution()
499     m_g, m_g2=u[1], m_g
500     p_g, p_g2=u[2], p_g
501     suction, suction2=u[0], suction
502    
503     norm_p_g=Lsup(p_g)
504     norm_suction=Lsup(suction)
505     norm_dp_g = Lsup(p_g-p_g2)
506     norm_dsuction = Lsup(suction- suction2)
507     print inf(u[0])
508     print sup(u[0])
509    
510     print " m_g correction %e/%e"%(Lsup(m_g-m_g2), Lsup(m_g))
511     print " p_g correction %e/%e"%(norm_dp_g, norm_p_g)
512     print " suction correction %e/%e"%(norm_dsuction, norm_suction)
513    
514     niter+=1
515     if niter >1: 1/0
516    
517    
518     print "Time step %d completed."%n
519     n+=1
520     t+=dt
521    
522     1/0
523    
524    
525    
526     #===============================
527     def meanValue(arg):
528     return integrate(arg,Function(arg.getDomain()))/(H*L)
529    
530     #Boundary conditions:
531     # axial loading: they applied a stress inversely proportional to the acoustic emission rate. We could have the axial forcing a stress or velocity inversely proportional to dalpha/dt (only when it is positive, and with the applied forcing rate going to zero when damage accumulation rate goes to a value we can determine in a test run with constant forcing). If this is to challenging or time consuming we could have a constant axial strain rate with very short time steps (at least when alpha increases above 0.3).
532    
533     #Variables calculated and written to an output file:
534     # time
535     # differential stress (S_33-S_11)
536     # deviatoric stress (S_33 - p)
537     # Axial and transverse strain
538     # damage and damage rate
539    
540    
541    
542     # # material parameters
543     C_V = 2.e-5/(U.Mega*U.Pa)
544     C_D = 3/U.sec
545     C_1 = 1e-10/U.sec #1e-12/U.sec
546     C_2 = 0.02 #0.03
547     XI_0 = -0.56
548     LAME_0 = 29e9*U.Pa
549     MU_0 = 19e9*U.Pa
550     ALPHA_0 = 0.
551     ALPHA_pert = 0.6
552     RHO = 2800*U.kg/U.m**3
553     G = 10*U.m/U.sec**2 *0
554     # # BC parameters
555     SIGMA_N=+50.*U.Mega*U.Pa #confining -50 , +50 gives tension -> high damage levels
556     DIM=2
557     VMAX=-1.*U.m/U.sec/800000. #*2
558     DT_MAX=50000.*U.sec
559     DT=400.*U.sec/10000.
560     #T_END=30*DT
561     if DIM==2:
562     xc=[L/2,H/2]
563     else:
564     xc=[L/2,L/2,H/2]
565     WWW=min(H,L)*0.02 # for one-element with a_0=pert+ALPHA_0 should be ~element length (=L/NE_L)
566    
567     VERBOSE=True
568     DT_VIS=T_END/100 # time distane between two visulaization files
569     DN_VIS=20 # maximum counter increment between two visulaization files
570     VIS_DIR="damage-localization" # name of the directory for vis files <<<<<<<<<<<<<<<<<<<
571     DIAGNOSE_FN="diagnose_damage.csv"
572    
573     ODE_TOL=0.1
574     ODE_ITER_TOL=1.e-6
575     ODE_ITER_MAX=15
576     DEPS_MAX=0.01
577     DU_ITER_TOL=1e-4
578     DU_ITER_MAX=20
579    
580     diagnose=FileWriter(DIAGNOSE_FN,append=False)
581     #===================================
582     S=0.5*XI_0*((2.*MU_0+3.*LAME_0)/(3.-XI_0**2) + LAME_0)
583     GAMMA_M=S + sqrt(S**2+2.*MU_0*(2.*MU_0+3.*LAME_0)/(3.-XI_0**2))
584    
585     class DivergenceError(Exception):
586     pass
587    
588     def solveODE(u0, a, b, dt):
589     """
590     solves du/dt=a*exp(b*u) u(t=0)=u0 and return approximation at t=dt
591    
592     we use backwards Euler by solving u-u0 = dt * a*exp(b*u) alpha=solveODE(alpha_old, a,b, dt)
593     with newton scheme u <- u - (u-u0 - dt * a*exp(b*u)) / (1-dt*a*b*exp(b*u))
594     """
595     u=u0.copy()
596     norm_du=1.
597     norm_u=0.
598     n=0
599     while norm_du > ODE_ITER_TOL * norm_u:
600     H=-dt*a*exp(b*u)
601     du=-(u-u0+H)/(1+b*H)
602     u+=du
603     norm_du = Lsup(du)
604     norm_u = Lsup(u)
605     n+=1
606     if n>ODE_ITER_MAX:
607     print "*** ODE iteration has given up after %s steps with correction %e"%(n, norm_du)
608     raise DivergenceError,"ODE iteration failed after %s steps."%n
609     print " ODE iteration completed after %d steps with correction %e."%(n,norm_du)
610     return u
611    
612     #======================
613     t=0 # time stamp
614     n=0 # time step counter
615     dt=DT # current time step size
616     t_vis=0
617     n_vis=0
618     counter_vis=0
619     mkDir(VIS_DIR)
620     #=========================
621     #
622     # ........set up domain
623     #
624    
625     BBOX=boundingBox(dom)
626     DIM=dom.getDim()
627     x=dom.getX()
628     #
629     # initial values:
630     #
631     sigma=Tensor(0.,Function(dom))
632     eps_e=Tensor(0.,Function(dom))
633     #
634     # nucleation site - alpha perturbation:
635    
636     alpha=ALPHA_0+ALPHA_pert*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
637     #MU_0 = MU_0 + MU_0/100*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
638    
639     #saveVTK("alpha0_test-cd_0-mu0_pert.vtu",alpha=alpha) # map initial damage.vtu <<<<<<<<<<<<<<<<<<
640    
641    
642     # pde.getSolverOptions().setPreconditioner(pde.getSolverOptions().AMG)
643     pde.getSolverOptions().setTolerance(DU_ITER_TOL**2)
644    
645    
646    
647     fixed_v_mask=Vector(0,Solution(dom))
648     v0=Vector(0.,ContinuousFunction(dom))
649    
650     for d in range(DIM):
651     fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)
652     if d == DIM-1:
653     fixed_v_mask+=whereZero(x[d]-BBOX[d][1])*unitVector(d,DIM)
654     v0[d]=(x[d]-BBOX[d][0])/(BBOX[d][1]-BBOX[d][0])*VMAX
655    
656     pde.setValue(Y=-G*RHO*kronecker(DIM)[DIM-1], q=fixed_v_mask)
657     du=Vector(0.,Solution(dom))
658     u=Vector(0.,Solution(dom))
659     norm_du=0.
660     deps=Tensor(0,Function(dom))
661     i_eta=0
662     #
663     # ........let the show begin:
664     #
665     k3=kronecker(DIM)
666     k3Xk3=outer(k3,k3)
667     alpha_old=alpha
668     dt_old=None
669     diagnose.write("t, -e22, e11, s00-s22, alpha_max, alpha_av\n")
670    
671     while t<T_END:
672    
673     print "===== start time step %d ===== "%(n+1,)
674    
675     eps_e_old = eps_e
676     sigma_old = sigma
677     alpha_old, alpha_oold = alpha, alpha_old
678     du_old=du
679     converged=False
680     # start the iteration for deps on a time step: deps from the last time step is used as an initial guess:
681     while not converged:
682     iter_du=0
683     if dt_old!=None:
684     du=du_old*(dt/dt_old)
685     else:
686     du=du_old
687     norm_ddu=Lsup(du)
688     norm_u = Lsup(u)
689     deps=symmetric(grad(du))
690     sigma=sigma_old
691     # start iteration :
692     try:
693     while norm_ddu > DU_ITER_TOL * norm_u or iter_du == 0 :
694     print " start iteration step %d at time step %d:"%(iter_du,n+1)
695    
696     eps_e = eps_e_old + deps-(dt/2)*i_eta*deviatoric(sigma)
697    
698     I1=trace(eps_e)
699     sqI2=length(eps_e)
700     xi=safeDiv(I1,sqI2)
701     # ......update damage parameter:
702     m=wherePositive(xi-XI_0)
703     a=sqI2**2*(xi-XI_0)*(m*C_D + (1-m)* C_1)
704     b=(1-m)*(1./C_2)
705    
706     alpha=solveODE(alpha_old, a,b, dt)
707     if sup(alpha) > 1.:
708     print "*** damage parameter %e > 1"%(sup(alpha), )
709     raise DivergenceError,"damage parameter %e > 1"%(sup(alpha), )
710    
711     alpha_dot=(alpha-alpha_old)/dt
712     i_eta = clip(2*C_V*alpha_dot,minval=0.)
713    
714     gamma=alpha*GAMMA_M
715     lame=LAME_0
716     mu=MU_0*(1-alpha)
717    
718     print " alpha = [ %e, %e]"%(inf(alpha),sup(alpha))
719     print " xi = [ %e, %e]"%(inf(xi),sup(xi))
720     print " gamma = [ %e, %e]"%(inf(gamma),sup(gamma))
721    
722     sigma = lame * I1 * k3 + 2* mu * eps_e - gamma * ( sqI2 * k3 + xi * eps_e )
723    
724     pde.setValue(A = mu * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame*k3Xk3)
725     pde.setValue(X=-sigma, y=SIGMA_N*dom.getNormal(), r=dt*v0-du)
726    
727    
728     ddu=pde.getSolution()
729     deps+=symmetric(grad(ddu))
730     du=du+ddu
731     norm_ddu=Lsup(ddu)
732     norm_u=Lsup(u+du)
733     print " displacement change update = %e of %e"%(norm_ddu, norm_u)
734     iter_du+=1
735     if iter_du > DU_ITER_MAX:
736     print "*** displacement iteration has given up after %s steps with rel. correction %e"%(iter_du, norm_ddu/norm_u)
737     raise DivergenceError,"displacement iteration failed after %s steps."%iter_du
738     converged=True
739     print "displacement iteration converged after %d steps (rel. increment = %e)."%(iter_du, norm_ddu/norm_u)
740    
741     except DivergenceError:
742     converged=False
743     dt*=0.5
744     print "*** iteration is resumed with new time step size = %e"%(dt,)
745    
746     u+=du
747     n+=1
748     t+=dt
749     diagnose.write(("%e,"*6+"\n")%(t, meanValue(-symmetric(grad(u))[DIM-1,DIM-1]),
750     meanValue(symmetric(grad(u))[0,0]),
751     meanValue(sigma[0,0]-sigma[DIM-1,DIM-1]),
752     sup(alpha),
753     meanValue(alpha))) #meanValue(alpha_dot)))
754     print "time step %s (t=%e) completed."%(n,t)
755     #
756     # ............visualization:
757     #
758     # dom.setX(dom.getX()+du)
759     if t>=t_vis or n>n_vis:
760     saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),u=u, alpha=alpha, I1=trace(eps_e), eps_12=(symmetric(grad(u))[0,1]), xi=safeDiv(trace(eps_e),length(eps_e)), alpha_dot=alpha_dot, sqI2=sqI2) # tau=sigma[0,1], I2=length(eps_e)**2
761     print "visualization file %d for time step %e generated."%(counter_vis,t)
762     counter_vis+=1
763     t_vis+=DT_VIS
764     n_vis+=DN_VIS
765     #
766     # ............control time step size:
767     #
768     # ss=sup(length(deps))
769     # if ss>0:
770     # dt_new=DEPS_MAX/ss*dt
771     # print " time step size to control strain increment %s."%(dt_new,)
772     # else:
773     dt_new=dt
774     if dt_old != None:
775     dd_alpha=2.*dt_old*(alpha-alpha_old)+(alpha_oold-alpha_old)*dt/(dt*dt_old*(dt_old+dt))
776     norm_alpha=Lsup(alpha)
777     fac=Lsup(dd_alpha)
778     if norm_alpha > 0: fac*=1./norm_alpha
779     if fac>0:
780     error=fac*0.5*dt**2
781     print " estimated local error for time step size %e is %e"%(dt,error)
782     dt_new=sqrt(ODE_TOL*2/fac)
783     print " new time step size to maintain error level = %e."%(dt_new,)
784     if n<10:
785     fac=10.
786     else:
787     fac=1.3
788     dt_new=min(max(dt_new,dt/fac),dt*fac) # aviod rapid changes
789     print "new time step size %e"%dt_new
790     dt, dt_old=dt_new, dt
791     print "#### Time integartion completed #########"

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26