/[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 2439 by gross, Tue May 26 02:26:52 2009 UTC revision 2440 by gross, Wed May 27 08:45:55 2009 UTC
# Line 12  Line 12 
12  """  """
13  this is a convection simulation over a domain [0,L] X [0,L] x [0,H]  this is a convection simulation over a domain [0,L] X [0,L] x [0,H]
14    
15    It is solved in dimensionless form
16    
17  """  """
18  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2008 by University of Queensland
19  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
# Line 22  http://www.opensource.org/licenses/osl-3 Line 24  http://www.opensource.org/licenses/osl-3
24  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
25    
26  from esys.escript import *  from esys.escript import *
27  from esys.escript.models import TemperatureCartesian, StokesProblemCartesian  from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian
28  from esys.finley import Rectangle, Brick, LoadMesh  from esys.finley import Rectangle, Brick, LoadMesh
29  from optparse import OptionParser  from optparse import OptionParser
30  from math import pi, ceil  from math import pi, ceil
# Line 32  import time Line 34  import time
34  # ======================= Default Values ==================================================  # ======================= Default Values ==================================================
35  DIM=2                           # spatial dimension  DIM=2                           # spatial dimension
36  H=1.                            # height  H=1.                            # height
37  L=4*H                           # length  L=H                           # length
38  NE=10                           # number of elements in H-direction.  NE=50                           # number of elements in H-direction.
39  PERTURBATION=0.1                # initial temperature perturbation  PERT=0.05               # initial temperature perturbation
40  DT=1.e-4                        # initial time step size  DT=1.e-7                        # initial time step size
41  USE_BACKWARD_EULER=False        # use backward Euler in time integartion of temperature otherwise Crank-Nicholson is used  CREATE_TOPOGRAPHY=False          # create topgraphy
42  CREATE_TOPGRAPHY=False          # create topgraphy  TOPOGRAPHY_LOAD=False        # considers the loading of the topography in the flow
43  TOL=1.e-4                       # tolerance  DT_MIN=1.e-10                    # minumum time step size
 DT_MIN=1.e-5                    # minumum time step size  
44  T_END=0.1                       # end time  T_END=0.1                       # end time
45    
46  DT_OUT=T_END/500  RHO_0=1.                        # surface density
47  Dn_OUT=2  G=1.                            # gravitational constant
48    ALPHA_0=4.e3                      # thermal expansion coefficient
49    T_0=0.                          # surface temperature
50    T_1=1.                          # bottom temperature
51    C_P=1                           # heat capacity
52    K=1.                            # thermal conductivity
53    CHI=0.                          # Taylor Quinny coefficient
54    MUE=None                        # elastic shear modulus
55    TAU_Y=3*10**(2.5)               # Drucker-Prager cohesion factor
56    BETA=0                          # Drucker-Prager friction factor
57    TAU_0=2*10**(2.5)               # transition stress
58    N=3                             # power for power law
59    
60    E=23                            # activation energy
61    V=18*0                        # activation volume
62    T_OFFSET=1                      # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0)
63    R=1                             # gas constant
64    ETA_N0=1.                       # viscosity at surface
65    
66    T_TOL=1.e-4                     # tolerance temperature transport
67    TOL2=0.1                        # tolerance for flow update at a timestep. (large value will do only one correction step)
68    FLOW_TOL=1.e-5                  # tolerance for inconcompressible flow solver
69    FLOW_SUB_TOL=1.e-10             # sub-tolerance for inconcompressible flow solver
70    NUSSELT_FN="nusselt.csv"
71  VERBOSE=True  VERBOSE=True
72  CREATE_RESTARTFILES_EVERY_STEP=10  DT_VIS=T_END/500                # time distane between two visulaization files
73  if True:  DN_VIS=5                        # maximum counter increment between two visulaization files
74     # this is a simple linear Stokes model:  VIS_DIR="results"               # name of the director for vis files
75     RA=1.e6 # Rayleigh number  DN_RESTART=1000000              # create a restart file every DN_RESTART step
76     A=22 # Arenious number  PREFIX_RESTART="restart"        # name prefix for restart directories
    DI = 0.  # dissipation number  
    MU=None  
    ETA0=1.  
    TAU0=None  
    N=None  
    NPL=None  
    ETAP0=ETA0  
    TAUY=None  
    useJAUMANNSTRESS=False  
    # this is a simple linear Stokes model:  
    RA=1.e5 # Rayleigh number  
    A=11 # Arenious number  
    DI = 0.  # dissipation number  
    MU=None  
    ETA0=1.  
    TAU0=250.  
    N=None  
    NPL=None  
    ETAP0=ETA0  
    TAUY=TAU0  
    useJAUMANNSTRESS=False  
 else:  
    RA=1.e4 # Rayleigh number  
    A=22 # Arenious number  
    DI = 0.  # dissipation number  
    MU=1.e4  
    ETA0=1.  
    TAU0=0.866*10**2.5  
    N=3  
    NPL=14  
    TAUY=TAU0  
    ETAP0=ETA0  
    useJAUMANNSTRESS=True  
   
   
77  # =========================================================================================  # =========================================================================================
78    
79  def removeRestartDirectory(dir_name):  def removeRestartDirectory(dir_name):
# Line 95  def removeRestartDirectory(dir_name): Line 84  def removeRestartDirectory(dir_name):
84         os.rmdir(dir_name)         os.rmdir(dir_name)
85         print "Restart files %s have been removed."%dir_name         print "Restart files %s have been removed."%dir_name
86     dom.MPIBarrier()     dom.MPIBarrier()
87    
88  # =========================================================================================  # =========================================================================================
89  #  #
90  # read options:  # read options:
# Line 108  restart=options.restart or (options.rest Line 98  restart=options.restart or (options.rest
98  #  #
99  #  overwrite the default options:  #  overwrite the default options:
100  #  #
101  print "Execution started ",time.asctime()  print "<%s> Execution started."%time.asctime()
102  if options.param !=None:  if options.param !=None:
103       exec(open(options.param,'r'))       exec(open(options.param,'r'))
104       print "Parameters are imported from file ",options.param       print "Parameters are imported from file ",options.param
105  print "Parameters:"  
106  print "\tDimension\tDIM=\t",DIM  CREATE_TOPOGRAPHY=CREATE_TOPOGRAPHY or TOPOGRAPHY_LOAD
107  print "\tHeight:\tH=\t",H  
108  print "\tLength:\tL=\t",L  print "Input Parameters:"
109  print "\telements in H:\tNE=\t",NE  print "\tDimension                     DIM\t\t=\t",DIM
110  print "\ttotal #element\t\t=\t",NE**DIM*int(L/H)**(DIM-1)  print "\tHeight                        H\t\t\t=\t",H
111  print "\ttolerance\tTOL=\t",TOL  print "\tLength                        L\t\t\t=\t",L
112  print "\ttemperature perturbation\tPERTURBATION=\t",PERTURBATION  print "\telements in H                 NE\t\t=\t",NE
113  print "\tinitial time step size\tDT=\t",DT  print "\ttemperature perturbation      PERT\t\t=\t",PERT
114  print "\tminimum time step size\tDT_MIN=\t",tDT_MIN  print "\tinitial time step size        DT\t\t=\t",DT
115  print "\tend time\tT_END=\t",T_END  print "\tminimum time step size        DT_MIN\t\t=\t",DT_MIN
116  print "\tbackward Euler?\tUSE_BACKWARD_EULER=\t",USE_BACKWARD_EULER  print "\tend time                      T_END\t\t=\t",T_END
117  print "\ttopgraphy?\tCREATE_TOPOGRAPHY=\t",CREATE_TOPOGRAPHY  print "\tcreate topgraphy              CREATE_TOPOGRAPHY\t=\t",CREATE_TOPOGRAPHY
118    print "\tloading by topography         TOPOGRAPHY_LOAD\t=\t",TOPOGRAPHY_LOAD
119    print "\tsurface density               RHO_0\t\t=\t",RHO_0
120    print "\tgravitational constant        G\t\t\t=\t",G
121    print "\tthermal expansion coefficient ALPHA_0\t\t=\t",ALPHA_0
122    print "\tsurface temperature           T_0\t\t=\t",T_0
123    print "\tbottom temperature            T_1\t\t=\t",T_1
124    print "\theat capacity                 C_P\t\t=\t",C_P
125    print "\tthermal conductivity          K\t\t\t=\t",K
126    print "\tTaylor-Qinny coefficient      CHI\t\t=\t",CHI
127    print "\telastic shear modulus         MUE\t\t=\t",MUE
128    print "\tcohesion factor               TAU_Y\t\t=\t",TAU_Y
129    print "\tfriction factor               BETA\t\t=\t",BETA
130    print "\ttransition stress             TAU_0\t\t=\t",TAU_0
131    print "\tpower for power law           N\t\t\t=\t",N
132    print "\tviscosity at surface          ETA_N0\t\t=\t",ETA_N0
133    print "\tactivation energy             E\t\t\t=\t",E
134    print "\tactivation volume             V\t\t\t=\t",V
135    print "\ttemperature offset            T_OFFSET\t\t=\t",T_OFFSET
136    print "\tgas constant                  R\t\t\t=\t",R
137    
138    print "\ttransport tolerance           T_TOL\t\t=\t",T_TOL
139    print "\ttolerance for flow updates    TOL2\t\t=\t",TOL2
140    print "\tflow tolerance                FLOW_TOL\t\t=\t",FLOW_TOL
141    print "\tflow sub-tolerance            FLOW_SUB_TOL\t=\t",FLOW_SUB_TOL
142    print "\tfile for Nusselt number       NUSSELT_FN\t=\t",NUSSELT_FN
143    print "\tmin. time incr. for vis file  DT_VIS\t\t=\t",DT_VIS
144    print "\tmin. count incr. for vis file DN_VIS\t\t=\t",DN_VIS
145    print "\tvisualization dir             VIS_DIR\t\t=\t",VIS_DIR
146    print "\trestart counter incerement    DN_RESTART\t\t=\t",DN_RESTART
147    print "\tprefix for restart dirs       PREFIX_RESTART\t\t=\t",PREFIX_RESTART
148    print "\tverbosity                     VERBOSE\t\t=\t",VERBOSE
149    
150    print "Control Parameters:"
151    t_REF=RHO_0*C_P*H**2/K
152    P_REF=ETA_N0/t_REF
153    Ar=E/R/(T_1-T_0)
154    if E>0 and V>0:
155       V_REF=P_REF*V/E
156    else:
157      V_REF=0
158    T_OFFSET_REF=T_OFFSET/(T_1-T_0)
159    Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF
160    Di=ALPHA_0*G*H/C_P
161    CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2)
162    if MUE == None:
163      De=None
164    else:
165      De=ETA_N0/MUE/t_REF
166    ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0
167    print "\ttotal #element                \t\t\t=\t",NE**DIM*int(L/H)**(DIM-1)
168    print "\treference time                t_REF\t\t=\t%e"%t_REF
169    print "\treference pressure            P_REF\t\t=\t%e"%P_REF
170    print "\treference Taylor-Qinny        CHI_REF\t\t=\t%e"%CHI_REF
171    print "\tDissipation number            DI\t\t=\t%e"%Di
172    print "\tRayleigh number surface       Ra\t\t=\t%e"%Ra
173    if MUE == None:
174       print "\tDebora number surface         De\t\t=\t",None
175    else:
176       print "\tDebora number surface         De\t\t=\t%e"%De
177        
178    print "\tbottom viscosity              \t\t\t=\t%e"%ETA_BOT
179    print "\tRayleigh number bottom        \t\t\t=\t%e"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT)
180    if MUE == None:
181       print "\tDebora number bottom          \t\t\t=\t",None
182    else:
183       print "\tDebora number bottom          \t\t\t=\t%d"%(ETA_BOT/MUE/t_REF)
184    print "\tArrhenius                     Ar\t\t=\t%e"%Ar
185    print "\tscaled activation volume      V_REF\t\t=\t%e"%V_REF
186  #  some control variables (will be overwritten in  case of a restart:  #  some control variables (will be overwritten in  case of a restart:
187  #  #
188  t=0         # time stamp  t=0         # time stamp
189  n=0         # time step counter  n=0         # time step counter
190  dt=DT       # current time step size  dt=DT       # current time step size
191    t_vis=0    
192  t_out=0     #  n_vis=0
193  n_out=0  counter_vis=0
194  out_count=0  if getMPIRankWorld()==0:
195  a=None      if not os.path.isdir(VIS_DIR): os.mkdir(VIS_DIR)
196  dt_a=None  #=========================
197  #  #
198  #   set up domain or read restart file  #   set up domain or read restart file
199  #  #
# Line 143  if restart: Line 201  if restart:
201     if options.restart_dir ==None:     if options.restart_dir ==None:
202        restart_files=[]        restart_files=[]
203        for f in os.listdir("."):        for f in os.listdir("."):
204            if f.startswith("restart"): restart_files.append(f)            if f.startswith(PREFIX_RESTART): restart_files.append(f)
205        if len(restart_files)==0:        if len(restart_files)==0:
206            raise IOError,"no restart files"            raise IOError,"no restart files"
207        restart_files.sort()        restart_files.sort()
208        f=restart_files[-1]        f=restart_files[-1]
209     else:     else:
210         f=options.restart_dir         f=options.restart_dir
    print ">>>Restart from directory ",f  
