PyMCT and PyCPL: Refactoring CCSM Using Python Michael Tobis (1) , - - PowerPoint PPT Presentation

pymct and pycpl refactoring ccsm using python
SMART_READER_LITE
LIVE PREVIEW

PyMCT and PyCPL: Refactoring CCSM Using Python Michael Tobis (1) , - - PowerPoint PPT Presentation

PyMCT and PyCPL: Refactoring CCSM Using Python Michael Tobis (1) , Michael Steder (1) , Robert L. Jacob (2,3) , Raymond T. Pierrehumbert (1) , Everest T. Ong (4) , & J. Walter Larson (2,3,5,6) Affiliations: (1) Department of Geosciences,


slide-1
SLIDE 1

PyMCT and PyCPL: Refactoring CCSM Using Python

Michael Tobis(1), Michael Steder(1), Robert L. Jacob(2,3), Raymond T. Pierrehumbert(1), Everest T. Ong(4), & J. Walter Larson(2,3,5,6) Affiliations: (1) Department of Geosciences, University of Chicago; (2) Mathematics and Computer Science Division, Argonne National Laboratory (3) Computation Institute, University of Chicago (4) Department of Atmospheric and Oceanic Sciences, University of Wisconsin (5) ANU Supercomputer Facility, The Australian National University (6) Australian Partnership for Advanced Computing (APAC) Presented at the Python Workshop at the Biennial Computational Techniques and Applications (CTAC ‘06) Townsville, Queensland, Australia 6 July 2006 1
slide-2
SLIDE 2

CTAC ’06 Python Workshop

Overview

  • Climate Models and the Parallel Coupling
Problem
  • More specifically, CCSM, MCT, and CPL6
  • Enter Python
  • More than you probably want to know
about MCT
  • Python Bindings for MCT, a.k.a. pyMCT
(plus example)
  • Re-implementation of CPL6, a.k.a. pyCPL
  • Conclusions + Future Work
2
slide-3
SLIDE 3

CTAC ’06 Python Workshop

ACT I: Parallel Coupled Models and CCSM

3
slide-4
SLIDE 4

CTAC ’06 Python Workshop

Schematic (Directed Graph) for the Climate System

ATM OCN LAND SEA- ICE

Nodes represent subsystem models Arcs represent data to be delivered from source to target (states and fluxes) 4
slide-5
SLIDE 5

CTAC ’06 Python Workshop

Complexity Barriers

The traditional approach was: Model the individual subsystems (atmosphere, ocean, sea-ice, and land) in isolation
  • Idealize interactions with outside world through prescription (e.g.,
climatological or time-series data) or simplified physics (e.g., bucket hydrology, swamp or mixed-layer ocean)
  • Why? Three complexity barriers to overcome:
  • Knowledge, a consequence of specialization
  • overcome through interdisciplinary teams
  • Computational, i.e., getting all the math done
  • overcome by faster processors, better algorithms, and parallel
computing
  • Software: build system, language barriers, interactions between
physics packages (very important here) 5
slide-6
SLIDE 6

CTAC ’06 Python Workshop

One of Life’s Little Ironies...

  • The solution to one problem can create a new

problem

  • In this case, parallel computing in the form of

message-passing parallelism (let’s face it--MPI is pretty much the standard) helps surmount the computational complexity barrier

  • Distributed-memory parallelism complicates

intercourse between coupled models due to the traffic in distributed data

  • This is the Parallel Coupling Problem
6
slide-7
SLIDE 7

CTAC ’06 Python Workshop

Parallel Coupling Problem

Given: N mutually interacting models (C1,C2,...CN), each of which may employ message-passing parallelism Goal: Build an efficient parallel coupled model Aspects of the problem:
  • Architecture
  • Parallel data processing
  • Environment
  • Language barriers
  • Build issues
Often viewed simplistically as the “M x N” data transfer problem 7
slide-8
SLIDE 8

CTAC ’06 Python Workshop

Architectural Aspects

Two sources that shape parallel coupled models:

  • Science of the system under study:
  • Connectivity--who talks to whom?
  • Coupling event scheduling (e.g., periodic?)
  • Domain overlap--lower-dimensional vs. colocation
  • Timescale separation/interaction & domain overlap
  • Tightness
  • Implementation choices:
  • Resource allocation
  • Scheduling of model execution
  • Number of executable images
  • Mechanism
