MMSP User's Guide
The page is a meant to provide both a tutorial-style introduction to MMSP as well as a reasonably complete programming reference. Links to the actual MMSP source code are given where appropriate.
Contents
1.1 The MMSP concept
1.2 Prerequisites
1.3 MMSP components
1.3.1 Point data structures
1.3.2 Grid objects
1.3.2 Mesh objects
1.3.3 The update() function
2.1 Single processor simulations
2.2 Parallel processor simulations
2.3 Microstructure generation
2.4 Microstructure analysis
2.4 File formats
2.5 Visualization
3.1 Point data structures
3.2 Grid objects
3.3 The update() function
4.1 MCgrid2D
4.2 MCgrid3D
4.3 PFgrid2D
4.4 PFgrid3D
4.5 sparsePF2D
4.6 sparsePF3D
4.7 CAgrid2D
4.8 CAgrid3D
4.9 FDgrid2D
4.10 FDgrid3D
1 Introduction to MMSP
The Mesoscale Microstructure Simulation Package (MMSP) is nothing more than a collection of C++ header files that declare a number of grid or mesh objects (classes) and define how most of their methods (member functions) are implemented. The following sections introduce the concepts behind MMSP's design, knowledge prerequisites for using and extending the MMSP package, and a brief description of MMSP's core components.
1.1 The MMSP concept
The design of MMSP is based on several observations about how mesoscale simulations are used by materials scientists:
- Most mesoscale simulations discretize the spatial domain using a rectilinear grid or a finite element mesh. A data structure is associated with each grid or mesh node that has a particular size (scalar, vector, etc.) and value type (integer, floating point) depending on the simulation method.
- Most simulations update the data structure at each node in a way that falls under some common methodology (Monte Carlo, phase field, etc.) but has features unique to each given physical process.
What this roughly means is that most mesoscale simulations use a common spatial discretization, but we usually need to tweak the details of how we represent spatial data and how we update it. The unfortunate truth is that, typically, a researcher wishing to model a particular physical process produces code with a focus mainly on the process itself (the "tweak"), largely ignoring the problem of how to design reusable data structures. What happens when we decide that we should try a different simulation method, or when we realize that we need to use a parallel implementation? A top-down design doesn't help us much here.
MMSP is meant to help keep its users from reinventing the grid and mesh data structures, file I/O, parallelization, boundary condition handling, etc. while retaining enough flexibility to model a large number of physical processes. The goal of MMSP is to provide a simple, consistent, and extensible programming interface for all grid–and mesh-based microstructure evolution methods. Simple means that the package has a very small learning curve, and for most routine simulations, only a minimal amount of code must be written. By consistent we mean, for example, that 2D simulation code is nearly identical to that for 3D simulations, single processor programs are easily parallelized, and fundamentally different methods like Monte Carlo or phase field have the same "look and feel." Finally, extensible means that it’s straightforward to add new grid types or physical behaviors to the package. Other considerations include efficiency, and to a lesser extent, portability (MMSP is written entirely in ISO compliant C++). The rest of this user's guide will help you write programs that benefit from the MMSP design.
1.2 Prerequisites
MMSP is not the kind of software that reads a few parameters or a file specified by the user and cranks out some generic computation. In fact, MMSP relies on the user to provide code for all of the real physics that a simulation is intended to capture. This isn't as scary as it sounds, as MMSP was designed to make this procedure as simple as possible. See the section titled "The update() function" and the [worked MMSP example] for more.
1.3 MMSP components
The core of MMSP is its collection of point data structures, grid and mesh objects implemented as C++ template classes. This section introduces the basic components (classes) defined within MMSP.
1.3.1 Point data structures
Point data structures are template classes which represent the data associated with each node of a rectilinear grid or finite element mesh (see the first observation above). There are currently three such classes: scalar, vector, and sparse. The first two probably agree with your intuition - a scalar is a single value and a vector is an array of values. The sparse data structure is used in cases where we might otherwise need a vector data structure with a large number of elements, except that the value of only a small number of components is relevant. Each data structure contains data of a single fundamental numeric type (int, float, double, etc.), which is the template parameter. See MMSP.data.hpp for the implementation of each point data structure.
1.3.2 Grid objects
The MMSP base class grid is defined in the header MMSP.grid.hpp. The derived template classes grid1D, grid2D, and grid3D provide the base for the high level grid objects associated with various simulation methods. The template parameters for each grid object are a point data structure and a fundamental numeric type. The grid objects contain a large number of inherited methods (member functions) that perform common tasks such as file I/O, data storage and retrieval, etc. Member data for the grid classes determine the grid size, boundary conditions and point spacing along each coordinate axis.
High level grid objects usually have a name that indicates both their intended use as well as their dimensionality, e.g. PFgrid3D and MCgrid2D. These are the grid objects that are used in performing simulations or other computations. Detailed information about each existing grid class is given in the reference section of this page. See also the source distribution.
MMSP grid classes are designed so that programming with them is as intuitive as possible. For example, data can be accessed by subscript operators in exactly the same way we would use data stored in a subscripted array
// First create a grid object. PFgrid2D grid(100,100,2); // Assign a value to the phase field grid. grid[10][20][1] = 1.0; // Use a value from the grid in a computation. float value = grid[10][20][1];
Because of this, C/C++ code that performs computations using "hard-coded" arrays automatically works for the corresponding MMSP grid object. This makes it trivial to insert MMSP grid objects into existing grid generation, simulation, or analysis code.
Each grid object also provides a neighbor() function that returns a grid data structure using coordinates relative to a given position (x,y[,z])
// First create a grid object. MCgrid2D grid(100,100); // Assign a value to the Monte Carlo grid. grid[10][20] = 5; // Access the value using the neighbor() function. int spin = grid.neighbor(9,22,1,-2);
The neighbor function automatically accounts for the grid boundary conditions set by the user. See the next subsection for more details on programming with the grid classes.
Finally, let's consider using grid objects in a parallel computing setting. In single processor programs, we typically create a single grid object within the main() function. This grid object is considered to be the entire simulation domain. In parallel programs, however, the grid object declared within main() is a subgrid of a larger, global simulation domain. MMSP can automatically decompose (or recompose) a global grid object into local subgrid "slabs" with cuts perpendicular to the x axis, a method which has a consistent (and simple) implementation regardless of grid dimension. MMSP allows the user to chose the number of ghost cells retained in each subgrid, and the boundary conditions of the global grid object are satisfied when the neighbor() function is used. The implementation of parallel grid functionality is contained within a separate header MMSP.grid.parallel.hpp for the convenience of users without a local MPI distribution.
1.3.3 Mesh objects
1.3.4 The update() function
In MMSP, the update() function is the one and only function that determines how the microstructure evolves. It is a member function of each grid or mesh object, but its implementation is always written by the user. The update() function is overloaded for single and parallel simulations. For every grid object, the update() member functions are declared as
update(int steps = 1); update(int steps, int id, int np, int ng = 1);
The first form is for single processor programs. Without any arguments, update() performs a single time step (however that might be defined), and with the argument steps, update() performs a given number of time steps. The second form is for parallel simulations. The arguments id and np refer to the process rank and total number of processes, respectively. The optional argument ng is the number of ghost cell layers that separates each subgrid. Calling update without ng implies a single ghost cell layer.
In most cases, changing the physical process that a simulation models only involves changing the update() function. This is usually quite simple. As an example, let's consider the FDgrid3D class, which is used for general-purpose 3D finite difference computations. A generic update() function for the FDgrid3D class in a single processor implementation looks like (see FDgrid.update.hpp)
void FDgrid3D::update(int steps) { FDgrid3D& grid = *this; FDgrid3D update(nx[0],nx[1],nx[2]); for (int step=0; step<steps; step++) { for (int x=0; x<nx[0]; x++) for (int y=0; y<nx[1]; y++) for (int z=0; z<nx[2]; z++) { // Compute update here ... update[x][y][z] = grid[x][y][z]; } grid.swap(update); } }
Depending on your own point of view, this is about as simple as we could possibly ask for. Obviously, this function as it's written does nothing but replace the old grid with a copy of itself, but hopefully you get the idea. In this example, we've introduced an additional member function, swap(), that helps facilitate replacing the old grid with the new. Now, for a parallel implementation, update() looks something like (see FDgrid.update.parallel.hpp)
void FDgrid3D::update(int steps, int id, int np, int ng) { FDgrid3D& grid = *this; FDgrid3D update(nx[0],nx[1],nx[2]); for (int step=0; step<steps; step++) { for (int x=ng; x<nx[0]-ng; x++) for (int y=0; y<nx[1]; y++) for (int z=0; z<nx[2]; z++) { // Compute update here ... update[x][y][z] = grid[x][y][z]; } grid.swap(update); grid.ghostswap(id,np,ng); } }
Notice that, in the nested loops, only the limits in x have changed. This is because MMSP subdivides the global grid object into slabs along the x direction, and the update on the subgrid only needs to be computed for non-ghost cells. We've also used the member function ghostswap() to trade ghost cells between processes. It should be clear that to correctly swap ghost cells, the update() function needs to know at least what processor it's running on, how many other processes are currently running, and the number of cells to swap, explaining why all three are function parameters for the parallel version. Also note that there's very little difference between single and parallel implementations, which is one of the goals of MMSP design. With MMSP, changing a single processor code to a parallel processor code typically involves only a copy and paste operation on the code inserted under "Compute update here". In fact, we've included files with generic update() functions for each grid or mesh class in the source distribution so that the copy and paste operation is probably all you'll ever need to do to parallelize your code. Additionally, several files are included with example update() functions that actually do something useful.
2 Using MMSP
This section is a guide to using MMSP to write new programs, with a focus on simple microstructure evolution codes. We assume that the user will be using only existing components; later on we'll describe how to add new components to the package.
2.1 Single processor simulations
Let's first look at an example of a minimal single processor simulation code using MMSP. Suppose we want to perform a 3D simulation of isotropic grain growth using the phase field method, and we already have a valid grid data file named file1.PF. The file graingrowth.hpp, is assumed to contain the update() function we need. We create a new source file, graingrowth.cpp, and enter
// graingrowth.cpp // A minimal phase field code using MMSP #include"graingrowth.hpp" int main(int argc, char* argv[]) { PFgrid3D grid("file1.PF"); grid.update(100); grid.output("file2.PF"); }
The first line in the main() function declares the grid object of type PFgrid3D and reads the data contained in file1.PF into the empty grid data structure. The update() function performs 100 phase field time steps, and the output() function writes the new, updated grid to file2.PF.
Note that there is one line of code for each operation on the grid. This is a typical feature of the MMSP design. In any case, we now have a 3D phase field code that will read a grid from a file, perform a given number of time steps, and write the result to another file.
All that's left is to compile and run the program. These steps are, of course, platform-dependent. As an example, to compile the above program on a Linux machine with most versions of gcc, type
$ g++ -O3 -ansi graingrowth.cpp -o phase_field $ ./phase_field
The user is encouraged to use appropriate (safe) optimizations for the target machine. The user is also strongly encouraged to force the compiler to accept only ISO compliant C++ source code. While doing so should not be necessary, this prevents any possible name conflicts between MMSP and non-standard compiler extensions.
In the example given above, we declared an object of type PFgrid3D. As you might guess, the way to recast this program for a 2D grid is to change the object declaration to PFgrid2D. Remember, though, that when the program is compiled and run, the actual grid that is input from the file must be a valid 2D or 3D grid, depending on the object declaration. Such I/O errors are easy to catch, as they terminate the program and write a message to stderr indicating a grid dimension error.
We've shown one example of using the MMSP to write a simple grain growth program using the phase field method. Now suppose we want to use a different computational technique, say the sparse phase field (sparsePF) method. In this case, we only need to replace the PFgrid3D declaration with sparsePF3D, the standard 3D grid object for sparsePF simulations, and choose the appropriate header file to include. Obviously, the input file must now be a valid 3D sparsePF grid or an I/O error occurs.
2.2 Parallel processor simulations
There are only a few significant changes that must be made to the minimal single processor program described above to turn it into a parallel program, provided that the update() function has been parallelized. The minimal parallel program is
// PF.graingrowth.parallel.cpp // A minimal phase field code (parallel implementation) #include"PFgrid.graingrowth.h" int main(int argc, char* argv[]) { MPI::Init(argc,argv); int id = MPI::COMM_WORLD.Get_rank(); int np = MPI::COMM_WORLD.Get_size(); PFgrid3D grid("file1.PF",id,np); grid.update(100,id,np); grid.output("file2.PF",id,np); MPI::Finalize(); }
The four lines
MPI::Init(argc,argv); int id = MPI::COMM_WORLD.Get_rank(); int np = MPI::COMM_WORLD.Get_size(); ... MPI::Finalize();
are required of any MPI program (or at least something like them). The other changes are, of course, that the PFgrid3D constructor, output(), and update() member functions now require the processor rank and the total number of processors as additional parameters. These overloaded versions of the single processor functions are smart enough to subdivide the domain and perform the necessary interprocess communication during the simulation. Thus no external codes to decompose or recombine the grid are necessary (we can use the same grid as we did in the single processor case), and only the minimum number of MPI subroutine calls must be made by the user.
Compiling and running parallel programs depends on both the user's compiler and local MPI distribution. Most MPI distributions will use commands that looks something like
$ mpic++ -O3 -ansi graingrowth.parallel.cpp -o phase_field_parallel $ mpirun -np 4 ./phase_field_parallel
The mpirun command, of course, can only be used after MPI has been successfully booted on some number of remote hosts.
2.3 Microstructure generation
While MMSP is mainly intended to be used as a simulation package, there's nothing at all to prevent us from using the grid objects to create new microstructures. Grid objects can be declared in several ways, but perhaps the simplest is to declare the object and indicate its extents, i.e.
// Declare a 100x100 PFgrid2D object // with 10 order parameters. PFgrid2D grid1(100,100,10); ... // Declare a 100x200 CAgrid2D object. CAgrid2D grid2(100,200); ... // Declare a 100x100x100 MCgrid3D object. MCgrid3D grid3(100,100,100); // etc.
You may notice that the PFgrid2D object takes three, not two, parameters, the last telling us how large to make each node vector data structure. Now we can enter the data at each grid point based on some scheme or desired geometry, e.g.
// Declare a 100x100 MCgrid2D object ... MCgrid2D grid(100,100); // initialize the data ... for (int x=0; x<100; x++) for (int y=0; y<100; y++) grid[x][y] = rand()%1000; // and output the result to a file. grid.output("random.MC");
Optionally, we can set the attributes (boundary conditions, node spacing) of a grid object using set functions
// Change the attributes of // a MCgrid2D object "grid" // Make "grid" periodic in x but finite in y grid.set_boundary(0) = true; grid.set_boundary(1) = false; // Set y spacing to twice that in x direction grid.set_spacing(0) = 1.0; grid.set_spacing(1) = 2.0*grid.get_spacing(0);
The defaults for all grid types are finite boundary conditions and grid spacing of 1.0 in each direction. Grid attributes can be obtained from a grid by using get functions with the same syntax.
A common method for generating microstructures is to initialize a grid with random noise or some other easy to code initial state, then allow the microstructure to evolve. Since MMSP is built especially for microstructure simulation, using MMSP is quite convenient for this purpose. Future releases of MMSP (or perhaps an add-on package) may include higher level microstructure generation methods.
2.3 Microstructure analysis
At present, MMSP does not include any microstructure analysis code. However, since the grid objects in MMSP have a user interface that mimics subscripted arrays, it's often possible to adapt existing code written for analysis of "hard-coded" data structures. For example, suppose we have some code used to compute the volume (3D) or area (2D) of a given spin in an array-based Monte Carlo data structure,
// An array-based MC data structure // Note that the array size must be constant const int nx = 100; const int ny = 100; const int nz = 100; int grid[nx][ny][nz]; ... // Compute the volume (number of voxels) // associated with some given spin number int spin = 123; int volume = 0; for (int x=0; x<nx; x++) for (int y=0; y<ny; y++) for (int z=0; z<nz; z++) volume += (grid[x][y][z]==spin);
This particular piece of code, starting with "Compute the volume...", could be used verbatim to compute the same quantity in a MMSP grid object such as MCgrid3D or even CAgrid3D. In fact, if we replace the array declaration with
MCgrid3D grid(nx,ny,nz);
then we've successfully adapted the entire code to using the MCgrid3D object. Remember - any code that works for subscripted arrays also works for their MMSP equivalent. As with microstructure generation, future releases of MMSP might include more direct methods for microstructure analysis.
2.4 File formats
The file format for each high level grid type object in MMSP is largely the same. Each grid file is written in binary format with a header that is consistent for each grid type. The header contains, in order, the grid dimension (of type int), global number of fields/components (int), the grid extent (two or three ints), the grid boundary conditions (two or three bools), and the grid spacing (two or three floats). The grid data is then written for each point, looping through x, y, and, if applicable, z coordinates. The data is written in a way that depends on whether the node data structure is a scalar, vector, or sparse. For more detail, see the file MMSP.data.hpp, which contains the implementation of the node data structure buffer I/O member functions.
2.5 Visualization
Each high level grid object in MMSP has a corresponding set of utility programs to facilitate visualization. These programs use headers found in the MMSP source distribution to convert grid types to VTK file formats. Such files can be viewed with a number of freely available programs such as ParaView or VisIt. At present, only single processor visualization is supported, although MMSP will soon include code to convert grids to partitioned VTK file formats for parallel visualization.
3 Writing new MMSP components
One of the best design features of MMSP is that it's simple to add new components, whether they be new simulation methods (grid objects) or new physical processes (update() functions).
3.1 Point data structures
MMSP's existing point data structures are suitable for most applications. In the case that a new data structure is needed, the following steps should be taken. First, declare the new point data structure class as a template class within the data namespace in header MMSP.data.hpp. The class should accept a single template parameter (regardless of whether this parameter is meaningful). The class must contain a constructor that accepts an arbitrary number of arguments, i.e. the last parameter in the constructor declaration must be the ellipsis token, although only the first argument should be used. The class must also contain a function that returns the number of data fields and buffer I/O functions declared as
int fields() const; int buffer_size() const; void to_buffer() const; void from_buffer();
See the existing point data structure classes scalar, vector, and sparse for examples.
3.2 Grid objects
New grid objects can be introduced into MMSP with only a little effort (some knowledge of C++ classes and inheritance is required). For instance, consider the MCgrid2D class definition, provided in MCgrid.hpp
class MCgrid2D : public grid::grid2D<data::scalar,int> { public: // constructors MCgrid2D(int x, int y) : grid::grid2D<data::scalar,int>(x,y) {} MCgrid2D(const char* filename) : grid::grid2D<data::scalar,int>(filename) {} MCgrid2D(const char* filename, int id, int np, int ng = 1) : grid::grid2D<data::scalar,int>(filename,id,np,ng) {} // the update function void update(int steps = 1); void update(int steps, int id, int np, int ng = 1); };
MCgrid2D inherits all the member functions and member data of a grid2D (declared within the namespace "grid") with node data structure scalar and data type int. The C++ class model requires explicit declaration of constructors for all new objects, hence the member functions labeled constructors. But these are really just calls to functions already implemented by the grid2D class. Likewise, the update() functions must be declared here, since their implementation is always given by the user. We now have a working grid object MCgrid2D, having used a minimum of new source code. Sometimes it's also useful to declare utility functions, e.g. in PFgrid.hpp we've included a function to calculate the Laplacian operator at a point.
The basic procedure presented here for creating new grid classes has been used in creating all of the high level grid objects in the MMSP source distribution. Users interested in writing their own grid object(s) are encouraged to study their declarations.
3.3 The update() function
As described above, any newly created grid class must declare update() functions before it can be used in microstructure simulations. We've already seen some examples of generic update() functions in an earlier section. In large part, users can base any new update() function on the existing generic functions provided with the source distribution. Remember, when writing new update() functions, use the new class name and scope resolution operator to ensure that each is a member function of the new class.
4 MMSP reference
4.1 MCgrid2D
4.2 MCgrid3D
4.3 PFgrid2D
4.4 PFgrid3D
4.5 sparsePF2D
4.6 sparsePF3D
4.7 CAgrid2D
4.8 CAgrid3D
4.9 FDgrid2D
4.10 FDgrid3D
