/[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 2440 - (show annotations)
Wed May 27 08:45:55 2009 UTC (11 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 15094 byte(s)
convection works now 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=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
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-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 # =========================================================================================
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 x=dom.getX()
254 #
255 # set up heat problem:
256 #
257 heat=TemperatureCartesian(dom,useBackwardEuler=False)
258 print "<%s> Temperature transport has been set up."%time.asctime()
259 heat.setTolerance(T_TOL)
260 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
261 heat.setInitialTemperature(clip(T,T_0))
262 heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
263 #
264 # velocity constraints:
265 #
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 if TOPOGRAPHY_LOAD and d==DIM-1:
274 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 #
278 # set up velovity problem
279 #
280 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 # let the show begin:
292 #
293 t1 = time.time()
294 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
295 while t<T_END:
296 #======= 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 for d in range(DIM):
304 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
305 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 n+=1
313 t+=dt
314 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 # =========================
340 #
341 # create restart files:
342 #
343 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
349 if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
350 dom.MPIBarrier()
351 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 v.dump(os.path.join(new_restart_dir,"v.nc"))
353 p.dump(os.path.join(new_restart_dir,"p.nc"))
354 T.dump(os.path.join(new_restart_dir,"T.nc"))
355 removeRestartDirectory(old_restart_dir)
356 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

  ViewVC Help
Powered by ViewVC 1.1.26