This presents a large and difficult-to-analyze decision space (see Larson’s talk from CTAC ‘06) 8
slide-9
SLIDE 9

CTAC ’06 Python Workshop

Parallel Data Processing

  • Description of data to be exchanged during coupling
  • Physical fields/variables
  • Mesh or representation associated with the data
  • Domain decomposition
  • Transfer of data--a.k.a. the MxN problem
  • Transformation of data
  • Intermesh interpolation/transformation between
representations (and associated conservation issues)
  • Time transformation
  • Diagnostic/variable transformations
  • Merging of data from multiple sources
9
slide-10
SLIDE 10

CTAC ’06 Python Workshop

CCSM

  • CCSM := Community Climate System
Model
  • Four physical components:
Atmosphere, Ocean, Sea-Ice, and Land- Surface models
  • One coupler component
  • Together, they comprise a hub and
spokes architecture (see box at right)
  • As implemented, CCSM is an
application-specific software framework 10
slide-11
SLIDE 11

CTAC ’06 Python Workshop

CCSM’s Coupler CPL6

Model Substitution Modification of fields under exchange via “Configurable Plugs”

CPL6 enables one to do--with aplomb-- some neat things to CCSM. 11
slide-12
SLIDE 12

CTAC ’06 Python Workshop

This is Great...But...

  • Alteration of CCSM’s architecture is more difficult,

