Vector 2 expression

Vector expressions

This example show how to use expressions to do complex calculation on particles

Initialization

He define some useful constants

 
constexpr int A = 0;
constexpr int B = 1;
constexpr int C = 2;
constexpr int D = 3;
 

Kernel function

We define a kernel function. In general a kernel is a function like \( K(x_p,x_q,A_p,A_q) \). Where \( x_p \) and \( x_q \) are positions of the 2 particles p and q. \( A_p \) and \( A_q \) are the properties on p and q.

A second version more General is

\( K(p,q,A_p,A_q,particles) \)

Where \( p \) in the index of the particle and \( q \) are the index of the particles produced by the cell-List. and particles is the vector containing the particle p and q (In multi-phase SPH this can also be not true)

Note

On multi-phase this is not true

In this case we are defining an exponential-like function \( A_q e^{\frac{|x_p-x_q|^2}{\sigma}} \)

 
// Kernel are structures
struct exp_kernel
{
    // variance of the exponential
    float sigma;
 
    // Constructor require to define the variance
    exp_kernel(float sigma)
    :sigma(sigma)
    {}
 
    // The kernel function itself K(x_p,x_q,A_p,A_q). The name MUST be value and require 4 arguments.
    // p,q,A_p,A_q Position of p, position of q, value of the property A on p, value of the property A on q
    // and it return the same type of the property A
    inline Point<3,double> value(const Point<3,float> & p, const Point<3,float> & q,Point<3,double> & Pp,Point<3,double> & Pq)
    {
        // calculate the distance between p and q
        float dist = norm(p-q);
 
        // Calculate the exponential
        return Pq * exp(- dist * dist / sigma);
    }
 
    // The kernel function itself K(p,q,A_p,A_q,particles). The name MUST be value and require 5 arguments.
    // p,q,A_p,A_q,particles size_t index of p, size_t index of q, value of the property A on p, value of the property A on q,
    // original vector of the particles p, and it return the same type of the property A
    inline Point<3,double> value(size_t p, size_t q, Point<3,double> & Pp, Point<3,double> & Pq, const vector_dist<3,double, aggregate<double,double,Point<3,double>,Point<3,double>> > & v)
    {
        // calculate the distance between p and q
        float dist = norm(p-q);
 
        // Calculate the exponential
        return Pq * exp(- dist * dist / sigma);
    }
};
 

Initialization

Here we Initialize the library, we create a Box that define our domain, boundary conditions, and ghost

See also

Initialization

 
    // initialize the library
    openfpm_init(&argc,&argv);
 
    // Here we define our domain a 2D box with intervals from 0 to 1.0
    Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
 
    // Here we define the boundary conditions of our problem
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    // extended boundary around the domain, and the processor domain
    Ghost<3,double> g(0.01);
 
    // Delta t
    double dt = 0.01;
 

%Vector create

We create 2 vectors one with 4 properties 2 scalar and 2 vectors. and one with one scalar

See also

Vector instantiation

 
    // distributed vectors
    vector_dist<3,double, aggregate<double,double,Point<3,double>,Point<3,double>> > vd(4096,domain,bc,g);
    vector_dist<3,double, aggregate<double> > vd2(4096,domain,bc,g);
 

Loop over particles

We assign random position to particles for both vectors

See also

Assign position

 
    // Assign random position to the vector dist
    auto it = vd.getDomainIterator();
 
    while (it.isNext())
    {
        auto p = it.get();
 
        vd.getPos(p)[0]  = (double)rand() / RAND_MAX;
        vd.getPos(p)[1]  = (double)rand() / RAND_MAX;
        vd.getPos(p)[2]  = (double)rand() / RAND_MAX;
 
        vd2.getPos(p)[0] = (double)rand() / RAND_MAX;
        vd2.getPos(p)[1]  = (double)rand() / RAND_MAX;
        vd2.getPos(p)[2]  = (double)rand() / RAND_MAX;
 
        ++it;
    }
 

Mapping particles

Redistribute the particles according to the decomposition

See also

Mapping particles

 
    vd.map();
    vd.map();
 

Getting expression properties

Before using expression we have to get alias for each properties we want to use for the expression

 
    // vA is an alias for the property A of the vector vd
    // vB is an alias for the property V of the vector vd
    // vC is ...
    auto vA = getV<A>(vd);
    auto vB = getV<B>(vd);
    auto vC = getV<C>(vd);
    auto vD = getV<D>(vd);
 
    // This indicate the particle position
    auto vPOS = getV<PROP_POS>(vd);
 
    // same concept for the vector v2
    auto v2A = getV<0>(vd);
    auto v2POS = getV<PROP_POS>(vd);
 

Expressions

Once we have the expression properties we can use them to compose expressions

 
    // Assign 1 to the property A
    vA = 1;
 
    // Equal the property A to B
    vB = vA;
 
    // Assign to the property C and for each component 4.0
    vC = 4;
 
    // Assign the position of the particle to the property D
    vD = vPOS;
 

All these expression are applied point wise for each particles

Note

