user/examples/diffusion.py

# $Id$
from mytools import Helmholtz
from esys.escript import Lsup,length,whereNegative,saveVTK
from esys.finley import Rectangle
#... set some parameters ...
xc=[0.02,0.002]
r=0.001
qc=50.e6
Tref=0.
rhocp=2.6e6
eta=75.
kappa=240.
tend=5.
# ... time, time step size and counter ...
t=0
h=0.1
i=0
#... generate domain ...
mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
#... open PDE ...
mypde=Helmholtz(mydomain)
# ... set heat source: ....
x=mydomain.getX()
q=qc*whereNegative(length(x-xc)-r)
# ... set initial temperature ....
T=Tref
# ... start iteration:
while t<tend:
      i+=1
      t+=h
      print "time step :",t
      mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)
      T=mypde.getSolution()
      saveVTK("T%d.xml"%i,temp=T)
1	jgs	102	# $Id$
2	jgs	108	from mytools import Helmholtz
3	gross	327	from esys.escript import Lsup,length,whereNegative,saveVTK
4	jgs	108	from esys.finley import Rectangle
5	jgs	102	#... set some parameters ...
6	jgs	108	xc=[0.02,0.002]
7	jgs	102	r=0.001
8	jgs	108	qc=50.e6
9	jgs	102	Tref=0.
10			rhocp=2.6e6
11			eta=75.
12			kappa=240.
13	jgs	108	tend=5.
14			# ... time, time step size and counter ...
15			t=0
16	jgs	102	h=0.1
17			i=0
18			#... generate domain ...
19	jgs	108	mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
20	jgs	102	#... open PDE ...
21			mypde=Helmholtz(mydomain)
22			# ... set heat source: ....
23			x=mydomain.getX()
24	gross	327	q=qc*whereNegative(length(x-xc)-r)
25	jgs	102	# ... set initial temperature ....
26			T=Tref
27			# ... start iteration:
28	jgs	108	while t<tend:
29	jgs	102	i+=1
30			t+=h
31			print "time step :",t
32			mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/hT,eta=eta,g=etaTref)
33			T=mypde.getSolution()
34	gross	327	saveVTK("T%d.xml"%i,temp=T)
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision