/[escript]/trunk/dudley/test/python/convection.py
ViewVC logotype

Contents of /trunk/dudley/test/python/convection.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3774 - (show annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 15847 byte(s)
dudley, pasowrap, pycad

1 ########################################################
2 #
3 # Copyright (c) 2003-2010 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-2010 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 import DataManager as DM
28 from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian, Mountains, SubSteppingException
29 from esys.dudley import Rectangle, Brick, LoadMesh
30 from optparse import OptionParser
31 from math import pi, ceil
32 import sys
33 import time
34
35 # ============================= Default Values ===============================
36
37 DIM=2 # spatial dimension
38 H=1. # height
39 L=2*H # length
40 NE=30 # number of elements in H-direction
41 PERT=0.15 # initial temperature perturbation
42 DT=1.e-4 # initial time step size
43 CREATE_TOPO=False # create topography
44 DT_MIN=1.e-10 # minimum time step size
45 T_END=10. # end time
46
47 RHO_0=100. # surface density (lf ~ RHO_0**2)
48 G=1. # gravitational constant
49 ALPHA_0=0.1 # thermal expansion coefficient (Ra ~ RHO_0**2 * ALPHA_0 = lf * ALPHA_0)
50 T_0=0. # surface temperature
51 T_1=1. # bottom temperature
52 C_P=1 # heat capacity
53 K=1. # thermal conductivity
54 CHI=0. # Taylor-Quinny coefficient
55 MUE=None # elastic shear modulus
56 TAU_Y=5*10**(2.5) # Drucker-Prager cohesion factor
57 BETA=0 # Drucker-Prager friction factor
58 TAU_0=2*10**(2.5) # transition stress
59 N=3 # power for power law
60
61 E=23*0 # activation energy
62 V=18*0 # activation volume
63 T_OFFSET=1 # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0)
64 R=1 # gas constant
65 ETA_N0=1. # viscosity at surface
66
67 TOPO_SMOOTH=1e-5 # smoothing factor for extrapolation of surface velocity to interior
68 T_TOL=1.e-4 # tolerance temperature transport
69 FLOW_TOL=1.e-3 # tolerance for inconcompressible flow solver
70 TOPO_TOL=0.1 # tolerance for update of topography
71 DIAGNOSTICS_FN="diag.csv" # filename for diagnostics
72 VERBOSE=False # output more messages from solvers
73 DT_VIS=T_END/500 # time difference between visualization files
74 DN_VIS=5 # max. number of steps between visualization files
75 DN_RESTART=1000 # create restart files every DN_RESTART steps
76 PREFIX_RESTART="restart" # name prefix for restart directories
77 TOPO_ITER_MAX=20 # max. number of iteration steps to update topography
78
79 # ============================================================================
80
81 #
82 # read options:
83 #
84 parser = OptionParser(usage="%prog [Options]")
85 parser.add_option("-r", "--restart", dest="restart",
86 help="restart from latest checkpoint directory. It will be deleted after new data is exported.", default=False, action="store_true")
87 parser.add_option("-d", "--dir", dest="restart_dir",
88 help="locate/create restart directories under DIR.", metavar="DIR", default='.')
89 parser.add_option("-p", "--param", dest="param",
90 help="name of file to be imported ", metavar="PARAM", default=None)
91 (options, args) = parser.parse_args()
92 restart=options.restart
93
94 #
95 # overwrite the default options:
96 #
97 print(("<%s> Execution started."%time.asctime()))
98 if options.param !=None:
99 exec(open(options.param,'r'))
100 print(("Parameters imported from file ",options.param))
101
102 print("Input Parameters:")
103 print(("\tDimension DIM\t\t= %d"%DIM))
104 print(("\tHeight H\t\t\t= %s"%H))
105 print(("\tLength L\t\t\t= %s"%L))
106 print(("\tElements in H NE\t\t= %d"%NE))
107 print(("\tTemperature perturbation PERT\t\t= %s"%PERT))
108 print(("\tInitial time step size DT\t\t= %s"%DT))
109 print(("\tMinimum time step size DT_MIN\t\t= %s"%DT_MIN))
110 print(("\tEnd time T_END\t\t= %s"%T_END))
111 print(("\tCreate topography CREATE_TOPO\t= %s"%CREATE_TOPO))
112 print(("\tSurface density RHO_0\t\t= %s"%RHO_0))
113 print(("\tGravitational constant G\t\t\t= %s"%G))
114 print(("\tThermal expansion coefficient ALPHA_0\t\t= %s"%ALPHA_0))
115 print(("\tSurface temperature T_0\t\t= %s"%T_0))
116 print(("\tBottom temperature T_1\t\t= %s"%T_1))
117 print(("\tHeat capacity C_P\t\t= %s"%C_P))
118 print(("\tThermal conductivity K\t\t\t= %s"%K))
119 print(("\tTaylor-Quinny coefficient CHI\t\t= %s"%CHI))
120 print(("\tElastic shear modulus MUE\t\t= %s"%MUE))
121 print(("\tCohesion factor TAU_Y\t\t= %s"%TAU_Y))
122 print(("\tFriction factor BETA\t\t= %s"%BETA))
123 print(("\tTransition stress TAU_0\t\t= %s"%TAU_0))
124 print(("\tPower for power law N\t\t\t= %s"%N))
125 print(("\tViscosity at surface ETA_N0\t\t= %s"%ETA_N0))
126 print(("\tActivation energy E\t\t\t= %s"%E))
127 print(("\tActivation volume V\t\t\t= %s"%V))
128 print(("\tTemperature offset T_OFFSET\t\t= %s"%T_OFFSET))
129 print(("\tGas constant R\t\t\t= %s"%R))
130 print(("\tTopography smoothing TOPO_SMOOTH\t= %s"%TOPO_SMOOTH))
131 print(("\tTolerance for topography TOPO_TOL\t\t= %s"%TOPO_TOL))
132 print(("\tTransport tolerance T_TOL\t\t= %s"%T_TOL))
133 print(("\tFlow tolerance FLOW_TOL\t\t= %s"%FLOW_TOL))
134 #print("\tFile for diagnostics DIAGNOSTICS_FN\t= %s"%DIAGNOSTICS_FN)
135 print(("\tRestart counter increment DN_RESTART\t= %d"%DN_RESTART))
136 print(("\tPrefix for restart dirs PREFIX_RESTART\t= %s"%PREFIX_RESTART))
137 print(("\tVerbosity VERBOSE\t\t= %s"%VERBOSE))
138
139 print("Control Parameters:")
140 t_REF=RHO_0*C_P*H**2/K
141 P_REF=ETA_N0/t_REF
142 Ar=E/R/(T_1-T_0)
143 if E>0 and V>0:
144 V_REF=P_REF*V/E
145 else:
146 V_REF=0
147 T_OFFSET_REF=T_OFFSET/(T_1-T_0)
148 Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF
149 Di=ALPHA_0*G*H/C_P
150 CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2)
151 if CREATE_TOPO:
152 SURFACE_LOAD=RHO_0*G*H/P_REF
153 else:
154 SURFACE_LOAD=0.
155 if MUE == None:
156 De=None
157 else:
158 De=ETA_N0/MUE/t_REF
159 ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0
160 print(("\tTotal #elements \t\t\t= %d"%(NE**DIM*int(L/H)**(DIM-1))))
161 print(("\tReference time t_REF\t\t= %s"%t_REF))
162 print(("\tReference pressure P_REF\t\t= %s"%P_REF))
163 print(("\tReference Taylor-Quinny CHI_REF\t\t= %s"%CHI_REF))
164 print(("\tDissipation number DI\t\t= %s"%Di))
165 print(("\tRayleigh number surface Ra\t\t= %s"%Ra))
166 print(("\tDebora number surface De\t\t= %s"%De))
167 print(("\tBottom viscosity \t\t\t= %s"%ETA_BOT))
168 print(("\tRayleigh number bottom \t\t\t= %s"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT)))
169 if MUE == None:
170 print("\tDebora number bottom \t\t\t= None")
171 else:
172 print(("\tDebora number bottom \t\t\t= %s"%(ETA_BOT/MUE/t_REF)))
173 print(("\tArrhenius Ar\t\t= %s"%Ar))
174 print(("\tSurface load factor SURFACE_LOAD\t= %s"%SURFACE_LOAD))
175 print(("\tScaled activation volume V_REF\t\t= %s"%V_REF))
176 print()
177
178 # some control variables (will be overwritten in case of a restart:
179 t=0 # time stamp
180 n=0 # time step counter
181 t_vis=DT_VIS # time of next visualization file export
182 dt=DT # current time step size
183 #=========================
184 #
185 # set up domain and data from scratch or from restart files
186 #
187 dataMgr=DM(formats=[DM.RESTART,DM.VTK], work_dir=options.restart_dir, restart_prefix=PREFIX_RESTART, do_restart=restart)
188 dataMgr.setCheckpointFrequency(DN_RESTART/DN_VIS)
189 if dataMgr.hasData():
190 dom=dataMgr.getDomain()
191 t=dataMgr.getValue('t')
192 t_vis=dataMgr.getValue('t_vis')
193 n=dataMgr.getValue('n')
194 dt=dataMgr.getValue('dt')
195 stress=dataMgr.getValue('stress')
196 v=dataMgr.getValue('v')
197 p=dataMgr.getValue('p')
198 T=dataMgr.getValue('T')
199 if CREATE_TOPO:
200 topography=dataMgr.getValue('topography')
201
202 #diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)
203 print(("<%s> Restart at time step %d (t=%e) completed."%(time.asctime(),n,t)))
204 else:
205 if DIM==2:
206 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=-1,optimize=True)
207 else:
208 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)
209 x=dom.getX()
210 T=Scalar(1,Solution(dom))
211 for d in range(DIM):
212 if d == DIM-1:
213 T*=sin(x[d]/H*pi)
214 else:
215 T*=cos(x[d]/L*pi)
216
217 T=(1.-x[DIM-1])+PERT*T
218 v=Vector(0,Solution(dom))
219 stress=Tensor(0,Function(dom))
220 x2=ReducedSolution(dom).getX()
221 p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)
222
223 if CREATE_TOPO:
224 topography=Scalar(0.,Solution(dom))
225 #diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=False)
226 #diagnostics_file.write("Ra = %e Lambda= %e\n"%(Ra, SURFACE_LOAD))
227
228 p_last=p
229 x=dom.getX()
230 #
231 # set up heat problem:
232 #
233 heat=TemperatureCartesian(dom,useBackwardEuler=False)
234 print(("<%s> Temperature transport has been set up."%time.asctime()))
235 heat.getSolverOptions().setTolerance(T_TOL)
236 heat.getSolverOptions().setVerbosity(VERBOSE)
237 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
238 heat.setInitialTemperature(clip(T,T_0))
239 heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
240 #
241 # velocity constraints:
242 #
243 fixed_v_mask=Vector(0,Solution(dom))
244 faces=Scalar(0.,Solution(dom))
245 for d in range(DIM):
246 if d == DIM-1:
247 ll = H
248 else:
249 ll = L
250 if CREATE_TOPO and d==DIM-1:
251 fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
252 else:
253 s=whereZero(x[d])+whereZero(x[d]-ll)
254 faces+=s
255 fixed_v_mask+=s*unitVector(d,DIM)
256 #
257 # set up velocity problem
258 #
259 flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE)
260 flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1]))
261
262 flow.setElasticShearModulus(MUE)
263 flow.setTolerance(FLOW_TOL)
264 flow.setEtaTolerance(FLOW_TOL)
265 flow.setExternals(fixed_v_mask=fixed_v_mask)
266 print(("<%s> Flow solver has been set up."%time.asctime()))
267 #
268 # topography setup
269 #
270 boundary=FunctionOnBoundary(dom).getX()[DIM-1]
271 top_boundary_mask=whereZero(boundary-sup(boundary))
272 surface_area=integrate(top_boundary_mask)
273 if CREATE_TOPO:
274 mts=Mountains(dom,eps=TOPO_SMOOTH)
275 mts.setTopography(topography)
276 print(("<%s> topography has been set up."%time.asctime()))
277
278 #
279 # let the show begin:
280 #
281 t1 = time.time()
282 print(("<%s> Start time step %d (t=%s)."%(time.asctime(),n,t)))
283 while t<T_END:
284 if CREATE_TOPO: topography_old=topography
285 v_old, p_old, stress_old=v, p, stress
286 T_old=T
287 #========================= solve for velocity ============================
288 eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
289 #print("viscosity range :", inf(eta_N), sup(eta_N))
290 flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0], [1,N])
291 flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM))
292 # if dt<=0 or not CREATE_TOPO:
293 if not CREATE_TOPO:
294 flow.update(dt, iter_max=100, verbose=False)
295 else:
296 topography_last=topography
297 Topo_norm, error_Topo=1,1
298 i=0
299 while error_Topo > TOPO_TOL * Topo_norm:
300 flow.setStatus(t, v_old, p_old, stress_old)
301 flow.setExternals(f=-SURFACE_LOAD*(topography-dt*v)*unitVector(DIM-1,DIM)*top_boundary_mask, restoration_factor=SURFACE_LOAD*dt*top_boundary_mask)
302 flow.update(dt, iter_max=100, verbose=False)
303 v=flow.getVelocity()
304 mts.setTopography(topography_old)
305 mts.setVelocity(v)
306 topography_last, topography=topography, mts.update(dt, allow_substeps=True)
307 error_Topo=sqrt(integrate(((topography-topography_last)*top_boundary_mask)**2))
308 Topo_norm=sqrt(integrate((topography*top_boundary_mask)**2))
309 #print("topography update step %d error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v))
310 i+=1
311 if i > TOPO_ITER_MAX:
312 raise RuntimeError("topography did not converge after %d steps."%TOPO_ITER_MAX)
313 v=flow.getVelocity()
314 #for d in range(DIM):
315 #print("range %d-velocity "%d,inf(v[d]),sup(v[d]))
316 #print("Courant = ",inf(dom.getSize()/(length(v)+1e-19)), inf(dom.getSize()**2))
317 #print("<%s> flow solver completed."%time.asctime())
318 n+=1
319 t+=dt
320 #print("influx= ",integrate(inner(v,dom.getNormal())), sqrt(integrate(length(v)**2,FunctionOnBoundary(dom))), integrate(1., FunctionOnBoundary(dom)))
321 #print("<%s> Time step %d (t=%e) completed."%(time.asctime(),n,t))
322
323 #====================== setup temperature problem ========================
324 heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
325 dt=heat.getSafeTimeStepSize()
326 #print("<%s> New time step size is %e"%(time.asctime(),dt))
327 #=========================== setup topography ============================
328 if CREATE_TOPO:
329 dt=min(mts.getSafeTimeStepSize()*0.5,dt)
330 #print("<%s> New time step size is %e"%(time.asctime(),dt))
331 #print("<%s> Start time step %d (t=%e)."%(time.asctime(),n+1,t+dt))
332 #
333 # solve temperature:
334 #
335 T=heat.getSolution(dt)
336 #print("Temperature range ",inf(T),sup(T))
337 #print("<%s> temperature update completed."%time.asctime())
338 #============================== analysis =================================
339 #
340 # .... Nusselt number
341 #
342 dTdz=grad(T)[DIM-1]
343 Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
344 eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
345 Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
346 #print("nusselt number = ",Nu)
347 #print("avg. eta = ",eta_bar)
348 #print("effective Rayleigh number = ",Ra_eff)
349 if CREATE_TOPO:
350 topo_level=integrate(topography*top_boundary_mask)/surface_area
351 valleys_deep=inf(topography)
352 mountains_heigh=sup(topography)
353 #print("topography level = ",topo_level)
354 #print("valleys deep = ",valleys_deep)
355 #print("mountains_heigh = ",mountains_heigh)
356 #diagnostics_file.write("%e %e %e %e %e %e %e\n"%(t,Nu, topo_level, valleys_deep, mountains_heigh, eta_bar, Ra_eff))
357 #else:
358 #diagnostics_file.write("%e %e %e %e\n"%(t,Nu, eta_bar, Ra_eff))
359 # ========================================================================
360 #
361 # create restart/visualization files:
362 #
363 if t>=t_vis or n%DN_VIS==0:
364 dataMgr.setTime(t)
365 t_vis+=DT_VIS
366 dataMgr.addData(t=t,n=n,t_vis=t_vis,dt=dt,T=T,v=v,eta=flow.getCurrentEtaEff(),stress=stress,p=p)
367 if CREATE_TOPO: dataMgr.addData(topography=topography)
368 dataMgr.export()
369 print(("<%s> Cycle %d (time %s) exported."%(time.asctime(),dataMgr.getCycle(),t)))
370
371 print(("<%s> Calculation finalized after %s seconds."%(time.asctime(),time.time()-t1)))
372

  ViewVC Help
Powered by ViewVC 1.1.26