User Plugins

Overview

Users can implement custom funcionality through a C++ API which provides access to internal simulation data. This is done by writing C++ code conforming to the API, compiling the code to a shared library, and providing the path to this shared library in a simulation.

Each simulation only loads one shared library. Provide the path to this shared library in the GUI under “Simulation” -> “Advanced - Plugin” -> “Plugin Path”. If multiple plugin features are used then they must be compiled together into this one shared library. Plugin features can be implemented in different C++ source files as long as they are ultimately compiled into the same shared library.

Setup

User plugins are supported on both Windows and Linux. Your system needs to be setup with a compatible compiler in order to create a shared library from the C++ code you write.

Windows

On Windows the M-Star CFD code is compiled with the Microsoft Visual C++ Compiler. While you are free to use any compiler that is binary compatible, here, we will detail the steps for installing the Microsoft Visual C++ compiler which is freely available and easy to install. This process will install command-line tools but not the Visual Studio user interface. If you prefer you could instead choose to use a Visual Studio installion with the user interface.

  1. Download the installer for the Visual Studio 2017 Build Tools. As of the time this was written, you can download it here: https://www.visualstudio.com/downloads/#build-tools-for-visual-studio-2017. Download “Build Tools for Visual Studio 2017” as shown below.

_images/visual-studio-download.PNG
  1. Run the installer and select to install “Visual C++ build tools”. And under “Optional” keep the Windows SDK option selected while deselecting the CMake option. Note that this will require a few GB of hard disk space.

_images/visual-studio-tools-install.PNG
  1. Restart you computer if a message appears prompting you to do so.

The compiler should now be installed on your system. To compile a plugin C++ file follow these steps:

  1. Launch the Visual Studio Build Tools command prompt. Start Menu -> Visual Studio 2017 -> x64 Native Tools Command Prompt for VS 2017

  2. Change directories (using the cd command) to the directory containing your code file.

  3. Compile using the following command, adjusting for the path to your M-Star CFD installation and the name of your code file:

    cl /LD /O2 -I"C:\Program Files\M-Star CFD" plugin.cpp
    

Linux

On Linux the M-Star CFD code is compiled with GCC 5.5. Again, any binary compatible compiler will work. Older linux installations based on GCC 5.0 or lower will not be compatible. A newer version of GCC can be compiled from source or you can use a package manager such as Conda which can quickly install it without root access. Once you have a compatible compiler you can compile a plugin C++ file with this command:

g++ -fPIC -shared -O3 -I/path/to/mstar/share plugin.cpp -o plugin.so

How to Write Plugins

The plugin API is defined by the MStarPluginAPI.hpp header file which is located the the main installation directory on Windows and in the share directory on Linux. This file contains C++ classes and functions provided to simplify plugin development. To create the plugin shared library include this header while providing the relevant functions from the list of available user defined functions below.

Cell Data

The M-Star CFD code operates on a mesh comprised of non-overlapping, rectangular Blocks with uniform spacing. These Blocks are distributed and communicate with neighbors using MPI. Ghost zones are used to aid in this communication process.

Each Block holds data for each Variable in the simulation, where a Variable refers to the velocity, pressure, species concentration, etc at each cell. Variable data is stored in the form of four dimensional arrays with three components for space and one for Variable component. The CellData class is a four dimensional array that gives you to access the raw Variable data on a Block. CellData is templated based on the type of data at each cell (floating point number, integer, etc) and the number of components in the Variable.

// Example for retrieving the velocity CellData
// Given: Block block
CellData<real_t,3>& velocityData = block.getCellData<real_t,3>(VarNames::Velocity);

Breaking down the code above, real_t refers the the floating point precision that the code was compiled with. The 3 indicates the number of components in the velocity vector at each cell. And the VarNames namespace provides a list of valid Variable string names.

Note

It is recommended to use the type real_t as listed in the Variables table below as opposed to float or double. This enables your code to be compatible with both single or double precision versions of the code. Keep in mind, however, that literal decimals in C++ are double so for code with large numbers of numeric literals it is recommended that you perform all your calculations with double while reading and writing data to real_t templated CellData objects.

For functions that pass a box representing the cell indices that should be updated during the function call, iterate over the three dimensions of this box like this:

// This code will loop over a given index box and update the species concentration at each cell
// Given: Block block, Box<int> box
CellData<real_t,1>& speciesData = block.getCellData<real_t,1>("SpeciesName");
for (int ix = box.xmin; ix <= box.xmax; ix++) {
    for (int iy = box.ymin; iy <= box.ymax; iy++) {
        for (int iz = box.zmin; iz <= box.zmax; iz++) {
            real_t s = speciesData(ix,iy,iz);
            real_t newS = ...; // perform some calculations here
            speciesData(ix,iy,iz) = newS;
        }
    }
}

For Variables with only one component you can use the syntax above with only the x, y, and z indices: data(ix,iy,iz). With Variables with more components you can use this syntax: data(ix,iy,iz,2), where 2 would refer to the third component.

Important

Unless otherwise stated in the function description, do not access cells outside the provided index box. These cells may be outside the block or need to be updated with ghost communication.

Remember that you are directly accessing data used by the code and, therefore, the units will not be in physical units many times. For this reason, the API provides the Units class for converting variables from and to physical units as needed.

// Convert velocity from code to real units and then back to code units
// Given: Block block, Box<int> box, Units units
CellData<real_t,3>& velocityData = block.getCellData<real_t,3>(VarNames::Velocity);
for (int ix = box.xmin; ix <= box.xmax; ix++) {
    for (int iy = box.ymin; iy <= box.ymax; iy++) {
        for (int iz = box.zmin; iz <= box.zmax; iz++) {
            real_t vy = units.velToReal( velocityData(ix,iy,iz,1) );
            real_t newVy = ...; // perform some calculations here
            velocityData(ix,iy,iz,1) = units.velToLB(newVy);
        }
    }
}

Each cell is assigned boundary conditions for each model requiring them. This means that there are boundary conditions corresponding to the fluid solver at each cell and also separate boundary conditions for scalar fields. The appropriate boundary conditions should be used based on the function being implemented. Using the boundary conditions CellData with type char, retrieve the boundary condition at the current cell. Then assess which calculations, if any, should be done at each cell accordingly. For example, many calculations are only valid at fluid cells as opposed to solid wall or void cells. This can be determined by calling the isFluid functions in the FluidBCs and ScalarBCs classes.

// Update species only at fluid cells
// Given: Block block, Box<int> box
CellData<real_t,1>& bcData = block.getCellData<real_t,1>(VarNames::ScalarBoundaryCondition);
CellData<real_t,1>& speciesData = block.getCellData<real_t,1>("SpeciesName");
for (int ix = box.xmin; ix <= box.xmax; ix++) {
    for (int iy = box.ymin; iy <= box.ymax; iy++) {
        for (int iz = box.zmin; iz <= box.zmax; iz++) {

            char bc = bcData(ix,iy,iz);
            if (!ScalarBCs::isFluid(bc)) continue;

            real_t s = speciesData(ix,iy,iz);
            real_t newS = ...; // perform some calculations here
            speciesData(ix,iy,iz) = newS;
        }
    }
}
Cell Data Variables

Variable

Data Type

# Components

Units

Velocity

real_t

3

Lattice Boltzmann

Strain Rate

real_t

1

Lattice Boltzmann

Viscosity

real_t

1

Lattice Boltzmann

Body Acceleration

real_t

3

Lattice Boltzmann

Volume Fraction

real_t

1

N/A

Temperature

real_t

1

Non-dimensionalized

Mean Age

real_t

1

seconds

Species Concentration

real_t

1

mol/L

Particle Volume Fraction

real_t

1

N/A

Particle Surface Area

real_t

1

Lattice Boltzmann

Fluid Boundary Conditions

char

1

See FluidBCs class

Scalar Boundary Conditions

char

1

See ScalarBCs class

Note that species’ names are used to access their CellData from a block in the getCellData function. And for particle volumes and surface areas the variable name should be appended by the origin index (for example, ParticleVolumeFraction0 or ParticleSurfaceArea0).

Particles

API funcionality for particles:

  • Iteration over particles

  • Iteration over a particle’s neighbors

  • Get/set particle variables

  • Get/set cell data at particle locations

  • Addition/removal of particles

  • Set custom scalar and vector variables for each particle