211     try:     try:
212        dom=LoadMesh("mesh.nc")        dom=LoadMesh("mesh.nc")
213     except:     except:
214        pass        pass
215     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(";")
216     t=float(FF[0])          # time stamp     t=float(FF[0])          # time stamp
217     t_out=float(FF[1])      #     n=int(FF[3])            # time step counter
218     n_out=int(FF[2])     t_vis=float(FF[1])      #
219     n=int(FF[3])     n_vis=int(FF[2])
220     out_count=int(FF[4])     counter_vis=int(FF[4])
221     dt=float(FF[5])     dt=float(FF[5])
222     stress=load(os.path.join(f,"stress.nc"),dom)     stress=load(os.path.join(f,"stress.nc"),dom)
223     v=load(os.path.join(f,"v.nc"),dom)     v=load(os.path.join(f,"v.nc"),dom)
224     p=load(os.path.join(f,"p.nc"),dom)     p=load(os.path.join(f,"p.nc"),dom)
225     T=load(os.path.join(f,"T.nc"),dom)     T=load(os.path.join(f,"T.nc"),dom)
226     if n>1:    
227        dt_a=float(FF[6])     nusselt_file=FileWriter(NUSSELT_FN,append=True)
228        a=load(os.path.join(f,"a.nc"),dom)     print "<%s> Restart from file %s at time step %s (t=%s) completed."%(time.asctime(),f,t,n)
    else:  
       dt_a=None  
       a=None  
    if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","a")  
