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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2440 - (hide annotations)
Wed May 27 08:45:55 2009 UTC (11 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 15094 byte(s)
convection works now again.
1 ksteube 1811 ########################################################
2     #
3     # Copyright (c) 2003-2008 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 gross 2438 """
13     this is a convection simulation over a domain [0,L] X [0,L] x [0,H]
14 ksteube 1811
15 gross 2440 It is solved in dimensionless form
16    
17 gross 2438 """
18 ksteube 1811 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
19     Earth Systems Science Computational Center (ESSCC)
20     http://www.uq.edu.au/esscc
21     Primary Business: Queensland, Australia"""
22     __license__="""Licensed under the Open Software License version 3.0
23     http://www.opensource.org/licenses/osl-3.0.php"""
24 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
25 ksteube 1811
26 gross 1417 from esys.escript import *
27 gross 2440 from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian
28 gross 1474 from esys.finley import Rectangle, Brick, LoadMesh
29 gross 1472 from optparse import OptionParser
30 gross 1417 from math import pi, ceil
31 artak 1483 import sys
32     import time
33    
34 gross 2438 # ======================= Default Values ==================================================
35     DIM=2 # spatial dimension
36     H=1. # height
37 gross 2440 L=H # length
38     NE=50 # number of elements in H-direction.
39     PERT=0.05 # initial temperature perturbation
40     DT=1.e-7 # initial time step size
41     CREATE_TOPOGRAPHY=False # create topgraphy
42     TOPOGRAPHY_LOAD=False # considers the loading of the topography in the flow
43     DT_MIN=1.e-10 # minumum time step size
44 gross 2438 T_END=0.1 # end time
45 artak 1483
46 gross 2440 RHO_0=1. # surface density
47     G=1. # gravitational constant
48     ALPHA_0=4.e3 # thermal expansion coefficient
49     T_0=0. # surface temperature
50     T_1=1. # bottom temperature
51     C_P=1 # heat capacity
52     K=1. # thermal conductivity
53     CHI=0. # Taylor Quinny coefficient
54     MUE=None # elastic shear modulus
55     TAU_Y=3*10**(2.5) # Drucker-Prager cohesion factor
56     BETA=0 # Drucker-Prager friction factor
57     TAU_0=2*10**(2.5) # transition stress
58     N=3 # power for power law
59 gross 1474
60 gross 2440 E=23 # activation energy
61     V=18*0 # activation volume
62     T_OFFSET=1 # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0)
63     R=1 # gas constant
64     ETA_N0=1. # viscosity at surface
65 gross 1472
66 gross 2440 T_TOL=1.e-4 # tolerance temperature transport
67     TOL2=0.1 # tolerance for flow update at a timestep. (large value will do only one correction step)
68     FLOW_TOL=1.e-5 # tolerance for inconcompressible flow solver
69     FLOW_SUB_TOL=1.e-10 # sub-tolerance for inconcompressible flow solver
70     NUSSELT_FN="nusselt.csv"
71     VERBOSE=True
72     DT_VIS=T_END/500 # time distane between two visulaization files
73     DN_VIS=5 # maximum counter increment between two visulaization files
74     VIS_DIR="results" # name of the director for vis files
75     DN_RESTART=1000000 # create a restart file every DN_RESTART step
76     PREFIX_RESTART="restart" # name prefix for restart directories
77 gross 2438 # =========================================================================================
78    
79     def removeRestartDirectory(dir_name):
80     if dom.onMasterProcessor() and os.path.isdir(dir_name):
81     for root, dirs, files in os.walk(dir_name, topdown=False):
82     for name in files: os.remove(os.path.join(root,name))
83     for name in dirs: os.remove(os.path.join(root,name))
84     os.rmdir(dir_name)
85     print "Restart files %s have been removed."%dir_name
86     dom.MPIBarrier()
87 gross 2440
88 gross 2438 # =========================================================================================
89 gross 1417 #
90 gross 2438 # read options:
91 gross 1417 #
92 gross 2438 parser = OptionParser(usage="%prog [Options]")
93     parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory. It will be deleted after a new directory has been created.", default=False, action="store_true")
94     parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR. The directory will not be deleted but new restart directories are created.",metavar="DIR", default=None)
95     parser.add_option("-p", "--param", dest="param", help="name of file to be imported ",metavar="PARAM", default=None)
96     (options, args) = parser.parse_args()
97     restart=options.restart or (options.restart_dir !=None)
98     #
99     # overwrite the default options:
100     #
101 gross 2440 print "<%s> Execution started."%time.asctime()
102 gross 2438 if options.param !=None:
103     exec(open(options.param,'r'))
104     print "Parameters are imported from file ",options.param
105 gross 2440
106     CREATE_TOPOGRAPHY=CREATE_TOPOGRAPHY or TOPOGRAPHY_LOAD
107    
108     print "Input Parameters:"
109     print "\tDimension DIM\t\t=\t",DIM
110     print "\tHeight H\t\t\t=\t",H
111     print "\tLength L\t\t\t=\t",L
112     print "\telements in H NE\t\t=\t",NE
113     print "\ttemperature perturbation PERT\t\t=\t",PERT
114     print "\tinitial time step size DT\t\t=\t",DT
115     print "\tminimum time step size DT_MIN\t\t=\t",DT_MIN
116     print "\tend time T_END\t\t=\t",T_END
117     print "\tcreate topgraphy CREATE_TOPOGRAPHY\t=\t",CREATE_TOPOGRAPHY
118     print "\tloading by topography TOPOGRAPHY_LOAD\t=\t",TOPOGRAPHY_LOAD
119     print "\tsurface density RHO_0\t\t=\t",RHO_0
120     print "\tgravitational constant G\t\t\t=\t",G
121     print "\tthermal expansion coefficient ALPHA_0\t\t=\t",ALPHA_0
122     print "\tsurface temperature T_0\t\t=\t",T_0
123     print "\tbottom temperature T_1\t\t=\t",T_1
124     print "\theat capacity C_P\t\t=\t",C_P
125     print "\tthermal conductivity K\t\t\t=\t",K
126     print "\tTaylor-Qinny coefficient CHI\t\t=\t",CHI
127     print "\telastic shear modulus MUE\t\t=\t",MUE
128     print "\tcohesion factor TAU_Y\t\t=\t",TAU_Y
129     print "\tfriction factor BETA\t\t=\t",BETA
130     print "\ttransition stress TAU_0\t\t=\t",TAU_0
131     print "\tpower for power law N\t\t\t=\t",N
132     print "\tviscosity at surface ETA_N0\t\t=\t",ETA_N0
133     print "\tactivation energy E\t\t\t=\t",E
134     print "\tactivation volume V\t\t\t=\t",V
135     print "\ttemperature offset T_OFFSET\t\t=\t",T_OFFSET
136     print "\tgas constant R\t\t\t=\t",R
137    
138     print "\ttransport tolerance T_TOL\t\t=\t",T_TOL
139     print "\ttolerance for flow updates TOL2\t\t=\t",TOL2
140     print "\tflow tolerance FLOW_TOL\t\t=\t",FLOW_TOL
141     print "\tflow sub-tolerance FLOW_SUB_TOL\t=\t",FLOW_SUB_TOL
142     print "\tfile for Nusselt number NUSSELT_FN\t=\t",NUSSELT_FN
143     print "\tmin. time incr. for vis file DT_VIS\t\t=\t",DT_VIS
144     print "\tmin. count incr. for vis file DN_VIS\t\t=\t",DN_VIS
145     print "\tvisualization dir VIS_DIR\t\t=\t",VIS_DIR
146     print "\trestart counter incerement DN_RESTART\t\t=\t",DN_RESTART
147     print "\tprefix for restart dirs PREFIX_RESTART\t\t=\t",PREFIX_RESTART
148     print "\tverbosity VERBOSE\t\t=\t",VERBOSE
149    
150     print "Control Parameters:"
151     t_REF=RHO_0*C_P*H**2/K
152     P_REF=ETA_N0/t_REF
153     Ar=E/R/(T_1-T_0)
154     if E>0 and V>0:
155     V_REF=P_REF*V/E
156     else:
157     V_REF=0
158     T_OFFSET_REF=T_OFFSET/(T_1-T_0)
159     Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF
160     Di=ALPHA_0*G*H/C_P
161     CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2)
162     if MUE == None:
163     De=None
164     else:
165     De=ETA_N0/MUE/t_REF
166     ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0
167     print "\ttotal #element \t\t\t=\t",NE**DIM*int(L/H)**(DIM-1)
168     print "\treference time t_REF\t\t=\t%e"%t_REF
169     print "\treference pressure P_REF\t\t=\t%e"%P_REF
170     print "\treference Taylor-Qinny CHI_REF\t\t=\t%e"%CHI_REF
171     print "\tDissipation number DI\t\t=\t%e"%Di
172     print "\tRayleigh number surface Ra\t\t=\t%e"%Ra
173     if MUE == None:
174     print "\tDebora number surface De\t\t=\t",None
175     else:
176     print "\tDebora number surface De\t\t=\t%e"%De
177    
178     print "\tbottom viscosity \t\t\t=\t%e"%ETA_BOT
179     print "\tRayleigh number bottom \t\t\t=\t%e"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT)
180     if MUE == None:
181     print "\tDebora number bottom \t\t\t=\t",None
182     else:
183     print "\tDebora number bottom \t\t\t=\t%d"%(ETA_BOT/MUE/t_REF)
184     print "\tArrhenius Ar\t\t=\t%e"%Ar
185     print "\tscaled activation volume V_REF\t\t=\t%e"%V_REF
186 gross 2438 # some control variables (will be overwritten in case of a restart:
187     #
188     t=0 # time stamp
189     n=0 # time step counter
190     dt=DT # current time step size
191 gross 2440 t_vis=0
192     n_vis=0
193     counter_vis=0
194     if getMPIRankWorld()==0:
195     if not os.path.isdir(VIS_DIR): os.mkdir(VIS_DIR)
196     #=========================
197 gross 2438 #
198     # set up domain or read restart file
199     #
200 gross 1472 if restart:
201     if options.restart_dir ==None:
202     restart_files=[]
203     for f in os.listdir("."):
204 gross 2440 if f.startswith(PREFIX_RESTART): restart_files.append(f)
205 gross 1472 if len(restart_files)==0:
206 gross 1474 raise IOError,"no restart files"
207 gross 1472 restart_files.sort()
208     f=restart_files[-1]
209     else:
210     f=options.restart_dir
211 artak 1483 try:
212     dom=LoadMesh("mesh.nc")
213     except:
214     pass
215 gross 1639 FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
216 gross 2438 t=float(FF[0]) # time stamp
217 gross 2440 n=int(FF[3]) # time step counter
218     t_vis=float(FF[1]) #
219     n_vis=int(FF[2])
220     counter_vis=int(FF[4])
221 gross 1639 dt=float(FF[5])
222     stress=load(os.path.join(f,"stress.nc"),dom)
223 gross 1474 v=load(os.path.join(f,"v.nc"),dom)
224     p=load(os.path.join(f,"p.nc"),dom)
225     T=load(os.path.join(f,"T.nc"),dom)
226 gross 2440
227     nusselt_file=FileWriter(NUSSELT_FN,append=True)
228     print "<%s> Restart from file %s at time step %s (t=%s) completed."%(time.asctime(),f,t,n)
229 gross 1417 else:
230 gross 1472 if DIM==2:
231 gross 2440 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=2, useFullElementOrder=True,optimize=True)
232 gross 1472 else:
233 gross 2440 dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L/H,l1=L/H,l2=1,order=2, useFullElementOrder=True,optimize=True)
234 artak 1483 try:
235     dom.dump("mesh.nc")
236     except:
237     pass
238 gross 1472 x=dom.getX()
239     T=Scalar(1,ReducedSolution(dom))
240     for d in range(DIM):
241     if d == DIM-1:
242     T*=sin(x[d]/H*pi)
243     else:
244 gross 1552 T*=cos(x[d]/L*pi)
245 gross 1417
246 gross 2440 T=(1.-x[DIM-1])+PERT*T
247 gross 1472 v=Vector(0,Solution(dom))
248 gross 1639 stress=Tensor(0,Function(dom))
249     x2=ReducedSolution(dom).getX()
250 gross 2440 p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)
251     nusselt_file=FileWriter(NUSSELT_FN,append=False)
252 gross 1472
253 gross 1417 x=dom.getX()
254     #
255     # set up heat problem:
256     #
257 gross 2440 heat=TemperatureCartesian(dom,useBackwardEuler=False)
258     print "<%s> Temperature transport has been set up."%time.asctime()
259     heat.setTolerance(T_TOL)
260 gross 1417 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
261 gross 2440 heat.setInitialTemperature(clip(T,T_0))
262     heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
263 gross 1417 #
264 gross 1639 # velocity constraints:
265 gross 1417 #
266     x2=ReducedSolution(dom).getX()
267     fixed_v_mask=Vector(0,Solution(dom))
268     for d in range(DIM):
269     if d == DIM-1:
270     ll = H
271     else:
272     ll = L
273 gross 2440 if TOPOGRAPHY_LOAD and d==DIM-1:
274 gross 2438 fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
275     else:
276     fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
277 gross 1417 #
278 gross 1639 # set up velovity problem
279     #
280 gross 2440 flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE)
281     flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1]))
282    
283     flow.setElasticShearModulus(MUE)
284     flow.setTolerance(TOL2)
285     flow.setFlowTolerance(FLOW_TOL)
286     flow.setFlowSubTolerance(FLOW_SUB_TOL)
287     flow.setEtaTolerance(FLOW_TOL)
288     flow.setExternals(fixed_v_mask=fixed_v_mask)
289     print "<%s> Flow solver has been set up."%time.asctime()
290     #
291 gross 1417 # let the show begin:
292     #
293 gross 2438 t1 = time.time()
294 gross 2440 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
295 gross 1417 while t<T_END:
296 gross 2440 #======= solve for velovity ====================================================================
297     eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
298     print "viscosity range :", inf(eta_N), sup(eta_N)
299     flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0], [1,N])
300     flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM)) #f= add the mountanins here
301     flow.update(dt, iter_max=100, inner_iter_max=200, verbose=False)
302     v=flow.getVelocity()
303 gross 1417 for d in range(DIM):
304     print "range %d-velocity"%d,inf(v[d]),sup(v[d])
305 gross 2440 print "<%s> flow solver completed."%time.asctime()
306     if t>=t_vis or n>n_vis:
307     saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())
308     print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)
309     counter_vis+=1
310     t_vis+=DT_VIS
311     n_vis+=DN_VIS
312 gross 1417 n+=1
313     t+=dt
314 gross 2440 print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)
315     #======= set Temperature problem ====================================================================
316     #
317     heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
318     dt=heat.getSafeTimeStepSize()
319     print "<%s> New time step size is %e"%(time.asctime(),dt)
320     print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)
321     #
322     # solve temperature:
323     #
324     T=heat.getSolution(dt, verbose=VERBOSE)
325     print "Temperature range ",inf(T),sup(T)
326     print "<%s> temperature update completed."%time.asctime()
327     #======= anaysis ==================================================================================
328     #
329     # .... Nusselt number
330     #
331     dTdz=grad(T)[DIM-1]
332     Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
333     print "nusselt number = ",Nu
334     nusselt_file.write("%e %e\n"%(t,Nu))
335     eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
336     Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
337     print "av. eta = %e"%eta_bar
338     print "effective Rayleigh number = %e"%Ra_eff
339 gross 1474 # =========================
340     #
341     # create restart files:
342     #
343 gross 2440 if DN_RESTART>0:
344     if (n-1)%DN_RESTART == 0:
345     c=(n-1)/DN_RESTART
346     old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)
347     new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)
348 gross 1474
349 ksteube 1877 if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
350     dom.MPIBarrier()
351 gross 2440 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_vis, n_vis, n, counter_vis, dt))
352 gross 2100 v.dump(os.path.join(new_restart_dir,"v.nc"))
353     p.dump(os.path.join(new_restart_dir,"p.nc"))
354 gross 1474 T.dump(os.path.join(new_restart_dir,"T.nc"))
355     removeRestartDirectory(old_restart_dir)
356 gross 2440 print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)
357     print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)
358 gross 1474

  ViewVC Help
Powered by ViewVC 1.1.26