/[escript]/trunk/esys2/doc/user/examples/diffusion.py
ViewVC logotype

Diff of /trunk/esys2/doc/user/examples/diffusion.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 102 by jgs, Wed Dec 15 07:08:39 2004 UTC revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2  from mytools import *  from mytools import Helmholtz
3  from esys.escript import *  from esys.escript import Lsup
4  import esys.finley  from esys.finley import Rectangle
5  #... set some parameters ...  #... set some parameters ...
6  x_c=[0.02,0.002]  xc=[0.02,0.002]
7  r=0.001  r=0.001
8  q0=50.e6  qc=50.e6
9  Tref=0.  Tref=0.
10  rhocp=2.6e6  rhocp=2.6e6
11  eta=75.  eta=75.
12  kappa=240.  kappa=240.
13  t_end=5.  tend=5.
14  # ...time step size and counter ...  # ... time, time step size and counter ...
15    t=0
16  h=0.1  h=0.1
17  i=0  i=0
 t=0  
18  #... generate domain ...  #... generate domain ...
19  mydomain = esys.finley.Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)  mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
20  #... open PDE ...  #... open PDE ...
21  mypde=Helmholtz(mydomain)  mypde=Helmholtz(mydomain)
22  # ... set heat source: ....  # ... set heat source: ....
23  x=mydomain.getX()  x=mydomain.getX()
24  q=q0*(length(x-x_c)-r).whereNegative()  q=qc*(length(x-xc)-r).whereNegative()
25  # ... set initial temperature ....  # ... set initial temperature ....
26  T=Tref  T=Tref
27  # ... start iteration:  # ... start iteration:
28  while t<t_end:  while t<tend:
29        i+=1        i+=1
30        t+=h        t+=h
31        print "time step :",t        print "time step :",t
32        mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)        mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)
33        T=mypde.getSolution()        T=mypde.getSolution()
34        T.saveDX("T%d.dx"%i)        T.saveDX("T%d.dx"%i)
Legend:
Removed from v.102  
changed lines
  Added in v.108
  ViewVC Help
Powered by ViewVC 1.1.26