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

  ViewVC Help
Powered by ViewVC 1.1.26