/[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 2619 by jfenwick, Mon Jul 20 06:43:47 2009 UTC revision 2620 by gross, Thu Aug 20 06:24:00 2009 UTC
# Line 24  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, IncompressibleIsotropicFlowCartesian  from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian, Mountains, SubSteppingException
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 35  import time Line 35  import time
35  DIM=2                           # spatial dimension  DIM=2                           # spatial dimension
36  H=1.                            # height  H=1.                            # height
37  L=H                           # length  L=H                           # length
38  NE=50                           # number of elements in H-direction.  NE=40                           # number of elements in H-direction.
39  PERT=0.05               # initial temperature perturbation  PERT=0.05               # initial temperature perturbation
40  DT=1.e-7                        # initial time step size  DT=1.e-7                        # initial time step size
41  CREATE_TOPOGRAPHY=False          # create topgraphy  CREATE_TOPOGRAPHY=True         # create topgraphy
 TOPOGRAPHY_LOAD=False        # considers the loading of the topography in the flow  
42  DT_MIN=1.e-10                    # minumum time step size  DT_MIN=1.e-10                    # minumum time step size
43  T_END=0.1                       # end time  T_END=0.1                       # end time
44    
45  RHO_0=1.                        # surface density  RHO_0=100.                     # surface density  (lf ~ RHO_0**2)
46  G=1.                            # gravitational constant  G=1.                            # gravitational constant
47  ALPHA_0=4.e3                      # thermal expansion coefficient  ALPHA_0=0.1                    # thermal expansion coefficient (Ra ~ RHO_0**2 * ALPHA_0 = lf * ALPHA_0)
48  T_0=0.                          # surface temperature  T_0=0.                          # surface temperature
49  T_1=1.                          # bottom temperature  T_1=1.                          # bottom temperature
50  C_P=1                           # heat capacity  C_P=1                           # heat capacity
# Line 57  BETA=0                          # Drucke Line 56  BETA=0                          # Drucke
56  TAU_0=2*10**(2.5)               # transition stress  TAU_0=2*10**(2.5)               # transition stress
57  N=3                             # power for power law  N=3                             # power for power law
58    
59  E=23                            # activation energy  E=23*0                        # activation energy
60  V=18*0                        # activation volume  V=18*0                        # activation volume
61  T_OFFSET=1                      # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0)  T_OFFSET=1                      # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0)
62  R=1                             # gas constant  R=1                             # gas constant
63  ETA_N0=1.                       # viscosity at surface  ETA_N0=1.                       # viscosity at surface
64    
65    TOPO_SMOOTH=1e-5                # smoothing factor of extrapolation of surface velocity to interior
66    
67  T_TOL=1.e-4                     # tolerance temperature transport  T_TOL=1.e-4                     # tolerance temperature transport
68  TOL2=0.1                        # tolerance for flow update at a timestep. (large value will do only one correction step)  TOL2=0.1                        # tolerance for flow update at a timestep. (large value will do only one correction step)
69  FLOW_TOL=1.e-4                  # tolerance for inconcompressible flow solver  FLOW_TOL=1.e-6                  # tolerance for inconcompressible flow solver
70  FLOW_SUB_TOL=1.e-8             # sub-tolerance for inconcompressible flow solver  FLOW_SUB_TOL=1.e-12             # sub-tolerance for inconcompressible flow solver
71  NUSSELT_FN="nusselt.csv"  TOPO_TOL=1.e-6                 # tolerance for update of topography
72    DIAGNOSTICS_FN="diagnostics.csv"
73  VERBOSE=True  VERBOSE=True
74  DT_VIS=T_END/500                # time distane between two visulaization files  DT_VIS=T_END/500                # time distane between two visulaization files
75  DN_VIS=5                        # maximum counter increment between two visulaization files  DN_VIS=5                        # maximum counter increment between two visulaization files
76  VIS_DIR="results"               # name of the director for vis files  VIS_DIR="results"               # name of the director for vis files
77  DN_RESTART=1000000              # create a restart file every DN_RESTART step  DN_RESTART=1000000              # create a restart file every DN_RESTART step
78  PREFIX_RESTART="restart"        # name prefix for restart directories  PREFIX_RESTART="restart"        # name prefix for restart directories
79    TOPO_ITER_MAX=20                # maximum number of iteration steps to update topography
80  # =========================================================================================  # =========================================================================================
81    
82  def removeRestartDirectory(dir_name):  def removeRestartDirectory(dir_name):
# Line 103  if options.param !=None: Line 106  if options.param !=None:
106       exec(open(options.param,'r'))       exec(open(options.param,'r'))
107       print "Parameters are imported from file ",options.param       print "Parameters are imported from file ",options.param
108    
 CREATE_TOPOGRAPHY=CREATE_TOPOGRAPHY or TOPOGRAPHY_LOAD  
   
109  print "Input Parameters:"  print "Input Parameters:"
110  print "\tDimension                     DIM\t\t=\t",DIM  print "\tDimension                     DIM\t\t=\t",DIM
111  print "\tHeight                        H\t\t\t=\t",H  print "\tHeight                        H\t\t\t=\t",H
# Line 115  print "\tinitial time step size        D Line 116  print "\tinitial time step size        D
116  print "\tminimum time step size        DT_MIN\t\t=\t",DT_MIN  print "\tminimum time step size        DT_MIN\t\t=\t",DT_MIN
117  print "\tend time                      T_END\t\t=\t",T_END  print "\tend time                      T_END\t\t=\t",T_END
118  print "\tcreate topgraphy              CREATE_TOPOGRAPHY\t=\t",CREATE_TOPOGRAPHY  print "\tcreate topgraphy              CREATE_TOPOGRAPHY\t=\t",CREATE_TOPOGRAPHY
 print "\tloading by topography         TOPOGRAPHY_LOAD\t=\t",TOPOGRAPHY_LOAD  
119  print "\tsurface density               RHO_0\t\t=\t",RHO_0  print "\tsurface density               RHO_0\t\t=\t",RHO_0
120  print "\tgravitational constant        G\t\t\t=\t",G  print "\tgravitational constant        G\t\t\t=\t",G
121  print "\tthermal expansion coefficient ALPHA_0\t\t=\t",ALPHA_0  print "\tthermal expansion coefficient ALPHA_0\t\t=\t",ALPHA_0
# Line 134  print "\tactivation energy             E Line 134  print "\tactivation energy             E
134  print "\tactivation volume             V\t\t\t=\t",V  print "\tactivation volume             V\t\t\t=\t",V
135  print "\ttemperature offset            T_OFFSET\t\t=\t",T_OFFSET  print "\ttemperature offset            T_OFFSET\t\t=\t",T_OFFSET
136  print "\tgas constant                  R\t\t\t=\t",R  print "\tgas constant                  R\t\t\t=\t",R
137    print "\ttopography smoothing          TOPO_SMOOTH\t=\t",TOPO_SMOOTH
138    
139    print "\ttolerance for topography      TOPO_TOL\t\t=\t",TOPO_TOL
140  print "\ttransport tolerance           T_TOL\t\t=\t",T_TOL  print "\ttransport tolerance           T_TOL\t\t=\t",T_TOL
141  print "\ttolerance for flow updates    TOL2\t\t=\t",TOL2  print "\ttolerance for flow updates    TOL2\t\t=\t",TOL2
142  print "\tflow tolerance                FLOW_TOL\t\t=\t",FLOW_TOL  print "\tflow tolerance                FLOW_TOL\t\t=\t",FLOW_TOL
143  print "\tflow sub-tolerance            FLOW_SUB_TOL\t=\t",FLOW_SUB_TOL  print "\tflow sub-tolerance            FLOW_SUB_TOL\t=\t",FLOW_SUB_TOL
144  print "\tfile for Nusselt number       NUSSELT_FN\t=\t",NUSSELT_FN  print "\tfile for diagnostics          DIAGNOSTICS_FN\t=\t",DIAGNOSTICS_FN
145  print "\tmin. time incr. for vis file  DT_VIS\t\t=\t",DT_VIS  print "\tmin. time incr. for vis file  DT_VIS\t\t=\t",DT_VIS
146  print "\tmin. count incr. for vis file DN_VIS\t\t=\t",DN_VIS  print "\tmin. count incr. for vis file DN_VIS\t\t=\t",DN_VIS
147  print "\tvisualization dir             VIS_DIR\t\t=\t",VIS_DIR  print "\tvisualization dir             VIS_DIR\t\t=\t",VIS_DIR
148  print "\trestart counter incerement    DN_RESTART\t\t=\t",DN_RESTART  print "\trestart counter incerement    DN_RESTART\t=\t",DN_RESTART
149  print "\tprefix for restart dirs       PREFIX_RESTART\t\t=\t",PREFIX_RESTART  print "\tprefix for restart dirs       PREFIX_RESTART\t=\t",PREFIX_RESTART
150  print "\tverbosity                     VERBOSE\t\t=\t",VERBOSE  print "\tverbosity                     VERBOSE\t\t=\t",VERBOSE
151    
152  print "Control Parameters:"  print "Control Parameters:"
# Line 159  T_OFFSET_REF=T_OFFSET/(T_1-T_0) Line 161  T_OFFSET_REF=T_OFFSET/(T_1-T_0)
161  Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF  Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF
162  Di=ALPHA_0*G*H/C_P  Di=ALPHA_0*G*H/C_P
163  CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2)  CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2)
164    if CREATE_TOPOGRAPHY:
165       SURFACE_LOAD=RHO_0*G*H/P_REF
166    else:
167       SURFACE_LOAD=0.
168  if MUE == None:  if MUE == None:
169    De=None    De=None
170  else:  else:
# Line 175  if MUE == None: Line 181  if MUE == None:
181  else:  else:
182     print "\tDebora number surface         De\t\t=\t%e"%De     print "\tDebora number surface         De\t\t=\t%e"%De
183            
184  print "\tbottom viscosity              \t\t\t=\t%e"%ETA_BOT  print "\tBottom viscosity              \t\t\t=\t%e"%ETA_BOT
185  print "\tRayleigh number bottom        \t\t\t=\t%e"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT)  print "\tRayleigh number bottom        \t\t\t=\t%e"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT)
186  if MUE == None:  if MUE == None:
187     print "\tDebora number bottom          \t\t\t=\t",None     print "\tDebora number bottom          \t\t\t=\t",None
188  else:  else:
189     print "\tDebora number bottom          \t\t\t=\t%d"%(ETA_BOT/MUE/t_REF)     print "\tDebora number bottom          \t\t\t=\t%d"%(ETA_BOT/MUE/t_REF)
190  print "\tArrhenius                     Ar\t\t=\t%e"%Ar  print "\tArrhenius                     Ar\t\t=\t%e"%Ar
191    print "\tsurface load factor           SURFACE_LOAD\t=\t%e"%SURFACE_LOAD
192  print "\tscaled activation volume      V_REF\t\t=\t%e"%V_REF  print "\tscaled activation volume      V_REF\t\t=\t%e"%V_REF
193  #  some control variables (will be overwritten in  case of a restart:  #  some control variables (will be overwritten in  case of a restart:
194  #  #
# Line 191  dt=DT       # current time step size Line 198  dt=DT       # current time step size
198  t_vis=0      t_vis=0    
199  n_vis=0  n_vis=0
200  counter_vis=0  counter_vis=0
201  if getMPIRankWorld()==0:  mkDir(VIS_DIR)
     if not os.path.isdir(VIS_DIR): os.mkdir(VIS_DIR)  
