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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2445 - (show annotations)
Fri May 29 03:23:25 2009 UTC (11 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 15102 byte(s)
some improvement on convection which ios now working again
1 ########################################################
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 """
13 this is a convection simulation over a domain [0,L] X [0,L] x [0,H]
14
15 It is solved in dimensionless form
16
17 """
18 __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 __url__="https://launchpad.net/escript-finley"
25
26 from esys.escript import *
27 from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian
28 from esys.finley import Rectangle, Brick, LoadMesh
29 from optparse import OptionParser
30 from math import pi, ceil
31 import sys
32 import time
33
34 # ======================= Default Values ==================================================
35 DIM=2 # spatial dimension
36 H=1. # height
37 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 T_END=0.1 # end time
45
46 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=5*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
60 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
66 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-4 # tolerance for inconcompressible flow solver
69 FLOW_SUB_TOL=1.e-8 # 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 # =========================================================================================
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
88 # =========================================================================================
89 #
90 # read options:
91 #
92 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 print "<%s> Execution started."%time.asctime()
102 if options.param !=None:
103 exec(open(options.param,'r'))
104 print "Parameters are imported from file ",options.param
105
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 # 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 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 #
198 # set up domain or read restart file
199 #
200 if restart:
201 if options.restart_dir ==None:
202 restart_files=[]
203 for f in os.listdir("."):
204 if f.startswith(PREFIX_RESTART): restart_files.append(f)
205 if len(restart_files)==0:
206 raise IOError,"no restart files"
207 restart_files.sort()
208 f=restart_files[-1]
209 else:
210 f=options.restart_dir
211 try:
212 dom=LoadMesh("mesh.nc")
213 except:
214 pass
215 FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
216 t=float(FF[0]) # time stamp
217 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 dt=float(FF[5])
222 stress=load(os.path.join(f,"stress.nc"),dom)
223 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
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 else:
230 if DIM==2:
231 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=2, useFullElementOrder=True,optimize=True)
232 else:
233 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 try:
235 dom.dump("mesh.nc")
236 except:
237 pass
238 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 T*=cos(x[d]/L*pi)
245
246 T=(1.-x[DIM-1])+PERT*T
247 v=Vector(0,Solution(dom))
248 stress=Tensor(0,Function(dom))
249 x2=ReducedSolution(dom).getX()
250 p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)
251 nusselt_file=FileWriter(NUSSELT_FN,append=False)
252
253 p_last=p
254 x=dom.getX()
255 #
256 # set up heat problem:
257 #
258 heat=TemperatureCartesian(dom,useBackwardEuler=False)
259 print "<%s> Temperature transport has been set up."%time.asctime()
260 heat.setTolerance(T_TOL)
261 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
262 heat.setInitialTemperature(clip(T,T_0))
263 heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
264 #
265 # velocity constraints:
266 #
267 x2=ReducedSolution(dom).getX()
268 fixed_v_mask=Vector(0,Solution(dom))
269 for d in range(DIM):
270 if d == DIM-1:
271 ll = H
272 else:
273 ll = L
274 if TOPOGRAPHY_LOAD and d==DIM-1:
275 fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
276 else:
277 fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
278 #
279 # set up velovity problem
280 #
281 flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE)
282 flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1]))
283
284 flow.setElasticShearModulus(MUE)
285 flow.setTolerance(TOL2)
286 flow.setFlowTolerance(FLOW_TOL)
287 flow.setFlowSubTolerance(FLOW_SUB_TOL)
288 flow.setEtaTolerance(FLOW_TOL)
289 flow.setExternals(fixed_v_mask=fixed_v_mask)
290 print "<%s> Flow solver has been set up."%time.asctime()
291 #
292 # let the show begin:
293 #
294 t1 = time.time()
295 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
296 while t<T_END:
297 #======= solve for velovity ====================================================================
298 eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
299 print "viscosity range :", inf(eta_N), sup(eta_N)
300 flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0], [1,N])
301 flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM)) #f= add the mountanins here
302 flow.update(dt, iter_max=100, inner_iter_max=200, verbose=False)
303 v=flow.getVelocity()
304 for d in range(DIM):
305 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
306 print "<%s> flow solver completed."%time.asctime()
307 if t>=t_vis or n>n_vis:
308 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())
309 print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)
310 counter_vis+=1
311 t_vis+=DT_VIS
312 n_vis+=DN_VIS
313 n+=1
314 t+=dt
315 print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)
316 #======= set Temperature problem ====================================================================
317 #
318 heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
319 dt=heat.getSafeTimeStepSize()
320 print "<%s> New time step size is %e"%(time.asctime(),dt)
321 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)
322 #
323 # solve temperature:
324 #
325 T=heat.getSolution(dt, verbose=VERBOSE)
326 print "Temperature range ",inf(T),sup(T)
327 print "<%s> temperature update completed."%time.asctime()
328 #======= anaysis ==================================================================================
329 #
330 # .... Nusselt number
331 #
332 dTdz=grad(T)[DIM-1]
333 Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
334 print "nusselt number = ",Nu
335 nusselt_file.write("%e %e\n"%(t,Nu))
336 eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
337 Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
338 print "av. eta = %e"%eta_bar
339 print "effective Rayleigh number = %e"%Ra_eff
340 # =========================
341 #
342 # create restart files:
343 #
344 if DN_RESTART>0:
345 if (n-1)%DN_RESTART == 0:
346 c=(n-1)/DN_RESTART
347 old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)
348 new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)
349
350 if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
351 dom.MPIBarrier()
352 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))
353 v.dump(os.path.join(new_restart_dir,"v.nc"))
354 p.dump(os.path.join(new_restart_dir,"p.nc"))
355 T.dump(os.path.join(new_restart_dir,"T.nc"))
356 removeRestartDirectory(old_restart_dir)
357 print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)
358 print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)
359

  ViewVC Help
Powered by ViewVC 1.1.26