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

Diff of /branches/windows_from_1456_trunk_1490_merged_in/finley/test/python/convection.py

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

revision 1521 by trankine, Mon Apr 14 06:43:51 2008 UTC revision 1522 by phornby, Tue Apr 22 06:13:48 2008 UTC
# Line 31  if solver!='PCG': Line 31  if solver!='PCG':
31  else:  else:
32   extratol=1   extratol=1
33    
34  DIM=2  DIM=3
35  H=1.  H=1.
36  L=4*H  L=2*H
37  THETA=0.5  THETA=0.5
38  TOL=1.e-3  TOL=1.e-3
39  PERTURBATION=0.1  PERTURBATION=0.1
# Line 41  T_END=0.3 Line 41  T_END=0.3
41  DT_OUT=T_END/500  DT_OUT=T_END/500
42  VERBOSE=False  VERBOSE=False
43  RA=1.e5 # Rayleigh number  RA=1.e5 # Rayleigh number
44  A=22.  # Arenious number  A=0.  # Arenious number
45  DI = 0.  # dissipation number  DI = 0.  # dissipation number
46  SUPG=False  SUPG=False
47  create_restartfiles_every_step=10  create_restartfiles_every_step=10
# Line 152  for d in range(DIM): Line 152  for d in range(DIM):
152  while t<T_END:  while t<T_END:
153      v_last=v*1      v_last=v*1
154      print "============== solve for v ========================"      print "============== solve for v ========================"
155      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.))      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
156      print "viscosity range :", inf(viscosity), sup(viscosity)      print "viscosity range :", inf(viscosity), sup(viscosity)
157      sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)      sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)
158      #v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG')      #v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG')
# Line 164  while t<T_END: Line 164  while t<T_END:
164           print "range %d-velocity"%d,inf(v[d]),sup(v[d])           print "range %d-velocity"%d,inf(v[d]),sup(v[d])
165    
166      if t>=t_out:      if t>=t_out:
167        saveVTK("state.%d.vtu"%n_out,T=T,v=v,p=p)        saveVTK("state.%d.vtu"%n_out,T=T,v=v)
168        print "visualization file %d for time step %e generated."%(n_out,t)        print "visualization file %d for time step %e generated."%(n_out,t)
169        n_out+=1        n_out+=1
170        t_out+=DT_OUT        t_out+=DT_OUT

Legend:
Removed from v.1521  
changed lines
  Added in v.1522

  ViewVC Help
Powered by ViewVC 1.1.26