202  #=========================  #=========================
203  #  #
204  #   set up domain or read restart file  #   set up domain or read restart file
# Line 223  if restart: Line 229  if restart:
229     v=load(os.path.join(f,"v.nc"),dom)     v=load(os.path.join(f,"v.nc"),dom)
230     p=load(os.path.join(f,"p.nc"),dom)     p=load(os.path.join(f,"p.nc"),dom)
231     T=load(os.path.join(f,"T.nc"),dom)     T=load(os.path.join(f,"T.nc"),dom)
232       if CREATE_TOPOGRAPHY:
233           topography=load(os.path.join(f,"topo.nc"),dom)
234        
235     nusselt_file=FileWriter(NUSSELT_FN,append=True)     diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)
236     print "<%s> Restart from file %s at time step %s (t=%s) completed."%(time.asctime(),f,t,n)     print "<%s> Restart from file %s at time step %s (t=%s) completed."%(time.asctime(),f,t,n)
237  else:  else:
238    if DIM==2:    if DIM==2:
239      dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=2, useFullElementOrder=True,optimize=True)      dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=2, useFullElementOrder=False,optimize=True)
240    else:    else:
241      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)      dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L/H,l1=L/H,l2=1,order=2, useFullElementOrder=False,optimize=True)
242    try:    try:
243       dom.dump("mesh.nc")       dom.dump("mesh.nc")
244    except:    except:
# Line 248  else: Line 256  else:
256    stress=Tensor(0,Function(dom))    stress=Tensor(0,Function(dom))
257    x2=ReducedSolution(dom).getX()    x2=ReducedSolution(dom).getX()
258    p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)    p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)
259    nusselt_file=FileWriter(NUSSELT_FN,append=False)  
260      if CREATE_TOPOGRAPHY:
261          topography=Scalar(0.,ReducedSolution(dom))
262      diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=False)
263      diagnostics_file.write("Ra = %e Lambda= %e\n"%(Ra, SURFACE_LOAD))
264    
265  p_last=p  p_last=p
266  x=dom.getX()  x=dom.getX()
# Line 267  heat.setValue(rhocp=1,k=1,given_T_mask=f Line 279  heat.setValue(rhocp=1,k=1,given_T_mask=f
279  #  #
280  x2=ReducedSolution(dom).getX()  x2=ReducedSolution(dom).getX()
281  fixed_v_mask=Vector(0,Solution(dom))  fixed_v_mask=Vector(0,Solution(dom))
282    faces=Scalar(0.,Solution(dom))
283  for d in range(DIM):  for d in range(DIM):
284      if d == DIM-1:      if d == DIM-1:
285         ll = H         ll = H
286      else:      else:
287         ll = L         ll = L
288      if TOPOGRAPHY_LOAD and d==DIM-1:      if CREATE_TOPOGRAPHY and d==DIM-1:
289         fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)         fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
290      else:      else:
291         fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)         s=whereZero(x[d])+whereZero(x[d]-ll)
292           faces+=s
293           fixed_v_mask+=s*unitVector(d,DIM)
294  #  #
295  #   set up velovity problem  #   set up velovity problem
296  #  #
# Line 289  flow.setEtaTolerance(FLOW_TOL) Line 304  flow.setEtaTolerance(FLOW_TOL)
304  flow.setExternals(fixed_v_mask=fixed_v_mask)  flow.setExternals(fixed_v_mask=fixed_v_mask)
305  print "<%s> Flow solver has been set up."%time.asctime()  print "<%s> Flow solver has been set up."%time.asctime()
306  #  #
307    # topography set-up
308    #
309    boundary=FunctionOnBoundary(dom).getX()[DIM-1]
310    top_boundary_mask=whereZero(boundary-sup(boundary))
311    surface_area=integrate(top_boundary_mask)
312    if CREATE_TOPOGRAPHY:
313        mts=Mountains(dom,eps=TOPO_SMOOTH,reduced=True)
314        mts.setTopography(topography)
315        print "<%s> topography has been set up."%time.asctime()
316    #
317  #  let the show begin:  #  let the show begin:
318  #  #
319  t1 = time.time()  t1 = time.time()
320  print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)  print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
321  while t<T_END:  while t<T_END:
322        topography_old=topography
323        v_old, p_old, stress_old=v, p, stress
324        T_old=T
325      #======= solve for velovity ====================================================================      #======= solve for velovity ====================================================================
326      eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))      eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
327      print "viscosity range :", inf(eta_N), sup(eta_N)      print "viscosity range :", inf(eta_N), sup(eta_N)
328      flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0],  [1,N])      flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0],  [1,N])
329      flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM)) #f= add the mountanins here      flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM))
330      flow.update(dt, iter_max=100, inner_iter_max=200, verbose=False)      # if dt<=0 or not CREATE_TOPOGRAPHY:
331        if not CREATE_TOPOGRAPHY:
332                flow.setExternals(f=-SURFACE_LOAD*topography*unitVector(DIM-1,DIM))
333                flow.update(dt, iter_max=100, inner_iter_max=200, verbose=False)
334        else:
335            topography_last=topography
336            Topo_norm, error_Topo=1,1
337            i=0
338            print "DDDDD : ====",dt
339            while error_Topo > TOPO_TOL * Topo_norm:
340                flow.setStatus(t, v_old, p_old, stress_old)
341                flow.setExternals(f=-SURFACE_LOAD*(topography-dt*v)*unitVector(DIM-1,DIM)*top_boundary_mask, restoration_factor=SURFACE_LOAD*dt*top_boundary_mask)
342                flow.update(dt, iter_max=100, inner_iter_max=200, verbose=False)
343                v=flow.getVelocity()
344                mts.setTopography(topography_old)
345                mts.setVelocity(v)
346                topography_last, topography=topography, mts.update(dt, allow_substeps=True)
347                error_Topo=sqrt(integrate(((topography-topography_last)*top_boundary_mask)**2))
348                Topo_norm=sqrt(integrate((topography*top_boundary_mask)**2))
349                print "DDDDD :", "input=",integrate(v*top_boundary_mask)[DIM-1],"output=", integrate(topography*top_boundary_mask)/Lsup(topography), error_Topo, Topo_norm
350                print "topography update step %s error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v)
351                i+=1
352                if i > TOPO_ITER_MAX:
353                   raise RuntimeError,"topography did not converge after %s steps."%TOPO_ITER_MAX
354      v=flow.getVelocity()      v=flow.getVelocity()
355      for d in range(DIM):      for d in range(DIM):
356           print "range %d-velocity"%d,inf(v[d]),sup(v[d])           print "range %d-velocity"%d,inf(v[d]),sup(v[d])
357      print "<%s> flow solver completed."%time.asctime()      print "<%s> flow solver completed."%time.asctime()
     if t>=t_vis or n>n_vis:  
       saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())  
       print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)  
       counter_vis+=1  
       t_vis+=DT_VIS  
       n_vis+=DN_VIS  
