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

  ViewVC Help
Powered by ViewVC 1.1.26