and one must re-code the coupler MAIN in terms of the bits and pieces from the CPL6 toolkit+library--in Fortran :-(

  • Number of components
  • Changes in scientific details of the couplings beyond

the most simple alterations

  • e.g., currently take-it-or-leave-it temporal advance

coupling, not anything nearly as sophisticated as predictor/corrector But this is exactly what many scientists want to do!

12
slide-13
SLIDE 13

CTAC ’06 Python Workshop

CCSM’s Fortran Software Stack

13
slide-14
SLIDE 14

CTAC ’06 Python Workshop

Attempted Object- Oriented Programming in Fortran...or, Don’t Try This at Home...

14
slide-15
SLIDE 15

CTAC ’06 Python Workshop

OO Programming in Fortran-- Where’s the Swindle?

  • The introduction of Derived Types and Explicit Interfaces in the
Fortran90 standard have allowed one to implement some OO concepts (see Decyk et al., 1996 and 1997)
  • Encapsulation and Information Hiding
  • Inheritance
  • Polymorphism
  • But, one must implement them--they are not available as part and
parcel of the language, and the technique used is through a design pattern called delegation
  • The 2003 Fortran Standard introduces the notion of a Class, but
when this feature will be ubiquitous in compilers is unclear 15
slide-16
SLIDE 16

CTAC ’06 Python Workshop

MCT is a Collection of Fortran Datatypes...

...and a comprehensive set of support routines (a.k.a methods)

Data Transformation Data Description Data Transfer

KEY

AttrVect GlobalSegMap GeneralGrid Accumulator SparseMatrix SparseMatrixPlus MCTWorld Router Rearranger 16
slide-17
SLIDE 17

CTAC ’06 Python Workshop

Seven the Hard Way

  • MCT has seven classes that are linked via a class

hierarchy

  • MCT is an attempt to follow OO discipline as best as

possible in a Fortran context; i.e., provide a comprehensive set of support methods for each class

  • Inheritance was implemented the hard way through

hard-coding of these relationships (i.e., delegation)

  • This is the primary reason that MCT is such a small set
  • f classes, and why the datatypes we will encounter in

CPL6 are not in MCT

17
slide-18
SLIDE 18

CTAC ’06 Python Workshop

CPL6 Fortran Datatypes

Looks Like a Job for Python!

18
slide-19
SLIDE 19

CTAC ’06 Python Workshop

The Case for Python

  • Python offers true OO programming, thus CPL6 class

implementations are far easier to implement

  • Most of the complexity of the coupler is at initialization

time: scripting penalty is irrelevant

  • Most of the coupler loop is idle; calculations can be

farmed out to a compiled language ('Numeric' or 'numpy' modules)

  • Coupler itself is approximately 2kLOC of Fortran, re-

implementation in Python should be doable

  • Realization of long-term goal of roll-your-own/disposable

couplers

19
slide-20
SLIDE 20

CTAC ’06 Python Workshop

PyCCSM/PyCPL/PyMCT Parts List

  • Python bindings for MCT, a.k.a. pyMCT
  • Python-based reimplementation of Multi-Process

Handshaking utilities (MPH)

  • Implementation of CPL6 datatypes in Python (as

true classes), a.k.a. pyCPL

  • Python implementation of a specific coupler +

legacy component models = pyCCSM

20
slide-21
SLIDE 21

CTAC ’06 Python Workshop

ACT II: MCT

21
slide-22
SLIDE 22

CTAC ’06 Python Workshop

  • Q. Why go into so much detail about

MCT?

  • A. Because PyMCT is the clear non-

climate spin-off product from this project

22
slide-23
SLIDE 23

CTAC ’06 Python Workshop

MCT’s Universe of Discourse

  • Support coupling of MPI-based MPP models
  • Data transfer using a peer commmunication model
  • Description of physical meshes and associated field data
through linearization
  • Data transfer and transformation are viewed as multi-field,
pointwise operations
  • We leave numerous, high-level operations to the user’s
discretion (e.g., choice of linearization and interpolation schemes), while concentrating on automation of complex (but important!) low-level operations 23
slide-24
SLIDE 24

CTAC ’06 Python Workshop

Coupling as Peer Communication

  • MCT’s organizing principle is the component model,
  • r component (not same as CORBA, CCA, ESMF,

JavaBeans)

  • An MCT component is merely a model that is part
  • f the larger system and participates in coupling
  • In MCT, components interact directly as peers
  • The user codes these connections into the model

source

24
slide-25
SLIDE 25

CTAC ’06 Python Workshop

Linearization of Multi- Dimensional Space

  • Linearization (first used in MxN schemes by UMD’s

METACHAOS) is the mapping from an n-tuple index space to a single global location index

  • This approach allows for a single representation of

grids/arrays of aribtrary dimension

25
slide-26
SLIDE 26

CTAC ’06 Python Workshop

Index Space

  • Definition: An M-dimensional index space is

subset of ZM , each element of which can be uniquely identified by an M-tuple of integers (i1, i2,...,ik,..., iM)

  • This is what we normally use to describe

Cartesian meshes and associated field data stored in multidimensional arrays

26
slide-27
SLIDE 27

CTAC ’06 Python Workshop

Example: 2-D Cartesian

y-index j 1 2 3 1 2 3 4 x-index i

27
slide-28
SLIDE 28

CTAC ’06 Python Workshop

Index Space: Ordered Pairs

x y 1 2 3 1 2 3 4 (0,1) (1,1) (1,2) (1,3) (2,1) (2,2) (2,3) (3,1) (3,2) (3,3) (4,3) (4,2) (4,1) (0,0) (0,2) (0,3) (1,0) (2,0) (3,0) (4,0)

28
slide-29
SLIDE 29

CTAC ’06 Python Workshop

Virtual Linearization of Multidimensional Index

Definition: A virtual linearization f of an index space A is a mapping f : A —›B

  • where A is a subset of ZM
  • B is a subset of N the set of natural numbers
  • f is one-to-one and onto
  • That is, f maps an M-tuple of integers to a unique

natural number

29
slide-30
SLIDE 30

CTAC ’06 Python Workshop

Simple Example Virtual Linearizations

Valid for i = {0,1,2,3} j = {0, 1, 2, 3, 4}

30
slide-31
SLIDE 31

CTAC ’06 Python Workshop

f(i,j) = i + 5j + 1

x y 1 2 3 1 2 3 4 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1

31
slide-32
SLIDE 32

CTAC ’06 Python Workshop

f(i,j) = 4i + j + 1

x y 1 2 3 1 2 3 4 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1

32
slide-33
SLIDE 33

CTAC ’06 Python Workshop

Pointwise Operations

  • FACT: Most coupling between components

involves multivariate exchanges (e.g., ~10-12 fields between atmosphere and ocean in a climate model)

  • Consolidating operations on data to be

exchanged can result in performance boosts due to better cache re-use and lowered MPI latency charges

  • In MCT, implementation choices (primarily

storage order) have been made to exploit this situation

33
slide-34
SLIDE 34

CTAC ’06 Python Workshop

Automating (Most of) the Tough Stuff

  • We leave coupled model architecture decisions to the user
  • We aim to be minimally invasive--user still owns main(), still calls
MPI_Init(), et cetera...
  • The user decides how and when to couple
  • The user must describe application grids and decompositions using
MCT’s GeneralGrid and GlobalSegMap classes
  • The user must pack coupling data in AttrVect form
  • MCT’s library routines do the heavy lifting
34
slide-35
SLIDE 35

CTAC ’06 Python Workshop

Data Description: Domain Decomposition

  • Embodied in the GlobalSegMap class
  • Linearization creates a single unique index for each
location, and runs of consecutive indices--or segments--are stored by start index, segment length, and processor ID on which it resides (on the component’s communicator)
  • Capable of supporting arbitrary decompositions and haloed
decompositions
  • Support for global-to-local and local-to-global index
translation 35
slide-36
SLIDE 36

CTAC ’06 Python Workshop

Data Description: Physical Meshes

  • Embodied in the GeneralGrid class
  • Linearization implies this class stores real and integer
attributes of individual mesh points
  • Coordinates, length, cross-sectional area, and volume
elements, integer and real masks, grid indices
  • Able to describe meshes of arbitrary dimensionality and
also unstructured meshes
  • Method support for sorting points in lexicographic order
by coordinates and/or other attributes 36
slide-37
SLIDE 37

CTAC ’06 Python Workshop

Data Description: Field Data

  • Embodied in MCT AttrVect class, which is similar to Trilinos
multi-vector
  • Allows storage of both integer and real attributes
  • In implementation, major index is field index, keeping in line with
pointwise approach
  • Attributes are accessed via use of string tokens
  • Lists of attributes are set at initialization, which allows run-time
binding of what fields are stored
  • Numerous query and manipulation methods including importing and
exporting attributes and MergeSort keyed by attributes 37
slide-38
SLIDE 38

CTAC ’06 Python Workshop

Data Transfer: Registration

  • Singleton class MCTWorld holds a lightweight

component registry

  • Components are registered with a component ID
  • MCTWorld contains a rank translation table that

allows a component to message another remote component by using information from its local domain decomposition descriptor (i.e., we do not construct nor use intercommunicators)

38
slide-39
SLIDE 39

CTAC ’06 Python Workshop

Data Transfer: Communications Scheduling

  • Scheduling for one-way parallel data

transfers are handled by the Router class

  • Scheduling for two-way parallel data

transfers and parallel data redistributions are handled by the Rearranger class, which comprises two Routers

39
slide-40
SLIDE 40

CTAC ’06 Python Workshop

Data Transfer: Execution

  • MxN Communications between components are carried out by library
routines MCT_Send() and MCT_Recv(), both of which have (overall) blocking and non-blocking versions
  • Overall blocking--Prevents concurrently scheduled components from
  • utrunning each other
  • Overall nonblocking--allows for same communications approach to be
used for sequentially scheduled components
  • N.B.: non-blocking MPI operations are used for the individual point-to-
point messages within the MxN transfer operation
  • MxM redistributions are handled by the MCT library routine Rearrange()
  • Contains logic to prevent self-messaging (does a copy instead)
40
slide-41
SLIDE 41

CTAC ’06 Python Workshop

Data Transformation: Intermesh Interpolation

  • MCT’s linearization-based worldview casts interpolation as a linear transform and thus
implementable as a matrix-vector multiply
  • In practice, these interpolation matrices are quite sparse
  • MCT’s SparseMatrix class provides storage for nonzero matrix elements in COO
format
  • Method support includes query methods, computation of sparsity, sorting methods to aid
matrix decomposition and to boost performance
  • Library function performs multiplication y = Mx, where x and y are of AttrVect type, and
this function performs automatic token-based matching of attributes
  • Main implementation targets commodity cache-based processor platforms, and MCT’s
pointwise operation view makes good re-use of cache
  • Modifications for vector platforms also included
41
slide-42
SLIDE 42

CTAC ’06 Python Workshop

MCT’s Parallel Linear Transformation

  • MCT provides the SparseMatrixPlus class, which
encapsulates everything needed for a sparse linear transform
  • This class includes Rearrangers to schedule communications
and a SparseMatrix to store matrix elements
  • Library routine to compute y = Mx in parallel, again where x
and y are of AttrVect type, and automatic token-based matching of attributes is performed
  • Support for both commodity and vector processors
42
slide-43
SLIDE 43

CTAC ’06 Python Workshop

Data Transformation: Time Accumulation

target t = tc source target t = tc source target t = tc source target t = tc source target (e) c source (a) (b) (c) (d) t = t
  • Time evolution of multi-component systems can
involve data exchanges of instantaneous or/or integrated data
  • For instantaneous exchanges, use of one or more
AttrVects is sufficient
  • For integrated data exchanges, MCT offers
accumulation registers in the form of the Accumulator class, and the accumulate() library routine
  • The Accumulator provides registers for time
integration and averaging of data, and keeps track of progress over an accumulation cycle
  • The accumulate() library routine works with
AttrVect and Accumulator arguments, and automatically cross-indexes and accumulates attributes with matching tokens 43
slide-44
SLIDE 44

CTAC ’06 Python Workshop

Other Services--Spatial Integration and Averaging

  • MCT provides routines for computing spatial integrals and
averages
  • These routines are included to diagnose (and if wished, even
enforce) conservation of integrals in the interpolation process
  • Library routines operate on AttrVect objects
  • Spatial weight elements and masks can be provided either in
GeneralGrid or array form
  • Paired integral/average routines to compute integrals
simultaneously on source and target grids to minimize global sum latency costs 44
slide-45
SLIDE 45

CTAC ’06 Python Workshop

Other Services--Merging of Data from Multiple Sources

  • Often one must combine outputs from multiple components
for use as input for another component (e.g., fluxes/states from
  • cean, land, and sea-ice for use by atmosphere)
  • MCT provides a Merge facility, which is a set of routines to
combine multiple AttrVect data streams into a single AttrVect result
  • Real and Integer masks for the merge can be supplied either in
GeneralGrid or array form
  • Automatic token-based attribute matching
45
slide-46
SLIDE 46

CTAC ’06 Python Workshop

MCT’s Fortran Programming Model

Based on Fortran90, but applicable to

  • ther languages (including Python).
  • Use modules for access to MCT

classes and methods

  • Declare variables of MCT datatypes
  • Invoke MCT library routines to

accomplish parallel coupling operations

46
slide-47
SLIDE 47

CTAC ’06 Python Workshop

ACT III: pyMCT, pyCPL, and pyCCSM

47
slide-48
SLIDE 48

CTAC ’06 Python Workshop

OK...So Where’s the Python?

  • Note so far that there hasn’t been a single line of
Python code.
  • Don’t worry, we’ll get there; please be patient:-)
  • MCT of course is Fortran90
  • Interfacing Fortran 90/95 to Python is much harder
