/[escript]/trunk/doc/examples/cookbook/example01a.py
ViewVC logotype

Diff of /trunk/doc/examples/cookbook/example01a.py

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

revision 2977 by ahallam, Tue Mar 9 00:33:28 2010 UTC revision 3892 by jfenwick, Tue Apr 10 08:57:23 2012 UTC
# Line 1  Line 1 
1    
2    from __future__ import print_function
3  ########################################################  ########################################################
4  #  #
5  # Copyright (c) 2009 by University of Queensland  # Copyright (c) 2009 by University of Queensland
# Line 32  Author: Antony Hallam antony.hallam@uqco Line 33  Author: Antony Hallam antony.hallam@uqco
33  from esys.escript import * # This imports everything from the escript library  from esys.escript import * # This imports everything from the escript library
34  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
35  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
36  from esys.finley import Rectangle # This imports the rectangle domain function from finley  from esys.finley import Rectangle # This imports the rectangle domain function
37    
38  #################################################ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
39  #Domain related.  #Domain related.
# Line 42  ndx = 100 # mesh steps in x direction Line 43  ndx = 100 # mesh steps in x direction
43  ndy = 1 # mesh steps in y direction - one dimension means one element  ndy = 1 # mesh steps in y direction - one dimension means one element
44  boundloc = mx/2 # location of boundary between the two blocks  boundloc = mx/2 # location of boundary between the two blocks
45  #PDE related  #PDE related
46  rho = 7874. *kg/m**3 #kg/m^{3} density of iron  rho = 2750. *kg/m**3 #kg/m{3} density of iron
47  cp = 449.*J/(kg*K) # J/Kg.K thermal capacity  cp = 790.*J/(kg*K) # J/Kg.K thermal capacity
48  rhocp = rho*cp  rhocp = rho*cp
49  kappa = 80.*W/m/K   # watts/m.Kthermal conductivity  kappa = 2.2*W/m/K # watts/m.Kthermal conductivity
50  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source  qH=0 * J/(sec*m**3) # J/(sec.m{3}) no heat source
51  T1=20 * Celsius # initial temperature at Block 1  T1=20 * Celsius # initial temperature at Block 1
52  T2=2273. * Celsius # initial temperature at Block 2  T2=2273. * Celsius # base temperature at Block 2
53    
54    
55  ################################################ESTABLISHING PARAMETERS  ################################################ESTABLISHING PARAMETERS
56  t=0 * day  # our start time, usually zero  t=0 * day  # our start time, usually zero
# Line 56  tend=50 * yr # - time to end simulation Line 58  tend=50 * yr # - time to end simulation
58  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
59  h=(tend-t)/outputs #size of time step  h=(tend-t)/outputs #size of time step
60  #user warning statement  #user warning statement
61  print "Expected Number of time outputs is: ", (tend-t)/h  print("Expected Number of time outputs is: ", (tend-t)/h)
62  i=0 #loop counter  i=0 #loop counter
63  #the folder to put our outputs in, leave blank "" for script path  #the folder to put our outputs in, leave blank "" for script path
64  save_path= os.path.join("data","example01")  save_path= os.path.join("data","example01")
# Line 67  mkDir(save_path, os.path.join(save_path, Line 69  mkDir(save_path, os.path.join(save_path,
69  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
70    
71  ###############################################ESCRIPT PDE CONSTRUCTION  ###############################################ESCRIPT PDE CONSTRUCTION
72    #... open PDE and set coefficients ...
73    mypde=LinearPDE(blocks)
74  mypde.setSymmetryOn()  mypde.setSymmetryOn()
75  A=zeros((2,2))  A=zeros((2,2))
76  A[0,0]=kappa  A[0,0]=kappa
# Line 82  while t<tend: Line 86  while t<tend:
86        mypde.setValue(Y=qH+rhocp/h*T)        mypde.setValue(Y=qH+rhocp/h*T)
87        T=mypde.getSolution()        T=mypde.getSolution()
88        totE=integrate(rhocp*T)        totE=integrate(rhocp*T)
89        print "time step %s at t=%e days completed. total energy = %e."%(i,t/day,totE)        print("time step %s at t=%e days completed. total energy = %e."%(i,t/day,totE))

Legend:
Removed from v.2977  
changed lines
  Added in v.3892

  ViewVC Help
Powered by ViewVC 1.1.26