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

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

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

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:
104        dom=LoadMesh("mesh.nc")        dom=LoadMesh("mesh.nc")
105     except:     except:
106        pass        pass
107     ff=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")     FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
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])
114       stress=load(os.path.join(f,"stress.nc"),dom)
115     v=load(os.path.join(f,"v.nc"),dom)     v=load(os.path.join(f,"v.nc"),dom)
116     p=load(os.path.join(f,"p.nc"),dom)     p=load(os.path.join(f,"p.nc"),dom)
117     T=load(os.path.join(f,"T.nc"),dom)     T=load(os.path.join(f,"T.nc"),dom)
118     if n>1:     if n>1:
119        dt_a=float(ff[6])        dt_a=float(FF[6])
120        a=load(os.path.join(f,"a.nc"),dom)        a=load(os.path.join(f,"a.nc"),dom)
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)
169  heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)  heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
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  
   
174  fixed_v_mask=Vector(0,Solution(dom))  fixed_v_mask=Vector(0,Solution(dom))
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
180      fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)      fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
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
196      sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)      sp.initialize(eta_N=ETA0*FF, eta_0=ETA0*FF, F=T*(RA*unitVector(DIM-1,DIM)))
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))
212      heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*viscosity*length(symmetric(grad(v)))**2)      heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)
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

  ViewVC Help
Powered by ViewVC 1.1.26