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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3773 by caltinay, Fri Nov 12 01:19:02 2010 UTC revision 3774 by jfenwick, Wed Jan 18 06:29:34 2012 UTC
# Line 94  restart=options.restart Line 94  restart=options.restart
94  #  #
95  #  overwrite the default options:  #  overwrite the default options:
96  #  #
97  print("<%s> Execution started."%time.asctime())  print(("<%s> Execution started."%time.asctime()))
98  if options.param !=None:  if options.param !=None:
99      exec(open(options.param,'r'))      exec(open(options.param,'r'))
100      print("Parameters imported from file ",options.param)      print(("Parameters imported from file ",options.param))
101    
102  print("Input Parameters:")  print("Input Parameters:")
103  print("\tDimension                     DIM\t\t= %d"%DIM)  print(("\tDimension                     DIM\t\t= %d"%DIM))
104  print("\tHeight                        H\t\t\t= %s"%H)  print(("\tHeight                        H\t\t\t= %s"%H))
105  print("\tLength                        L\t\t\t= %s"%L)  print(("\tLength                        L\t\t\t= %s"%L))
106  print("\tElements in H                 NE\t\t= %d"%NE)  print(("\tElements in H                 NE\t\t= %d"%NE))
107  print("\tTemperature perturbation      PERT\t\t= %s"%PERT)  print(("\tTemperature perturbation      PERT\t\t= %s"%PERT))
108  print("\tInitial time step size        DT\t\t= %s"%DT)  print(("\tInitial time step size        DT\t\t= %s"%DT))
109  print("\tMinimum time step size        DT_MIN\t\t= %s"%DT_MIN)  print(("\tMinimum time step size        DT_MIN\t\t= %s"%DT_MIN))
110  print("\tEnd time                      T_END\t\t= %s"%T_END)  print(("\tEnd time                      T_END\t\t= %s"%T_END))
111  print("\tCreate topography             CREATE_TOPO\t= %s"%CREATE_TOPO)  print(("\tCreate topography             CREATE_TOPO\t= %s"%CREATE_TOPO))
112  print("\tSurface density               RHO_0\t\t= %s"%RHO_0)  print(("\tSurface density               RHO_0\t\t= %s"%RHO_0))
113  print("\tGravitational constant        G\t\t\t= %s"%G)  print(("\tGravitational constant        G\t\t\t= %s"%G))
114  print("\tThermal expansion coefficient ALPHA_0\t\t= %s"%ALPHA_0)  print(("\tThermal expansion coefficient ALPHA_0\t\t= %s"%ALPHA_0))
115  print("\tSurface temperature           T_0\t\t= %s"%T_0)  print(("\tSurface temperature           T_0\t\t= %s"%T_0))
116  print("\tBottom temperature            T_1\t\t= %s"%T_1)  print(("\tBottom temperature            T_1\t\t= %s"%T_1))
117  print("\tHeat capacity                 C_P\t\t= %s"%C_P)  print(("\tHeat capacity                 C_P\t\t= %s"%C_P))
118  print("\tThermal conductivity          K\t\t\t= %s"%K)  print(("\tThermal conductivity          K\t\t\t= %s"%K))
119  print("\tTaylor-Quinny coefficient     CHI\t\t= %s"%CHI)  print(("\tTaylor-Quinny coefficient     CHI\t\t= %s"%CHI))
120  print("\tElastic shear modulus         MUE\t\t= %s"%MUE)  print(("\tElastic shear modulus         MUE\t\t= %s"%MUE))
121  print("\tCohesion factor               TAU_Y\t\t= %s"%TAU_Y)  print(("\tCohesion factor               TAU_Y\t\t= %s"%TAU_Y))
122  print("\tFriction factor               BETA\t\t= %s"%BETA)  print(("\tFriction factor               BETA\t\t= %s"%BETA))
123  print("\tTransition stress             TAU_0\t\t= %s"%TAU_0)  print(("\tTransition stress             TAU_0\t\t= %s"%TAU_0))
124  print("\tPower for power law           N\t\t\t= %s"%N)  print(("\tPower for power law           N\t\t\t= %s"%N))
125  print("\tViscosity at surface          ETA_N0\t\t= %s"%ETA_N0)  print(("\tViscosity at surface          ETA_N0\t\t= %s"%ETA_N0))
126  print("\tActivation energy             E\t\t\t= %s"%E)  print(("\tActivation energy             E\t\t\t= %s"%E))
127  print("\tActivation volume             V\t\t\t= %s"%V)  print(("\tActivation volume             V\t\t\t= %s"%V))
128  print("\tTemperature offset            T_OFFSET\t\t= %s"%T_OFFSET)  print(("\tTemperature offset            T_OFFSET\t\t= %s"%T_OFFSET))
129  print("\tGas constant                  R\t\t\t= %s"%R)  print(("\tGas constant                  R\t\t\t= %s"%R))
130  print("\tTopography smoothing          TOPO_SMOOTH\t= %s"%TOPO_SMOOTH)  print(("\tTopography smoothing          TOPO_SMOOTH\t= %s"%TOPO_SMOOTH))
131  print("\tTolerance for topography      TOPO_TOL\t\t= %s"%TOPO_TOL)  print(("\tTolerance for topography      TOPO_TOL\t\t= %s"%TOPO_TOL))
132  print("\tTransport tolerance           T_TOL\t\t= %s"%T_TOL)  print(("\tTransport tolerance           T_TOL\t\t= %s"%T_TOL))
133  print("\tFlow tolerance                FLOW_TOL\t\t= %s"%FLOW_TOL)  print(("\tFlow tolerance                FLOW_TOL\t\t= %s"%FLOW_TOL))
134  #print("\tFile for diagnostics          DIAGNOSTICS_FN\t= %s"%DIAGNOSTICS_FN)  #print("\tFile for diagnostics          DIAGNOSTICS_FN\t= %s"%DIAGNOSTICS_FN)
135  print("\tRestart counter increment     DN_RESTART\t= %d"%DN_RESTART)  print(("\tRestart counter increment     DN_RESTART\t= %d"%DN_RESTART))
136  print("\tPrefix for restart dirs       PREFIX_RESTART\t= %s"%PREFIX_RESTART)  print(("\tPrefix for restart dirs       PREFIX_RESTART\t= %s"%PREFIX_RESTART))
137  print("\tVerbosity                     VERBOSE\t\t= %s"%VERBOSE)  print(("\tVerbosity                     VERBOSE\t\t= %s"%VERBOSE))
138    
139  print("Control Parameters:")  print("Control Parameters:")
140  t_REF=RHO_0*C_P*H**2/K  t_REF=RHO_0*C_P*H**2/K
# Line 157  if MUE == None: Line 157  if MUE == None:
157  else:  else:
158      De=ETA_N0/MUE/t_REF      De=ETA_N0/MUE/t_REF
159  ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0  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)))  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)  print(("\tReference time                t_REF\t\t= %s"%t_REF))
162  print("\tReference pressure            P_REF\t\t= %s"%P_REF)  print(("\tReference pressure            P_REF\t\t= %s"%P_REF))
163  print("\tReference Taylor-Quinny       CHI_REF\t\t= %s"%CHI_REF)  print(("\tReference Taylor-Quinny       CHI_REF\t\t= %s"%CHI_REF))
164  print("\tDissipation number            DI\t\t= %s"%Di)  print(("\tDissipation number            DI\t\t= %s"%Di))
165  print("\tRayleigh number surface       Ra\t\t= %s"%Ra)  print(("\tRayleigh number surface       Ra\t\t= %s"%Ra))
166  print("\tDebora number surface         De\t\t= %s"%De)  print(("\tDebora number surface         De\t\t= %s"%De))
167  print("\tBottom viscosity              \t\t\t= %s"%ETA_BOT)  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))  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:  if MUE == None:
170     print("\tDebora number bottom          \t\t\t= None")     print("\tDebora number bottom          \t\t\t= None")
171  else:  else:
172     print("\tDebora number bottom          \t\t\t= %s"%(ETA_BOT/MUE/t_REF))     print(("\tDebora number bottom          \t\t\t= %s"%(ETA_BOT/MUE/t_REF)))
173  print("\tArrhenius                     Ar\t\t= %s"%Ar)  print(("\tArrhenius                     Ar\t\t= %s"%Ar))
174  print("\tSurface load factor           SURFACE_LOAD\t= %s"%SURFACE_LOAD)  print(("\tSurface load factor           SURFACE_LOAD\t= %s"%SURFACE_LOAD))
175  print("\tScaled activation volume      V_REF\t\t= %s"%V_REF)  print(("\tScaled activation volume      V_REF\t\t= %s"%V_REF))
176  print  print()
177    
178  # some control variables (will be overwritten in case of a restart:  # some control variables (will be overwritten in case of a restart:
179  t=0          # time stamp  t=0          # time stamp
# Line 200  if dataMgr.hasData(): Line 200  if dataMgr.hasData():
200          topography=dataMgr.getValue('topography')          topography=dataMgr.getValue('topography')
201        
202      #diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)      #diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)
203      print("<%s> Restart at time step %d (t=%e) completed."%(time.asctime(),n,t))      print(("<%s> Restart at time step %d (t=%e) completed."%(time.asctime(),n,t)))
204  else:  else:
205      if DIM==2:      if DIM==2:
206          dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=-1,optimize=True)          dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=-1,optimize=True)
# Line 231  x=dom.getX() Line 231  x=dom.getX()
231  #   set up heat problem:  #   set up heat problem:
232  #  #
233  heat=TemperatureCartesian(dom,useBackwardEuler=False)  heat=TemperatureCartesian(dom,useBackwardEuler=False)
234  print("<%s> Temperature transport has been set up."%time.asctime())  print(("<%s> Temperature transport has been set up."%time.asctime()))
235  heat.getSolverOptions().setTolerance(T_TOL)  heat.getSolverOptions().setTolerance(T_TOL)
236  heat.getSolverOptions().setVerbosity(VERBOSE)  heat.getSolverOptions().setVerbosity(VERBOSE)
237  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
# Line 263  flow.setElasticShearModulus(MUE) Line 263  flow.setElasticShearModulus(MUE)
263  flow.setTolerance(FLOW_TOL)  flow.setTolerance(FLOW_TOL)
264  flow.setEtaTolerance(FLOW_TOL)  flow.setEtaTolerance(FLOW_TOL)
265  flow.setExternals(fixed_v_mask=fixed_v_mask)  flow.setExternals(fixed_v_mask=fixed_v_mask)
266  print("<%s> Flow solver has been set up."%time.asctime())  print(("<%s> Flow solver has been set up."%time.asctime()))
267  #  #
268  #   topography setup  #   topography setup
269  #  #
# Line 273  surface_area=integrate(top_boundary_mask Line 273  surface_area=integrate(top_boundary_mask
273  if CREATE_TOPO:  if CREATE_TOPO:
274      mts=Mountains(dom,eps=TOPO_SMOOTH)      mts=Mountains(dom,eps=TOPO_SMOOTH)
275      mts.setTopography(topography)      mts.setTopography(topography)
276      print("<%s> topography has been set up."%time.asctime())      print(("<%s> topography has been set up."%time.asctime()))
277    
278  #  #
279  #   let the show begin:  #   let the show begin:
280  #  #
281  t1 = time.time()  t1 = time.time()
282  print("<%s> Start time step %d (t=%s)."%(time.asctime(),n,t))  print(("<%s> Start time step %d (t=%s)."%(time.asctime(),n,t)))
283  while t<T_END:  while t<T_END:
284      if CREATE_TOPO: topography_old=topography      if CREATE_TOPO: topography_old=topography
285      v_old, p_old, stress_old=v, p, stress      v_old, p_old, stress_old=v, p, stress
# Line 366  while t<T_END: Line 366  while t<T_END:
366          dataMgr.addData(t=t,n=n,t_vis=t_vis,dt=dt,T=T,v=v,eta=flow.getCurrentEtaEff(),stress=stress,p=p)          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)          if CREATE_TOPO: dataMgr.addData(topography=topography)
368          dataMgr.export()          dataMgr.export()
369          print("<%s> Cycle %d (time %s) exported."%(time.asctime(),dataMgr.getCycle(),t))          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))  print(("<%s> Calculation finalized after %s seconds."%(time.asctime(),time.time()-t1)))
372    

Legend:
Removed from v.3773  
changed lines
  Added in v.3774

  ViewVC Help
Powered by ViewVC 1.1.26