229  else:  else:
230    if DIM==2:    if DIM==2:
231      dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)      dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=2, useFullElementOrder=True,optimize=True)
232    else:    else:
233      dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True,optimize=True)      dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L/H,l1=L/H,l2=1,order=2, useFullElementOrder=True,optimize=True)
234    try:    try:
235       dom.dump("mesh.nc")       dom.dump("mesh.nc")
236    except:    except:
# Line 190  else: Line 243  else:
243        else:        else:
244           T*=cos(x[d]/L*pi)           T*=cos(x[d]/L*pi)
245    
246    T=1.-x[DIM-1]+PERTURBATION*T    T=(1.-x[DIM-1])+PERT*T
247    v=Vector(0,Solution(dom))    v=Vector(0,Solution(dom))
248    stress=Tensor(0,Function(dom))    stress=Tensor(0,Function(dom))
249    x2=ReducedSolution(dom).getX()    x2=ReducedSolution(dom).getX()
250    p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)    p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)
251    if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","w")    nusselt_file=FileWriter(NUSSELT_FN,append=False)
252    
 vol=integrate(1.,Function(dom))  
 p-=integrate(p)/vol  
253  x=dom.getX()  x=dom.getX()
254  #  #
255  #   set up heat problem:  #   set up heat problem:
256  #  #
257  heat=TemperatureCartesian(dom,useBackwardEuler=USE_BACKWARD_EULER)  heat=TemperatureCartesian(dom,useBackwardEuler=False)
258  heat.setTolerance(TOL)  print "<%s> Temperature transport has been set up."%time.asctime()
259    heat.setTolerance(T_TOL)
260  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
261  print "initial Temperature range ",inf(T),sup(T)  heat.setInitialTemperature(clip(T,T_0))
262  heat.setInitialTemperature(clip(T,0))  heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
 heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)  