the * is a scalar product in case of vectors simple multiplication in case of scalar

 
    // In vA we set the distance from the origin for each particle
    vA = norm(vPOS);
 
    //
    // For each particle p calculate the expression under
    //
    // NOTE sin(2.0 * vD)  and exp(5.0 * vD) are both vector
    //
    //              so this is a scalar product
    //                      |
    //                      V
    vA = vA + sin(2.0 * vD) * exp(1.2 * vD);
    //       |_____________________________|
    //                      |
    //                      V
    //              and this is a number
    //
    //
 
 
    // pmul indicate component-wise multiplication the same as .* in Matlab
    vC = pmul(vD,vD) * dt;
 
    // Normalization of the vector vD
    vD = vD / sqrt( vD * vD );
 

Visualization, write VTK files

With this function we output the particles in VTK format.

See also

Visualization, write VTK files

 
    // we write the result on file
    vd.write("output");
 

Lazy expression and convert object into expressions

If we want to use an object like Point into an expression we have to first convert into an expression with the function getVExpr. Every expression does not produce any code execution until we do not assin this expression to some property

 
    // Constant point
    Point<3,double> p0({0.5,0.5,0.5});
 
    // we cannot use p0 directly we have to create an expression from it
    auto p0_e = getVExpr(p0);
 
    // here we are creating an expression, expr1 collect the full expression but does not produce any
    // code execution
    auto expr1 = 1.0/2.0/sqrt(M_PI)*exp(-(v2POS-p0_e)*(v2POS-p0_e)/2.0);
 
    // here the expression is executed on all particles of v2 and the result assigned to the property A of v2
    v2A = expr1;
 

Each operator= produce an iteration across particles. It is possible to use the function assign to execute multiple assignment in one cycle. Assign support a variable number of expression

 
    assign(vA,1.0,               // vA = 1
           vB,expr1,          // vB = 1.0/2.0/sqrt(M_PI)*exp(-(v2POS-p0_e)*(v2POS-p0_e)/2.0)
           vC,vD/norm(vD),    // vC = vD/norm(vD)
           vD,2.0*vC);        // vD = 2.0*vC // here vC = vD/norm(vD)
 

Apply kernel

we shown the simple usage of simple point-wise function like sin and exp. Here we show the function applyKernel that encapsulate the operation

\(\sum_{q = Neighborhood(p)} Kernel(x_p,x_q,A_p,A_q) \)

where \( x_p \) is the position of the particle p, \( x_q \) is the position of the particle q, \( A_p \) is the value of the property carried by the particle p, \( A_q \) is the value of the property carried by the particle q

Note

On particle method we saw the following formulas

\( \sum_{q = Neighborhood(p)} A_q [D^{\beta}ker](x_p,x_q) V_q \) in SPH

\( \sum_{q = Neighborhood(p)} (A_p-A_q)ker(x_p,x_q) V_q \) in PSE

These 2 formulas and others can be embedded into the previous one where

\( Kernel(x_p,x_q,A_p,A_q) = A_q [D^{\beta}ker](x_p,x_q) V_q \)

\( Kernel(x_p,x_q,A_p,A_q) = (A_p-A_q)*ker(x_p,x_q) V_q \)

We create some useful object to use apply kernel plus we synchronize the ghost for the property 3 because apply kernel require to read the ghost part

 
    // Create a kernel with sigma 0.5
    exp_kernel ker(0.5);
 
    // Get a Cell list with 0.1 of cutoff-radius
    auto cl = vd.getCellList(0.1);
 
    // We are going to do some calculation that require ghost, so sync it 3 == vD
    vd.ghost_get<3>();
 

Than we construct an expression with applyKernel

 
    vC = applyKernel_in( vD / norm(vD) ,vd,cl,ker) + vD;
 

A renormalizaion of the vector vD is calculated on fly for each particle q neighborhood of p. The the following expression is calculated for each particle p

\( (\sum_{q = Neighborhood p} Ker(x_p,x_q,(\frac{vD}{norm(vD)})_p,(\frac{vD}{norm(vD)})_q)) + vD \)

Note

Cerefull here ApplyKernel is a special function if it act on property vD we cannot write on property vD ( vD = ... is an error and produce unpredictable result )

Exist also a second version suitable for more General cases

 
    // Second version in this case the kernel is called with the particles id for p
    // and
    vC = applyKernel_in_gen( vD / norm(vD) ,vd,cl,ker) + vD;
 

In the case of General-Finite-Differences, DC-PSE and Multi-Phase SPH. \( K(x_p,x_q,A_p,A_q)\) is limitative and we have to use the following \( K(p,q,A_p,A_q,particles) \). In this case the position \( x_p \) and \( x_q \) are substituted by the id of p and id of q comming from the cell list + the vector on which we are operating.

Note

in case of Multi-Phase multiple vectors are needed. But such vector can passed in the Kernel structure

 
    // Second version in this case the kernel is called with the particles id for p
    // and
    vC = applyKernel_in_gen( vD / norm(vD) ,vd,cl,ker) + vD;
 

If the expression is point-wise the rule the propety on the right cannot appear on the left does not apply

VTK and visualization

Write a VTK file and visualize. Before write the particles also delete Ghost (VTK writer in general output also the ghost particles)

See also

Visualization, write VTK files

 
    openfpm_finalize();
 

Finalize

At the very end of the program we have always de-initialize the library

 
    openfpm_finalize();