From 20b2b64d19b44b01815e7755a8f29578837086f8 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <i-bird@localhost.localdomain>
Date: Wed, 20 Jul 2016 18:59:54 +0200
Subject: [PATCH] Small changes + example refactor
---
example/Numerics/PSE/1_Diffusion_1D/Makefile | 2 +-
example/Numerics/PSE/1_Diffusion_1D/main.cpp | 487 +++++++++++-------
example/Vector/2_molecular_dynamic/Makefile | 21 -
example/Vector/2_molecular_dynamic/config.cfg | 2 -
example/Vector/2_molecular_dynamic/main.cpp | 379 --------------
openfpm_numerics | 2 +-
src/Grid/grid_dist_id.hpp | 4 +-
src/Vector/vector_dist.hpp | 17 +
8 files changed, 334 insertions(+), 580 deletions(-)
delete mode 100644 example/Vector/2_molecular_dynamic/Makefile
delete mode 100644 example/Vector/2_molecular_dynamic/config.cfg
delete mode 100644 example/Vector/2_molecular_dynamic/main.cpp
diff --git a/example/Numerics/PSE/1_Diffusion_1D/Makefile b/example/Numerics/PSE/1_Diffusion_1D/Makefile
index a4360521..4cf11983 100644
--- a/example/Numerics/PSE/1_Diffusion_1D/Makefile
+++ b/example/Numerics/PSE/1_Diffusion_1D/Makefile
@@ -7,7 +7,7 @@ LDIR =
OBJ = main.o
%.o: %.cpp
- $(CC) -O0 -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
+ $(CC) -O3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
diff_1d: $(OBJ)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
diff --git a/example/Numerics/PSE/1_Diffusion_1D/main.cpp b/example/Numerics/PSE/1_Diffusion_1D/main.cpp
index 0fb0d3db..93299113 100644
--- a/example/Numerics/PSE/1_Diffusion_1D/main.cpp
+++ b/example/Numerics/PSE/1_Diffusion_1D/main.cpp
@@ -3,7 +3,11 @@
*
* ## Simple example
*
- * In this example we show 1D PSE derivative function approximation
+ * In this example we show The 1D PSE Diffusion equation
+ *
+ * in particular we integrate the following equation
+ *
+ * \$ \frac{u}{t} = D \laplacian u \$
*
* ### WIKI END ###
*
@@ -21,7 +25,7 @@
/*
* ### WIKI 2 ###
*
- * Here we define the function x*e^(-x*x) and its
+ * Here we define the function x*e^(-x*x) the initial condition
*
*/
@@ -35,6 +39,7 @@ double f_xex2(Point<1,double> & x)
return x.get(0)*exp(-x.get(0)*x.get(0));
}
+
/*
*
* ### WIKI 3 ###
@@ -52,15 +57,15 @@ double f_xex2(Point<1,double> & x)
* \param cl CellList
*
*/
-template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key, vector_dist<1,double, aggregate<double>, CartDecomposition<1,double> > & vd, double eps, double spacing, CellL & cl)
+template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key, vector_dist<1,double, aggregate<double,double> > & vd, double eps, double spacing, CellL & cl)
{
- // Laplacian PSE kernel 1 dimension, on double, second order
+ // Laplacian PSE kernel<1,double,2> = 1 dimension, on double, second order
Lap_PSE<1,double,2> lker(eps);
// double PSE integration accumulator
double pse = 0;
- // Get f(x) at the position of the particle
+ // f(x_p)
double prp_x = vd.template getProp<0>(key);
// Get the neighborhood of the particles
@@ -69,17 +74,16 @@ template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key,
{
auto nnp = NN.get();
- // Calculate contribution given by the kernel value at position p,
- // given by the Near particle
+ // Calculate contribution given by the particle q near to p
if (nnp != key.getKey())
{
- // W(x-y)
+ // W(x_p-x_q) x = position of p, y = position of q
double ker = lker.value(p,vd.getPos(nnp));
- // f(y)
+ // f(x_q)
double prp_y = vd.template getProp<0>(nnp);
- // 1.0/(eps)^2 [f(y)-f(x)]*W(x,y)*V_q
+ // 1.0/(eps)^2 [f(x_q)-f(x_p)]*W(x_p-x_q)*V_q
double prp = 1.0/eps/eps * (prp_y - prp_x) * spacing;
pse += prp * ker;
}
@@ -88,9 +92,72 @@ template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key,
++NN;
}
+ // return the calculated laplacian
return pse;
}
+/*
+ *
+ * ### WIKI 4 ###
+ *
+ * It mirror the particles at the border of the domain, in particular on the left we do the operation
+ * show in figure
+ *
+ * \verbatim
+ *
+ * -0.6 -0.5 0.5 0.6 Strength
+ +----+----*----*----*-
+ 0.0 Position
+
+ with * = real particle
+ + = mirror particle
+
+ \endverbatim
+ *
+ * We create particle on the negative part of the domain and we assign a strength opposite to the particle in the
+ * positive part
+ *
+ * On the right instead we prolong the initial condition (Ideally should remain almost fixed)
+ *
+ * \param vd vector of particles
+ * \param key particle
+ * \param m_pad border box(interval) on the left
+ * \param m_pad2 border box(interval) on the right
+ * \param box domain
+ *
+ */
+inline void mirror(vector_dist<1,double, aggregate<double,double> > & vd, vect_dist_key_dx & key, const Box<1,double> & m_pad, Box<1,double> & m_pad2, const Box<1,double> & box)
+{
+ // If the particle is inside the box (interval), it mean that is at the left border
+ if (m_pad.isInsideNB(vd.getPos(key)) == true)
+ {
+ // Add a new particle
+ vd.add();
+
+ //Set the added particle position, symmetrically positioned on the negative part of the domain
+ vd.getLastPos()[0] = - vd.getPos(key)[0];
+
+ // Set the property of the particle to be negative
+ vd.template getLastProp<0>() = - vd.template getProp<0>(key);
+ }
+
+ // If the particle is inside the box (interval), it mean that is at the right border
+ if (m_pad2.isInsideNB(vd.getPos(key)) == true)
+ {
+ // Add a new particle
+ vd.add();
+
+ //Set the added particle position, symmetrically positioned around x = 4.0
+ // 4.0 Added
+ // *------|------*
+ //
+ vd.getLastPos()[0] = 2.0 * box.getHigh(0) - vd.getPos(key)[0];
+
+ // Prolongate the initial condition
+ vd.template getLastProp<0>() = f_xex2(vd.getLastPos()[0]);
+ }
+}
+
/*
*
* ### WIKI END ###
@@ -100,28 +167,24 @@ template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key,
int main(int argc, char* argv[])
{
//
- // ### WIKI 3 ###
- //
- // Some useful parameters. Like
+ // ### WIKI 5 ###
//
- // * Number of particles
- // * Minimum number of padding particles
- // * The computational domain
- // * The spacing
- // * The mollification length
- // * Second order Laplacian kernel in 1D
+ // Some useful parameters.
//
// Number of particles
const size_t Npart = 1000;
+
+ // A different way to write Npart
+ size_t NpartA[1] = {Npart+1};
- // Number of step
+ // Number of steps
const size_t Nstep = 1000;
// Delta t
const double dt = 10.0 / Nstep;
- // Number of padding particles (At least)
+ // Number of padding particles, particles outside the domain
const size_t Npad = 20;
// The domain
@@ -133,8 +196,15 @@ int main(int argc, char* argv[])
// The mollification length
const double eps = 2*spacing;
+ // Diffusion constant
+ double D = 1e-4;
+
+ // 2 constants
+ constexpr int U = 0;
+ constexpr int Lap_U = 1;
+
//
- // ### WIKI 2 ###
+ // ### WIKI 6 ###
//
// Here we Initialize the library and we define Ghost size
// and non-periodic boundary conditions
@@ -142,62 +212,64 @@ int main(int argc, char* argv[])
openfpm_init(&argc,&argv);
Vcluster & v_cl = create_vcluster();
- size_t bc[1]={NON_PERIODIC};
+ size_t bc[1]={NON_PERIODIC};
Ghost<1,double> g(12*eps);
//
- // ### WIKI 3 ###
- //
- // Here we are creating a distributed vector defined by the following parameters
+ // ### WIKI 7 ###
//
- // we create a set of N+1 particles to have a fully covered domain of particles between 0.0 and 4.0
- // Suppose we have a spacing given by 1.0 you need 4 +1 particles to cover your domain
+ // Here we are creating a distributed vector defined by the following parameters
+ //
+ // <1,double, aggregate<double,double> > = 1 dimension, space with type double, 2 properties double,double (aggregate is like tuples, a way to encapsulate types). The first property is the field value u, the second is the Laplacian of u
+ // with Npart+1 particles, on the box domain [0.0,4.0], with bc boundary condition NON_PERIODIC, and an extended ghost part of size g
//
- vector_dist<1,double, aggregate<double>, CartDecomposition<1,double> > vd(Npart+1,box,bc,g);
+ vector_dist<1,double, aggregate<double,double> > vd(Npart+1,box,bc,g);
//
- // ### WIKI 4 ###
- //
- // We assign the position to the particles, the scalar property is set to
- // the function x*e^(-x*x) value.
- // Each processor has parts of the particles and fill part of the space, the
- // position is assigned independently from the decomposition.
- //
- // In this case, if there are 1001 particles and 3 processors the in the
- // domain from 0.0 to 4.0
+ // ### WIKI 8 ###
//
- // * processor 0 place particles from 0.0 to 1.332 (334 particles)
- // * processor 1 place particles from 1.336 to 2.668 (334 particles)
- // * processor 2 place particles from 2.672 to 4.0 (333 particles)
+ // We assign the position to the particles in order to be equally spaced in a grid like position, the scalar property is set to
+ // the function x*e^(-x*x) value the initial condition
+ //
+ // P.S. In a Multi processor context getGridIterator adapt based on the underline space division across processors
//
+
+ // Create a grid like iterator
+ auto itp = vd.getGridIterator(NpartA);
- // It return how many particles the processors with id < rank has in total
- size_t base = vd.init_size_accum(Npart+1);
- auto it2 = vd.getIterator();
-
- while (it2.isNext())
+ // Counter starting from 0
+ size_t p = 0;
+
+ // For each point in the grid
+ while (itp.isNext())
{
- auto key = it2.get();
-
- // set the position of the particles
- vd.getPos(key)[0] = (key.getKey() + base) * spacing;
- //set the property of the particles
- vd.template getProp<0>(key) = f_xex2((key.getKey() + base) * spacing);
-
- ++it2;
+ // Get the grid index
+ auto key = itp.get();
+
+ // Assign the particle position
+ vd.getPos(p)[0] = key.get(0) * spacing;
+
+ // Assign the property 0 based on the initial condition
+ vd.template getProp<U>(p) = f_xex2(key.get(0) * spacing);
+
+ // Next particle
+ ++p;
+
+ // Next grid point
+ ++itp;
}
//
- // ### WIKI 5 ###
+ // ### WIKI 9 ###
//
// Once defined the position, we distribute them across the processors
- // following the decomposition, finally we get the ghost part
+ // following the decomposition, finally we update the ghost part for the vield U
//
vd.map();
- vd.ghost_get<0>();
+ vd.ghost_get<U>();
//
- // ### WIKI 6 ###
+ // ### WIKI 10 ###
//
// near the boundary we have two options, or we use one sided kernels,
// or we add additional particles, it is required that such particles
@@ -219,11 +291,19 @@ int main(int argc, char* argv[])
// \endverbatim
//
//
+
+ // The particle that we have to mirror are ...
+
+ // ... on the box m_pad for the laft part
Box<1,double> m_pad({0.0},{0.1});
+
+ // ... on the box m_pad2 for the right part
Box<1,double> m_pad2({3.9},{4.0});
- double enlarge = 0.1;
+
+ // It contain how much we enlarge the domain
+ double enlarge = m_pad.getHigh(0) - m_pad.getLow(0);
- // Create a box
+ // In general 0.1 should be a sufficent padding, but in case it is not we recalculate the left and right boxes
if (Npad * spacing > 0.1)
{
m_pad.setHigh(0,Npad * spacing);
@@ -231,144 +311,205 @@ int main(int argc, char* argv[])
enlarge = Npad * spacing;
}
+
+ //
+ // ### WIKI 11 ###
+ //
+ // Here we mirror the particles if needed
+ //
+ //
+
auto it = vd.getDomainIterator();
+ size_t n_part = vd.size_local();
while (it.isNext())
{
auto key = it.get();
- // set the position of the particles
- if (m_pad.isInsideNB(vd.getPos(key)) == true)
- {
- vd.add();
- vd.getLastPos()[0] = - vd.getPos(key)[0];
- vd.template getLastProp<0>() = - vd.template getProp<0>(key);
- }
-
- // set the position of the particles
- if (m_pad2.isInsideNB(vd.getPos(key)) == true)
- {
- vd.add();
- vd.getLastPos()[0] = 2.0 * box.getHigh(0) - vd.getPos(key)[0];
- vd.template getLastProp<0>() = f_xex2(vd.getLastPos()[0]);
- }
+ mirror(vd,key,m_pad,m_pad2,box);
++it;
}
//
- // ### WIKI 6 ###
- //
- // We create a CellList with cell spacing 12 sigma
+ // ### WIKI 12 ###
//
-
- // get and construct the Cell list
+ // We create a CellList with cell spacing 8 epsilon, on the part of space domain + ghost + padding
+ // The padding area, can be bigger than the ghost part. In practice there is no reason to do it, but
+ // keeping ghost size and padding area unrelated give us the possibility to show how to create a CellList
+ // on a area bigger than the domain + ghost
+
Ghost<1,double> gp(enlarge);
- auto cl = vd.getCellList(8*eps,gp);
+
+ // Create a Cell list with Cell spaping 8*epsilon, the CellList is created on a domain+ghost+enlarge space
+ auto cl = vd.getCellList(8*eps,gp);
- // Maximum infinity norm
- double linf = 0.0;
+ // Maximum infinity norm
+ double linf = 0.0;
//
- // ### WIKI 7 ###
+ // ### WIKI 13 ###
//
- // Euler time integration
- //
- //
+ // Euler time integration until t = 10.0
+ //
+ //
- for (double t = 0; t <= 10 ; t += dt)
- {
- auto it_p = vd.getDomainIterator();
- while (it_p.isNext())
- {
- // key
- vect_dist_key_dx key = it_p.get();
+ for (double t = 0; t <= 10.0 ; t += dt)
+ {
+
+ // Iterate all the particles, including padding because are real particles
+ auto it_p = vd.getDomainIterator();
+ while (it_p.isNext())
+ {
+ // Get particle p
+ auto p = it_p.get();
- Point<1,double> p = vd.getPos(key);
+ // Get the position of the particle p
+ Point<1,double> x_p = vd.getPos(p);
- // We are not interested in calculating out the domain
- // note added padding particle are considered domain particles
- if (p.get(0) < 0.0 || p.get(0) >= 4.0)
+ // We are not interested in calculating out the domain (So on the padded particles)
+ if (x_p.get(0) < 0.0 || x_p.get(0) >= 4.0)
{
++it_p;
continue;
}
- double pse = calcLap(p,key,vd,eps,spacing,cl);
-
- // Euler update step
- vd.template getProp<0>(key) = vd.template getProp<0>(key) + pse * dt;
+ // Here we calculate the Laplacian of u on position x_p and we store the result on the particle p
+ vd.template getProp<Lap_U>(p) = calcLap(x_p,p,vd,eps,spacing,cl);
+ // Next particle
++it_p;
- }
- }
-
- // Once we have the value of the Laplacian we update with the eulerian step
-
-
- //
- // ### WIKI 8 ###
- //
- // Collect the solution in one processor
- //
-
- // Resize a vector with the local size of the vector
- openfpm::vector<double> v;
- v.resize(vd.size_local());
-
- // y contain 2 vectors the first is the calculated solution
- // the second contain the analytical solution
- openfpm::vector<openfpm::vector<double>> y;
-
- // Copy the property we want to Plot
- for (size_t i = 0 ; i < vd.size_local() ; i++)
- y.get(i) = vd.template getProp<0>(i);
-
- // Gather the solution on master
- v_cl.SGather(v,y.get(0),0);
-
- //
- // ### WIKI 9 ###
- //
- // Here we plot
- //
-
- // If there are too many points terminate
- if (y.get(0).size() > 16000)
- {
- std::cerr << "Too many points for plotting " << std::endl;
- return 1;
- }
-
- openfpm::vector<double> x;
-
- // It fill the vector x with vg.size() points with interval values between 0.0 and 4.0
- // and fill the vector y;
-
- // Total points to plot
- size_t tot_p = y.get(0).size();
-
- Fill1D(0.0,4.0,tot_p,x);
-
- // It fill y.get(1) with the Analytical solution
- Fill1D(0.0,4.0,tot_p,y.get(1),f_xex2);
-
- // Plot with GooglePlot
- // Google charts options
- GCoptions options;
-
- options.yAxis = std::string("Y Axis");
- options.xAxis = std::string("X Axis");
- options.lineWidth = 1.0;
-
- GoogleChart cg;
- cg.AddLinesGraph(x,y,options);
- cg.write("PSE_plot.html");
-
- //
- // ### WIKI 8 ###
- //
- // Deinitialize the library
- //
- openfpm_finalize();
+ }
+
+ // Once we calculated the Laplacian we do not need enymore the padding particle eliminate them
+ // n_part is the original size of the vector without padding particles
+ vd.resize(n_part);
+
+ // Iterate again on all particles
+ auto it_p2 = vd.getDomainIterator();
+ while (it_p2.isNext())
+ {
+ // Get particle p
+ auto p = it_p2.get();
+
+ // Get the Laplacian
+ double pse = vd.template getProp<Lap_U>(p);
+
+ // Integrate with Euler step
+ vd.template getProp<0>(p) += D * pse * dt;
+
+ // mirror again the particles if needed (this for the next step)
+ mirror(vd,p,m_pad,m_pad2,box);
+
+ // Next particle
+ ++it_p2;
+ }
+
+ // Update the ghost part
+ vd.ghost_get<U>();
+ }
+
+ //
+ // ### WIKI 14 ###
+ //
+ // U now contain the solution scattered across processors. Because it is not possible plot on multiple
+ // processors (Google chart does not support it) we collect the solution on one processor
+ //
+
+ struct pos_val
+ {
+ // position of the particle
+ double pos;
+
+ // value of U of the particle
+ double val;
+
+ // usefull to sort the particle by position
+ bool operator<(const pos_val & p)
+ {
+ return pos < p.pos;
+ }
+ };
+
+ // Resize a vector with the local size of the vector
+ openfpm::vector<pos_val> v(vd.size_local());
+
+ // This will contail all the particles on the master processor
+ openfpm::vector<pos_val> v_tot;
+
+
+ // Copy the particle position and the field U, into the vector v
+ for (size_t i = 0 ; i < vd.size_local() ; i++)
+ {
+ // particle position
+ v.get(i).pos = vd.getPos(i)[0];
+
+ // Field U
+ v.get(i).val = vd.template getProp<U>(i);
+ }
+
+ // from the partial solution vector v we Gather the total solution v_tot on master (processor=0)
+ v_cl.SGather(v,v_tot,0);
+
+ //
+ // ### WIKI 15 ###
+ //
+ // Once we gather the solution on master we plot it. Only processor 0 do this
+ //
+
+ // Only processor = 0
+ if (v_cl.getProcessUnitID() == 0)
+ {
+ // sort the solution by particle position
+ v_tot.sort();
+
+ // y contain 2 vectors the first is the calculated solution
+ // the second contain the initial condition
+ openfpm::vector<openfpm::vector<double>> y(2);
+
+ // vector with the particle position for plotting
+ openfpm::vector<double> x;
+
+ // We fill the plotting variables
+ for (size_t i = 0 ; i < v_tot.size() ; i++)
+ {
+ // fill x with the particle position
+ x.add(v_tot.get(i).pos);
+
+ // fill y with the value of the field U contained by the particles
+ y.get(0).add(v_tot.get(i).val);
+
+ // fill y with the initial condition
+ y.get(1).add(f_xex2(v_tot.get(i).pos));
+ }
+
+ // Plot with GooglePlot
+ // Google charts options
+ GCoptions options;
+
+ // Assign a name for the Y axis
+ options.yAxis = std::string("U Axis");
+
+ // Assign a name for the X Axis
+ options.xAxis = std::string("X Axis");
+
+ // width of the line
+ options.lineWidth = 1.0;
+
+ // Create a google plot object
+ GoogleChart cg;
+
+ // add a line plot
+ cg.AddLines(x,y,options);
+
+ // write the plot file
+ cg.write("PSE_plot.html");
+ }
+
+ //
+ // ### WIKI 16 ###
+ //
+ // Deinitialize the library
+ //
+ openfpm_finalize();
}
diff --git a/example/Vector/2_molecular_dynamic/Makefile b/example/Vector/2_molecular_dynamic/Makefile
deleted file mode 100644
index c4a841ff..00000000
--- a/example/Vector/2_molecular_dynamic/Makefile
+++ /dev/null
@@ -1,21 +0,0 @@
-include ../../example.mk
-
-CC=mpic++
-
-LDIR =
-
-OBJ = main.o
-
-%.o: %.cpp
- $(CC) -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
-
-md_dyn: $(OBJ)
- $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
-
-all: md_dyn
-
-.PHONY: clean all
-
-clean:
- rm -f *.o *~ core md_dyn
-
diff --git a/example/Vector/2_molecular_dynamic/config.cfg b/example/Vector/2_molecular_dynamic/config.cfg
deleted file mode 100644
index 1eecbac3..00000000
--- a/example/Vector/2_molecular_dynamic/config.cfg
+++ /dev/null
@@ -1,2 +0,0 @@
-[pack]
-files = main.cpp Makefile
diff --git a/example/Vector/2_molecular_dynamic/main.cpp b/example/Vector/2_molecular_dynamic/main.cpp
deleted file mode 100644
index b8ae1efe..00000000
--- a/example/Vector/2_molecular_dynamic/main.cpp
+++ /dev/null
@@ -1,379 +0,0 @@
-
-#include "Vector/vector_dist.hpp"
-#include "Decomposition/CartDecomposition.hpp"
-#include "data_type/aggregate.hpp"
-#include "Plot/GoogleChart.hpp"
-#include "Plot/util.hpp"
-
-/*
- * ### WIKI 1 ###
- *
- * ## Molecular Dynamic with Lennard-Jones potential
- *
- * This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime
- *
- * ### WIKI END ###
- *
- */
-
-constexpr int velocity = 0;
-constexpr int force = 1;
-
-/* ### WIKI 11 ###
- *
- * The function to calculate the forces between particles. It require the vector of particles
- * Cell list and scaling factor.
- *
- */
-void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double L)
-{
-
- // ### WIKI 12 ###
- //
- // Update the cell list from the actual particle configuration
- //
- //
- vd.updateCellList(NN);
-
- // ### WIKI 13 ###
- //
- // Calculate the forces
- //
- auto it2 = vd.getDomainIterator();
- while (it2.isNext())
- {
- auto p = it2.get();
-
- Point<3,double> xp = vd.getPos(p);
-
- vd.template getProp<force>(p)[0] = 0.0;
- vd.template getProp<force>(p)[1] = 0.0;
- vd.template getProp<force>(p)[2] = 0.0;
-
- // For each neighborhood particle
- auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
-
- while (Np.isNext())
- {
- // Neighborhood particle q
- auto q = Np.get();
-
- if (q == p.getKey()) {++Np; continue;};
-
- // repulsive
- Point<3,double> xq = vd.getPos(q);
- Point<3,double> r = xp - xq;
-
- // take the norm, normalize
- float rn = r.norm();
- r /= rn;
- rn *= L;
-
- // Calculate the force, using pow is slower
- Point<3,double> f = 24.0*(2.0 / (rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn) - 1.0 / (rn*rn*rn*rn*rn*rn*rn)) * r;
-
- // we sum the force produced by q on p
- vd.template getProp<force>(p)[0] += f.get(0);
- vd.template getProp<force>(p)[1] += f.get(1);
- vd.template getProp<force>(p)[2] += f.get(2);
-
- ++Np;
- }
-
- ++it2;
- }
-}
-
-/* ### WIKI 14 ###
- *
- * The function to calculate the total energy. It require the same parameter as calculate forces
- *
- */
-double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double L)
-{
- // ### WIKI 15 ###
- //
- // Reset the counter for the energy and
- // update the cell list from the actual particle configuration
- //
- //
- double E = 0.0;
- vd.updateCellList(NN);
-
- // ### WIKI 16 ###
- //
- // Calculate the forces
- //
- auto it2 = vd.getDomainIterator();
- while (it2.isNext())
- {
- auto p = it2.get();
-
- Point<3,double> xp = vd.getPos(p);
-
- vd.template getProp<force>(p)[0] = 0.0;
- vd.template getProp<force>(p)[1] = 0.0;
- vd.template getProp<force>(p)[2] = 0.0;
-
- // For each neighborhood particle
- auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
-
- while (Np.isNext())
- {
- // Neighborhood particle q
- auto q = Np.get();
-
- if (q == p.getKey()) {++Np; continue;};
-
- Point<3,double> xq = vd.getPos(q);
- Point<3,double> r = xp - xq;
-
- // take the normalized direction
- float rn = r.norm();
- r /= rn;
-
- rn *= L;
-
- // potential energy (using pow is slower)
- E += 4.0 * ( 1.0 / (rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn*rn) + 1.0 / ( rn*rn*rn*rn*rn*rn) );
-
- // Kinetic energy
- E += sqrt(vd.template getProp<force>(p)[0]*vd.template getProp<force>(p)[0] +
- vd.template getProp<force>(p)[1]*vd.template getProp<force>(p)[1] +
- vd.template getProp<force>(p)[2]*vd.template getProp<force>(p)[2]);
-
- ++Np;
- }
-
- ++it2;
- }
-
- return E;
-}
-
-int main(int argc, char* argv[])
-{
- //
- // ### WIKI 2 ###
- //
- // Here we define important parameters or the simulation, time step integration,
- // size of the box, and cut-off radius of the interaction
- //
- double dt = 0.005;
- double L = 10.0;
- float r_cut = 0.3;
-
- openfpm::vector<double> x;
- openfpm::vector<openfpm::vector<double>> y;
-
- //
- // ### WIKI 2 ###
- //
- // Here we Initialize the library, we create a Box that define our domain, boundary conditions, ghost
- // and the grid size
- //
- openfpm_init(&argc,&argv);
- Vcluster & v_cl = create_vcluster();
-
- // we create a 10x10x10 Grid iterator
- size_t sz[3] = {10,10,10};
-
- Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
-
- // Boundary conditions
- size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
-
- // ghost, big enough to contain the interaction radius
- Ghost<3,float> ghost(r_cut);
-
- //
- // ### WIKI 3 ###
- //
- // Here we define a distributed vector in 3D, containing 3 properties, a
- // scalar double, a vector double[3], and a tensor or rank 2 double[3][3].
- // In this case the vector contain 0 particles in total
- //
- vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost);
-
- //
- // ### WIKI 4 ###
- //
- // We define a grid iterator, to create particles on a grid like way.
- // An important note is that the grid iterator, iterate only on the
- // local nodes for each processor for example suppose to have a domain like
- // the one in figure
- //
- // +---------+
- // |* * *|* *|
- // | 2 | |
- // |* * *|* *|
- // | --- |
- // |* *|* * *|
- // | | |
- // |* *|* * *|
- // | | 1 |
- // |* *|* * *|
- // +---------+
- //
- // divided in 2 processors, the processor 1 will iterate only on the points
- // inside the portion of space marked with one. Also the grid iterator follow the
- // boundary condition specified in vector. For a perdiodic 2D 5x5 grid we have
- //
- // +---------+
- // * * * * * |
- // | |
- // * * * * * |
- // | |
- // * * * * * |
- // | |
- // * * * * * |
- // | |
- // *-*-*-*-*-+
- //
- // Because the right border is equivalent to the left border, while for a non periodic we have the
- // following distribution of points
- //
- // *-*-*-*-*
- // | |
- // * * * * *
- // | |
- // * * * * *
- // | |
- // * * * * *
- // | |
- // *-*-*-*-*
- //
- // The loop will place particles on a grid like on each processor
- //
- auto it = vd.getGridIterator(sz);
-
- while (it.isNext())
- {
- vd.add();
-
- auto key = it.get();
-
- vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
- vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
- vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
-
- vd.template getLastProp<velocity>()[0] = 0.0;
- vd.template getLastProp<velocity>()[1] = 0.0;
- vd.template getLastProp<velocity>()[2] = 0.0;
-
- vd.template getLastProp<force>()[0] = 0.0;
- vd.template getLastProp<force>()[1] = 0.0;
- vd.template getLastProp<force>()[2] = 0.0;
-
- ++it;
- }
-
- vd.map();
-
- //
- // ### WIKI 5 ###
- //
- // we get the cell list to compute the neighborhood of the particles
- auto NN = vd.getCellList(r_cut);
-
- // calculate forces a(tn)
- calc_forces(vd,NN,L);
- unsigned long int f = 0;
-
- //
- // ### WIKI 6 ###
- //
- // Here we do 100 MD steps using verlet integrator
- //
- // $$ \vec{v}(t_{n+1/2}) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) $$
- // $$ \vec{x}(t_{n}) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) $$
- //
- // calculate the forces $$ \vec{a} (t_{n}) $$ from $$ \vec{x} (t_{n}) $$
- //
- // $$ \vec{v}(t_{n+1}) = \vec{v}_p(t_n+1/2) + \frac{1}{2} \delta t \vec{a}(t_n+1) $$
- //
- //
- for (size_t i = 0; i < 10000 ; i++)
- {
- auto it3 = vd.getDomainIterator();
-
- while (it3.isNext())
- {
- auto p = it3.get();
-
- // here we calculate v(tn + 0.5)
- vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
- vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
- vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
-
- // here we calculate x(tn + 1)
- vd.getPos(p)[0] += vd.template getProp<velocity>(p)[0]*dt;
- vd.getPos(p)[1] += vd.template getProp<velocity>(p)[1]*dt;
- vd.getPos(p)[2] += vd.template getProp<velocity>(p)[2]*dt;
-
- ++it3;
- }
-
- vd.map();
- vd.template ghost_get<>();
-
- // calculate forces or a(tn + 1)
- calc_forces(vd,NN,L);
-
- auto it4 = vd.getDomainIterator();
-
- while (it4.isNext())
- {
- auto p = it4.get();
-
- // here we calculate v(tn + 1)
- vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
- vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
- vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
-
- ++it4;
- }
-
- if (i % 100 == 0)
- {
- vd.deleteGhost();
- vd.write("particles_",f);
- double energy = calc_energy(vd,NN,L);
- auto & vcl = create_vcluster();
- vcl.sum(energy);
- vcl.execute();
-
- vd.ghost_get<>();
-
- x.add(i);
- y.add({energy});
- if (vcl.getProcessUnitID() == 0)
- std::cout << "Energy: " << energy << std::endl;
-
- f++;
- }
- }
-
- // Google charts options
- GCoptions options;
-
- options.title = std::string("Energy with time");
- options.yAxis = std::string("Energy");
- options.xAxis = std::string("iteration");
- options.lineWidth = 1.0;
-
- GoogleChart cg;
- cg.AddLinesGraph(x,y,options);
- cg.write("gc_plot2_out.html");
-
- //
- // ### WIKI 10 ###
- //
- // Deinitialize the library
- //
- openfpm_finalize();
-}
-
-
-
-
diff --git a/openfpm_numerics b/openfpm_numerics
index ce33014d..df19a7de 160000
--- a/openfpm_numerics
+++ b/openfpm_numerics
@@ -1 +1 @@
-Subproject commit ce33014d668526a01d02cc1bb28e1ee91e6ab0bb
+Subproject commit df19a7de48ebbbea0c3c2a1f8ebca60bdfef6618
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index 57929aea..674eb447 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -1405,9 +1405,7 @@ public:
Point<dim,St> offset = getOffset(i);
vtk_g.add(loc_grid.get(i),offset,cd_sm.getCellBox().getP2(),gdb_ext.get(i).Dbox);
}
- vtk_g.write(output + "_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk");
-
- write_ie_boxes(output);
+ vtk_g.write(output + "_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk");
return true;
}
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 14ce55bf..e2e9a195 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -1600,6 +1600,23 @@ public:
v_prp.resize(g_m);
}
+ /*! \brief Resize the vector (locally)
+ *
+ * \warning It automaticallt delete the ghosts
+ *
+ * \param rs
+ *
+ */
+ void resize(size_t rs)
+ {
+ deleteGhost();
+
+ v_pos.resize(rs);
+ v_prp.resize(rs);
+
+ g_m = rs;
+ }
+
/*! \brief Output particle position and properties
*
* \param out output
--
GitLab