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

revision 2950 by gross, Thu Feb 25 07:33:16 2010 UTC revision 2977 by ahallam, Tue Mar 9 00:33:28 2010 UTC
# Line 22  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22  """  """
23  Author: Antony Hallam antony.hallam@uqconnect.edu.au  Author: Antony Hallam antony.hallam@uqconnect.edu.au
24  """  """
26    # example01b.py
27    # Model temperature diffusion between two granite blocks of unequal
28    # initial temperature. Solve for total energy in the system. Use
29    # matplotlib to visualise the answer.
30
31    #######################################################EXTERNAL MODULES
32  # 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.
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 *
# Line 35  import pylab as pl #Plotting package. Line 41  import pylab as pl #Plotting package.
41  import numpy as np #Array package.  import numpy as np #Array package.
42  import os, sys #This package is necessary to handle saving our data.  import os, sys #This package is necessary to handle saving our data.
43
44  ##ESTABLISHING VARIABLES  #################################################ESTABLISHING VARIABLES
45  #Domain related.  #Domain related.
46  mx = 500*m #meters - model length  mx = 500*m #meters - model length
47  my = 100*m #meters - model width  my = 100*m #meters - model width
# Line 51  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no h Line 57  qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no h
57  T1=20 * Celsius # initial temperature at Block 1  T1=20 * Celsius # initial temperature at Block 1
58  T2=2273. * Celsius # initial temperature at Block 2  T2=2273. * Celsius # initial temperature at Block 2
59
60    ################################################ESTABLISHING PARAMETERS
61  t=0 * day  # our start time, usually zero  t=0 * day  # our start time, usually zero
62  tend=50 * yr # - time to end simulation  tend=50 * yr # - time to end simulation
63  outputs = 200 # number of time steps required.  outputs = 200 # number of time steps required.
# Line 63  save_path= os.path.join("data","example0 Line 70  save_path= os.path.join("data","example0
70  #ensure the dir exists  #ensure the dir exists
71  mkDir(save_path, os.path.join(save_path,"tempT"))  mkDir(save_path, os.path.join(save_path,"tempT"))
72
73  #... generate domain ...  ####################################################DOMAIN CONSTRUCTION
74  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)  blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
75
76    ###############################################ESCRIPT PDE CONSTRUCTION
77  #... open PDE and set coefficients ...  #... open PDE and set coefficients ...
78  mypde=LinearPDE(blocks)  mypde=LinearPDE(blocks)
79  mypde.setSymmetryOn()  mypde.setSymmetryOn()
# Line 78  T= T1*whereNegative(x[0]-boundloc)+T2*(1 Line 87  T= T1*whereNegative(x[0]-boundloc)+T2*(1
87  # ... open a collector for the time marks and corresponding total energy  # ... open a collector for the time marks and corresponding total energy
88  t_list=[]  t_list=[]
89  E_list=[]  E_list=[]
90  # ... start iteration:  ########################################################START ITERATION
91  while t<tend:  while t<tend:
92        i+=1        i+=1
93        t+=h        t+=h
# Line 89  while t<tend: Line 98  while t<tend:
98        t_list.append(t)        t_list.append(t)
99        E_list.append(totE)        E_list.append(totE)
100
101    ###############################################################PLOTTING
102  # plot the total energy over time:  # plot the total energy over time:
103  if getMPIRankWorld() == 0:  if getMPIRankWorld() == 0:
104      pl.plot(t_list,E_list)      pl.plot(t_list,E_list)

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