/[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 2862 - (show annotations)
Thu Jan 21 04:45:39 2010 UTC (9 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 18133 byte(s)
more review in the cookbook
1 ########################################################
2 #
3 # Copyright (c) 2003-2009 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-2009 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, Mountains, SubSteppingException
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=2*H # length
38 NE=30 # number of elements in H-direction.
39 PERT=0.15 # initial temperature perturbation
40 DT=1.e-4 # initial time step size
41 CREATE_TOPOGRAPHY=True # create topgraphy
42 DT_MIN=1.e-10 # minumum time step size
43 T_END=10. # end time
44
45 RHO_0=100. # surface density (lf ~ RHO_0**2)
46 G=1. # gravitational constant
47 ALPHA_0=0.1 # thermal expansion coefficient (Ra ~ RHO_0**2 * ALPHA_0 = lf * ALPHA_0)
48 T_0=0. # surface temperature
49 T_1=1. # bottom temperature
50 C_P=1 # heat capacity
51 K=1. # thermal conductivity
52 CHI=0. # Taylor Quinny coefficient
53 MUE=None # elastic shear modulus
54 TAU_Y=5*10**(2.5) # Drucker-Prager cohesion factor
55 BETA=0 # Drucker-Prager friction factor
56 TAU_0=2*10**(2.5) # transition stress
57 N=3 # power for power law
58
59 E=23*0 # activation energy
60 V=18*0 # activation volume
61 T_OFFSET=1 # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0)
62 R=1 # gas constant
63 ETA_N0=1. # viscosity at surface
64
65 TOPO_SMOOTH=1e-5 # smoothing factor of extrapolation of surface velocity to interior
66
67 T_TOL=1.e-4 # tolerance temperature transport
68 FLOW_TOL=1.e-3 # tolerance for inconcompressible flow solver
69 TOPO_TOL=0.1 # tolerance for update of topography
70 DIAGNOSTICS_FN="diagnostics.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 TOPO_ITER_MAX=20 # maximum number of iteration steps to update topography
78 # =========================================================================================
79
80 def removeRestartDirectory(dir_name):
81 if dom.onMasterProcessor() and os.path.isdir(dir_name):
82 for root, dirs, files in os.walk(dir_name, topdown=False):
83 for name in files: os.remove(os.path.join(root,name))
84 for name in dirs: os.remove(os.path.join(root,name))
85 os.rmdir(dir_name)
86 print "Restart files %s have been removed."%dir_name
87 dom.MPIBarrier()
88
89 # =========================================================================================
90 #
91 # read options:
92 #
93 parser = OptionParser(usage="%prog [Options]")
94 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")
95 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)
96 parser.add_option("-p", "--param", dest="param", help="name of file to be imported ",metavar="PARAM", default=None)
97 (options, args) = parser.parse_args()
98 restart=options.restart or (options.restart_dir !=None)
99 #
100 # overwrite the default options:
101 #
102 print "<%s> Execution started."%time.asctime()
103 if options.param !=None:
104 exec(open(options.param,'r'))
105 print "Parameters are imported from file ",options.param
106
107 print "Input Parameters:"
108 print "\tDimension DIM\t\t=\t",DIM
109 print "\tHeight H\t\t\t=\t",H
110 print "\tLength L\t\t\t=\t",L
111 print "\telements in H NE\t\t=\t",NE
112 print "\ttemperature perturbation PERT\t\t=\t",PERT
113 print "\tinitial time step size DT\t\t=\t",DT
114 print "\tminimum time step size DT_MIN\t\t=\t",DT_MIN
115 print "\tend time T_END\t\t=\t",T_END
116 print "\tcreate topgraphy CREATE_TOPOGRAPHY\t=\t",CREATE_TOPOGRAPHY
117 print "\tsurface density RHO_0\t\t=\t",RHO_0
118 print "\tgravitational constant G\t\t\t=\t",G
119 print "\tthermal expansion coefficient ALPHA_0\t\t=\t",ALPHA_0
120 print "\tsurface temperature T_0\t\t=\t",T_0
121 print "\tbottom temperature T_1\t\t=\t",T_1
122 print "\theat capacity C_P\t\t=\t",C_P
123 print "\tthermal conductivity K\t\t\t=\t",K
124 print "\tTaylor-Qinny coefficient CHI\t\t=\t",CHI
125 print "\telastic shear modulus MUE\t\t=\t",MUE
126 print "\tcohesion factor TAU_Y\t\t=\t",TAU_Y
127 print "\tfriction factor BETA\t\t=\t",BETA
128 print "\ttransition stress TAU_0\t\t=\t",TAU_0
129 print "\tpower for power law N\t\t\t=\t",N
130 print "\tviscosity at surface ETA_N0\t\t=\t",ETA_N0
131 print "\tactivation energy E\t\t\t=\t",E
132 print "\tactivation volume V\t\t\t=\t",V
133 print "\ttemperature offset T_OFFSET\t\t=\t",T_OFFSET
134 print "\tgas constant R\t\t\t=\t",R
135 print "\ttopography smoothing TOPO_SMOOTH\t=\t",TOPO_SMOOTH
136
137 print "\ttolerance for topography TOPO_TOL\t\t=\t",TOPO_TOL
138 print "\ttransport tolerance T_TOL\t\t=\t",T_TOL
139 print "\tflow tolerance FLOW_TOL\t\t=\t",FLOW_TOL
140 print "\tfile for diagnostics DIAGNOSTICS_FN\t=\t",DIAGNOSTICS_FN
141 print "\tmin. time incr. for vis file DT_VIS\t\t=\t",DT_VIS
142 print "\tmin. count incr. for vis file DN_VIS\t\t=\t",DN_VIS
143 print "\tvisualization dir VIS_DIR\t\t=\t",VIS_DIR
144 print "\trestart counter incerement DN_RESTART\t=\t",DN_RESTART
145 print "\tprefix for restart dirs PREFIX_RESTART\t=\t",PREFIX_RESTART
146 print "\tverbosity VERBOSE\t\t=\t",VERBOSE
147
148 print "Control Parameters:"
149 t_REF=RHO_0*C_P*H**2/K
150 P_REF=ETA_N0/t_REF
151 Ar=E/R/(T_1-T_0)
152 if E>0 and V>0:
153 V_REF=P_REF*V/E
154 else:
155 V_REF=0
156 T_OFFSET_REF=T_OFFSET/(T_1-T_0)
157 Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF
158 Di=ALPHA_0*G*H/C_P
159 CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2)
160 if CREATE_TOPOGRAPHY:
161 SURFACE_LOAD=RHO_0*G*H/P_REF
162 else:
163 SURFACE_LOAD=0.
164 if MUE == None:
165 De=None
166 else:
167 De=ETA_N0/MUE/t_REF
168 ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0
169 print "\ttotal #element \t\t\t=\t",NE**DIM*int(L/H)**(DIM-1)
170 print "\treference time t_REF\t\t=\t%e"%t_REF
171 print "\treference pressure P_REF\t\t=\t%e"%P_REF
172 print "\treference Taylor-Qinny CHI_REF\t\t=\t%e"%CHI_REF
173 print "\tDissipation number DI\t\t=\t%e"%Di
174 print "\tRayleigh number surface Ra\t\t=\t%e"%Ra
175 if MUE == None:
176 print "\tDebora number surface De\t\t=\t",None
177 else:
178 print "\tDebora number surface De\t\t=\t%e"%De
179
180 print "\tBottom viscosity \t\t\t=\t%e"%ETA_BOT
181 print "\tRayleigh number bottom \t\t\t=\t%e"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT)
182 if MUE == None:
183 print "\tDebora number bottom \t\t\t=\t",None
184 else:
185 print "\tDebora number bottom \t\t\t=\t%d"%(ETA_BOT/MUE/t_REF)
186 print "\tArrhenius Ar\t\t=\t%e"%Ar
187 print "\tsurface load factor SURFACE_LOAD\t=\t%e"%SURFACE_LOAD
188 print "\tscaled activation volume V_REF\t\t=\t%e"%V_REF
189 # some control variables (will be overwritten in case of a restart:
190 #
191 t=0 # time stamp
192 n=0 # time step counter
193 dt=DT # current time step size
194 t_vis=0
195 n_vis=0
196 counter_vis=0
197 mkDir(VIS_DIR)
198 #=========================
199 #
200 # set up domain or read restart file
201 #
202 if restart:
203 if options.restart_dir ==None:
204 restart_files=[]
205 for f in os.listdir("."):
206 if f.startswith(PREFIX_RESTART): restart_files.append(f)
207 if len(restart_files)==0:
208 raise IOError,"no restart files"
209 restart_files.sort()
210 f=restart_files[-1]
211 else:
212 f=options.restart_dir
213 try:
214 dom=LoadMesh("mesh.nc")
215 except:
216 pass
217 FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
218 t=float(FF[0]) # time stamp
219 n=int(FF[3]) # time step counter
220 t_vis=float(FF[1]) #
221 n_vis=int(FF[2])
222 counter_vis=int(FF[4])
223 dt=float(FF[5])
224 stress=load(os.path.join(f,"stress.nc"),dom)
225 v=load(os.path.join(f,"v.nc"),dom)
226 p=load(os.path.join(f,"p.nc"),dom)
227 T=load(os.path.join(f,"T.nc"),dom)
228 if CREATE_TOPOGRAPHY:
229 topography=load(os.path.join(f,"topo.nc"),dom)
230
231 diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)
232 print "<%s> Restart from file %s at time step %s (t=%s) completed."%(time.asctime(),f,t,n)
233 else:
234 if DIM==2:
235 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=-1,optimize=True)
236 else:
237 dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L/H,l1=L/H,l2=1,order=-1,optimize=True)
238 try:
239 dom.dump("mesh.nc")
240 except:
241 pass
242 x=dom.getX()
243 T=Scalar(1,Solution(dom))
244 for d in range(DIM):
245 if d == DIM-1:
246 T*=sin(x[d]/H*pi)
247 else:
248 T*=cos(x[d]/L*pi)
249
250 T=(1.-x[DIM-1])+PERT*T
251 v=Vector(0,Solution(dom))
252 stress=Tensor(0,Function(dom))
253 x2=ReducedSolution(dom).getX()
254 p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)
255
256 if CREATE_TOPOGRAPHY:
257 topography=Scalar(0.,Solution(dom))
258 diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=False)
259 diagnostics_file.write("Ra = %e Lambda= %e\n"%(Ra, SURFACE_LOAD))
260
261 p_last=p
262 x=dom.getX()
263 #
264 # set up heat problem:
265 #
266 heat=TemperatureCartesian(dom,useBackwardEuler=False)
267 print "<%s> Temperature transport has been set up."%time.asctime()
268 heat.getSolverOptions().setTolerance(T_TOL)
269 heat.getSolverOptions().setVerbosity(VERBOSE)
270 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
271 heat.setInitialTemperature(clip(T,T_0))
272 heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
273 #
274 # velocity constraints:
275 #
276 fixed_v_mask=Vector(0,Solution(dom))
277 faces=Scalar(0.,Solution(dom))
278 for d in range(DIM):
279 if d == DIM-1:
280 ll = H
281 else:
282 ll = L
283 if CREATE_TOPOGRAPHY and d==DIM-1:
284 fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
285 else:
286 s=whereZero(x[d])+whereZero(x[d]-ll)
287 faces+=s
288 fixed_v_mask+=s*unitVector(d,DIM)
289 #
290 # set up velovity problem
291 #
292 flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE)
293 flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1]))
294
295 flow.setElasticShearModulus(MUE)
296 flow.setTolerance(FLOW_TOL)
297 flow.setEtaTolerance(FLOW_TOL)
298 flow.setExternals(fixed_v_mask=fixed_v_mask)
299 print "<%s> Flow solver has been set up."%time.asctime()
300 #
301 # topography set-up
302 #
303 boundary=FunctionOnBoundary(dom).getX()[DIM-1]
304 top_boundary_mask=whereZero(boundary-sup(boundary))
305 surface_area=integrate(top_boundary_mask)
306 if CREATE_TOPOGRAPHY:
307 mts=Mountains(dom,eps=TOPO_SMOOTH)
308 mts.setTopography(topography)
309 print "<%s> topography has been set up."%time.asctime()
310 #
311 # let the show begin:
312 #
313 t1 = time.time()
314 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
315 while t<T_END:
316 if CREATE_TOPOGRAPHY: topography_old=topography
317 v_old, p_old, stress_old=v, p, stress
318 T_old=T
319 #======= solve for velovity ====================================================================
320 eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
321 print "viscosity range :", inf(eta_N), sup(eta_N)
322 flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0], [1,N])
323 flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM))
324 # if dt<=0 or not CREATE_TOPOGRAPHY:
325 if not CREATE_TOPOGRAPHY:
326 flow.update(dt, iter_max=100, verbose=False)
327 else:
328 topography_last=topography
329 Topo_norm, error_Topo=1,1
330 i=0
331 print "DDDDD : ====",dt
332 while error_Topo > TOPO_TOL * Topo_norm:
333 flow.setStatus(t, v_old, p_old, stress_old)
334 flow.setExternals(f=-SURFACE_LOAD*(topography-dt*v)*unitVector(DIM-1,DIM)*top_boundary_mask, restoration_factor=SURFACE_LOAD*dt*top_boundary_mask)
335 flow.update(dt, iter_max=100, verbose=False)
336 v=flow.getVelocity()
337 mts.setTopography(topography_old)
338 mts.setVelocity(v)
339 topography_last, topography=topography, mts.update(dt, allow_substeps=True)
340 error_Topo=sqrt(integrate(((topography-topography_last)*top_boundary_mask)**2))
341 Topo_norm=sqrt(integrate((topography*top_boundary_mask)**2))
342 print "DDDDD :", "input=",integrate(v*top_boundary_mask)[DIM-1],"output=", integrate(topography*top_boundary_mask)/Lsup(topography), error_Topo, Topo_norm
343 print "topography update step %s error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v)
344 i+=1
345 if i > TOPO_ITER_MAX:
346 raise RuntimeError,"topography did not converge after %s steps."%TOPO_ITER_MAX
347 v=flow.getVelocity()
348 for d in range(DIM):
349 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
350 print "<%s> flow solver completed."%time.asctime()
351 n+=1
352 t+=dt
353 print "influx= ",integrate(inner(v,dom.getNormal())), sqrt(integrate(length(v)**2,FunctionOnBoundary(dom))), integrate(1., FunctionOnBoundary(dom))
354 print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)
355 #======= setup Temperature problem ====================================================================
356 #
357 heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
358 dt=heat.getSafeTimeStepSize()
359 print "<%s> New time step size is %e"%(time.asctime(),dt)
360 #======= set-up topography ==================================================================================
361 if CREATE_TOPOGRAPHY:
362 dt=min(mts.getSafeTimeStepSize()*0.5,dt)
363 print "<%s> New time step size is %e"%(time.asctime(),dt)
364 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)
365 #
366 # solve temperature:
367 #
368 T=heat.getSolution(dt)
369 print "Temperature range ",inf(T),sup(T)
370 print "<%s> temperature update completed."%time.asctime()
371 #======= anaysis ==================================================================================
372 #
373 # .... Nusselt number
374 #
375 dTdz=grad(T)[DIM-1]
376 Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
377 eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
378 Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
379 print "nusselt number = %e"%Nu
380 print "av. eta = %e"%eta_bar
381 print "effective Rayleigh number = %e"%Ra_eff
382 if CREATE_TOPOGRAPHY:
383 topo_level=integrate(topography*top_boundary_mask)/surface_area
384 valleys_deep=inf(topography)
385 mountains_heigh=sup(topography)
386 print "topography level = ",topo_level
387 print "valleys deep = ", valleys_deep
388 print "mountains_heigh = ", mountains_heigh
389 diagnostics_file.write("%e %e %e %e %e %e %e\n"%(t,Nu, topo_level, valleys_deep, mountains_heigh, eta_bar, Ra_eff))
390 else:
391 diagnostics_file.write("%e %e %e %e\n"%(t,Nu, eta_bar, Ra_eff))
392 #
393 # .... visualization
394 #
395 if t>=t_vis or n>n_vis:
396 if CREATE_TOPOGRAPHY:
397 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff(), topography=topography_old, vex=mts.getVelocity())
398 else:
399 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())
400 print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)
401 counter_vis+=1
402 t_vis+=DT_VIS
403 n_vis+=DN_VIS
404 # =========================
405 #
406 # create restart files:
407 #
408 if DN_RESTART>0:
409 if (n-1)%DN_RESTART == 0:
410 c=(n-1)/DN_RESTART
411 old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)
412 new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)
413 mkDir(new_restart_dir)
414 dom.MPIBarrier()
415 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))
416 v.dump(os.path.join(new_restart_dir,"v.nc"))
417 p.dump(os.path.join(new_restart_dir,"p.nc"))
418 T.dump(os.path.join(new_restart_dir,"T.nc"))
419 if CREATE_TOPOGRAPHY: topography.dump(os.path.join(new_restart_dir,"topo.nc"))
420 removeRestartDirectory(old_restart_dir)
421 print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)
422 print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)
423

  ViewVC Help
Powered by ViewVC 1.1.26