358      n+=1      n+=1
359      t+=dt      t+=dt
360        print "influx= ",integrate(inner(v,dom.getNormal())), sqrt(integrate(length(v)**2,FunctionOnBoundary(dom))), integrate(1., FunctionOnBoundary(dom))
361      print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)      print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)
362      #======= set Temperature problem ====================================================================      #======= setup Temperature problem ====================================================================
363      #      #
364      heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())      heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
365      dt=heat.getSafeTimeStepSize()      dt=heat.getSafeTimeStepSize()
366      print "<%s> New time step size is %e"%(time.asctime(),dt)      print "<%s> New time step size is %e"%(time.asctime(),dt)
367        #======= set-up topography ==================================================================================
368        if CREATE_TOPOGRAPHY:
369            dt=min(mts.getSafeTimeStepSize()*0.5,dt)
370            print "<%s> New time step size is %e"%(time.asctime(),dt)
371      print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)      print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)
372      #      #
373      #  solve temperature:      #  solve temperature:
# Line 331  while t<T_END: Line 381  while t<T_END:
381      #      #
382      dTdz=grad(T)[DIM-1]      dTdz=grad(T)[DIM-1]
383      Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)      Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
     print "nusselt number = ",Nu  
     nusselt_file.write("%e %e\n"%(t,Nu))  