Since Blocks are the objects into which cell data is decomposed and distributed, particles are treated similarly. Each mesh Block has ownership over the particles inside the space occupied by the Block. In addition, each Block will have ghost particles belonging to its neighboring Blocks. All of these particles are sorted into bins which enables efficient searching for neighbor particles.

The M-Star CFD code allows particles with different properties and origins to coexist in the same simulation. To support this ability particles are organized into separate ParticleSets each acting independently and owning separate data structures. The ParticleSet class is foundational for the particle funcionality in this API so you will need to access the ParticleSets in a simulation which is done by iterating through the “particleSets” vector member of Block.

// Accessing ParticleSets
// Given: Block block
for (int i = 0; i < block.particleSets.size(); i++) {
    ParticleSet& pset = block.particleSets[i];
}

After retrieving a ParticleSet you may want to iterate over its particles. This is done by calling a ParticleSet’s “loopOverParticles” function which takes a function object as an argument. This function object is where you will place the code that is run for each particle. It’s operator() member function must have the form “void operator()(ParticleSet& pset, Particle& particle)”.

// Iterate over a ParticleSet's particles

struct ParticleIteration
{
    void operator()(ParticleSet& pset, Particle& particle)
    {
        // code per particle
    }
};

void someFunction(Block& block)
{
    for (int i = 0; i < block.particleSets.size(); i++) {
        ParticleSet& pset = block.particleSets[i];
        ParticleIteration iteration;
        pset.loopOverParticles(iteration);
    }
}

To get/set variable data for a particle, use the provided functions as shown here:

// Get/set particle data

struct ParticleIteration
{
    void operator()(ParticleSet& pset, Particle& particle)
    {
        // get particle data
        int state = pset.getState(particle);
        Array<real_t,3> pos = pset.getPosition(particle);
        Array<real_t,3> vel = pset.getVelocity(particle);
        Array<real_t,3> accel = pset.getAccel(particle);
        int id = pset.getID(particle);
        real_t timeAdded = pset.getTimeAdded(particle);
        real_t diam = pset.getDiameter(particle);
        real_t volume = pset.getVolume(particle);
        real_t customScalar = pset.getCustomScalar(particle);
        Array<real_t,3> customVector = pset.getCustomVector(particle);

        // set particle data
        pset.setVelocity(particle, vel);
        pset.setAccel(particle, accel);
        pset.setDiameter(particle, diam);
        pset.setVolume(particle, volume);
        pset.setCustomScalar(particle, customScalar);
        pset.setCustomVector(particle, customVector);
    }
};

void someFunction(Block& block)
{
    for (int i = 0; i < block.particleSets.size(); i++) {
        ParticleSet& pset = block.particleSets[i];
        ParticleIteration iteration;
        pset.loopOverParticles(iteration);
    }
}

It is also possible to get/set cell data at the location of a particle by converting the particle’s position to a 3D cell index.

// Get/set cell data at particle location

struct ParticleIteration
{
    Block& block;
    CellData<real_t,1>& varData;

    ParticleIteration(Block& block_, CellData<real_t,1>& varData_):
        block(block_), varData(varData_) {}

    void operator()(ParticleSet& pset, Particle& particle)
    {
        Array<real_t,3> pos = pset.getPosition(particle);
        Array<int,3> cellIdx = block.coords2Idx(pos);
        real_t varValue = varData(cellIdx[0],cellIdx[1],cellIdx[2]);
        varData(cellIdx[0],cellIdx[1],cellIdx[2]) = 3.7;
    }
};

void someFunction(Block& block)
{
    CellData<real_t,1>& varData = block.getCellData<real_t,1>(varName);
    for (int i = 0; i < block.particleSets.size(); i++) {
        ParticleSet& pset = block.particleSets[i];
        // need to pass block and varData to iteration function object here
        ParticleIteration iteration(block, varData);
        pset.loopOverParticles(iteration);
    }
}

Adding and removing particles is done by calling the appropriate functions on a ParticleSet (“addParticle” and “removeParticle”). These functions can be called at any time (during iteration or not). Particle removals are effective immediately meaning further iteration will ignore the particle. Particle additions are delayed and will only be performed after the user defined function is complete and returns to the main code.

