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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5706 - (hide annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 3 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 jfenwick 3981 ##############################################################################
2 ksteube 1811 #
3 jfenwick 5593 # Copyright (c) 2003-2015 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 ksteube 1811 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 gross 2438 """
16     this is a convection simulation over a domain [0,L] X [0,L] x [0,H]
17 ksteube 1811
18 gross 2440 It is solved in dimensionless form
19    
20 gross 2438 """
21 sshaw 5706
22     from __future__ import print_function, division
23    
24 jfenwick 5593 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
25 jfenwick 3981 http://www.uq.edu.au
26 ksteube 1811 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
30 ksteube 1811
31 gross 1417 from esys.escript import *
32 caltinay 3346 from esys.escript import DataManager as DM
33 gross 2620 from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian, Mountains, SubSteppingException
34 jfenwick 3087 from esys.dudley import Rectangle, Brick, LoadMesh
35 gross 1472 from optparse import OptionParser
36 gross 1417 from math import pi, ceil
37 artak 1483 import sys
38     import time
39    
40 caltinay 3346 # ============================= Default Values ===============================
41 artak 1483
42 caltinay 3346 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 gross 1474
52 caltinay 3346 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 gross 1472
66 caltinay 3346 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 gross 2620
72 caltinay 3346 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 gross 2438
84 caltinay 3346 # ============================================================================
85 gross 2440
86 gross 1417 #
87 gross 2438 # read options:
88 gross 1417 #
89 gross 2438 parser = OptionParser(usage="%prog [Options]")
90 caltinay 3346 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 gross 2438 (options, args) = parser.parse_args()
97 caltinay 3346 restart=options.restart
98    
99 gross 2438 #
100     # overwrite the default options:
101     #
102 jfenwick 3774 print(("<%s> Execution started."%time.asctime()))
103 gross 2438 if options.param !=None:
104 caltinay 3346 exec(open(options.param,'r'))
105 jfenwick 3774 print(("Parameters imported from file ",options.param))
106 gross 2440
107 caltinay 3346 print("Input Parameters:")
108 jfenwick 3774 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 caltinay 3346 #print("\tFile for diagnostics DIAGNOSTICS_FN\t= %s"%DIAGNOSTICS_FN)
140 jfenwick 3774 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 gross 2440
144 caltinay 3346 print("Control Parameters:")
145 gross 2440 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 caltinay 3346 V_REF=P_REF*V/E
150 gross 2440 else:
151 caltinay 3346 V_REF=0
152 gross 2440 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 caltinay 3346 if CREATE_TOPO:
157     SURFACE_LOAD=RHO_0*G*H/P_REF
158 gross 2620 else:
159 caltinay 3346 SURFACE_LOAD=0.
160 gross 2440 if MUE == None:
161 caltinay 3346 De=None
162 gross 2440 else:
163 caltinay 3346 De=ETA_N0/MUE/t_REF
164 gross 2440 ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0
165 jfenwick 3774 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 gross 2440 if MUE == None:
175 caltinay 3346 print("\tDebora number bottom \t\t\t= None")
176 gross 2440 else:
177 jfenwick 3774 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 caltinay 3346
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 gross 2440 #=========================
189 gross 2438 #
190 caltinay 3346 # set up domain and data from scratch or from restart files
191 gross 2438 #
192 caltinay 3346 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 gross 2440
207 caltinay 3346 #diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)
208 jfenwick 3774 print(("<%s> Restart at time step %d (t=%e) completed."%(time.asctime(),n,t)))
209 gross 1417 else:
210 caltinay 3346 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 gross 1417
222 caltinay 3346 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 gross 1472
228 caltinay 3346 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 gross 2620
233 gross 2445 p_last=p
234 gross 1417 x=dom.getX()
235     #
236     # set up heat problem:
237     #
238 gross 3793 heat=TemperatureCartesian(dom)
239 jfenwick 3774 print(("<%s> Temperature transport has been set up."%time.asctime()))
240 gross 2474 heat.getSolverOptions().setTolerance(T_TOL)
241     heat.getSolverOptions().setVerbosity(VERBOSE)
242 gross 1417 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
243 gross 2440 heat.setInitialTemperature(clip(T,T_0))
244     heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
245 gross 1417 #
246 gross 1639 # velocity constraints:
247 gross 1417 #
248     fixed_v_mask=Vector(0,Solution(dom))
249 gross 2620 faces=Scalar(0.,Solution(dom))
250 gross 1417 for d in range(DIM):
251     if d == DIM-1:
252 caltinay 3346 ll = H
253 gross 1417 else:
254 caltinay 3346 ll = L
255     if CREATE_TOPO and d==DIM-1:
256     fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
257 gross 2438 else:
258 caltinay 3346 s=whereZero(x[d])+whereZero(x[d]-ll)
259     faces+=s
260     fixed_v_mask+=s*unitVector(d,DIM)
261 gross 1417 #
262 caltinay 3346 # set up velocity problem
263 gross 1639 #
264 gross 2440 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 gross 2862 flow.setTolerance(FLOW_TOL)
269 gross 2440 flow.setEtaTolerance(FLOW_TOL)
270     flow.setExternals(fixed_v_mask=fixed_v_mask)
271 jfenwick 3774 print(("<%s> Flow solver has been set up."%time.asctime()))
272 gross 2440 #
273 caltinay 3346 # topography setup
274 gross 2620 #
275     boundary=FunctionOnBoundary(dom).getX()[DIM-1]
276     top_boundary_mask=whereZero(boundary-sup(boundary))
277     surface_area=integrate(top_boundary_mask)
278 caltinay 3346 if CREATE_TOPO:
279 gross 2857 mts=Mountains(dom,eps=TOPO_SMOOTH)
280 gross 2620 mts.setTopography(topography)
281 jfenwick 3774 print(("<%s> topography has been set up."%time.asctime()))
282 caltinay 3346
283 gross 2620 #
284 caltinay 3346 # let the show begin:
285 gross 1417 #
286 gross 2438 t1 = time.time()
287 jfenwick 3774 print(("<%s> Start time step %d (t=%s)."%(time.asctime(),n,t)))
288 gross 1417 while t<T_END:
289 caltinay 3346 if CREATE_TOPO: topography_old=topography
290 gross 2620 v_old, p_old, stress_old=v, p, stress
291     T_old=T
292 caltinay 3346 #========================= solve for velocity ============================
293 gross 2440 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 caltinay 3346 #print("viscosity range :", inf(eta_N), sup(eta_N))
295 gross 2440 flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0], [1,N])
296 gross 2620 flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM))
297 caltinay 3346 # if dt<=0 or not CREATE_TOPO:
298     if not CREATE_TOPO:
299     flow.update(dt, iter_max=100, verbose=False)
300 gross 2620 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 gross 2862 flow.update(dt, iter_max=100, verbose=False)
308 gross 2620 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 caltinay 3346 #print("topography update step %d error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v))
315 gross 2620 i+=1
316     if i > TOPO_ITER_MAX:
317 caltinay 3346 raise RuntimeError("topography did not converge after %d steps."%TOPO_ITER_MAX)
318 gross 2440 v=flow.getVelocity()
319 caltinay 3346 #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 gross 1417 n+=1
324     t+=dt
325 caltinay 3346 #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 gross 2440 heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
330     dt=heat.getSafeTimeStepSize()
331 caltinay 3346 #print("<%s> New time step size is %e"%(time.asctime(),dt))
332     #=========================== setup topography ============================
333     if CREATE_TOPO:
334 gross 2620 dt=min(mts.getSafeTimeStepSize()*0.5,dt)
335 caltinay 3346 #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 gross 2440 #
338     # solve temperature:
339     #
340 gross 2474 T=heat.getSolution(dt)
341 caltinay 3346 #print("Temperature range ",inf(T),sup(T))
342     #print("<%s> temperature update completed."%time.asctime())
343     #============================== analysis =================================
344 gross 2440 #
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 caltinay 3346 #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 gross 2620 #
366 caltinay 3346 # create restart/visualization files:
367 gross 2620 #
368 caltinay 3346 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 jfenwick 3774 print(("<%s> Cycle %d (time %s) exported."%(time.asctime(),dataMgr.getCycle(),t)))
375 gross 1474
376 jfenwick 3774 print(("<%s> Calculation finalized after %s seconds."%(time.asctime(),time.time()-t1)))
377 caltinay 3346

  ViewVC Help
Powered by ViewVC 1.1.26