384      eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())      eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
385      Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar      Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
386        print "nusselt number = %e"%Nu
387      print "av. eta = %e"%eta_bar      print "av. eta = %e"%eta_bar
388      print "effective Rayleigh number = %e"%Ra_eff      print "effective Rayleigh number = %e"%Ra_eff
389        if CREATE_TOPOGRAPHY:
390           topo_level=integrate(topography*top_boundary_mask)/surface_area
391           valleys_deep=inf(topography)
392           mountains_heigh=sup(topography)
393           print "topography level = ",topo_level
394           print "valleys deep = ", valleys_deep
395           print "mountains_heigh = ", mountains_heigh
396           diagnostics_file.write("%e %e %e %e %e %e %e\n"%(t,Nu, topo_level, valleys_deep, mountains_heigh, eta_bar, Ra_eff))
397        else:
398           diagnostics_file.write("%e %e %e %e\n"%(t,Nu, eta_bar, Ra_eff))
399        #
400        #  .... visualization
401        #
402        if t>=t_vis or n>n_vis:
403          if CREATE_TOPOGRAPHY:
404             saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff(), topography=topography_old, vex=mts.getVelocity())
405          else:
406             saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())
407          print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)
408          counter_vis+=1
409          t_vis+=DT_VIS
410          n_vis+=DN_VIS
411      # =========================      # =========================
412      #      #
413      #    create restart files:      #    create restart files:
# Line 346  while t<T_END: Line 417  while t<T_END:
417           c=(n-1)/DN_RESTART           c=(n-1)/DN_RESTART
418           old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)           old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)
419           new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)           new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)
420             mkDir(new_restart_dir)
          if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)  
421       dom.MPIBarrier()       dom.MPIBarrier()
422           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))           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))
423           v.dump(os.path.join(new_restart_dir,"v.nc"))           v.dump(os.path.join(new_restart_dir,"v.nc"))
424           p.dump(os.path.join(new_restart_dir,"p.nc"))           p.dump(os.path.join(new_restart_dir,"p.nc"))
425           T.dump(os.path.join(new_restart_dir,"T.nc"))           T.dump(os.path.join(new_restart_dir,"T.nc"))
426             if CREATE_TOPOGRAPHY: topography.dump(os.path.join(new_restart_dir,"topo.nc"))
427           removeRestartDirectory(old_restart_dir)           removeRestartDirectory(old_restart_dir)
428           print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)           print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)
429  print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)  print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)

Legend:
Removed from v.2619  
changed lines
  Added in v.2620

  ViewVC Help
Powered by ViewVC 1.1.26