// Add/remove particles
// Given: ParticleSet pset, Particle particle, new particle properties
NewParticle newParticle(position, velocity, acceleration, id, timeAdded, diameter, volume, customScalar, customVector);
pset.addParticle(newParticle);
pset.removeParticle(particle);

For models involving iteractions between particles, ParticleSets have another looping function for iterating over a particle’s neighbors “loopOverNeighbors”. Again, this function requires a function object as an argument. It also requires a bin stencil size which determines how many bins away to search relative to the current particle’s bin. Block has a convenience function “getNbrBinStencil” which will determine this argument based on a distance in meters. An example of neighbor iteration is given here:

// Iterate over all particles and each particle's neighbors

struct NeighborIteration
{
    void operator()(ParticleSet& pset, Particle& particle)
    {
        // logic per neighbor here
    }
};

struct ParticleIteration
{
    int nbrBinStencil;

    ParticleIteration(int nbrBinStencil_):
        nbrBinStencil(nbrBinStencil_) {}

    void operator()(ParticleSet& pset, Particle& particle)
    {
        NeighborIteration iteration;
        pset.loopOverNeighbors(particle, nbrBinStencil, iteration);
    }
};

void someFunction(Block& block, Units& units)
{
    int nbrBinStencil = block.getNbrBinStencil(0.005, units);
    for (int i = 0; i < block.particleSets.size(); i++) {
        ParticleSet& pset = block.particleSets[i];
        ParticleIteration iteration(nbrBinStencil);
        pset.loopOverParticles(iteration);
    }
}
Particle Data Variables

Variable

Data Type

# Components

Units

Position

real_t

3

Lattice Boltzmann (coarsest level)

Velocity

real_t

3

Lattice Boltzmann

Acceleration

real_t

3

Lattice Boltzmann

Diameter

real_t

1

meters

Volume

real_t

1

Lattice Boltzmann (coarsest level)

ID

int

1

N/A

Time Added

real_t

1

seconds

State

char

1

See ParticleState class

Custom Scalar

real_t

1

User Defined

Custom Vector

real_t

3

User Defined

User Defined Functions

Reaction Chemistry

Description: Given the species and temperature, update species concentrations to model a chemical reaction
How to Use: If the “Reaction Implementation Type” property of any one of the scalar fields is set to the plugin option then the solver will use the plugin funcionality and ignore all other reaction chemistry property inputs.
Cell Variable Access:
- All Species: read and write, stencil size 0
- Temperature: read only, stencil size 0
- Viscosity: read only, stencil size 0
- Velocity: read only, stencil size 1
- Strain Rate: read only, stencil size 0
- Volume Fraction: read only, stencil size 1
- Mean Age: read only, stencil size 0
- Particle Volume Fraction: read only, stencil size 0
- Particle Surface Area: read only, stencil size 0
- Scalar Boundary Conditions: read only, stencil size 1
- Fluid Boundary Conditions: read only, stencil size 1
Function Declaration:
DLLEXPORT void scalarReactionsUpdate(Block& block, Box<int> box, Units units, real_t simTime, real_t dt);
Arguments:
- Block& block: current Block being processed in the mesh
- Box<int> box: min and max indices for each dimension over which to iterate
- Units units: for converting between code and real units
- real_t simTime: simulation time in seconds
- real_t dt: timestep in seconds
Example:
#include "MStarPluginAPI.hpp"

using namespace MStarAPI;

DLLEXPORT void scalarReactionsUpdate(Block& block, Box<int> box, Units units, real_t simTime, real_t dt)
{
    CellData<char,1>& bcData = block.getCellData<char,1>(VarNames::ScalarBoundaryCondition);
    CellData<real_t,1>& species1Data = block.getCellData<real_t,1>("Species1Name");
    CellData<real_t,1>& species2Data = block.getCellData<real_t,1>("Species2Name");
    CellData<real_t,1>& tempData = block.getCellData<real_t,1>(VarNames::Temperature);

    for (int ix = box.xmin; ix <= box.xmax; ix++) {
        for (int iy = box.ymin; iy <= box.ymax; iy++) {
            for (int iz = box.zmin; iz <= box.zmax; iz++) {

                char bc = bcData(ix,iy,iz);
                if (!ScalarBCs::isFluid(bc)) continue;

                double s1 = species1Data(ix,iy,iz);
                double s2 = species2Data(ix,iy,iz);
                double temp = units.tempToReal(tempData(ix,iy,iz));
                double newS1 = ...; // perform some calculations here
                double newS2 = ...; // perform some calculations here
                species1Data(ix,iy,iz) = newS1;
                species2Data(ix,iy,iz) = newS2;
            }
        }
    }
}

