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

revision 1552 by gross, Thu May 8 08:52:41 2008 UTC revision 1639 by gross, Mon Jul 14 08:55:25 2008 UTC
# Line 1  Line 1
1  from esys.escript import *  from esys.escript import *
2  from esys.escript.models import TemperatureCartesian, StokesProblemCartesian  from esys.escript.models import TemperatureCartesian, PlateMantelModel
3  from esys.finley import Rectangle, Brick, LoadMesh  from esys.finley import Rectangle, Brick, LoadMesh
4  from optparse import OptionParser  from optparse import OptionParser
5  from math import pi, ceil  from math import pi, ceil
6
7
8  def removeRestartDirectory(dir_name):  def removeRestartDirectory(dir_name):
9     if os.path.isdir(dir_name):     if os.path.isdir(dir_name):
10         for root, dirs, files in os.walk(dir_name, topdown=False):         for root, dirs, files in os.walk(dir_name, topdown=False):
# Line 21  if (len(sys.argv)>=3): Line 23  if (len(sys.argv)>=3):
23  else:  else:
24   NE=20   NE=20
25  NE=64  NE=64
26    NE=20
27
28  if (len(sys.argv)>=2):  if (len(sys.argv)>=2):
29   solver=sys.argv[1]   solver=sys.argv[1]
# Line 33  else: Line 36  else:
36   extratol=1   extratol=1
37  extratol=0.1  extratol=0.1
38
39  DIM=3
40    DIM=2
41  H=1.  H=1.
42  L=2*H  L=4*H
43  THETA=0.5  THETA=0.5 # time stepping THETA=0.5 cranck nicolson
44  TOL=1.e-3  TOL=1.e-4
45  PERTURBATION=0.1  PERTURBATION=0.1
46    DT=1.e-4
47    DT_MIN=1.e-5
48  T_END=0.1  T_END=0.1
49  DT_OUT=T_END/500  DT_OUT=T_END/500
50  Dn_OUT=2  Dn_OUT=2
51  VERBOSE=False  VERBOSE=True
RA=1.e6 # Rayleigh number
A=22 # Arenious number
DI = 0.  # dissipation number
52  SUPG=False  SUPG=False
53  create_restartfiles_every_step=10  create_restartfiles_every_step=10
54    if True:
55       # this is a simple linear Stokes model:
56       RA=1.e6 # Rayleigh number
57       A=22 # Arenious number
58       DI = 0.  # dissipation number
59       MU=None
60       ETA0=1.
61       TAU0=None
62       N=None
63       NPL=None
64       ETAP0=ETA0
65       TAUY=None
66       useJAUMANNSTRESS=False
67    else:
68       RA=1.e4 # Rayleigh number
69       A=22 # Arenious number
70       DI = 0.  # dissipation number
71       MU=1.e4
72       ETA0=1.
73       TAU0=0.866*10**2.5
74       N=3
75       NPL=14
76       TAUY=TAU0
77       ETAP0=ETA0
78       useJAUMANNSTRESS=True
79
80  print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)  print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
81
# Line 76  if restart: Line 104  if restart:
105     except:     except:
106        pass        pass
108     t=float(ff[0])     t=float(FF[0])
109     t_out=float(ff[1])     t_out=float(FF[1])
110     n_out=int(ff[2])     n_out=int(FF[2])
111     n=int(ff[3])     n=int(FF[3])
112     out_count=int(ff[4])     out_count=int(FF[4])
113     dt=float(ff[5])     dt=float(FF[5])
118     if n>1:     if n>1:
119        dt_a=float(ff[6])        dt_a=float(FF[6])
121     else:     else:
122        dt_a=None        dt_a=None
# Line 112  else: Line 141  else:
141
142    T=1.-x[DIM-1]+PERTURBATION*T    T=1.-x[DIM-1]+PERTURBATION*T
143    v=Vector(0,Solution(dom))    v=Vector(0,Solution(dom))
144      stress=Tensor(0,Function(dom))
145      x2=ReducedSolution(dom).getX()
146      p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
147    if dom.getMPIRank() ==0: nusselt_file=open("nusselt.csv","w")    if dom.getMPIRank() ==0: nusselt_file=open("nusselt.csv","w")
148    t=0    t=0
149    t_out=0    t_out=0
150    n_out=0    n_out=0
151    n=0    n=0
152    out_count=0    out_count=0
153    dt=None    dt=DT
154    a=None    a=None
155    dt_a=None    dt_a=None
156
157  vol=integrate(1.,Function(dom))  vol=integrate(1.,Function(dom))
158    p-=integrate(p)/vol
159  x=dom.getX()  x=dom.getX()
160  #  #
161  #   set up heat problem:  #   set up heat problem:
# Line 135  heat.setInitialTemperature(T) Line 168  heat.setInitialTemperature(T)
168  print "initial Temperature range ",inf(T),sup(T)  print "initial Temperature range ",inf(T),sup(T)
170  #  #
171  #   set up velovity problem  #   velocity constraints:
172  #  #
sp=StokesProblemCartesian(dom)
sp.setTolerance(TOL*extratol)
sp.setToleranceReductionFactor(TOL)
173  x2=ReducedSolution(dom).getX()  x2=ReducedSolution(dom).getX()
p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
p-=integrate(p)/vol

