/[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

trunk/doc/cookbook/onedheatdiff.py revision 2597 by ahallam, Thu Aug 6 03:30:09 2009 UTC trunk/doc/examples/cookbook/example01a.py revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
1    
2  ########################################################  from __future__ import print_function
3    ##############################################################################
4  #  #
5  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2009-2012 by University of Queensland
6  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
7  #  #
8  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
9  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
10  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
11  #  #
12  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13    # Development since 2012 by School of Earth Sciences
14    #
15    ##############################################################################
16    
17  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2009-2012 by University of Queensland
18  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
19  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
20  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
21  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 22  __url__="https://launchpad.net/escript-f Line 24  __url__="https://launchpad.net/escript-f
24  """  """
25  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
26  """  """
27    ############################################################FILE HEADER
28    # example01a.py
29    # Model temperature diffusion between two granite blocks of unequal
30    # initial temperature. Solve for total energy in the system.
31    
32    #######################################################EXTERNAL MODULES
33  # To solve the problem it is necessary to import the modules we require.  # To solve the problem it is necessary to import the modules we require.
34  from esys.escript import * # This imports everything from the escript library  from esys.escript import * # This imports everything from the escript library
35    from esys.escript.unitsSI import *
36  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE  from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
37  from esys.finley import Rectangle # This imports the rectangle domain function from finley  from esys.finley import Rectangle # This imports the rectangle domain function
 import os #This package is necessary to handle saving our data.  
   
   
38    
39  ##ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
40    #Domain related.
41    mx = 500*m #meters - model length
42    my = 100*m #meters - model width
43    ndx = 100 # mesh steps in x direction
44    ndy = 1 # mesh steps in y direction - one dimension means one element
45    boundloc = mx/2 # location of boundary between the two blocks
46  #PDE related  #PDE related
47  q=50.e6 #our heat source temperature  rho = 2750. *kg/m**3 #kg/m{3} density of iron
48  Tref=0. #the starting temperature of our iron bar  cp = 790.*J/(kg*K) # J/Kg.K thermal capacity
49  rho=2.6e6  rhocp = rho*cp
50  eta=0#75.  kappa = 2.2*W/m/K # watts/m.Kthermal conductivity
51  kappa=240.  qH=0 * J/(sec*m**3) # J/(sec.m{3}) no heat source
52  #Script/Iteration Related  T1=20 * Celsius # initial temperature at Block 1
53  t=0 #our start time, usually zero  T2=2273. * Celsius # base temperature at Block 2
54  tend=5.#the time we want to end the simulation  
55  h=0.05 #size of time step  
56    ################################################ESTABLISHING PARAMETERS
57  print "Expected Number of Output Files is: ", (tend-t)/h  t=0 * day  # our start time, usually zero
58    tend=50 * yr # - time to end simulation
59  i=0 #loop counter  outputs = 200 # number of time steps required.
60  save_path = "data/onedheatdiff" #the folder to put our outputs in, leave blank "" for script path - note this folder path must exist to work  h=(tend-t)/outputs #size of time step
61    #user warning statement
62  #... generate domain ...  print("Expected Number of time outputs is: ", (tend-t)/h)
63  rod = Rectangle(l0=0.05,l1=.01,n0=500, n1=1)  i=0 #loop counter
64  # extract finite points  #the folder to put our outputs in, leave blank "" for script path
65  x=rod.getX()  save_path= os.path.join("data","example01")
66  #... open PDE ...  #ensure the dir exists
67  mypde=LinearPDE(rod)  mkDir(save_path, os.path.join(save_path,"tempT"))
68    
69    ####################################################DOMAIN CONSTRUCTION
70    blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
71    
72    ###############################################ESCRIPT PDE CONSTRUCTION
73    #... open PDE and set coefficients ...
74    mypde=LinearPDE(blocks)
75  mypde.setSymmetryOn()  mypde.setSymmetryOn()
76  mypde.setValue(A=kappa*kronecker(rod),D=rho/h,d=eta,y=eta*Tref)  A=zeros((2,2))
77  # ... set heat source: ....  A[0,0]=kappa
78    mypde.setValue(A=A,D=rhocp/h)
 qH=q*whereZero(x[0])  
79  # ... set initial temperature ....  # ... set initial temperature ....
80  T=Tref  x=Solution(blocks).getX()
81    T= T1*whereNegative(x[0]-boundloc)+T2*(1-whereNegative(x[0]-boundloc))
82    
83  # ... start iteration:  ########################################################START ITERATION
84  while t<=tend:  while t<tend:
85        i+=1        i+=1
86        t+=h        t+=h
87        mypde.setValue(Y=qH+rho/h*T)        mypde.setValue(Y=qH+rhocp/h*T)
88        T=mypde.getSolution()        T=mypde.getSolution()
89        print T        totE=integrate(rhocp*T)
90        saveVTK(os.path.join(save_path,"data%03d.vtu") %i,sol=T)        print("time step %s at t=%e days completed. total energy = %e."%(i,t/day,totE))
         
 #~ command = ('mencoder',  
            #~ 'mf://*.png',  
            #~ '-mf',  
            #~ 'type=png:w=800:h=600:fps=25',  
            #~ '-ovc',  
            #~ 'lavc',  
            #~ '-lavcopts',  
            #~ 'vcodec=mpeg4',  
            #~ '-oac',  
            #~ 'copy',  
            #~ '-o',  
            #~ 'output.avi')  
   

Legend:
Removed from v.2597  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26