/[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 2438 - (show annotations)
Tue May 26 02:26:52 2009 UTC (12 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 10675 byte(s)
new version of power law (blame Hans)
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 """
16 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
17 Earth Systems Science Computational Center (ESSCC)
18 http://www.uq.edu.au/esscc
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 from esys.escript import *
25 from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
26 from esys.finley import Rectangle, Brick, LoadMesh
27 from optparse import OptionParser
28 from math import pi, ceil
29 import sys
30 import time
31
32 # ======================= Default Values ==================================================
33 DIM=2 # spatial dimension
34 H=1. # height
35 L=4*H # length
36 NE=10 # number of elements in H-direction.
37 PERTURBATION=0.1 # initial temperature perturbation
38 DT=1.e-4 # initial time step size
39 USE_BACKWARD_EULER=False # use backward Euler in time integartion of temperature otherwise Crank-Nicholson is used
40 CREATE_TOPGRAPHY=False # create topgraphy
41 TOL=1.e-4 # tolerance
42 DT_MIN=1.e-5 # minumum time step size
43 T_END=0.1 # end time
44
45 DT_OUT=T_END/500
46 Dn_OUT=2
47 VERBOSE=True
48 CREATE_RESTARTFILES_EVERY_STEP=10
49 if True:
50 # this is a simple linear Stokes model:
51 RA=1.e6 # Rayleigh number
52 A=22 # Arenious number
53 DI = 0. # dissipation number
54 MU=None
55 ETA0=1.
56 TAU0=None
57 N=None
58 NPL=None
59 ETAP0=ETA0
60 TAUY=None
61 useJAUMANNSTRESS=False
62 # this is a simple linear Stokes model:
63 RA=1.e5 # Rayleigh number
64 A=11 # Arenious number
65 DI = 0. # dissipation number
66 MU=None
67 ETA0=1.
68 TAU0=250.
69 N=None
70 NPL=None
71 ETAP0=ETA0
72 TAUY=TAU0
73 useJAUMANNSTRESS=False
74 else:
75 RA=1.e4 # Rayleigh number
76 A=22 # Arenious number
77 DI = 0. # dissipation number
78 MU=1.e4
79 ETA0=1.
80 TAU0=0.866*10**2.5
81 N=3
82 NPL=14
83 TAUY=TAU0
84 ETAP0=ETA0
85 useJAUMANNSTRESS=True
86
87
88 # =========================================================================================
89
90 def removeRestartDirectory(dir_name):
91 if dom.onMasterProcessor() and os.path.isdir(dir_name):
92 for root, dirs, files in os.walk(dir_name, topdown=False):
93 for name in files: os.remove(os.path.join(root,name))
94 for name in dirs: os.remove(os.path.join(root,name))
95 os.rmdir(dir_name)
96 print "Restart files %s have been removed."%dir_name
97 dom.MPIBarrier()
98 # =========================================================================================
99 #
100 # read options:
101 #
102 parser = OptionParser(usage="%prog [Options]")
103 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")
104 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)
105 parser.add_option("-p", "--param", dest="param", help="name of file to be imported ",metavar="PARAM", default=None)
106 (options, args) = parser.parse_args()
107 restart=options.restart or (options.restart_dir !=None)
108 #
109 # overwrite the default options:
110 #
111 print "Execution started ",time.asctime()
112 if options.param !=None:
113 exec(open(options.param,'r'))
114 print "Parameters are imported from file ",options.param
115 print "Parameters:"
116 print "\tDimension\tDIM=\t",DIM
117 print "\tHeight:\tH=\t",H
118 print "\tLength:\tL=\t",L
119 print "\telements in H:\tNE=\t",NE
120 print "\ttotal #element\t\t=\t",NE**DIM*int(L/H)**(DIM-1)
121 print "\ttolerance\tTOL=\t",TOL
122 print "\ttemperature perturbation\tPERTURBATION=\t",PERTURBATION
123 print "\tinitial time step size\tDT=\t",DT
124 print "\tminimum time step size\tDT_MIN=\t",tDT_MIN
125 print "\tend time\tT_END=\t",T_END
126 print "\tbackward Euler?\tUSE_BACKWARD_EULER=\t",USE_BACKWARD_EULER
127 print "\ttopgraphy?\tCREATE_TOPOGRAPHY=\t",CREATE_TOPOGRAPHY
128 # some control variables (will be overwritten in case of a restart:
129 #
130 t=0 # time stamp
131 n=0 # time step counter
132 dt=DT # current time step size
133
134 t_out=0 #
135 n_out=0
136 out_count=0
137 a=None
138 dt_a=None
139 #
140 # set up domain or read restart file
141 #
142 if restart:
143 if options.restart_dir ==None:
144 restart_files=[]
145 for f in os.listdir("."):
146 if f.startswith("restart"): restart_files.append(f)
147 if len(restart_files)==0:
148 raise IOError,"no restart files"
149 restart_files.sort()
150 f=restart_files[-1]
151 else:
152 f=options.restart_dir
153 print ">>>Restart from directory ",f
154 try:
155 dom=LoadMesh("mesh.nc")
156 except:
157 pass
158 FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
159 t=float(FF[0]) # time stamp
160 t_out=float(FF[1]) #
161 n_out=int(FF[2])
162 n=int(FF[3])
163 out_count=int(FF[4])
164 dt=float(FF[5])
165 stress=load(os.path.join(f,"stress.nc"),dom)
166 v=load(os.path.join(f,"v.nc"),dom)
167 p=load(os.path.join(f,"p.nc"),dom)
168 T=load(os.path.join(f,"T.nc"),dom)
169 if n>1:
170 dt_a=float(FF[6])
171 a=load(os.path.join(f,"a.nc"),dom)
172 else:
173 dt_a=None
174 a=None
175 if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","a")
176 else:
177 if DIM==2:
178 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
179 else:
180 dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True,optimize=True)
181 try:
182 dom.dump("mesh.nc")
183 except:
184 pass
185 x=dom.getX()
186 T=Scalar(1,ReducedSolution(dom))
187 for d in range(DIM):
188 if d == DIM-1:
189 T*=sin(x[d]/H*pi)
190 else:
191 T*=cos(x[d]/L*pi)
192
193 T=1.-x[DIM-1]+PERTURBATION*T
194 v=Vector(0,Solution(dom))
195 stress=Tensor(0,Function(dom))
196 x2=ReducedSolution(dom).getX()
197 p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
198 if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","w")
199
200 vol=integrate(1.,Function(dom))
201 p-=integrate(p)/vol
202 x=dom.getX()
203 #
204 # set up heat problem:
205 #
206 heat=TemperatureCartesian(dom,useBackwardEuler=USE_BACKWARD_EULER)
207 heat.setTolerance(TOL)
208 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
209 print "initial Temperature range ",inf(T),sup(T)
210 heat.setInitialTemperature(clip(T,0))
211 heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
212 #
213 # velocity constraints:
214 #
215 x2=ReducedSolution(dom).getX()
216 fixed_v_mask=Vector(0,Solution(dom))
217 for d in range(DIM):
218 if d == DIM-1:
219 ll = H
220 else:
221 ll = L
222 if CREATE_TOPOGRAPHY and d==DIM-1:
223 fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
224 else:
225 fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
226 1/0
227 #
228 # set up velovity problem
229 # ????????????????????
230 sp=StokesProblemCartesian(dom)
231 # ,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
232 # sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
233 # sp.initialize(mu=MU, tau_0=TAU0, n=N, eta_Y=ETAP0, tau_Y=TAU0, n_Y=NPL, q=fixed_v_mask)
234 sp.setTolerance(TOL*extratol)
235 # sp.setToleranceReductionFactor(TOL)
236 #
237 # let the show begin:
238 #
239 t1 = time.time()
240 while t<T_END:
241 v_last=v*1
242 print "============== solve for v ========================"
243 FF=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
244 print "viscosity range :", inf(FF)*ETA0, sup(FF)*ETA0
245 # sp.initialize(eta_N=ETA0*FF, eta_0=ETA0*FF, F=T*(RA*unitVector(DIM-1,DIM)))
246 sp.initialize(eta=ETA0*FF, f=T*(RA*unitVector(DIM-1,DIM)),fixed_u_mask=fixed_v_mask)
247 v,p=sp.solve(v,p,max_iter=20, verbose=VERBOSE, show_details=False)
248
249 for d in range(DIM):
250 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
251 if t>=t_out or n>n_out:
252 saveVTK("state.%d.vtu"%out_count,T=T,v=v)
253 print "visualization file %d for time step %e generated."%(out_count,t)
254 out_count+=1
255 t_out+=DT_OUT
256 n_out+=Dn_OUT
257 # calculation of nusselt number:
258 se=ETA0*FF*length(symmetric(grad(v)))**2 # this propably wroing!
259 print "Xse:",inf(se),sup(se)
260 Nu=1.+integrate(se)/(RA*vol)
261 if dom.onMasterProcessor(): nusselt_file.write("%e %e\n"%(t,Nu))
262 heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)
263 print "Xnusselt number = ",Nu, "dt =",dt
264 if n>0:
265 a,a_alt = (v_last-v)/dt, a
266 dt_a,dt_a_alt = dt, dt_a
267 if n>1:
268 z=(a-a_alt)/((dt_a+dt_a_alt)/2)
269 f=Lsup(z)/Lsup(v)
270 print "estimated error ",f*dt**2
271 # dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))
272 dt_new=sqrt(0.05/f)
273 dt=min(dt_new,heat.getSafeTimeStepSize())
274 else:
275 dt=heat.getSafeTimeStepSize()
276 dt=max(DT_MIN,dt)
277 print n,". time step t=",t," step size ",dt
278 print "============== solve for T ========================"
279 T=heat.getSolution(dt, verbose=VERBOSE)
280 print "Temperature range ",inf(T),sup(T)
281 n+=1
282 t+=dt
283 # =========================
284 #
285 # create restart files:
286 #
287 if create_restartfiles_every_step>0:
288 if (n-1)%create_restartfiles_every_step == 0:
289 c=(n-1)/create_restartfiles_every_step
290 old_restart_dir="restart_%s_"%(c-1)
291 new_restart_dir="restart_%s_"%c
292
293 print "Write restart files to ",new_restart_dir
294 if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
295 dom.MPIBarrier()
296 # sp.getStress().dump(os.path.join(new_restart_dir,"stress.nc"))
297 v.dump(os.path.join(new_restart_dir,"v.nc"))
298 p.dump(os.path.join(new_restart_dir,"p.nc"))
299 T.dump(os.path.join(new_restart_dir,"T.nc"))
300 if n>1:
301 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e; %e;\n"%(t, t_out, n_out, n, out_count, dt, dt_a))
302 a.dump(os.path.join(new_restart_dir,"a.nc"))
303 else:
304 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_out, n_out, n, out_count, dt))
305 removeRestartDirectory(old_restart_dir)
306 elapsed = time.time() - t1
307 print "plot","\t",NE,"\t",elapsed
308

  ViewVC Help
Powered by ViewVC 1.1.26