OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
No Matches
Vector 2 expression

Vector expressions

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


He define some useful constants

constexpr int A = 0;
constexpr int B = 1;
constexpr int C = 2;
constexpr int D = 3;
Derivative second order on h (spacing)

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)

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)
// 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
// 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);
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Distributed vector.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
__device__ __host__ float value(const Point< 3, float > &p, const Point< 3, float > &q, float pA, float pB)
Result of the exponential kernel.


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

See also
// initialize the library
// 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
// extended boundary around the domain, and the processor domain
Ghost<3,double> g(0.01);
// Delta t
double dt = 0.01;
This class represent an N-dimensional box.
Definition Box.hpp:61

%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> > 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;

Mapping particles

Redistribute the particles according to the decomposition

See also
Mapping particles;;

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);


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

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

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

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

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 \)

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.

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


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