263  #  #
264  #   velocity constraints:  #   velocity constraints:
265  #  #
# Line 219  for d in range(DIM): Line 270  for d in range(DIM):
270         ll = H         ll = H
271      else:      else:
272         ll = L         ll = L
273      if CREATE_TOPOGRAPHY and d==DIM-1:      if TOPOGRAPHY_LOAD and d==DIM-1:
274         fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)         fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
275      else:      else:
276         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)
 1/0  
277  #  #
278  #   set up velovity problem  #   set up velovity problem
279  # ????????????????????  #
280  sp=StokesProblemCartesian(dom)  flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE)
281  # ,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)  flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1]))
282  # sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)  
283  # sp.initialize(mu=MU, tau_0=TAU0, n=N, eta_Y=ETAP0, tau_Y=TAU0, n_Y=NPL, q=fixed_v_mask)  flow.setElasticShearModulus(MUE)
284  sp.setTolerance(TOL*extratol)  flow.setTolerance(TOL2)
285  # sp.setToleranceReductionFactor(TOL)  flow.setFlowTolerance(FLOW_TOL)
286    flow.setFlowSubTolerance(FLOW_SUB_TOL)
287    flow.setEtaTolerance(FLOW_TOL)
288    flow.setExternals(fixed_v_mask=fixed_v_mask)
289    print "<%s> Flow solver has been set up."%time.asctime()
290  #  #
291  #  let the show begin:  #  let the show begin:
292  #  #
293  t1 = time.time()  t1 = time.time()
294    print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
295  while t<T_END:  while t<T_END:
296      v_last=v*1      #======= solve for velovity ====================================================================
297      print "============== solve for v ========================"      eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
298      FF=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))      print "viscosity range :", inf(eta_N), sup(eta_N)
299      print "viscosity range :", inf(FF)*ETA0, sup(FF)*ETA0      flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0],  [1,N])
300      # sp.initialize(eta_N=ETA0*FF, eta_0=ETA0*FF, F=T*(RA*unitVector(DIM-1,DIM)))      flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM)) #f= add the mountanins here
301      sp.initialize(eta=ETA0*FF, f=T*(RA*unitVector(DIM-1,DIM)),fixed_u_mask=fixed_v_mask)      flow.update(dt, iter_max=100, inner_iter_max=200, verbose=False)
302      v,p=sp.solve(v,p,max_iter=20, verbose=VERBOSE, show_details=False)      v=flow.getVelocity()
   