Fluid Viscosity and Acceleration

Description: Set fluid viscosity and body acceleration at each point with custom expression
How to Use: Set “Expression Type” property under Custom Fluid to plugin option
Cell Variable Access:
- Viscosity: read and write, stencil size 0
- Body Acceleration: read and write, stencil size 0 (should increase/decrease acceleration rather than overwriting)
- All Species: read only, stencil size 0
- Temperature: read only, stencil size 0
- Velocity: read only, stencil size 1
- Strain Rate: read only, stencil size 0
- Volume Fraction: read only, stencil size 1
- Mean Age: read only, stencil size 0
- Particle Volume Fraction: read only, stencil size 0
- Particle Surface Area: read only, stencil size 0
- Scalar Boundary Conditions: read only, stencil size 1
- Fluid Boundary Conditions: read only, stencil size 1
Function Declaration:
DLLEXPORT void viscosityUpdate(Block& block, Box<int> box, Units units, real_t simTime, real_t dt);
Arguments:
- Block& block: current Block being processed in the mesh
- Box<int> box: min and max indices for each dimension over which to iterate
- Units units: for converting between code and real units
- real_t simTime: simulation time in seconds
- real_t dt: timestep in seconds
Example:
#include "MStarPluginAPI.hpp"

using namespace MStarAPI;

DLLEXPORT void viscosityUpdate(Block& block, Box<int> box, Units units, real_t simTime, real_t dt)
{
    CellData<char,1>& bcData = block.getCellData<char,1>(VarNames::FluidBoundaryCondition);
    CellData<real_t,1>& viscData = block.getCellData<real_t,1>(VarNames::Viscosity);
    CellData<real_t,1>& strainRateData = block.getCellData<real_t,1>(VarNames::StrainRate);

    for (int ix = box.xmin; ix <= box.xmax; ix++) {
        for (int iy = box.ymin; iy <= box.ymax; iy++) {
            for (int iz = box.zmin; iz <= box.zmax; iz++) {

                char bc = bcData(ix,iy,iz);
                if (!FluidBCs::isFluid(bc)) continue;

                real_t strainRate = units.strainRateToReal(strainRateData(ix,iy,iz));
                real_t viscosity = ...; // perform some calculations here
                viscData(ix,iy,iz) = units.viscosityToLB(viscosity);
            }
        }
    }
}

Particle Update

Description: Change particle variables, add/remove particles, change cell data at particles
How to Use: Turn on “Particle Plugin” property under Simulation -> Advanced - Plugin
Cell Variable Access:
- Viscosity: read and write, stencil size 0
- Body Acceleration: read and write, stencil size 0 (should increase/decrease acceleration rather than overwriting)
- All Species: read only, stencil size 1
- Temperature: read only, stencil size 1
- Velocity: read only, stencil size 1
- Strain Rate: read only, stencil size 1
- Volume Fraction: read only, stencil size 1
- Mean Age: read only, stencil size 0
- Particle Volume Fraction: read only, stencil size 0
- Particle Surface Area: read only, stencil size 0
- Scalar Boundary Conditions: read only, stencil size 1
- Fluid Boundary Conditions: read only, stencil size 1
Function Declaration:
DLLEXPORT void particleUpdate(Block& block, Units units, real_t simTime, real_t dt);
Arguments:
- Block& block: current Block being processed in the mesh
- Units units: for converting between code and real units
- real_t simTime: simulation time in seconds
- real_t dt: timestep in seconds
Example:
#include <cstdlib>
#include "MStarPluginAPI.hpp"
using namespace MStarAPI;