than binding Python to f77
  • This is entirely the fault of the ANSI Fortran J3 and
ISO WG5 Fortran committees
  • If this makes you angry, go to
http://fortranstatement.com
  • So, let’s first discuss this problem...
“Ever had the feeling you’ve been cheated?”
  • Johnny Rotten
48
slide-49
SLIDE 49

CTAC ’06 Python Workshop

Fortran “Dope Vectors”

  • The Fortran90 Standard introduced the notion of

allocatable arrays and pointers

  • The standard was utterly silent on how compiler

developers should implement array descriptors

  • Furthermore, the standard failed to offer a

standardized API to query them

  • This creates a BIG PROBLEM for language

interoperability

  • The Fortran2003 standard promises to solve this with

its BINDC feature

49
slide-50
SLIDE 50

CTAC ’06 Python Workshop

Enter Babel

  • You can wait until the Fortran 2003 standard is implemented and use
BIND C and go from there, but don’t hold your breath
  • The lag between the standard specification and its widespread
availability in compilers is usually 5-7 years
  • There is a solution available today: Babel
  • Babel is a general language interoperability tool that automates
construction of interlanguage glue to connect codes written in C, C++, f77, Fortran, Java, and Python
  • Babel accomplishes this by processing user-authored interface
descriptions in Scientific Interface Description Language (SIDL) 50
slide-51
SLIDE 51

