Gluing Stand-Alone Applications

This section follows closely the approach of Langtangen: Python Scripting for Computational Science and Engineering, chapter 2.3.

Scripting

Suppose you have a stand-alone code for solving the ODE:

and you wish to have simulation results for varoius choices of the parameters. For instance, the provided

oscillator.f inp run

(gfortran -O3 -o oscillator oscillator.f)

reads the following parameters from the standard input, in the listed order:

, name of f(y) function, A,

The values of the parameters may be placed in a file inp: and the program can be run as

oscillator < inp

The output consists of data points , a suitable format to plot, e.g. with gnuplot.

Under the link you’ll find fortran, C and python versions of this stand-alon code.

Our first goal is to simplify the user’s interaction with the program; with the script simviz1.py it is possible to adjust the parameters through command-line options:

python simviz1.py -m 3.3 -b 0.9 -dt 0.05

The result should be png-files and an otional plot on the screen. The files for each such experiment should go in a subdirectory.

#!/usr/bin/env python
import sys, math

# default values of input parameters:
m = 1.0; b = 0.7; c = 5.0; func = 'y'; A = 5.0; w = 2*math.pi
y0 = 0.2; tstop = 30.0; dt = 0.05; case = 'tmp1'; screenplot = 1

# read variables from the command line, one by one:
while len(sys.argv) > 1:
    option = sys.argv[1];           del sys.argv[1]
    if   option == '-m':
        m = float(sys.argv[1]);     del sys.argv[1]
    elif option == '-b':
        b = float(sys.argv[1]);     del sys.argv[1]
    elif option == '-c':
        c = float(sys.argv[1]);     del sys.argv[1]
    elif option == '-func':
        func = sys.argv[1];         del sys.argv[1]
    elif option == '-A':
        A = float(sys.argv[1]);     del sys.argv[1]
    elif option == '-w':
        w = float(sys.argv[1]);     del sys.argv[1]
    elif option == '-y0':
        y0 = float(sys.argv[1]);    del sys.argv[1]
    elif option == '-tstop':
        tstop = float(sys.argv[1]); del sys.argv[1]
    elif option == '-dt':
        dt = float(sys.argv[1]);    del sys.argv[1]
    elif option == '-noscreenplot':
        screenplot = 0
    elif option == '-case':
        case = sys.argv[1];         del sys.argv[1]
    else:
        print sys.argv[0],': invalid option',option
        sys.exit(1)

# create a subdirectory:
d = case              # name of subdirectory
import os, shutil
if os.path.isdir(d):  # does d exist?
    shutil.rmtree(d)  # yes, remove old directory
os.mkdir(d)           # make new directory d
os.chdir(d)           # move to new directory d

# make input file to the program:
f = open('%s.i' % case, 'w')
# write a multi-line (triple-quoted) string with
# variable interpolation:
f.write("""
        %(m)g
        %(b)g
        %(c)g
        %(func)s
        %(A)g
        %(w)g
        %(y0)g
        %(tstop)g
        %(dt)g
        """ % vars())
f.close()
# run simulator:
cmd = '../oscillator < %s.i' % case  # command to run
failure = os.system(cmd)
if failure:
    print 'running the oscillator code failed'; sys.exit(1)

# make file with gnuplot commands:
f = open(case + '.gnuplot', 'w')
f.write("""
set title '%s: m=%g b=%g c=%g f(y)=%s A=%g w=%g y0=%g dt=%g';
""" % (case, m, b, c, func, A, w, y0, dt))
if screenplot:
    f.write("plot 'sim.dat' title 'y(t)' with lines;\n")
f.write("""
set size ratio 0.3 1.5, 1.0;  
# define the postscript output format:
set term postscript eps monochrome dashed 'Times-Roman' 28;
# output file containing the plot:
set output '%s.ps';
# basic plot command:
plot 'sim.dat' title 'y(t)' with lines;
# make a plot in PNG format:
set term png small;
set output '%s.png';
plot 'sim.dat' title 'y(t)' with lines;
""" % (case, case))
f.close()
# make plot:
cmd = 'gnuplot -geometry 800x200 -persist ' + case + '.gnuplot'
failure = os.system(cmd)
if failure:
    print 'running gnuplot failed'; sys.exit(1)

One can easily generate an automatic LaTeX-report or HTML-report containing the values of the parameters and the corresponding pictures.

Conducting Numerical Experiments

Suppose we want to keep most of the parameters fixed and vary one parameter between a minimum and a maximum value; the results should be collected in a HTML-report:

python loop4simviz1.py m_min m_max m_increment [ simviz1.py options ]

#!/usr/bin/env python
"""calls simviz1.py with different m values (in a loop)"""
import sys, os
usage = 'Usage: %s m_min m_max m_increment [ simviz1.py options ]' \
        % sys.argv[0]

try:
    m_min = float(sys.argv[1])
    m_max = float(sys.argv[2])
    dm    = float(sys.argv[3])
except IndexError:
    print usage;  sys.exit(1)

simviz1_options = ' '.join(sys.argv[4:])

html = open('tmp_mruns.html', 'w')
html.write('<HTML><BODY BGCOLOR="white">\n')
psfiles = []  # plot files in PostScript format

m = m_min
while m <= m_max:
    case = 'tmp_m_%g' % m
    cmd = 'python simviz1.py %s -m %g -case %s' % \
          (simviz1_options, m, case)
    print 'running', cmd
    os.system(cmd)
    html.write('<H1>m=%g</H1> <IMG SRC="%s">\n' \
               % (m,os.path.join(case,case+'.png')))
    psfiles.append(os.path.join(case,case+'.ps'))
    m += dm
html.write('</BODY></HTML>\n')

or let vary any parameter:

loop4simviz2.py parameter min max increment [ simviz2.py options ]

#!/usr/bin/env python
"""
As loop4simviz1.py, but here we call simviz2.py, make movies,
and also allow any simviz2.py option to be varied in a loop.
"""
import sys, os
usage = 'Usage: %s parameter min max increment '\
        '[ simviz2.py options ]' % sys.argv[0]
try:
    option_name = sys.argv[1]
    min = float(sys.argv[2])
    max = float(sys.argv[3])
    incr = float(sys.argv[4])
except:
    print usage;  sys.exit(1)

simviz2_options = ' '.join(sys.argv[5:])

html = open('tmp_%s_runs.html' % option_name, 'w')
html.write('<HTML><BODY BGCOLOR="white">\n')
psfiles = []    # plot files in PostScript format
pngfiles = []   # plot files in PNG format

value = min
while value <= max:
    case = 'tmp_%s_%g' % (option_name,value)
    cmd = 'python simviz2.py %s -%s %g -case %s' % \
          (simviz2_options, option_name, value, case)
    print 'running', cmd
    os.system(cmd)
    psfile = os.path.join(case,case+'.ps')
    pngfile = os.path.join(case,case+'.png')
    html.write('<H1>%s=%g</H1> <IMG SRC="%s">\n' \
               % (option_name, value, pngfile))
    psfiles.append(psfile)
    pngfiles.append(pngfile)
    value += incr
cmd = 'convert -delay 50 -loop 1000 %s tmp_%s.gif' \
      % (' '.join(pngfiles), option_name)
print 'converting PNG files to animated GIF:\n', cmd
os.system(cmd)
html.write('<H1>Movie</H1> <IMG SRC="tmp_%s.gif">\n' % option_name)
html.write('</BODY></HTML>\n')
html.close()