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

revision 1462 by gross, Mon Feb 25 04:45:48 2008 UTC revision 1463 by gross, Tue Apr 1 05:20:57 2008 UTC
# Line 2  from esys.escript import * Line 2  from esys.escript import *
2  from esys.escript.models import TemperatureCartesian, StokesProblemCartesian  from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
3  from esys.finley import Rectangle, Brick  from esys.finley import Rectangle, Brick
4  from math import pi, ceil  from math import pi, ceil
5  NE=20  NE=40
6  DIM=2  DIM=3
7  H=1.  H=1.
8  L=4.  L=4*H
9  THETA=0.5  THETA=0.5
10  TOL=1.e-3  TOL=1.e-3
11  PERTURBATION=0.1  PERTURBATION=0.1
12  T_END=0.2  T_END=0.3
13  DT_OUT=T_END/500*1e99  DT_OUT=T_END/500
14  VERBOSE=False  VERBOSE=False
15  RA=1.e6 # Rayleigh number  RA=1.e5 # Rayleigh number
16  A=25.  # Arenious number  A=22.  # Arenious number
17  DI = 0.  # dissipation number  DI = 0.  # dissipation number
18  SUPG=True  SUPG=False
19    print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
20  #  #
21  #   set up domain:  #   set up domain:
22  #  #
# Line 38  for d in range(DIM): Line 39  for d in range(DIM):
39      if d == DIM-1:      if d == DIM-1:
40         T*=sin(x[d]/H*pi)         T*=sin(x[d]/H*pi)
41      else:      else:
42         T*=cos(x[d]/L*pi)         T*=cos(x[d]/H*pi)
43  T=1.-x[DIM-1]+PERTURBATION*T  T=1.-x[DIM-1]+PERTURBATION*T
44  heat.setInitialTemperature(T)  heat.setInitialTemperature(T)
45  print "initial Temperature range ",inf(T),sup(T)  print "initial Temperature range ",inf(T),sup(T)
# Line 72  dt=None Line 73  dt=None
73  dt_new=1.  dt_new=1.
74  a=0  a=0
75  if dom.getMPIRank() ==0:  if dom.getMPIRank() ==0:
76       if SUPG:       nusselt_file=open("nusselt.csv","w")
nusselt_file=open("nusselt_supg.csv","w")
else:
nusselt_file=open("nusselt.csv","w")
77  while t<T_END:  while t<T_END:
78      v_last=v*1.      v_last=v*1.
79      print "============== solve for v ========================"      print "============== solve for v ========================"
80      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-(2./3.)))      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-0.5))
81      print "viscosity range :", inf(viscosity), sup(viscosity)      print "viscosity range :", inf(viscosity), sup(viscosity)