CTAC ’06 Python Workshop

The Babel Recipe

To accomplish calling from language A to code written in language B,
  • ne must:
  • 1. Write a SIDL interface description for the code implemented in
language B
  • 2. Run babel to create
  • Skeletons in Language A (also called .IMPL files)
  • The Internal Object Representation (IOR), which is implemented
in C
  • 3. Paste into the skeletons either calls to the package in Language A, or
maybe the whole source
  • 4. Run babel again to create call-in stubs in Language B
  • 5. Use and enjoy the Language B stubs from your source code written
in Language B 51
slide-52
SLIDE 52

CTAC ’06 Python Workshop

But...Beware!

  • Babel builds are complex, and associated

language binding builds are confusing

  • We found the best environment for

confronting these problems was a cluster

  • ver which we had administrative control!
  • MCT’s restricted SIDL API had bugs
  • Fortran’s 1-based indexing vs. Python’s 0-

based indexing required modification to MCT

52
slide-53
SLIDE 53

CTAC ’06 Python Workshop

Python and MPI Issues

  • There exist numerous Python MPI binding

implementations (pyMPI, MPI-B, ANU’s pyPar...)

  • But, we have specific needs in messaging between

Python and non-Python codes, as well as the need to support multiple executable images

  • Required a custom re-wrapping of MPI to confront

