/[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 3458 - (hide annotations)
Mon Jan 31 07:06:42 2011 UTC (8 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 14195 byte(s)
more on AMG under MPI
1 gross 3458 #######################################################
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 import *
24     from esys.escript.linearPDEs import LinearPDESystem
25     from esys.finley import Rectangle, Brick
26     from esys.escript import unitsSI as U
27     from esys.weipa import saveVTK
28     from math import pi, ceil
29     import sys
30     import time
31    
32     # ======================= Default Values ==================================================
33     H=2.*U.m # height
34     L=H # length
35     NE_H=5 # number of elements in H-direction.
36     NE_L=int(ceil(L*NE_H/H))
37     DIM=2
38     T_END=40000.0*U.sec # end time at max strain = 0.5
39    
40     CHI=0.5 # =0 = explicit Euler, 0.5 = Crank-Nicholson, 1. = implicit Euler
41     G=-9.81*U.m/U.sec**2
42    
43    
44     class Porosity(object):
45     def __init__(self, phi0=1., alpha=0., Ks=1., eps_v=0, p0=0., eps_L=0, p_L=0):
46     self.phi0=phi0
47     self.alpha=alpha
48     self.eps_v=eps_v
49     self.p0=p0
50     self.p_L=p_L
51     self.eps_L = eps_L
52     self.Ks = Ks
53    
54     def getValue(self,p=0):
55     return self.phi0*+ self.alpha*(self.eps_v + (p-self.p0)/Ks + self.eps_L * self.p_L * (self.p0-p)/((self.p_L + self.p_0) * (p + self.p_L)))
56    
57     def __call__(self, p=0):
58     return self.getValue(p)
59    
60     def getDp(self, p=0):
61     return self.alpha * (1./self.Ks - self.eps_L * self.p_L/(p-self.p_L)**2 )
62    
63     class Density(object):
64     def __init__(self, p_a =0, rho_a=0, p0 =0, rho0=0):
65     self.p_a=p_a
66     self.rho_a=rho_a
67     self.p0=p0
68     self.rho0=rho0
69     def getValue(self,p=0):
70     return (self.rho_a - self.rho0)/(self.p_a - self.p0)*(p-self.p0) + self.rho0
71     def __call__(self, p=0):
72     return self.getValue(p)
73    
74     def getDp(self, p=0):
75     return (self.rho_a - self.rho0)/(self.p_a - self.p0)
76    
77    
78     class GasStorageCapacity(object):
79     def __init__(self, V_L=0, p_L=0.):
80     self.V_L=V_L
81     self.p_L=p_L
82    
83     def getValue(self,p=0):
84     return self.V_L*p/(self.p_L +p)
85     def __call__(self, p=0):
86     return self.getValue(p)
87    
88     def getDp(self, p=0):
89     return self.V_L*p/(self.p_L +p)**2
90    
91     class WaterVolumeFraction(object):
92     def __init__(self, S_res=0.1, m =3, p_d=1.)
93     self.S_res = S_res
94     self.m = m
95     self.p_d = p_d
96    
97     def getValue(self,p=0):
98     s_hat = (1+(abs(p)/self.p_d)**self.m)**(1./self.m-1)
99     return self.S_res+s_hat*(1-self.S_res)
100    
101     def __call__(self, p=0):
102     return self.getValue(p)
103    
104     def getDp(self, p=0):
105     tmp = (abs(p)/self.p_d)**self.m
106     return (1-self.S_res) * (self.m-1)/self.m*((1+tmp)**(1./self.m) * tmp*1/p
107    
108     if DIM==2:
109     dom=Rectangle(NE_L,NE_H,l0=L,l1=H,order=1,optimize=True)
110     else:
111     dom=Brick(NE_L,NE_L,NE_H,l0=L,l1=L,l2=H,order=1,optimize=True)
112    
113    
114    
115     Rho_g=Density(p_a = PA, rho_a= RHO_AG, p0 =0, rho0=0)
116     S_w=WaterVolumeFraction(S_res=S_RES, m =M, p_d=P_D)
117     C_c=GasStorageCapacity(V_L=V_L, p_L=P_L)
118    
119    
120     pde=LinearPDE(dom,numEquations=3)
121    
122     t=0
123     n=0
124    
125     z0=sup(dom.getX()[DIM-1])
126     p_w=RHO_AW*G*(z-z0)+P_A
127     p_g=RHO_AG*G*(z-z0)+P_A
128     omega=p_w-p_g
129     m_g=(1-s_w) * ( rho_g * phi + RHO_AG * RHO_COAL * c)
130    
131    
132     while :
133     s_w=S_w(omega)
134     rho_g=Rho_g(p_g)
135     c = C_c(p_g)
136    
137    
138     D=pde.getNewCoefficient("D")
139     A=pde.getNewCoefficient("A")
140     X=pde.getNewCoefficient("X")
141     Y=pde.getNewCoefficient("Y")
142    
143     D[0,0]=phi*DS_w*Dp
144     D[0,2]=s_w*DPhiDp
145    
146     D[1,1]=1.
147    
148     Hp =phi*DRho_gDp+rho_g*DPhiDp+RHO_AG*RHO_COAL*DC_cDp
149     Homega=Ds_wDo*(rho_g*phi+RHO_AG*RHO_COAL*c)
150     D[2,0]=-Homega
151     D[2,1]=-1.
152     D[2,2]= Hp
153    
154    
155    
156    
157     K00=chi*dt*k_w/ETA_W
158     K11=chi*dt*rho_g*D_G
159     K12=chi*dt*rho_g*k_g/ETA_G
160     for i in range(DIM):
161     A[0,i,0,i]=K00
162     A[0,i,2,i]=K00
163     A[1,i,1,i]=K11
164     A[1,i,2,i]=K12
165    
166    
167     Y[0]=phi*DS_w*Dp*p_g_old+s_w*DPhiDp*omega_old
168     Y[1]=m_g_old
169     Y[2]=h_p * p_g_old - m_g_old - Homega*omega
170    
171    
172     X[0,:] = dt * k_w/ETA_W * ( (1-chi) * (grad(omega_old) + grad(omega_g_old) ) - g * RHO_W * kronecker(DIM)[DIM-1] )
173     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 * k_g/eta_g * kronecker(DIM)[DIM-1]) )
174    
175     #===============================
176     def meanValue(arg):
177     return integrate(arg,Function(arg.getDomain()))/(H*L)
178    
179     #Boundary conditions:
180     # 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).
181    
182     #Variables calculated and written to an output file:
183     # time
184     # differential stress (S_33-S_11)
185     # deviatoric stress (S_33 - p)
186     # Axial and transverse strain
187     # damage and damage rate
188    
189    
190    
191     # # material parameters
192     C_V = 2.e-5/(U.Mega*U.Pa)
193     C_D = 3/U.sec
194     C_1 = 1e-10/U.sec #1e-12/U.sec
195     C_2 = 0.02 #0.03
196     XI_0 = -0.56
197     LAME_0 = 29e9*U.Pa
198     MU_0 = 19e9*U.Pa
199     ALPHA_0 = 0.
200     ALPHA_pert = 0.6
201     RHO = 2800*U.kg/U.m**3
202     G = 10*U.m/U.sec**2 *0
203     # # BC parameters
204     SIGMA_N=+50.*U.Mega*U.Pa #confining -50 , +50 gives tension -> high damage levels
205     DIM=2
206     VMAX=-1.*U.m/U.sec/800000. #*2
207     DT_MAX=50000.*U.sec
208     DT=400.*U.sec/10000.
209     #T_END=30*DT
210     if DIM==2:
211     xc=[L/2,H/2]
212     else:
213     xc=[L/2,L/2,H/2]
214     WWW=min(H,L)*0.02 # for one-element with a_0=pert+ALPHA_0 should be ~element length (=L/NE_L)
215    
216     VERBOSE=True
217     DT_VIS=T_END/100 # time distane between two visulaization files
218     DN_VIS=20 # maximum counter increment between two visulaization files
219     VIS_DIR="damage-localization" # name of the directory for vis files <<<<<<<<<<<<<<<<<<<
220     DIAGNOSE_FN="diagnose_damage.csv"
221    
222     ODE_TOL=0.1
223     ODE_ITER_TOL=1.e-6
224     ODE_ITER_MAX=15
225     DEPS_MAX=0.01
226     DU_ITER_TOL=1e-4
227     DU_ITER_MAX=20
228    
229     diagnose=FileWriter(DIAGNOSE_FN,append=False)
230     #===================================
231     S=0.5*XI_0*((2.*MU_0+3.*LAME_0)/(3.-XI_0**2) + LAME_0)
232     GAMMA_M=S + sqrt(S**2+2.*MU_0*(2.*MU_0+3.*LAME_0)/(3.-XI_0**2))
233    
234     class DivergenceError(Exception):
235     pass
236    
237     def solveODE(u0, a, b, dt):
238     """
239     solves du/dt=a*exp(b*u) u(t=0)=u0 and return approximation at t=dt
240    
241     we use backwards Euler by solving u-u0 = dt * a*exp(b*u) alpha=solveODE(alpha_old, a,b, dt)
242     with newton scheme u <- u - (u-u0 - dt * a*exp(b*u)) / (1-dt*a*b*exp(b*u))
243     """
244     u=u0.copy()
245     norm_du=1.
246     norm_u=0.
247     n=0
248     while norm_du > ODE_ITER_TOL * norm_u:
249     H=-dt*a*exp(b*u)
250     du=-(u-u0+H)/(1+b*H)
251     u+=du
252     norm_du = Lsup(du)
253     norm_u = Lsup(u)
254     n+=1
255     if n>ODE_ITER_MAX:
256     print "*** ODE iteration has given up after %s steps with correction %e"%(n, norm_du)
257     raise DivergenceError,"ODE iteration failed after %s steps."%n
258     print " ODE iteration completed after %d steps with correction %e."%(n,norm_du)
259     return u
260    
261     #======================
262     t=0 # time stamp
263     n=0 # time step counter
264     dt=DT # current time step size
265     t_vis=0
266     n_vis=0
267     counter_vis=0
268     mkDir(VIS_DIR)
269     #=========================
270     #
271     # ........set up domain
272     #
273    
274     BBOX=boundingBox(dom)
275     DIM=dom.getDim()
276     x=dom.getX()
277     #
278     # initial values:
279     #
280     sigma=Tensor(0.,Function(dom))
281     eps_e=Tensor(0.,Function(dom))
282     #
283     # nucleation site - alpha perturbation:
284    
285     alpha=ALPHA_0+ALPHA_pert*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
286     #MU_0 = MU_0 + MU_0/100*exp(-length(Function(dom).getX()-xc)**2/WWW**2)
287    
288     #saveVTK("alpha0_test-cd_0-mu0_pert.vtu",alpha=alpha) # map initial damage.vtu <<<<<<<<<<<<<<<<<<
289    
290    
291     # pde.getSolverOptions().setPreconditioner(pde.getSolverOptions().AMG)
292     pde.getSolverOptions().setTolerance(DU_ITER_TOL**2)
293    
294    
295    
296     fixed_v_mask=Vector(0,Solution(dom))
297     v0=Vector(0.,ContinuousFunction(dom))
298    
299     for d in range(DIM):
300     fixed_v_mask+=whereZero(x[d]-BBOX[d][0])*unitVector(d,DIM)
301     if d == DIM-1:
302     fixed_v_mask+=whereZero(x[d]-BBOX[d][1])*unitVector(d,DIM)
303     v0[d]=(x[d]-BBOX[d][0])/(BBOX[d][1]-BBOX[d][0])*VMAX
304    
305     pde.setValue(Y=-G*RHO*kronecker(DIM)[DIM-1], q=fixed_v_mask)
306     du=Vector(0.,Solution(dom))
307     u=Vector(0.,Solution(dom))
308     norm_du=0.
309     deps=Tensor(0,Function(dom))
310     i_eta=0
311     #
312     # ........let the show begin:
313     #
314     k3=kronecker(DIM)
315     k3Xk3=outer(k3,k3)
316     alpha_old=alpha
317     dt_old=None
318     diagnose.write("t, -e22, e11, s00-s22, alpha_max, alpha_av\n")
319    
320     while t<T_END:
321    
322     print "===== start time step %d ===== "%(n+1,)
323    
324     eps_e_old = eps_e
325     sigma_old = sigma
326     alpha_old, alpha_oold = alpha, alpha_old
327     du_old=du
328     converged=False
329     # start the iteration for deps on a time step: deps from the last time step is used as an initial guess:
330     while not converged:
331     iter_du=0
332     if dt_old!=None:
333     du=du_old*(dt/dt_old)
334     else:
335     du=du_old
336     norm_ddu=Lsup(du)
337     norm_u = Lsup(u)
338     deps=symmetric(grad(du))
339     sigma=sigma_old
340     # start iteration :
341     try:
342     while norm_ddu > DU_ITER_TOL * norm_u or iter_du == 0 :
343     print " start iteration step %d at time step %d:"%(iter_du,n+1)
344    
345     eps_e = eps_e_old + deps-(dt/2)*i_eta*deviatoric(sigma)
346    
347     I1=trace(eps_e)
348     sqI2=length(eps_e)
349     xi=safeDiv(I1,sqI2)
350     # ......update damage parameter:
351     m=wherePositive(xi-XI_0)
352     a=sqI2**2*(xi-XI_0)*(m*C_D + (1-m)* C_1)
353     b=(1-m)*(1./C_2)
354    
355     alpha=solveODE(alpha_old, a,b, dt)
356     if sup(alpha) > 1.:
357     print "*** damage parameter %e > 1"%(sup(alpha), )
358     raise DivergenceError,"damage parameter %e > 1"%(sup(alpha), )
359    
360     alpha_dot=(alpha-alpha_old)/dt
361     i_eta = clip(2*C_V*alpha_dot,minval=0.)
362    
363     gamma=alpha*GAMMA_M
364     lame=LAME_0
365     mu=MU_0*(1-alpha)
366    
367     print " alpha = [ %e, %e]"%(inf(alpha),sup(alpha))
368     print " xi = [ %e, %e]"%(inf(xi),sup(xi))
369     print " gamma = [ %e, %e]"%(inf(gamma),sup(gamma))
370    
371     sigma = lame * I1 * k3 + 2* mu * eps_e - gamma * ( sqI2 * k3 + xi * eps_e )
372    
373     pde.setValue(A = mu * ( swap_axes(k3Xk3,0,3)+swap_axes(k3Xk3,1,3) ) + lame*k3Xk3)
374     pde.setValue(X=-sigma, y=SIGMA_N*dom.getNormal(), r=dt*v0-du)
375    
376    
377     ddu=pde.getSolution()
378     deps+=symmetric(grad(ddu))
379     du=du+ddu
380     norm_ddu=Lsup(ddu)
381     norm_u=Lsup(u+du)
382     print " displacement change update = %e of %e"%(norm_ddu, norm_u)
383     iter_du+=1
384     if iter_du > DU_ITER_MAX:
385     print "*** displacement iteration has given up after %s steps with rel. correction %e"%(iter_du, norm_ddu/norm_u)
386     raise DivergenceError,"displacement iteration failed after %s steps."%iter_du
387     converged=True
388     print "displacement iteration converged after %d steps (rel. increment = %e)."%(iter_du, norm_ddu/norm_u)
389    
390     except DivergenceError:
391     converged=False
392     dt*=0.5
393     print "*** iteration is resumed with new time step size = %e"%(dt,)
394    
395     u+=du
396     n+=1
397     t+=dt
398     diagnose.write(("%e,"*6+"\n")%(t, meanValue(-symmetric(grad(u))[DIM-1,DIM-1]),
399     meanValue(symmetric(grad(u))[0,0]),
400     meanValue(sigma[0,0]-sigma[DIM-1,DIM-1]),
401     sup(alpha),
402     meanValue(alpha))) #meanValue(alpha_dot)))
403     print "time step %s (t=%e) completed."%(n,t)
404     #
405     # ............visualization:
406     #
407     # dom.setX(dom.getX()+du)
408     if t>=t_vis or n>n_vis:
409     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
410     print "visualization file %d for time step %e generated."%(counter_vis,t)
411     counter_vis+=1
412     t_vis+=DT_VIS
413     n_vis+=DN_VIS
414     #
415     # ............control time step size:
416     #
417     # ss=sup(length(deps))
418     # if ss>0:
419     # dt_new=DEPS_MAX/ss*dt
420     # print " time step size to control strain increment %s."%(dt_new,)
421     # else:
422     dt_new=dt
423     if dt_old != None:
424     dd_alpha=2.*dt_old*(alpha-alpha_old)+(alpha_oold-alpha_old)*dt/(dt*dt_old*(dt_old+dt))
425     norm_alpha=Lsup(alpha)
426     fac=Lsup(dd_alpha)
427     if norm_alpha > 0: fac*=1./norm_alpha
428     if fac>0:
429     error=fac*0.5*dt**2
430     print " estimated local error for time step size %e is %e"%(dt,error)
431     dt_new=sqrt(ODE_TOL*2/fac)
432     print " new time step size to maintain error level = %e."%(dt_new,)
433     if n<10:
434     fac=10.
435     else:
436     fac=1.3
437     dt_new=min(max(dt_new,dt/fac),dt*fac) # aviod rapid changes
438     print "new time step size %e"%dt_new
439     dt, dt_old=dt_new, dt
440     print "#### Time integartion completed #########"

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26