303      for d in range(DIM):      for d in range(DIM):
304           print "range %d-velocity"%d,inf(v[d]),sup(v[d])           print "range %d-velocity"%d,inf(v[d]),sup(v[d])
305      if t>=t_out or n>n_out:      print "<%s> flow solver completed."%time.asctime()
306        saveVTK("state.%d.vtu"%out_count,T=T,v=v)      if t>=t_vis or n>n_vis:
307        print "visualization file %d for time step %e generated."%(out_count,t)        saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())
308        out_count+=1        print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)
309        t_out+=DT_OUT        counter_vis+=1
310        n_out+=Dn_OUT        t_vis+=DT_VIS
311      # calculation of nusselt number:        n_vis+=DN_VIS
     se=ETA0*FF*length(symmetric(grad(v)))**2 # this propably wroing!  
     print "Xse:",inf(se),sup(se)  
     Nu=1.+integrate(se)/(RA*vol)  
     if dom.onMasterProcessor(): nusselt_file.write("%e %e\n"%(t,Nu))  
     heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)  
     print "Xnusselt number = ",Nu, "dt =",dt  
     if n>0:  
         a,a_alt = (v_last-v)/dt, a  
         dt_a,dt_a_alt = dt, dt_a  
     if n>1:  
        z=(a-a_alt)/((dt_a+dt_a_alt)/2)  
        f=Lsup(z)/Lsup(v)  
        print "estimated error ",f*dt**2  
        # dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))  
        dt_new=sqrt(0.05/f)  
        dt=min(dt_new,heat.getSafeTimeStepSize())  
     else:  
        dt=heat.getSafeTimeStepSize()  
     dt=max(DT_MIN,dt)  
     print n,". time step t=",t," step size ",dt  
     print "============== solve for T ========================"  
     T=heat.getSolution(dt, verbose=VERBOSE)  
     print "Temperature range ",inf(T),sup(T)  