language interoperability+MPI issues

  • called MMPI (http://www.penzilla.net/mmpi)
  • Also required an upgrade to MPI-2 (specifically, we

used mpich2)

53
slide-54
SLIDE 54

CTAC ’06 Python Workshop

(Finally!) Some Python in the form of a pyMCT Example

54
slide-55
SLIDE 55

CTAC ’06 Python Workshop

What It Is

  • Python implementation of a simple MCT example coupling
generic models to a hub coupler (yes, it uses MPI)
  • In this case, component scheduling is concurrent
  • Parts List:
  • 1. Driver (master.py)
  • 2. Generic Model (model.py)
  • 3. Coupler (coupler.py)
  • 4. Downloadable from MCT Web site:

http://www.mcs.anl.gov/mct

55
slide-56
SLIDE 56

CTAC ’06 Python Workshop

master.py (part 1)

# Check all necessary modules are available # Exits and prints an error if components aren't available import modulecheck # import standard modules import sys, traceback import Numeric # Courtesy of pyMPI import mpi # MCT from MCT import World,GlobalSegMap,AttrVect,Router # Example Codes import model import coupler 56
slide-57
SLIDE 57

CTAC ’06 Python Workshop

master.py -- main()

def main(): rank = mpi.rank size = mpi.size # Split MPI_COMM_WORLD into a communicator for each model if (rank == 0): print "Running PyMCT Example on %d Processors"%(size) ncomps = size color = (rank % ncomps) + 1 splitcomm = mpi.comm_create([mpi.WORLD[rank]]) # More complex test # Start up the models cpl = coupler.Coupler() pop = model.Model() csm = model.Model() csim = model.Model() clm = model.Model() if( color == 1 ): print "COUPLER COMPONENT STARTING ON PROCESSOR %d"%(rank) cpl.start(splitcomm, ncomps, color) elif( color == 2 ): print "COUPLER COMPONENT STARTING ON PROCESSOR %d"%(rank) pop.start(splitcomm, ncomps, color) elif( color == 3 ): print "COUPLER COMPONENT STARTING ON PROCESSOR %d"%(rank) csm.start(splitcomm, ncomps, color) elif( color == 4 ): print "COUPLER COMPONENT STARTING ON PROCESSOR %d"%(rank) csim.start(splitcomm, ncomps,color) elif( color == 5 ): print "COUPLER COMPONENT STARTING ON PROCESSOR %d"%(rank) clm.start(splitcomm, ncomps, color) else: print "! COLOR(%d) ERROR ON PROCESSOR %d"%(color,rank) """ """ while(1): atmdata = atmos.step() if (timestep % 5 == 0): if( mpi.rank == 0): # Im the coupler coupler.couple() if( mpi.rank == 6): atmos.couple() """ return if __name__=="__main__": main()

Delegates execution to appropriate model

57
slide-58
SLIDE 58

CTAC ’06 Python Workshop

model.py

""" model.py This model requires a square processor decomposition i.e., a 2x2, 3x3, 4x4 grid. """ # This checks that all modules are available # and prints errors about components. #import modulecheck # Actual imports import traceback import mpi from MCT import World,GlobalSegMap,AttrVect,Router import Numeric

Module Imports

58
slide-59
SLIDE 59

CTAC ’06 Python Workshop

model.py, Cont’d...

class Model: def __init__(self): return def start(self, comm, ncomps, compid): """ start( Communicator comm, int ncomps, int compid ) This method is called to start the model instance. """ mctWorld = World.World() mctGlobalSegMap = GlobalSegMap.GlobalSegMap() mctAttrVect = AttrVect.AttrVect() router = Router.Router() myrank, mysize = comm.comm_rank(),comm.comm_size() # Register self with MCT.World print "PROCESSOR:%d INITIALIZING AS COMPONENT %d of %d"%(mpi.WORLD.comm_rank(),compid,mpi.WORLD.comm_size()) mctWorld.initd0(ncomps, mpi.WORLD, comm, compid) print "COMPID %d INITIALIZED!"%(compid) # Initialize the global segmap print "COMPID:", compid, "INITIALIZING GLOBALSEGMAP" nx = 128 ny = 64 localsize = (nx*ny)/mysize 59
slide-60
SLIDE 60

CTAC ’06 Python Workshop

model.py, Cont’d...

# Create arrays for start and length # MCT expects Fortran arrays starting at 1 start = Numeric.zeros(2,Numeric.Int32) start[1] = ((myrank*localsize)+1) length = Numeric.zeros(2,Numeric.Int32) length[1] = (localsize) # Describe Decomposition with MCT GSMAP TYPE print "COMPID:",compid,"INITIALIZING GLOBALSEGMAP" mctGlobalSegMap.initd0(start, length, 0, comm, compid) # Initialize an ATTRIBUTE VECTOR print "COMPID:",compid,"INITIALIZNNG ATTRVECT" len = mctGlobalSegMap.lsize(comm) mctAttrVect.init("field1:field2","field1:field2",len) print "COMPID:",compid,"INITIALIZED ATTRVECT (SIZE=",len,")" # Fill AV with DATA using the import function # MCT expects fortran arrays starting at 1 print "COMPID",compid,"IMPORTING ATTRIBUTES" field1 = Numeric.zeros(len) field2 = Numeric.zeros(len) for i in xrange(len): field1[i] = 1.0 / len field2[i] = 0.5 mctAttrVect.importRAttr("field1",field1) mctAttrVect.importRAttr("field2",field2) 60
slide-61
SLIDE 61

CTAC ’06 Python Workshop

model.py, Cont’d...

# Initialize a ROUTER print "COMPID:",compid,"INITIALIZING ROUTER" router.init(1,mctGlobalSegMap,comm) # Send data over the router print "COMPID:",compid,"SENDING ATTRVECT TO COMP 1" tag = 888+compid print "COMPID:",compid,"TAG:",tag mctAttrVect.send(router, tag) # Deallocate and clean up: mctAttrVect.clean() mctGlobalSegMap.clean() router.clean() mctWorld.clean() # DONE print "COMPID:",compid,"FINISHED" return 61
slide-62
SLIDE 62

CTAC ’06 Python Workshop

coupler.py

""" Mike Steder 10/20/04 """ #First make sure all the necessary modules are available #import modulecheck # Now handle actual imports import sys, traceback,math import Numeric import mpi from MCT import World,GlobalSegMap,AttrVect,Router,SparseMatrix,SparseMatrixPlus 62
slide-63
SLIDE 63

CTAC ’06 Python Workshop

coupler.py, Cont’d...

class Coupler: def __init__(self): return def start(self,comm, ncomps, compid): # Get local rank and size myrank,mysize = comm.comm_rank(),comm.comm_size() # Initialize MCTWORLD print "PROCESSOR:%d INITIALIZING AS COMPONENT %d of %d"%(mpi.WORLD.comm_rank(), compid,mpi.WORLD.comm_size()) mctworld = World.World() mctworld.initd0(ncomps, mpi.WORLD, comm, compid) print "COMPID:%d WORLD INITIALIZED!"%(compid) # Set up latitude decomposition nx, ny = 128,64 localsize = (nx*ny)/mysize # Create start and length arrays start = Numeric.zeros(2,Numeric.Int32) start[1] = (myrank*localsize)+1 length = Numeric.zeros(2,Numeric.Int32) length[1] = localsize 63
slide-64
SLIDE 64

CTAC ’06 Python Workshop

coupler.py, Cont’d...

# Describe decomposition with MCT Global Seg Map print "COMPID",compid,"INITIALIZING ATMOSPHERE GLOBALSEGMAP" AtmGSMap = GlobalSegMap.GlobalSegMap() AtmGSMap.initd0(start,length, 0,comm, compid) print "COMPID",compid,"INITIALIZED ATMOS GLOBALSEGMAP" # Initialize a globalsegmap for the ocean # By setting up a checkboard grid decomposition nx,ny = 320,384 nxprocs = math.sqrt(mysize) nyprocs = math.sqrt(mysize) if(mysize==2): # Special Case Layout nxprocs=nyprocs=1 # Number of lat/lon's in *my* grid lonsize = nx / nxprocs latsize = ny / nyprocs rowindex = myrank/nxprocs colindex = myrank/nxprocs # Lower left vertex of grid box boxvertex = (colindex*lonsize+1) + (latsize*rowindex*nx) # Create more start/length arrays start = Numeric.zeros(latsize+1,Numeric.Int32) length = Numeric.zeros(latsize+1,Numeric.Int32) for i in range(0,int(latsize)): # Was (i-1) in C++, but the indexing is different here. start[i] = boxvertex+(i)*nx length[i] = lonsize print "COMPID",compid,"INITIALIZING OCEAN GLOBALSEGMAP" OcnGSMap = GlobalSegMap.GlobalSegMap() OcnGSMap.initd0(start, length, 0, comm, compid) 64
slide-65
SLIDE 65

CTAC ’06 Python Workshop

coupler.py, Cont’d...

# Initializng an attribute vector for the atmosphere # grid and ocean grid print "COMPID",compid,"INITIALIZING ATTRVECT" AtmAV = AttrVect.AttrVect() AtmAV.init("field1:field2","field1:field2",AtmGSMap.lsize(comm)) print "COMPID",compid,"INITIALIZED ATTRVECT (SIZE=",AtmAV.lsize(),")" # Setup a SPARSEMATRIX with some test data #print "COMPID",compid,"CREATING SPARSEMATRIX..." nRows = 320 * 384 nCols = 128 * 64 num_elements = 145548 rows = Numeric.zeros(num_elements,Numeric.Int32) columns = Numeric.zeros(num_elements, Numeric.Int32) weights = Numeric.zeros(num_elements, Numeric.Float32) 65
slide-66
SLIDE 66

CTAC ’06 Python Workshop

coupler.py, Cont’d...

print "COMPID",compid,"ENTERING RECV LOOP ..." for i in range( 0, ncomps-1 ): # Adjust from 0 index to component index i = i + 2 Rout = Router.Router() Rout.init(i, AtmGSMap, comm ) print "COMPID",compid,"RECEIVING ATTRVECT FROM COMP",i tag = 888 + i sum = False print "COMPID",compid,"TAG:",tag AtmAV.recv(Rout, tag, sum) # Remove for now due to buggy ness # Interpolate with a Parallel Sparsematrix Multiply #print "COMPID",compid,"INTERPOLATING DATA FROM COMPID",i #OcnAV.sMatAvMult_sMPlus( AtmAV, A2OMatPlus ) Rout.clean() print "COMPID",compid,"PRINTING FIELD2" field2 = AtmAV.exportRAttr("field2") #print field2 #sum = 0 #for i in range(AtmGSMap.lsize()): # if (field2.get(i) != 0.5): # print "CORRUPT DATA IN ATMAV =", field2.get(i) print "COMPID",compid,"FINISHED" return 66
slide-67
SLIDE 67

CTAC ’06 Python Workshop

About pyCCSM

  • The whole CPL6 toolkit and library layer was re-

implemented using Python

  • After that, the CPL6 main was re-written in Python,

and the legacy models (atmosphere, ocean, sea-ice, and land-surface) were left as-is

  • This all now works for CCSM’s “dead” and “data”

models

  • Live models now “close to working”
  • So far, the performance penalty is negligible
67
slide-68
SLIDE 68

CTAC ’06 Python Workshop

What We Have Learned

  • Control level of a big multiphysics application like CCSM should be
implemented as a script using a very high level of abstraction (unproven, but we are confident this will be the case)
  • MCT provides an excellent basis for a scalable and logically
decoupled middleware layer
  • Fortran 90/95 arrays still provide an considerable language
interoperability hurdle in spite of tools like Babel--either wait for BIND C in Fortran2003 or stick with f77 and pyFort or similar software if you can
  • CCSM’s initiative to go towards from a multiple to a single
executable image application is unnecessary and in fact may be counterproductive for lightweight coupling mechanisms 68
slide-69
SLIDE 69

CTAC ’06 Python Workshop

Where This is Going

  • The US NSF ITR grant that funded this work

will end later this year

  • We still hope to get a working version of

pyCCSM by the project’s end

  • This project has produced two useful spin-
  • ff software products--pyMCT and MMPI
  • We hope these packages will see public

release

69
slide-70
SLIDE 70

CTAC ’06 Python Workshop

FIN

70