test/python/convection.py

from esys.escript import *
from esys.escript.models import TemperatureCartesian, PlateMantelModel
from esys.finley import Rectangle, Brick, LoadMesh
from optparse import OptionParser
from math import pi, ceil


def removeRestartDirectory(dir_name):
   if os.path.isdir(dir_name):
       for root, dirs, files in os.walk(dir_name, topdown=False):
           for name in files: os.remove(os.path.join(root,name))
           for name in dirs: os.remove(os.path.join(root,name))
       os.rmdir(dir_name)
       print "Restart files %s have been removed."%dir_name


import sys
import time
t1 = time.time()

if (len(sys.argv)>=3):
 NE=int(sys.argv[2])
else:
 NE=20
NE=64
NE=20

if (len(sys.argv)>=2):
 solver=sys.argv[1]
else:
 solver='PCG'

if solver!='PCG':
 extratol=0.001
else:
 extratol=1
extratol=0.1


DIM=2
H=1.
L=4*H
THETA=0.5 # time stepping THETA=0.5 cranck nicolson
TOL=1.e-4
PERTURBATION=0.1
DT=1.e-4
DT_MIN=1.e-5
T_END=0.1
DT_OUT=T_END/500
Dn_OUT=2
VERBOSE=True
SUPG=False
create_restartfiles_every_step=10
if True:
   # this is a simple linear Stokes model:
   RA=1.e6 # Rayleigh number
   A=22 # Arenious number 
   DI = 0.  # dissipation number 
   MU=None
   ETA0=1.
   TAU0=None
   N=None
   NPL=None
   ETAP0=ETA0
   TAUY=None
   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

print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)

# read options:
parser = OptionParser(usage="%prog [-r [DIR]]")
parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory. It will be deleted after a new directory has been created.", default=False, action="store_true")
parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR. The directory will not be deleted but new restart directories are created.",metavar="DIR", default=None)
(options, args) = parser.parse_args()
restart=options.restart or (options.restart_dir !=None)
#
#   set up domain:
#
if restart:
   if options.restart_dir ==None:
      restart_files=[]
      for f in os.listdir("."):
          if f.startswith("restart"): restart_files.append(f)
      if len(restart_files)==0:
          raise IOError,"no restart files"
      restart_files.sort()
      f=restart_files[-1]
   else:
       f=options.restart_dir
   print "restart from directory ",f
   try:
      dom=LoadMesh("mesh.nc")
   except:
      pass
   FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
   t=float(FF[0])
   t_out=float(FF[1])
   n_out=int(FF[2])
   n=int(FF[3])
   out_count=int(FF[4])
   dt=float(FF[5])
   stress=load(os.path.join(f,"stress.nc"),dom)
   v=load(os.path.join(f,"v.nc"),dom)
   p=load(os.path.join(f,"p.nc"),dom)
   T=load(os.path.join(f,"T.nc"),dom)
   if n>1:
      dt_a=float(FF[6])
      a=load(os.path.join(f,"a.nc"),dom)
   else:
      dt_a=None
      a=None
   if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a")
else:
  if DIM==2:
    dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
  else:
    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)
  try:
     dom.dump("mesh.nc")
  except:
     pass
  x=dom.getX()
  T=Scalar(1,ReducedSolution(dom))
  for d in range(DIM):
      if d == DIM-1: 
         T*=sin(x[d]/H*pi)
      else:
         T*=cos(x[d]/L*pi)

  T=1.-x[DIM-1]+PERTURBATION*T
  v=Vector(0,Solution(dom))
  stress=Tensor(0,Function(dom))
  x2=ReducedSolution(dom).getX()
  p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
  if dom.getMPIRank() ==0: nusselt_file=open("nusselt.csv","w")
  t=0
  t_out=0
  n_out=0
  n=0
  out_count=0
  dt=DT
  a=None
  dt_a=None

vol=integrate(1.,Function(dom))
p-=integrate(p)/vol
x=dom.getX()
#
#   set up heat problem:
#
heat=TemperatureCartesian(dom,theta=THETA,useSUPG=SUPG)
heat.setTolerance(TOL*extratol)

fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
heat.setInitialTemperature(T)
print "initial Temperature range ",inf(T),sup(T)
heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
#
#   velocity constraints:
#
x2=ReducedSolution(dom).getX()
fixed_v_mask=Vector(0,Solution(dom))
for d in range(DIM):
    if d == DIM-1: 
       ll = H
    else:
       ll = L
    fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