312      n+=1      n+=1
313      t+=dt      t+=dt
314        print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)
315        #======= set Temperature problem ====================================================================
316        #
317        heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
318        dt=heat.getSafeTimeStepSize()
319        print "<%s> New time step size is %e"%(time.asctime(),dt)
320        print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)
321        #
322        #  solve temperature:
323        #
324        T=heat.getSolution(dt, verbose=VERBOSE)
325        print "Temperature range ",inf(T),sup(T)
326        print "<%s> temperature update completed."%time.asctime()
327        #======= anaysis ==================================================================================
328        #
329        #  .... Nusselt number
330        #
331        dTdz=grad(T)[DIM-1]
332        Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
333        print "nusselt number = ",Nu
334        nusselt_file.write("%e %e\n"%(t,Nu))
335        eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
336        Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
337        print "av. eta = %e"%eta_bar
338        print "effective Rayleigh number = %e"%Ra_eff
339      # =========================      # =========================
340      #      #
341      #    create restart files:      #    create restart files:
342      #      #
343      if create_restartfiles_every_step>0:      if DN_RESTART>0:
344         if (n-1)%create_restartfiles_every_step == 0:         if (n-1)%DN_RESTART == 0:
345           c=(n-1)/create_restartfiles_every_step           c=(n-1)/DN_RESTART
346           old_restart_dir="restart_%s_"%(c-1)           old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)
347           new_restart_dir="restart_%s_"%c           new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)
348    
          print "Write restart files to ",new_restart_dir  
349           if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)           if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
350       dom.MPIBarrier()       dom.MPIBarrier()
351           # sp.getStress().dump(os.path.join(new_restart_dir,"stress.nc"))           file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_vis, n_vis, n, counter_vis, dt))
352           v.dump(os.path.join(new_restart_dir,"v.nc"))           v.dump(os.path.join(new_restart_dir,"v.nc"))
353           p.dump(os.path.join(new_restart_dir,"p.nc"))           p.dump(os.path.join(new_restart_dir,"p.nc"))
354           T.dump(os.path.join(new_restart_dir,"T.nc"))           T.dump(os.path.join(new_restart_dir,"T.nc"))
          if n>1:  
              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))  
              a.dump(os.path.join(new_restart_dir,"a.nc"))  
          else:  
              file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_out, n_out, n, out_count, dt))  
355           removeRestartDirectory(old_restart_dir)           removeRestartDirectory(old_restart_dir)
356  elapsed = time.time() - t1           print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)
357  print "plot","\t",NE,"\t",elapsed  print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)
358    

Legend:
Removed from v.2439  
changed lines
  Added in v.2440

  ViewVC Help
Powered by ViewVC 1.1.26