/[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 3793 - (hide annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 15824 byte(s)
new implementation of FCT solver with some modifications to the python interface
1 ksteube 1811 ########################################################
2     #
3 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
4 ksteube 1811 # 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 gross 2438 """
13     this is a convection simulation over a domain [0,L] X [0,L] x [0,H]
14 ksteube 1811
15 gross 2440 It is solved in dimensionless form
16    
17 gross 2438 """
18 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
19 ksteube 1811 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
25 ksteube 1811
26 gross 1417 from esys.escript import *
27 caltinay 3346 from esys.escript import DataManager as DM
28 gross 2620 from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian, Mountains, SubSteppingException
29 jfenwick 3087 from esys.dudley import Rectangle, Brick, LoadMesh
30 gross 1472 from optparse import OptionParser
31 gross 1417 from math import pi, ceil
32 artak 1483 import sys
33     import time
34    
35 caltinay 3346 # ============================= Default Values ===============================
36 artak 1483
37 caltinay 3346 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 gross 1474
47 caltinay 3346 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 gross 1472
61 caltinay 3346 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 gross 2620
67 caltinay 3346 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 gross 2438
79 caltinay 3346 # ============================================================================
80 gross 2440
81 gross 1417 #
82 gross 2438 # read options:
83 gross 1417 #
84 gross 2438 parser = OptionParser(usage="%prog [Options]")
85 caltinay 3346 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 gross 2438 (options, args) = parser.parse_args()
92 caltinay 3346 restart=options.restart
93    
94 gross 2438 #
95     # overwrite the default options:
96     #
97 jfenwick 3774 print(("<%s> Execution started."%time.asctime()))
98 gross 2438 if options.param !=None:
99 caltinay 3346 exec(open(options.param,'r'))
100 jfenwick 3774 print(("Parameters imported from file ",options.param))
101 gross 2440
102 caltinay 3346 print("Input Parameters:")
103 jfenwick 3774 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 caltinay 3346 #print("\tFile for diagnostics DIAGNOSTICS_FN\t= %s"%DIAGNOSTICS_FN)
135 jfenwick 3774 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 gross 2440
139 caltinay 3346 print("Control Parameters:")
140 gross 2440 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 caltinay 3346 V_REF=P_REF*V/E
145 gross 2440 else:
146 caltinay 3346 V_REF=0
147 gross 2440 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 caltinay 3346 if CREATE_TOPO:
152     SURFACE_LOAD=RHO_0*G*H/P_REF
153 gross 2620 else:
154 caltinay 3346 SURFACE_LOAD=0.
155 gross 2440 if MUE == None:
156 caltinay 3346 De=None
157 gross 2440 else:
158 caltinay 3346 De=ETA_N0/MUE/t_REF
159 gross 2440 ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0
160 jfenwick 3774 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 gross 2440 if MUE == None:
170 caltinay 3346 print("\tDebora number bottom \t\t\t= None")
171 gross 2440 else:
172 jfenwick 3774 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 caltinay 3346
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 gross 2440 #=========================
184 gross 2438 #
185 caltinay 3346 # set up domain and data from scratch or from restart files
186 gross 2438 #
187 caltinay 3346 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 gross 2440
202 caltinay 3346 #diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)
203 jfenwick 3774 print(("<%s> Restart at time step %d (t=%e) completed."%(time.asctime(),n,t)))
204 gross 1417 else:
205 caltinay 3346 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 gross 1417
217 caltinay 3346 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 gross 1472
223 caltinay 3346 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 gross 2620
228 gross 2445 p_last=p
229 gross 1417 x=dom.getX()
230     #
231     # set up heat problem:
232     #
233 gross 3793 heat=TemperatureCartesian(dom)
234 jfenwick 3774 print(("<%s> Temperature transport has been set up."%time.asctime()))
235 gross 2474 heat.getSolverOptions().setTolerance(T_TOL)
236     heat.getSolverOptions().setVerbosity(VERBOSE)
237 gross 1417 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
238 gross 2440 heat.setInitialTemperature(clip(T,T_0))
239     heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
240 gross 1417 #
241 gross 1639 # velocity constraints:
242 gross 1417 #
243     fixed_v_mask=Vector(0,Solution(dom))
244 gross 2620 faces=Scalar(0.,Solution(dom))
245 gross 1417 for d in range(DIM):
246     if d == DIM-1:
247 caltinay 3346 ll = H
248 gross 1417 else:
249 caltinay 3346 ll = L
250     if CREATE_TOPO and d==DIM-1:
251     fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
252 gross 2438 else:
253 caltinay 3346 s=whereZero(x[d])+whereZero(x[d]-ll)
254     faces+=s
255     fixed_v_mask+=s*unitVector(d,DIM)
256 gross 1417 #
257 caltinay 3346 # set up velocity problem
258 gross 1639 #
259 gross 2440 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 gross 2862 flow.setTolerance(FLOW_TOL)
264 gross 2440 flow.setEtaTolerance(FLOW_TOL)
265     flow.setExternals(fixed_v_mask=fixed_v_mask)
266 jfenwick 3774 print(("<%s> Flow solver has been set up."%time.asctime()))
267 gross 2440 #
268 caltinay 3346 # topography setup
269 gross 2620 #
270     boundary=FunctionOnBoundary(dom).getX()[DIM-1]
271     top_boundary_mask=whereZero(boundary-sup(boundary))
272     surface_area=integrate(top_boundary_mask)
273 caltinay 3346 if CREATE_TOPO:
274 gross 2857 mts=Mountains(dom,eps=TOPO_SMOOTH)
275 gross 2620 mts.setTopography(topography)
276 jfenwick 3774 print(("<%s> topography has been set up."%time.asctime()))
277 caltinay 3346
278 gross 2620 #
279 caltinay 3346 # let the show begin:
280 gross 1417 #
281 gross 2438 t1 = time.time()
282 jfenwick 3774 print(("<%s> Start time step %d (t=%s)."%(time.asctime(),n,t)))
283 gross 1417 while t<T_END:
284 caltinay 3346 if CREATE_TOPO: topography_old=topography
285 gross 2620 v_old, p_old, stress_old=v, p, stress
286     T_old=T
287 caltinay 3346 #========================= solve for velocity ============================
288 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))
289 caltinay 3346 #print("viscosity range :", inf(eta_N), sup(eta_N))
290 gross 2440 flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0], [1,N])
291 gross 2620 flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM))
292 caltinay 3346 # if dt<=0 or not CREATE_TOPO:
293     if not CREATE_TOPO:
294     flow.update(dt, iter_max=100, verbose=False)
295 gross 2620 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 gross 2862 flow.update(dt, iter_max=100, verbose=False)
303 gross 2620 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 caltinay 3346 #print("topography update step %d error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v))
310 gross 2620 i+=1
311     if i > TOPO_ITER_MAX:
312 caltinay 3346 raise RuntimeError("topography did not converge after %d steps."%TOPO_ITER_MAX)
313 gross 2440 v=flow.getVelocity()
314 caltinay 3346 #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 gross 1417 n+=1
319     t+=dt
320 caltinay 3346 #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 gross 2440 heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
325     dt=heat.getSafeTimeStepSize()
326 caltinay 3346 #print("<%s> New time step size is %e"%(time.asctime(),dt))
327     #=========================== setup topography ============================
328     if CREATE_TOPO:
329 gross 2620 dt=min(mts.getSafeTimeStepSize()*0.5,dt)
330 caltinay 3346 #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 gross 2440 #
333     # solve temperature:
334     #
335 gross 2474 T=heat.getSolution(dt)
336 caltinay 3346 #print("Temperature range ",inf(T),sup(T))
337     #print("<%s> temperature update completed."%time.asctime())
338     #============================== analysis =================================
339 gross 2440 #
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 caltinay 3346 #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 gross 2620 #
361 caltinay 3346 # create restart/visualization files:
362 gross 2620 #
363 caltinay 3346 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 jfenwick 3774 print(("<%s> Cycle %d (time %s) exported."%(time.asctime(),dataMgr.getCycle(),t)))
370 gross 1474
371 jfenwick 3774 print(("<%s> Calculation finalized after %s seconds."%(time.asctime(),time.time()-t1)))
372 caltinay 3346

  ViewVC Help
Powered by ViewVC 1.1.26