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

  ViewVC Help
Powered by ViewVC 1.1.26