# Contents of /trunk/finley/test/python/convection.py

Revision 2856 - (show annotations)
Mon Jan 18 04:14:37 2010 UTC (9 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 18390 byte(s)
```FunctionSpaces provide now some information about their approximation order.
```
 1 ######################################################## 2 # 3 # Copyright (c) 2003-2009 by University of Queensland 4 # Earth Systems Science Computational Center (ESSCC) 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 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-2009 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 24 __url__= 25 26 from esys.escript import * 27 from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian, Mountains, SubSteppingException 28 from esys.finley import Rectangle, Brick, LoadMesh 29 from optparse import OptionParser 30 from math import pi, ceil 31 import sys 32 import time 33 34 # ======================= Default Values ================================================== 35 DIM=2 # spatial dimension 36 H=1. # height 37 L=H # length 38 NE=40 # number of elements in H-direction. 39 PERT=0.05 # initial temperature perturbation 40 DT=1.e-4 # initial time step size 41 CREATE_TOPOGRAPHY=False # create topgraphy 42 DT_MIN=1.e-10 # minumum time step size 43 T_END=10. # end time 44 45 RHO_0=100. # surface density (lf ~ RHO_0**2) 46 G=1. # gravitational constant 47 ALPHA_0=0.1 # thermal expansion coefficient (Ra ~ RHO_0**2 * ALPHA_0 = lf * ALPHA_0) 48 T_0=0. # surface temperature 49 T_1=1. # bottom temperature 50 C_P=1 # heat capacity 51 K=1. # thermal conductivity 52 CHI=0. # Taylor Quinny coefficient 53 MUE=None # elastic shear modulus 54 TAU_Y=5*10**(2.5) # Drucker-Prager cohesion factor 55 BETA=0 # Drucker-Prager friction factor 56 TAU_0=2*10**(2.5) # transition stress 57 N=3 # power for power law 58 59 E=23 # activation energy 60 V=18*0 # activation volume 61 T_OFFSET=1 # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0) 62 R=1 # gas constant 63 ETA_N0=1. # viscosity at surface 64 65 TOPO_SMOOTH=1e-5 # smoothing factor of extrapolation of surface velocity to interior 66 67 T_TOL=1.e-4 # tolerance temperature transport 68 TOL2=0.1e-3 # tolerance for flow update at a timestep. (large value will do only one correction step) 69 FLOW_TOL=1.e-6 # tolerance for inconcompressible flow solver 70 TOPO_TOL=1.e-6 # tolerance for update of topography 71 DIAGNOSTICS_FN="diagnostics.csv" 72 VERBOSE=True 73 DT_VIS=T_END/500 # time distane between two visulaization files 74 DN_VIS=5 # maximum counter increment between two visulaization files 75 VIS_DIR="results" # name of the director for vis files 76 DN_RESTART=1000000 # create a restart file every DN_RESTART step 77 PREFIX_RESTART="restart" # name prefix for restart directories 78 TOPO_ITER_MAX=20 # maximum number of iteration steps to update topography 79 # ========================================================================================= 80 81 def removeRestartDirectory(dir_name): 82 if dom.onMasterProcessor() and os.path.isdir(dir_name): 83 for root, dirs, files in os.walk(dir_name, topdown=False): 84 for name in files: os.remove(os.path.join(root,name)) 85 for name in dirs: os.remove(os.path.join(root,name)) 86 os.rmdir(dir_name) 87 print "Restart files %s have been removed."%dir_name 88 dom.MPIBarrier() 89 90 # ========================================================================================= 91 # 92 # read options: 93 # 94 parser = OptionParser(usage="%prog [Options]") 95 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") 96 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) 97 parser.add_option("-p", "--param", dest="param", help="name of file to be imported ",metavar="PARAM", default=None) 98 (options, args) = parser.parse_args() 99 restart=options.restart or (options.restart_dir !=None) 100 # 101 # overwrite the default options: 102 # 103 print "<%s> Execution started."%time.asctime() 104 if options.param !=None: 105 exec(open(options.param,'r')) 106 print "Parameters are imported from file ",options.param 107 108 print "Input Parameters:" 109 print "\tDimension DIM\t\t=\t",DIM 110 print "\tHeight H\t\t\t=\t",H 111 print "\tLength L\t\t\t=\t",L 112 print "\telements in H NE\t\t=\t",NE 113 print "\ttemperature perturbation PERT\t\t=\t",PERT 114 print "\tinitial time step size DT\t\t=\t",DT 115 print "\tminimum time step size DT_MIN\t\t=\t",DT_MIN 116 print "\tend time T_END\t\t=\t",T_END 117 print "\tcreate topgraphy CREATE_TOPOGRAPHY\t=\t",CREATE_TOPOGRAPHY 118 print "\tsurface density RHO_0\t\t=\t",RHO_0 119 print "\tgravitational constant G\t\t\t=\t",G 120 print "\tthermal expansion coefficient ALPHA_0\t\t=\t",ALPHA_0 121 print "\tsurface temperature T_0\t\t=\t",T_0 122 print "\tbottom temperature T_1\t\t=\t",T_1 123 print "\theat capacity C_P\t\t=\t",C_P 124 print "\tthermal conductivity K\t\t\t=\t",K 125 print "\tTaylor-Qinny coefficient CHI\t\t=\t",CHI 126 print "\telastic shear modulus MUE\t\t=\t",MUE 127 print "\tcohesion factor TAU_Y\t\t=\t",TAU_Y 128 print "\tfriction factor BETA\t\t=\t",BETA 129 print "\ttransition stress TAU_0\t\t=\t",TAU_0 130 print "\tpower for power law N\t\t\t=\t",N 131 print "\tviscosity at surface ETA_N0\t\t=\t",ETA_N0 132 print "\tactivation energy E\t\t\t=\t",E 133 print "\tactivation volume V\t\t\t=\t",V 134 print "\ttemperature offset T_OFFSET\t\t=\t",T_OFFSET 135 print "\tgas constant R\t\t\t=\t",R 136 print "\ttopography smoothing TOPO_SMOOTH\t=\t",TOPO_SMOOTH 137 138 print "\ttolerance for topography TOPO_TOL\t\t=\t",TOPO_TOL 139 print "\ttransport tolerance T_TOL\t\t=\t",T_TOL 140 print "\ttolerance for flow updates TOL2\t\t=\t",TOL2 141 print "\tflow tolerance FLOW_TOL\t\t=\t",FLOW_TOL 142 print "\tfile for diagnostics DIAGNOSTICS_FN\t=\t",DIAGNOSTICS_FN 143 print "\tmin. time incr. for vis file DT_VIS\t\t=\t",DT_VIS 144 print "\tmin. count incr. for vis file DN_VIS\t\t=\t",DN_VIS 145 print "\tvisualization dir VIS_DIR\t\t=\t",VIS_DIR 146 print "\trestart counter incerement DN_RESTART\t=\t",DN_RESTART 147 print "\tprefix for restart dirs PREFIX_RESTART\t=\t",PREFIX_RESTART 148 print "\tverbosity VERBOSE\t\t=\t",VERBOSE 149 150 print "Control Parameters:" 151 t_REF=RHO_0*C_P*H**2/K 152 P_REF=ETA_N0/t_REF 153 Ar=E/R/(T_1-T_0) 154 if E>0 and V>0: 155 V_REF=P_REF*V/E 156 else: 157 V_REF=0 158 T_OFFSET_REF=T_OFFSET/(T_1-T_0) 159 Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF 160 Di=ALPHA_0*G*H/C_P 161 CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2) 162 if CREATE_TOPOGRAPHY: 163 SURFACE_LOAD=RHO_0*G*H/P_REF 164 else: 165 SURFACE_LOAD=0. 166 if MUE == None: 167 De=None 168 else: 169 De=ETA_N0/MUE/t_REF 170 ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0 171 print "\ttotal #element \t\t\t=\t",NE**DIM*int(L/H)**(DIM-1) 172 print "\treference time t_REF\t\t=\t%e"%t_REF 173 print "\treference pressure P_REF\t\t=\t%e"%P_REF 174 print "\treference Taylor-Qinny CHI_REF\t\t=\t%e"%CHI_REF 175 print "\tDissipation number DI\t\t=\t%e"%Di 176 print "\tRayleigh number surface Ra\t\t=\t%e"%Ra 177 if MUE == None: 178 print "\tDebora number surface De\t\t=\t",None 179 else: 180 print "\tDebora number surface De\t\t=\t%e"%De 181 182 print "\tBottom viscosity \t\t\t=\t%e"%ETA_BOT 183 print "\tRayleigh number bottom \t\t\t=\t%e"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT) 184 if MUE == None: 185 print "\tDebora number bottom \t\t\t=\t",None 186 else: 187 print "\tDebora number bottom \t\t\t=\t%d"%(ETA_BOT/MUE/t_REF) 188 print "\tArrhenius Ar\t\t=\t%e"%Ar 189 print "\tsurface load factor SURFACE_LOAD\t=\t%e"%SURFACE_LOAD 190 print "\tscaled activation volume V_REF\t\t=\t%e"%V_REF 191 # some control variables (will be overwritten in case of a restart: 192 # 193 t=0 # time stamp 194 n=0 # time step counter 195 dt=DT # current time step size 196 t_vis=0 197 n_vis=0 198 counter_vis=0 199 mkDir(VIS_DIR) 200 #========================= 201 # 202 # set up domain or read restart file 203 # 204 if restart: 205 if options.restart_dir ==None: 206 restart_files=[] 207 for f in os.listdir("."): 208 if f.startswith(PREFIX_RESTART): restart_files.append(f) 209 if len(restart_files)==0: 210 raise IOError,"no restart files" 211 restart_files.sort() 212 f=restart_files[-1] 213 else: 214 f=options.restart_dir 215 try: 216 dom=LoadMesh("mesh.nc") 217 except: 218 pass 219 FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";") 220 t=float(FF[0]) # time stamp 221 n=int(FF[3]) # time step counter 222 t_vis=float(FF[1]) # 223 n_vis=int(FF[2]) 224 counter_vis=int(FF[4]) 225 dt=float(FF[5]) 226 stress=load(os.path.join(f,"stress.nc"),dom) 227 v=load(os.path.join(f,"v.nc"),dom) 228 p=load(os.path.join(f,"p.nc"),dom) 229 T=load(os.path.join(f,"T.nc"),dom) 230 if CREATE_TOPOGRAPHY: 231 topography=load(os.path.join(f,"topo.nc"),dom) 232 233 diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True) 234 print "<%s> Restart from file %s at time step %s (t=%s) completed."%(time.asctime(),f,t,n) 235 else: 236 if DIM==2: 237 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=-1,optimize=True) 238 else: 239 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) 240 try: 241 dom.dump("mesh.nc") 242 except: 243 pass 244 x=dom.getX() 245 T=Scalar(1,ReducedSolution(dom)) 246 for d in range(DIM): 247 if d == DIM-1: 248 T*=sin(x[d]/H*pi) 249 else: 250 T*=cos(x[d]/L*pi) 251 252 T=(1.-x[DIM-1])+PERT*T 253 v=Vector(0,Solution(dom)) 254 stress=Tensor(0,Function(dom)) 255 x2=ReducedSolution(dom).getX() 256 p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5) 257 258 if CREATE_TOPOGRAPHY: 259 topography=Scalar(0.,ReducedSolution(dom)) 260 diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=False) 261 diagnostics_file.write("Ra = %e Lambda= %e\n"%(Ra, SURFACE_LOAD)) 262 263 p_last=p 264 x=dom.getX() 265 # 266 # set up heat problem: 267 # 268 heat=TemperatureCartesian(dom,useBackwardEuler=False) 269 print "<%s> Temperature transport has been set up."%time.asctime() 270 heat.getSolverOptions().setTolerance(T_TOL) 271 heat.getSolverOptions().setVerbosity(VERBOSE) 272 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1]) 273 heat.setInitialTemperature(clip(T,T_0)) 274 heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at) 275 # 276 # velocity constraints: 277 # 278 x2=ReducedSolution(dom).getX() 279 fixed_v_mask=Vector(0,Solution(dom)) 280 faces=Scalar(0.,Solution(dom)) 281 for d in range(DIM): 282 if d == DIM-1: 283 ll = H 284 else: 285 ll = L 286 if CREATE_TOPOGRAPHY and d==DIM-1: 287 fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM) 288 else: 289 s=whereZero(x[d])+whereZero(x[d]-ll) 290 faces+=s 291 fixed_v_mask+=s*unitVector(d,DIM) 292 # 293 # set up velovity problem 294 # 295 flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE) 296 flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1])) 297 298 flow.setElasticShearModulus(MUE) 299 flow.setTolerance(TOL2) 300 flow.setEtaTolerance(FLOW_TOL) 301 flow.setExternals(fixed_v_mask=fixed_v_mask) 302 print "<%s> Flow solver has been set up."%time.asctime() 303 # 304 # topography set-up 305 # 306 boundary=FunctionOnBoundary(dom).getX()[DIM-1] 307 top_boundary_mask=whereZero(boundary-sup(boundary)) 308 surface_area=integrate(top_boundary_mask) 309 if CREATE_TOPOGRAPHY: 310 mts=Mountains(dom,eps=TOPO_SMOOTH,reduced=True) 311 mts.setTopography(topography) 312 print "<%s> topography has been set up."%time.asctime() 313 # 314 # let the show begin: 315 # 316 t1 = time.time() 317 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t) 318 while t TOPO_TOL * Topo_norm: 336 flow.setStatus(t, v_old, p_old, stress_old) 337 flow.setExternals(f=-SURFACE_LOAD*(topography-dt*v)*unitVector(DIM-1,DIM)*top_boundary_mask, restoration_factor=SURFACE_LOAD*dt*top_boundary_mask) 338 flow.update(dt, iter_max=100, inner_iter_max=200, verbose=False) 339 v=flow.getVelocity() 340 mts.setTopography(topography_old) 341 mts.setVelocity(v) 342 topography_last, topography=topography, mts.update(dt, allow_substeps=True) 343 error_Topo=sqrt(integrate(((topography-topography_last)*top_boundary_mask)**2)) 344 Topo_norm=sqrt(integrate((topography*top_boundary_mask)**2)) 345 print "DDDDD :", "input=",integrate(v*top_boundary_mask)[DIM-1],"output=", integrate(topography*top_boundary_mask)/Lsup(topography), error_Topo, Topo_norm 346 print "topography update step %s error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v) 347 i+=1 348 if i > TOPO_ITER_MAX: 349 raise RuntimeError,"topography did not converge after %s steps."%TOPO_ITER_MAX 350 v=flow.getVelocity() 351 for d in range(DIM): 352 print "range %d-velocity"%d,inf(v[d]),sup(v[d]) 353 print "<%s> flow solver completed."%time.asctime() 354 n+=1 355 t+=dt 356 print "influx= ",integrate(inner(v,dom.getNormal())), sqrt(integrate(length(v)**2,FunctionOnBoundary(dom))), integrate(1., FunctionOnBoundary(dom)) 357 print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t) 358 #======= setup Temperature problem ==================================================================== 359 # 360 heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff()) 361 dt=heat.getSafeTimeStepSize() 362 print "<%s> New time step size is %e"%(time.asctime(),dt) 363 #======= set-up topography ================================================================================== 364 if CREATE_TOPOGRAPHY: 365 dt=min(mts.getSafeTimeStepSize()*0.5,dt) 366 print "<%s> New time step size is %e"%(time.asctime(),dt) 367 print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt) 368 # 369 # solve temperature: 370 # 371 T=heat.getSolution(dt) 372 print "Temperature range ",inf(T),sup(T) 373 print "<%s> temperature update completed."%time.asctime() 374 #======= anaysis ================================================================================== 375 # 376 # .... Nusselt number 377 # 378 dTdz=grad(T)[DIM-1] 379 Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz) 380 eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff()) 381 Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar 382 print "nusselt number = %e"%Nu 383 print "av. eta = %e"%eta_bar 384 print "effective Rayleigh number = %e"%Ra_eff 385 if CREATE_TOPOGRAPHY: 386 topo_level=integrate(topography*top_boundary_mask)/surface_area 387 valleys_deep=inf(topography) 388 mountains_heigh=sup(topography) 389 print "topography level = ",topo_level 390 print "valleys deep = ", valleys_deep 391 print "mountains_heigh = ", mountains_heigh 392 diagnostics_file.write("%e %e %e %e %e %e %e\n"%(t,Nu, topo_level, valleys_deep, mountains_heigh, eta_bar, Ra_eff)) 393 else: 394 diagnostics_file.write("%e %e %e %e\n"%(t,Nu, eta_bar, Ra_eff)) 395 # 396 # .... visualization 397 # 398 if t>=t_vis or n>n_vis: 399 if CREATE_TOPOGRAPHY: 400 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff(), topography=topography_old, vex=mts.getVelocity()) 401 else: 402 saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff()) 403 print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t) 404 counter_vis+=1 405 t_vis+=DT_VIS 406 n_vis+=DN_VIS 407 # ========================= 408 # 409 # create restart files: 410 # 411 if DN_RESTART>0: 412 if (n-1)%DN_RESTART == 0: 413 c=(n-1)/DN_RESTART 414 old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1) 415 new_restart_dir="%s_%s_"%(PREFIX_RESTART,c) 416 mkDir(new_restart_dir) 417 dom.MPIBarrier() 418 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_vis, n_vis, n, counter_vis, dt)) 419 v.dump(os.path.join(new_restart_dir,"v.nc")) 420 p.dump(os.path.join(new_restart_dir,"p.nc")) 421 T.dump(os.path.join(new_restart_dir,"T.nc")) 422 if CREATE_TOPOGRAPHY: topography.dump(os.path.join(new_restart_dir,"topo.nc")) 423 removeRestartDirectory(old_restart_dir) 424 print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir) 425 print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1) 426

 ViewVC Help Powered by ViewVC 1.1.26