175  for d in range(DIM):  for d in range(DIM):
176      if d == DIM-1:      if d == DIM-1:
# Line 151  for d in range(DIM): Line 178  for d in range(DIM):
178      else:      else:
179         ll = L         ll = L
181    #
182    #   set up velovity problem
183    #
184    sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
185    sp.initialize(mu=MU, tau_0=TAU0, n=N, eta_Y=ETAP0, tau_Y=TAU0, n_Y=NPL, q=fixed_v_mask)
186    sp.setTolerance(TOL*extratol)
187    sp.setToleranceReductionFactor(TOL)
188  #  #
189  #  let the show begin:  #  let the show begin:
190  #  #
191  while t<T_END:  while t<T_END:
192      v_last=v*1      v_last=v*1
193      print "============== solve for v ========================"      print "============== solve for v ========================"
194      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))      FF=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
195      print "viscosity range :", inf(viscosity), sup(viscosity)      print "viscosity range :", inf(FF)*ETA0, sup(FF)*ETA0
197      v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG')      sp.update(dt,max_inner_iter=20, verbose=VERBOSE, show_details=False, tol=10., solver='PCG')
198      # v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='GMRES')      v=sp.getVelocity()
#v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='MINRES')
# v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver=solver)
199
200      for d in range(DIM):      for d in range(DIM):
201           print "range %d-velocity"%d,inf(v[d]),sup(v[d])           print "range %d-velocity"%d,inf(v[d]),sup(v[d])
# Line 174  while t<T_END: Line 205  while t<T_END:
205        out_count+=1        out_count+=1
206        t_out+=DT_OUT        t_out+=DT_OUT
207        n_out+=Dn_OUT        n_out+=Dn_OUT
208      Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)      # calculation of nusselt number:
209        se=sp.getStrainEnergy()
210        Nu=1.+integrate(se)/(RA*vol)
211      if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))      if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))
213      print "nusselt number = ",Nu,n      print "nusselt number = ",Nu
214      if n>0:      if n>0:
215          a,a_alt = (v_last-v)/dt, a          a,a_alt = (v_last-v)/dt, a
216          dt_a,dt_a_alt = dt, dt_a          dt_a,dt_a_alt = dt, dt_a
# Line 190  while t<T_END: Line 223  while t<T_END:
223         dt=min(dt_new,heat.getSafeTimeStepSize())         dt=min(dt_new,heat.getSafeTimeStepSize())
224      else:      else:
225         dt=heat.getSafeTimeStepSize()         dt=heat.getSafeTimeStepSize()
226        dt=max(DT_MIN,dt)
227      print n,". time step t=",t," step size ",dt      print n,". time step t=",t," step size ",dt
228      print "============== solve for T ========================"      print "============== solve for T ========================"
229      T=heat.solve(dt, verbose=VERBOSE)        T=heat.solve(dt, verbose=VERBOSE)
# Line 208  while t<T_END: Line 242  while t<T_END:
242
243           print "Write restart files to ",new_restart_dir           print "Write restart files to ",new_restart_dir
244           if not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)           if not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
245           v.dump(os.path.join(new_restart_dir,"v.nc"))           sp.getStress().dump(os.path.join(new_restart_dir,"stress.nc"))
246           p.dump(os.path.join(new_restart_dir,"p.nc"))           sp.getVelocity().dump(os.path.join(new_restart_dir,"v.nc"))
247             sp.getPressure().dump(os.path.join(new_restart_dir,"p.nc"))
248           T.dump(os.path.join(new_restart_dir,"T.nc"))           T.dump(os.path.join(new_restart_dir,"T.nc"))
249           if n>1:           if n>1:
250               file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e; %e;\n"%(t, t_out, n_out, n, out_count, dt, dt_a))               file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e; %e;\n"%(t, t_out, n_out, n, out_count, dt, dt_a))

Legend:
 Removed from v.1552 changed lines Added in v.1639