user/examples/diffusion.py

# $Id$
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
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=LinearPDE(mydomain)
mypde.setSymmetryOn()
mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
# ... 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(Y=q+rhocp/h*T)
      T=mypde.getSolution()
      saveVTK("T.%d.xml"%i,temp=T)
1	jgs	102	# $Id$
2	gross	569	from esys.escript import *
3			from esys.escript.linearPDEs import LinearPDE
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	gross	569	mypde=LinearPDE(mydomain)
22			mypde.setSymmetryOn()
23			mypde.setValue(A=kappakronecker(mydomain),D=rhocp/h,d=eta,y=etaTref)
24	jgs	102	# ... set heat source: ....
25			x=mydomain.getX()
26	gross	327	q=qc*whereNegative(length(x-xc)-r)
27	jgs	102	# ... set initial temperature ....
28			T=Tref
29			# ... start iteration:
30	jgs	108	while t<tend:
31	jgs	102	i+=1
32			t+=h
33			print "time step :",t
34	gross	569	mypde.setValue(Y=q+rhocp/h*T)
35	jgs	102	T=mypde.getSolution()
36	gross	569	saveVTK("T.%d.xml"%i,temp=T)
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision