# Diff of /trunk/doc/examples/cookbook/example01c.py

revision 2949 by gross, Thu Feb 25 05:23:11 2010 UTC revision 2977 by ahallam, Tue Mar 9 00:33:28 2010 UTC
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
25
27    # example01c.py
28    # Model temperature diffusion between two granite blocks of unequal
29    # initial temperature. Solve for the spatial distribution of temperature.
30
31  # 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.
32  from esys.escript import * # This imports everything from the escript library  from esys.escript import * # This imports everything from the escript library
33  from esys.escript.unitsSI import *  from esys.escript.unitsSI import *
# Line 35  import pylab as pl #Plotting package. Line 40  import pylab as pl #Plotting package.
40  import numpy as np #Array package.  import numpy as np #Array package.
41  import os, sys #This package is necessary to handle saving our data.  import os, sys #This package is necessary to handle saving our data.
42
43  # .. MPI WORLD CHECK  ########################################################MPI WORLD CHECK
44  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
45      import sys      import sys
46      print "This example will not run in an MPI world."      print "This example will not run in an MPI world."
47      sys.exit(0)      sys.exit(0)
48
49  ##ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
50  #Domain related.  #Domain related.
51  mx = 500*m #meters - model length  mx = 500*m #meters - model length
52  my = 100*m #meters - model width  my = 100*m #meters - model width
# Line 57  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no h Line 62  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no h
62  T1=20 * Celsius # initial temperature at Block 1  T1=20 * Celsius # initial temperature at Block 1
63  T2=2273. * Celsius # initial temperature at Block 2  T2=2273. * Celsius # initial temperature at Block 2
64
65    ################################################ESTABLISHING PARAMETERS
66  t=0 * day  # our start time, usually zero  t=0 * day  # our start time, usually zero
67  tend=50 * yr # - time to end simulation  tend=50 * yr # - time to end simulation
68  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
# Line 69  save_path= os.path.join("data","example0 Line 75  save_path= os.path.join("data","example0
75  #ensure the dir exists  #ensure the dir exists
76  mkDir(save_path, os.path.join(save_path,"tempT"))  mkDir(save_path, os.path.join(save_path,"tempT"))
77
78  #... generate domain ...  ####################################################DOMAIN CONSTRUCTION
79  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
80
81    ###############################################ESCRIPT PDE CONSTRUCTION
82  #... open PDE and set coefficients ...  #... open PDE and set coefficients ...
83  mypde=LinearPDE(blocks)  mypde=LinearPDE(blocks)
84  mypde.setSymmetryOn()  mypde.setSymmetryOn()
# Line 88  E_list=[] Line 96  E_list=[]
96  plx = x.toListOfTuples()  plx = x.toListOfTuples()
97  plx = np.array(plx) #convert to tuple to numpy array  plx = np.array(plx) #convert to tuple to numpy array
98  plx = plx[:,0] #extract x locations  plx = plx[:,0] #extract x locations
99  # ... start iteration:  ########################################################START ITERATION
100  while t<tend:  while t<tend:
101        i+=1        i+=1
102        t+=h        t+=h
# Line 109  while t<tend: Line 117  while t<tend:
117        #save figure to file        #save figure to file
118        pl.savefig(os.path.join(save_path,"tempT", "blockspyplot%03d.png"%i))        pl.savefig(os.path.join(save_path,"tempT", "blockspyplot%03d.png"%i))
119        pl.clf() #clear figure        pl.clf() #clear figure
120
121    ###############################################################PLOTTING
122  # plot the total energy over time:  # plot the total energy over time:
123  pl.figure(2)  pl.figure(2)
124  pl.plot(t_list,E_list)  pl.plot(t_list,E_list)
# Line 117  pl.axis([0,max(t_list),0,max(E_list)*1.1 Line 127  pl.axis([0,max(t_list),0,max(E_list)*1.1
127  pl.savefig(os.path.join(save_path,"totE.png"))  pl.savefig(os.path.join(save_path,"totE.png"))
128  pl.clf()  pl.clf()
129
130    ###########################################################MAKE A MOVIE
131  # compile the *.png files to create a*.avi video that show T change  # compile the *.png files to create a*.avi video that show T change
132  # with time. This opperation uses linux mencoder. For other operating  # with time. This opperation uses linux mencoder. For other operating
133  # systems it may be possible to use your favourite video compiler to  # systems it may be possible to use your favourite video compiler to

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