struct ParticleIteration
{
    Block& block;
    Units& units;
    CellData<real_t,1>& strainRateData;

    ParticleIteration(Block& block_, Units& units_, CellData<real_t,1>& strainRateData_):
        block(block_), units(units_), strainRateData(strainRateData_) {}

    void operator()(ParticleSet& pset, Particle& particle)
    {
        // get strain rate in real units at particle position
        Array<real_t,3> pos = pset.getPosition(particle);
        Array<int,3> cellIdx = block.coords2Idx(pos);
        real_t strainRate = units.strainRateToReal(strainRateData(cellIdx[0],cellIdx[1],cellIdx[2]));

        // randomly split particle into 2 particles if strain rate exceeds some paramerter
        if (strainRate > 10.0 && double(rand())/RAND_MAX < 0.01) {

            // only split particle if diameter is above minimum limit
            real_t diam = pset.getDiameter(particle);
            if (diam > 1.0e-5) {

                // new particles are near old particle but given a random kick
                Array<real_t,3> newPos1, newPos2;
                for (int i = 0; i < 3; i++) {
                    newPos1[i] = pos[i] + double(rand())/RAND_MAX * 0.3;
                    newPos2[i] = pos[i] + double(rand())/RAND_MAX * 0.3;
                }

                NewParticle newParticle;
                newParticle.vel = pset.getVelocity(particle);
                newParticle.accel = pset.getAccel(particle);
                newParticle.id = pset.getID(particle);
                newParticle.timeAdded = pset.getTimeAdded(particle);
                newParticle.diameter = diam * 0.7937; // diameter corresponding to halving volume
                newParticle.customScalar = pset.getCustomScalar(particle);
                newParticle.customVector = pset.getCustomVector(particle);

                // add two new smaller particles and remove existing larger particle
                newParticle.pos = newPos1;
                pset.addParticle(newParticle);
                newParticle.pos = newPos2;
                pset.addParticle(newParticle);
                pset.removeParticle(particle);
            }
        }
    }
};

DLLEXPORT void particleUpdate(Block& block, Units units, real_t simTime, real_t dt)
{
    CellData<real_t,1>& strainRateData = block.getCellData<real_t,1>(VarNames::StrainRate);
    for (int i = 0; i < block.particleSets.size(); i++) {
        ParticleSet& pset = block.particleSets[i];
        ParticleIteration iteration(block, units, strainRateData);
        pset.loopOverParticles(iteration);
    }
}

Additional Functionality

Array

This class is a high performance one dimensional array to simplify handling on Variables with multiple components. It is templated with the type of each item in the array and the array’s length. Note that there is also some builtin functionality for addition, subtraction, cross product, dot product, etc of these arrays.

Array<real_t,3> a(0,0,0);
Array<real_t,3> b = a + Array<real_t,3>(0, 1.3, 6.4);
real_t c = dot(a, b);
Array<real_t,3> v = b / norm(b);
Array<real_t,3> cp = crossProduct(a, b);

Units

The Units class provides functions for converting between physical units and the units used internally in the code.

real_t strainRateReal = units.strainRateToReal(srLB);
real_t strainRateLB = units.strainRateToLB(strainRateReal);

Gradient Calculation

The calcGradient functions perform simple finite difference calculations of CellData at the given x,y,z cell indices. Two functions are available - one for scalar Variables and one for three-component Variables.

DLLEXPORT void someFunction(Block& block, Box<int> box, Units units, real_t simTime, real_t dt)
{
    CellData<char,1>& bcData = block.getCellData<char,1>(VarNames::FluidBoundaryCondition);
    CellData<real_t,3>& velData = block.getCellData<real_t,3>(VarNames::Velocity);

    for (int ix = box.xmin; ix <= box.xmax; ix++) {
        for (int iy = box.ymin; iy <= box.ymax; iy++) {
            for (int iz = box.zmin; iz <= box.zmax; iz++) {

                char bc = bcData(ix,iy,iz);
                if (!FluidBCs::isFluid(bc)) continue;

                Matrix3x3 velGrad = calcGradient(velData,ix,iy,iz);
                // access components with velGrad[0][2] which is dvx/dz for example
            }
        }
    }
}