test/python/convection.py


########################################################
#
# Copyright (c) 2003-2008 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2008 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.uq.edu.au/esscc/escript-finley"

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 dom.onMasterProcessor() and 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
   dom.MPIBarrier()


import sys
import time
t1 = time.time()

extratol=1

# read options:
parser = OptionParser(usage="%prog [-r [DIR]] [-e [NE]] [-s [solver]]")
parser.add_option("-s", "--solver", dest="solver", help="solver to be used for saddle point problem. The possible options are PCG, GMRES, NewtonGMRES, MINRES and TFQMR.", metavar="solver", default="PCG")
parser.add_option("-e", "--elements", dest="NE", help="number of elements in one direction.",metavar="NE", default=16)

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)

solver=options.solver
NE=int(options.NE)

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
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
   # this is a simple linear Stokes model:
   RA=1.e5 # Rayleigh number
   A=0 # 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

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

#
#   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.onMasterProcessor(): 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.onMasterProcessor(): 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)
heat.setTolerance(TOL*extratol)

fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
print "initial Temperature range ",inf(T),sup(T)
heat.setInitialTemperature(clip(T,0))
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=solver)
    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.getMechanicalPower()
    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)
    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 dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
     dom.MPIBarrier()
         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	ksteube	1811
2			########################################################
3			#
4			# Copyright (c) 2003-2008 by University of Queensland
5			# Earth Systems Science Computational Center (ESSCC)
6			# http://www.uq.edu.au/esscc
7			#
8			# Primary Business: Queensland, Australia
9			# Licensed under the Open Software License version 3.0
10			# http://www.opensource.org/licenses/osl-3.0.php
11			#
12			########################################################
13
14			__copyright__="""Copyright (c) 2003-2008 by University of Queensland
15			Earth Systems Science Computational Center (ESSCC)
16			http://www.uq.edu.au/esscc
17			Primary Business: Queensland, Australia"""
18			__license__="""Licensed under the Open Software License version 3.0
19			http://www.opensource.org/licenses/osl-3.0.php"""
20			__url__="http://www.uq.edu.au/esscc/escript-finley"
21
22	gross	1417	from esys.escript import *
23	gross	1639	from esys.escript.models import TemperatureCartesian, PlateMantelModel
24	gross	1474	from esys.finley import Rectangle, Brick, LoadMesh
25	gross	1472	from optparse import OptionParser
26	gross	1417	from math import pi, ceil
27	gross	1639
28
29	gross	1474	def removeRestartDirectory(dir_name):
30	ksteube	1877	if dom.onMasterProcessor() and os.path.isdir(dir_name):
31	gross	1474	for root, dirs, files in os.walk(dir_name, topdown=False):
32			for name in files: os.remove(os.path.join(root,name))
33			for name in dirs: os.remove(os.path.join(root,name))
34			os.rmdir(dir_name)
35			print "Restart files %s have been removed."%dir_name
36	ksteube	1877	dom.MPIBarrier()
37	gross	1472
38	artak	1483
39			import sys
40			import time
41			t1 = time.time()
42
43	artak	1674	extratol=1
44	artak	1483
45	artak	1674	# read options:
46			parser = OptionParser(usage="%prog [-r [DIR]] [-e [NE]] [-s [solver]]")
47			parser.add_option("-s", "--solver", dest="solver", help="solver to be used for saddle point problem. The possible options are PCG, GMRES, NewtonGMRES, MINRES and TFQMR.", metavar="solver", default="PCG")
48			parser.add_option("-e", "--elements", dest="NE", help="number of elements in one direction.",metavar="NE", default=16)
49	artak	1483
50	artak	1674	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")
51			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)
52			(options, args) = parser.parse_args()
53			restart=options.restart or (options.restart_dir !=None)
54	artak	1483
55	artak	1674	solver=options.solver
56			NE=int(options.NE)
57	gross	1639
58			DIM=2
59	gross	1417	H=1.
60	gross	1639	L=4*H
61			THETA=0.5 # time stepping THETA=0.5 cranck nicolson
62			TOL=1.e-4
63	gross	1417	PERTURBATION=0.1
64	gross	1639	DT=1.e-4
65			DT_MIN=1.e-5
66	gross	1552	T_END=0.1
67	gross	1463	DT_OUT=T_END/500
68	gross	1552	Dn_OUT=2
69	gross	1639	VERBOSE=True
70	gross	1474	create_restartfiles_every_step=10
71	gross	1639	if True:
72			# this is a simple linear Stokes model:
73			RA=1.e6 # Rayleigh number
74			A=22 # Arenious number
75			DI = 0. # dissipation number
76			MU=None
77			ETA0=1.
78			TAU0=None
79			N=None
80			NPL=None
81			ETAP0=ETA0
82			TAUY=None
83			useJAUMANNSTRESS=False
84	gross	1659	# this is a simple linear Stokes model:
85			RA=1.e5 # Rayleigh number
86			A=0 # Arenious number
87			DI = 0. # dissipation number
88			MU=None
89			ETA0=1.
90			TAU0=250.
91			N=None
92			NPL=None
93			ETAP0=ETA0
94			TAUY=TAU0
95			useJAUMANNSTRESS=False
96	gross	1639	else:
97			RA=1.e4 # Rayleigh number
98			A=22 # Arenious number
99			DI = 0. # dissipation number
100			MU=1.e4
101			ETA0=1.
102			TAU0=0.86610*2.5
103			N=3
104			NPL=14
105			TAUY=TAU0
106			ETAP0=ETA0
107			useJAUMANNSTRESS=True
108	gross	1474
109	gross	1463	print "total number of elements = ",NE*DIMint(L/H)**(DIM-1)
110	gross	1472
111	gross	1417	#
112			# set up domain:
113			#
114	gross	1472	if restart:
115			if options.restart_dir ==None:
116			restart_files=[]
117			for f in os.listdir("."):
118			if f.startswith("restart"): restart_files.append(f)
119			if len(restart_files)==0:
120	gross	1474	raise IOError,"no restart files"
121	gross	1472	restart_files.sort()
122			f=restart_files[-1]
123			else:
124			f=options.restart_dir
125	gross	1474	print "restart from directory ",f
126	artak	1483	try:
127			dom=LoadMesh("mesh.nc")
128			except:
129			pass
130	gross	1639	FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
131			t=float(FF[0])
132			t_out=float(FF[1])
133			n_out=int(FF[2])
134			n=int(FF[3])
135			out_count=int(FF[4])
136			dt=float(FF[5])
137			stress=load(os.path.join(f,"stress.nc"),dom)
138	gross	1474	v=load(os.path.join(f,"v.nc"),dom)
139			p=load(os.path.join(f,"p.nc"),dom)
140			T=load(os.path.join(f,"T.nc"),dom)
141			if n>1:
142	gross	1639	dt_a=float(FF[6])
143	gross	1474	a=load(os.path.join(f,"a.nc"),dom)
144			else:
145			dt_a=None
146			a=None
147	ksteube	1877	if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","a")
148	gross	1417	else:
149	gross	1472	if DIM==2:
150			dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
151			else:
152			dom=Brick(int(ceil(LNE/H)),int(ceil(LNE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True,optimize=True)
153	artak	1483	try:
154			dom.dump("mesh.nc")
155			except:
156			pass
157	gross	1472	x=dom.getX()
158			T=Scalar(1,ReducedSolution(dom))
159			for d in range(DIM):
160			if d == DIM-1:
161			T=sin(x[d]/Hpi)
162			else:
163	gross	1552	T=cos(x[d]/Lpi)
164	gross	1417
165	gross	1472	T=1.-x[DIM-1]+PERTURBATION*T
166			v=Vector(0,Solution(dom))
167	gross	1639	stress=Tensor(0,Function(dom))
168			x2=ReducedSolution(dom).getX()
169			p=-RA(x2[DIM-1]-0.5x2[DIM-1]**2)
170	ksteube	1877	if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","w")
171	gross	1472	t=0
172			t_out=0
173	gross	1552	n_out=0
174	gross	1472	n=0
175	gross	1552	out_count=0
176	gross	1639	dt=DT
177	gross	1474	a=None
178			dt_a=None
179	gross	1472
180	gross	1417	vol=integrate(1.,Function(dom))
181	gross	1639	p-=integrate(p)/vol
182	gross	1417	x=dom.getX()
183			#
184			# set up heat problem:
185			#
186	gross	1841	heat=TemperatureCartesian(dom,theta=THETA)
187	gross	1552	heat.setTolerance(TOL*extratol)
188	gross	1417
189			fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
190			print "initial Temperature range ",inf(T),sup(T)
191	gross	1841	heat.setInitialTemperature(clip(T,0))
192	gross	1417	heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
193			#
194	gross	1639	# velocity constraints:
195	gross	1417	#
196			x2=ReducedSolution(dom).getX()
197			fixed_v_mask=Vector(0,Solution(dom))
198			for d in range(DIM):
199			if d == DIM-1:
200			ll = H
201			else:
202			ll = L
203			fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
204			#
205	gross	1639	# set up velovity problem
206			#
207			sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
208			sp.initialize(mu=MU, tau_0=TAU0, n=N, eta_Y=ETAP0, tau_Y=TAU0, n_Y=NPL, q=fixed_v_mask)
209			sp.setTolerance(TOL*extratol)
210			sp.setToleranceReductionFactor(TOL)
211			#
212	gross	1417	# let the show begin:
213			#
214			while t<T_END:
215	gross	1472	v_last=v*1
216	gross	1417	print "============== solve for v ========================"
217	gross	1639	FF=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
218			print "viscosity range :", inf(FF)ETA0, sup(FF)ETA0
219			sp.initialize(eta_N=ETA0FF, eta_0=ETA0FF, F=T(RAunitVector(DIM-1,DIM)))
220	artak	1674	sp.update(dt,max_inner_iter=20, verbose=VERBOSE, show_details=False, tol=10., solver=solver)
221	gross	1639	v=sp.getVelocity()
222	gross	1417
223			for d in range(DIM):
224			print "range %d-velocity"%d,inf(v[d]),sup(v[d])
225	gross	1552	if t>=t_out or n>n_out:
226			saveVTK("state.%d.vtu"%out_count,T=T,v=v)
227			print "visualization file %d for time step %e generated."%(out_count,t)
228			out_count+=1
229	gross	1417	t_out+=DT_OUT
230	gross	1552	n_out+=Dn_OUT
231	gross	1639	# calculation of nusselt number:
232	gross	1659	se=sp.getMechanicalPower()
233			print "Xse:",inf(se),sup(se)
234	gross	1639	Nu=1.+integrate(se)/(RA*vol)
235	ksteube	1877	if dom.onMasterProcessor(): nusselt_file.write("%e %e\n"%(t,Nu))
236	gross	1639	heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)
237	gross	1659	print "Xnusselt number = ",Nu, "dt =",dt
238	gross	1474	if n>0:
239			a,a_alt = (v_last-v)/dt, a
240			dt_a,dt_a_alt = dt, dt_a
241			if n>1:
242			z=(a-a_alt)/((dt_a+dt_a_alt)/2)
243			f=Lsup(z)/Lsup(v)
244			print "estimated error ",fdt*2
245	gross	1552	# dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))
246			dt_new=sqrt(0.05/f)
247	gross	1474	dt=min(dt_new,heat.getSafeTimeStepSize())
248	gross	1472	else:
249	gross	1474	dt=heat.getSafeTimeStepSize()
250	gross	1639	dt=max(DT_MIN,dt)
251	gross	1417	print n,". time step t=",t," step size ",dt
252			print "============== solve for T ========================"
253	gross	1841	T=heat.getSolution(dt, verbose=VERBOSE)
254	gross	1417	print "Temperature range ",inf(T),sup(T)
255			n+=1
256			t+=dt
257	gross	1474	# =========================
258			#
259			# create restart files:
260			#
261			if create_restartfiles_every_step>0:
262			if (n-1)%create_restartfiles_every_step == 0:
263			c=(n-1)/create_restartfiles_every_step
264			old_restart_dir="restart_%s_"%(c-1)
265			new_restart_dir="restart_%s_"%c
266
267			print "Write restart files to ",new_restart_dir
268	ksteube	1877	if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
269			dom.MPIBarrier()
270	gross	1639	sp.getStress().dump(os.path.join(new_restart_dir,"stress.nc"))
271			sp.getVelocity().dump(os.path.join(new_restart_dir,"v.nc"))
272			sp.getPressure().dump(os.path.join(new_restart_dir,"p.nc"))
273	gross	1474	T.dump(os.path.join(new_restart_dir,"T.nc"))
274			if n>1:
275	gross	1552	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))
276	gross	1474	a.dump(os.path.join(new_restart_dir,"a.nc"))
277			else:
278	gross	1552	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))
279	gross	1474	removeRestartDirectory(old_restart_dir)
280	artak	1483	elapsed = time.time() - t1
281			print "plot","\t",NE,"\t",elapsed
282	gross	1474