#
#   set up velovity problem
#
sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
sp.initialize(mu=MU, tau_0=TAU0, n=N, eta_Y=ETAP0, tau_Y=TAU0, n_Y=NPL, q=fixed_v_mask)
sp.setTolerance(TOL*extratol)
sp.setToleranceReductionFactor(TOL)
#
#  let the show begin:
#
while t<T_END:
    v_last=v*1
    print "============== solve for v ========================"
    FF=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
    print "viscosity range :", inf(FF)*ETA0, sup(FF)*ETA0
    sp.initialize(eta_N=ETA0*FF, eta_0=ETA0*FF, F=T*(RA*unitVector(DIM-1,DIM)))
    sp.update(dt,max_inner_iter=20, verbose=VERBOSE, show_details=False, tol=10., solver='PCG')
    v=sp.getVelocity()

    for d in range(DIM):
         print "range %d-velocity"%d,inf(v[d]),sup(v[d])
    if t>=t_out or n>n_out:
      saveVTK("state.%d.vtu"%out_count,T=T,v=v)
      print "visualization file %d for time step %e generated."%(out_count,t)
      out_count+=1
      t_out+=DT_OUT
      n_out+=Dn_OUT
    # calculation of nusselt number:
    se=sp.getStrainEnergy()
    Nu=1.+integrate(se)/(RA*vol)
    if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))
    heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)
    print "nusselt number = ",Nu
    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.solve(dt, verbose=VERBOSE)   
    print "Temperature range ",inf(T),sup(T)
    n+=1
    t+=dt
    # =========================
    #
    #    create restart files:
    #
    if create_restartfiles_every_step>0:
       if (n-1)%create_restartfiles_every_step == 0:
         c=(n-1)/create_restartfiles_every_step
         old_restart_dir="restart_%s_"%(c-1)
         new_restart_dir="restart_%s_"%c

         print "Write restart files to ",new_restart_dir
         if not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
         sp.getStress().dump(os.path.join(new_restart_dir,"stress.nc"))
         sp.getVelocity().dump(os.path.join(new_restart_dir,"v.nc"))
         sp.getPressure().dump(os.path.join(new_restart_dir,"p.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))
         removeRestartDirectory(old_restart_dir)
elapsed = time.time() - t1
print "plot","\t",NE,"\t",elapsed

1	from esys.escript import *
2	from esys.escript.models import TemperatureCartesian, PlateMantelModel
3	from esys.finley import Rectangle, Brick, LoadMesh
4	from optparse import OptionParser
5	from math import pi, ceil
6
7
8	def removeRestartDirectory(dir_name):
9	if os.path.isdir(dir_name):
10	for root, dirs, files in os.walk(dir_name, topdown=False):
11	for name in files: os.remove(os.path.join(root,name))
12	for name in dirs: os.remove(os.path.join(root,name))
13	os.rmdir(dir_name)
14	print "Restart files %s have been removed."%dir_name
15
16
17	import sys
18	import time
19	t1 = time.time()
20
21	if (len(sys.argv)>=3):
22	NE=int(sys.argv[2])
23	else:
24	NE=20
25	NE=64
26	NE=20
27
28	if (len(sys.argv)>=2):
29	solver=sys.argv[1]
30	else:
31	solver='PCG'
32
33	if solver!='PCG':
34	extratol=0.001
35	else:
36	extratol=1
37	extratol=0.1
38
39
40	DIM=2
41	H=1.
42	L=4*H
43	THETA=0.5 # time stepping THETA=0.5 cranck nicolson
44	TOL=1.e-4
45	PERTURBATION=0.1
46	DT=1.e-4
47	DT_MIN=1.e-5
48	T_END=0.1
49	DT_OUT=T_END/500
50	Dn_OUT=2
51	VERBOSE=True
52	SUPG=False
53	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.86610*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*DIMint(L/H)**(DIM-1)
81
82	# read options:
83	parser = OptionParser(usage="%prog [-r [DIR]]")
84	parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory. It will be deleted after a new directory has been created.", default=False, action="store_true")
85	parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR. The directory will not be deleted but new restart directories are created.",metavar="DIR", default=None)
86	(options, args) = parser.parse_args()
87	restart=options.restart or (options.restart_dir !=None)
88	#
89	# set up domain:
90	#
91	if restart:
92	if options.restart_dir ==None:
93	restart_files=[]
94	for f in os.listdir("."):
95	if f.startswith("restart"): restart_files.append(f)
96	if len(restart_files)==0:
97	raise IOError,"no restart files"
98	restart_files.sort()
99	f=restart_files[-1]
100	else:
101	f=options.restart_dir
102	print "restart from directory ",f
103	try:
104	dom=LoadMesh("mesh.nc")
105	except:
106	pass
107	FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
108	t=float(FF[0])
109	t_out=float(FF[1])
110	n_out=int(FF[2])
111	n=int(FF[3])
112	out_count=int(FF[4])
113	dt=float(FF[5])
114	stress=load(os.path.join(f,"stress.nc"),dom)
115	v=load(os.path.join(f,"v.nc"),dom)
116	p=load(os.path.join(f,"p.nc"),dom)
117	T=load(os.path.join(f,"T.nc"),dom)
118	if n>1:
119	dt_a=float(FF[6])
120	a=load(os.path.join(f,"a.nc"),dom)
121	else:
122	dt_a=None
123	a=None
124	if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a")
125	else:
126	if DIM==2:
127	dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
128	else:
129	dom=Brick(int(ceil(LNE/H)),int(ceil(LNE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True,optimize=True)
130	try:
131	dom.dump("mesh.nc")
132	except:
133	pass
134	x=dom.getX()
135	T=Scalar(1,ReducedSolution(dom))
136	for d in range(DIM):
137	if d == DIM-1:
138	T=sin(x[d]/Hpi)
139	else:
140	T=cos(x[d]/Lpi)
141
142	T=1.-x[DIM-1]+PERTURBATION*T
143	v=Vector(0,Solution(dom))
144	stress=Tensor(0,Function(dom))
145	x2=ReducedSolution(dom).getX()
146	p=-RA(x2[DIM-1]-0.5x2[DIM-1]**2)
147	if dom.getMPIRank() ==0: nusselt_file=open("nusselt.csv","w")
148	t=0
149	t_out=0
150	n_out=0
151	n=0
152	out_count=0
153	dt=DT
154	a=None
155	dt_a=None
156
157	vol=integrate(1.,Function(dom))
158	p-=integrate(p)/vol
159	x=dom.getX()
160	#
161	# set up heat problem:
162	#
163	heat=TemperatureCartesian(dom,theta=THETA,useSUPG=SUPG)
164	heat.setTolerance(TOL*extratol)
165
166	fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
167	heat.setInitialTemperature(T)
168	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)
170	#
171	# velocity constraints:
172	#
173	x2=ReducedSolution(dom).getX()
174	fixed_v_mask=Vector(0,Solution(dom))
175	for d in range(DIM):
176	if d == DIM-1:
177	ll = H
178	else:
179	ll = L
180	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:
190	#
191	while t<T_END:
192	v_last=v*1
193	print "============== solve for v ========================"
194	FF=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
195	print "viscosity range :", inf(FF)ETA0, sup(FF)ETA0
196	sp.initialize(eta_N=ETA0FF, eta_0=ETA0FF, F=T(RAunitVector(DIM-1,DIM)))
197	sp.update(dt,max_inner_iter=20, verbose=VERBOSE, show_details=False, tol=10., solver='PCG')
198	v=sp.getVelocity()
199
200	for d in range(DIM):
201	print "range %d-velocity"%d,inf(v[d]),sup(v[d])
202	if t>=t_out or n>n_out:
203	saveVTK("state.%d.vtu"%out_count,T=T,v=v)
204	print "visualization file %d for time step %e generated."%(out_count,t)
205	out_count+=1
206	t_out+=DT_OUT
207	n_out+=Dn_OUT
208	# 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))
212	heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)
213	print "nusselt number = ",Nu
214	if n>0:
215	a,a_alt = (v_last-v)/dt, a
216	dt_a,dt_a_alt = dt, dt_a
217	if n>1:
218	z=(a-a_alt)/((dt_a+dt_a_alt)/2)
219	f=Lsup(z)/Lsup(v)
220	print "estimated error ",fdt*2
221	# dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))
222	dt_new=sqrt(0.05/f)
223	dt=min(dt_new,heat.getSafeTimeStepSize())
224	else:
225	dt=heat.getSafeTimeStepSize()
226	dt=max(DT_MIN,dt)
227	print n,". time step t=",t," step size ",dt
228	print "============== solve for T ========================"
229	T=heat.solve(dt, verbose=VERBOSE)
230	print "Temperature range ",inf(T),sup(T)
231	n+=1
232	t+=dt
233	# =========================
234	#
235	# create restart files:
236	#
237	if create_restartfiles_every_step>0:
238	if (n-1)%create_restartfiles_every_step == 0:
239	c=(n-1)/create_restartfiles_every_step
240	old_restart_dir="restart_%s_"%(c-1)
241	new_restart_dir="restart_%s_"%c
242
243	print "Write restart files to ",new_restart_dir
244	if not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
245	sp.getStress().dump(os.path.join(new_restart_dir,"stress.nc"))
246	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"))
249	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))
251	a.dump(os.path.join(new_restart_dir,"a.nc"))
252	else:
253	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))
254	removeRestartDirectory(old_restart_dir)
255	elapsed = time.time() - t1
256	print "plot","\t",NE,"\t",elapsed
257