From ca16d911ecca49c27ba7663811c430da48725d05 Mon Sep 17 00:00:00 2001 From: Pietro Incardona <i-bird@localhost.localdomain> Date: Fri, 15 Apr 2016 01:55:38 +0200 Subject: [PATCH] Fixing getPos --- CHANGELOG.md | 4 +- .../PSE/0_Derivative_approx_1D/main.cpp | 16 +- .../main_float128.cpp | 16 +- example/Numerics/PSE/1_Diffusion_1D/main.cpp | 16 +- example/SE/0_classes/main.cpp | 4 +- example/Vector/0_simple/main.cpp | 6 +- example/Vector/1_celllist/main.cpp | 12 +- example/Vector/1_verlet/main.cpp | 10 +- example/Vector/2_molecular_dynamic/Makefile | 21 + example/Vector/2_molecular_dynamic/config.cfg | 2 + .../2_molecular_dynamic/gc_plot2_out.html | 129 ++++++ example/Vector/2_molecular_dynamic/main.cpp | 376 ++++++++++++++++++ images/vector.cpp | 12 +- src/Decomposition/nn_processor.hpp | 2 +- src/Vector/vector_dist.hpp | 53 ++- src/Vector/vector_dist_unit_test.hpp | 118 +++--- 16 files changed, 671 insertions(+), 126 deletions(-) create mode 100644 example/Vector/2_molecular_dynamic/Makefile create mode 100644 example/Vector/2_molecular_dynamic/config.cfg create mode 100644 example/Vector/2_molecular_dynamic/gc_plot2_out.html create mode 100644 example/Vector/2_molecular_dynamic/main.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 9995b933..5fcf57ef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@ All notable changes to this project will be documented in this file. ## [0.3.0] - ### Added -- Nothing to report +- Molacular Dynamic example ### Fixed - Nothing to report @@ -13,6 +13,8 @@ All notable changes to this project will be documented in this file. - Eliminated global_v_cluster, init_global_v_cluster, delete_global_v_cluster, substituted by create_vcluster, openfpm_init, openfpm_delete +- CartDecomposition parameter for the distributed structures is now optional +- template getPos<0>(), substituted by getPos() ## [0.2.1] - diff --git a/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp b/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp index 0116c627..06007039 100644 --- a/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp +++ b/example/Numerics/PSE/0_Derivative_approx_1D/main.cpp @@ -127,7 +127,7 @@ int main(int argc, char* argv[]) auto key = it2.get(); // set the position of the particles - vd.template getPos<0>(key)[0] = (key.getKey() + base) * spacing; + 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); @@ -185,19 +185,19 @@ int main(int argc, char* argv[]) auto key = it.get(); // set the position of the particles - if (m_pad.isInsideNB(vd.template getPos<0>(key)) == true) + if (m_pad.isInsideNB(vd.getPos(key)) == true) { vd.add(); - vd.template getLastPos<0>()[0] = - vd.template getPos<0>(key)[0]; + 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.template getPos<0>(key)) == true) + if (m_pad2.isInsideNB(vd.getPos(key)) == true) { vd.add(); - vd.template getLastPos<0>()[0] = 2.0 * box.getHigh(0) - vd.template getPos<0>(key)[0]; - vd.template getLastProp<0>() = f_xex2(vd.template getLastPos<0>()[0]); + vd.getLastPos()[0] = 2.0 * box.getHigh(0) - vd.getPos(key)[0]; + vd.template getLastProp<0>() = f_xex2(vd.getLastPos()[0]); } ++it; @@ -237,7 +237,7 @@ int main(int argc, char* argv[]) vect_dist_key_dx key = it_p.get(); // Get the position of the particles - Point<1,double> p = vd.template getPos<0>(key); + Point<1,double> p = vd.getPos(key); // We are not interested in calculating out the domain // note added padding particle are considered domain particles @@ -261,7 +261,7 @@ int main(int argc, char* argv[]) if (nnp != key.getKey()) { // W(x-y) - double ker = lker.value(p,vd.template getPos<0>(nnp)); + double ker = lker.value(p,vd.getPos(nnp)); // f(y) double prp_y = vd.template getProp<0>(nnp); diff --git a/example/Numerics/PSE/1_Derivative_approx_1D_mp/main_float128.cpp b/example/Numerics/PSE/1_Derivative_approx_1D_mp/main_float128.cpp index a5a8b2ad..b350abd8 100644 --- a/example/Numerics/PSE/1_Derivative_approx_1D_mp/main_float128.cpp +++ b/example/Numerics/PSE/1_Derivative_approx_1D_mp/main_float128.cpp @@ -130,7 +130,7 @@ int main(int argc, char* argv[]) auto key = it2.get(); // set the position of the particles - vd.template getPos<0>(key)[0] = (key.getKey() + base) * spacing; + 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); @@ -188,19 +188,19 @@ int main(int argc, char* argv[]) auto key = it.get(); // set the position of the particles - if (m_pad.isInsideNB(vd.template getPos<0>(key)) == true) + if (m_pad.isInsideNB(vd.getPos(key)) == true) { vd.add(); - vd.template getLastPos<0>()[0] = - vd.template getPos<0>(key)[0]; + 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.template getPos<0>(key)) == true) + if (m_pad2.isInsideNB(vd.getPos(key)) == true) { vd.add(); - vd.template getLastPos<0>()[0] = 2.0 * box.getHigh(0) - vd.template getPos<0>(key)[0]; - vd.template getLastProp<0>() = f_xex2(vd.template getLastPos<0>()[0]); + vd.getLastPos()[0] = 2.0 * box.getHigh(0) - vd.getPos(key)[0]; + vd.template getLastProp<0>() = f_xex2(vd.getLastPos()[0]); } ++it; @@ -241,7 +241,7 @@ int main(int argc, char* argv[]) vect_dist_key_dx key = it_p.get(); // Get the position of the particles - Point<1,float128> p = vd.template getPos<0>(key); + Point<1,float128> p = vd.getPos(key); // We are not interested in calculating out the domain // note added padding particle are considered domain particles @@ -265,7 +265,7 @@ int main(int argc, char* argv[]) if (nnp != key.getKey()) { // W(x-y) - float128 ker = lker.value(p,vd.template getPos<0>(nnp)); + float128 ker = lker.value(p,vd.getPos(nnp)); // f(y) float128 prp_y = vd.template getProp<0>(nnp); diff --git a/example/Numerics/PSE/1_Diffusion_1D/main.cpp b/example/Numerics/PSE/1_Diffusion_1D/main.cpp index 38937a37..0fb0d3db 100644 --- a/example/Numerics/PSE/1_Diffusion_1D/main.cpp +++ b/example/Numerics/PSE/1_Diffusion_1D/main.cpp @@ -74,7 +74,7 @@ template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key, if (nnp != key.getKey()) { // W(x-y) - double ker = lker.value(p,vd.template getPos<0>(nnp)); + double ker = lker.value(p,vd.getPos(nnp)); // f(y) double prp_y = vd.template getProp<0>(nnp); @@ -180,7 +180,7 @@ int main(int argc, char* argv[]) auto key = it2.get(); // set the position of the particles - vd.template getPos<0>(key)[0] = (key.getKey() + base) * spacing; + 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); @@ -238,19 +238,19 @@ int main(int argc, char* argv[]) auto key = it.get(); // set the position of the particles - if (m_pad.isInsideNB(vd.template getPos<0>(key)) == true) + if (m_pad.isInsideNB(vd.getPos(key)) == true) { vd.add(); - vd.template getLastPos<0>()[0] = - vd.template getPos<0>(key)[0]; + 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.template getPos<0>(key)) == true) + if (m_pad2.isInsideNB(vd.getPos(key)) == true) { vd.add(); - vd.template getLastPos<0>()[0] = 2.0 * box.getHigh(0) - vd.template getPos<0>(key)[0]; - vd.template getLastProp<0>() = f_xex2(vd.template getLastPos<0>()[0]); + vd.getLastPos()[0] = 2.0 * box.getHigh(0) - vd.getPos(key)[0]; + vd.template getLastProp<0>() = f_xex2(vd.getLastPos()[0]); } ++it; @@ -284,7 +284,7 @@ int main(int argc, char* argv[]) // key vect_dist_key_dx key = it_p.get(); - Point<1,double> p = vd.template getPos<0>(key); + Point<1,double> p = vd.getPos(key); // We are not interested in calculating out the domain // note added padding particle are considered domain particles diff --git a/example/SE/0_classes/main.cpp b/example/SE/0_classes/main.cpp index 952adf16..d476a76c 100644 --- a/example/SE/0_classes/main.cpp +++ b/example/SE/0_classes/main.cpp @@ -98,7 +98,7 @@ int main(int argc, char* argv[]) try { vect_dist_key_dx vt(5048); - auto it = vd.getPos<0>(vt); + auto it = vd.getPos(vt); } catch (size_t e) { @@ -185,7 +185,7 @@ int main(int argc, char* argv[]) try { vect_dist_key_dx vt(0); - auto it = vd1->getPos<0>(vt); + auto it = vd1->getPos(vt); } catch (size_t e) { diff --git a/example/Vector/0_simple/main.cpp b/example/Vector/0_simple/main.cpp index 9f26e5fc..85922d6f 100644 --- a/example/Vector/0_simple/main.cpp +++ b/example/Vector/0_simple/main.cpp @@ -91,8 +91,8 @@ int main(int argc, char* argv[]) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); ++it; } @@ -121,7 +121,7 @@ int main(int argc, char* argv[]) auto key = it.get(); // The template parameter is unuseful and will probably disappear - if (ct.isLocal(vd.template getPos<0>(key)) == false) + if (ct.isLocal(vd.getPos(key)) == false) std::cerr << "Error particle is not local" << "\n"; // set the all the properties to 0.0 diff --git a/example/Vector/1_celllist/main.cpp b/example/Vector/1_celllist/main.cpp index 45d1e56b..34ea6ede 100644 --- a/example/Vector/1_celllist/main.cpp +++ b/example/Vector/1_celllist/main.cpp @@ -103,9 +103,9 @@ int main(int argc, char* argv[]) auto key = it.get(); - vd.template getLastPos<0>()[0] = key.get(0) * it.getSpacing(0); - vd.template getLastPos<0>()[1] = key.get(1) * it.getSpacing(1); - vd.template getLastPos<0>()[2] = key.get(2) * it.getSpacing(2); + 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); ++it; } @@ -133,9 +133,9 @@ int main(int argc, char* argv[]) { auto p = it2.get(); - Point<3,float> xp = vd.getPos<0>(p); + Point<3,float> xp = vd.getPos(p); - auto Np = NN.getIterator(NN.getCell(vd.getPos<0>(p))); + auto Np = NN.getIterator(NN.getCell(vd.getPos(p))); while (Np.isNext()) { @@ -143,7 +143,7 @@ int main(int argc, char* argv[]) // repulsive - Point<3,float> xq = vd.getPos<0>(q); + Point<3,float> xq = vd.getPos(q); Point<3,float> f = (xp - xq); // we sum the distance of all the particles diff --git a/example/Vector/1_verlet/main.cpp b/example/Vector/1_verlet/main.cpp index b9dae2ad..1c8c44a3 100644 --- a/example/Vector/1_verlet/main.cpp +++ b/example/Vector/1_verlet/main.cpp @@ -102,9 +102,9 @@ int main(int argc, char* argv[]) auto key = it.get(); - vd.template getLastPos<0>()[0] = key.get(0) * it.getSpacing(0); - vd.template getLastPos<0>()[1] = key.get(1) * it.getSpacing(1); - vd.template getLastPos<0>()[2] = key.get(2) * it.getSpacing(2); + 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); ++it; } @@ -138,14 +138,14 @@ int main(int argc, char* argv[]) for (size_t i = 0 ; i < verlet.size() ; i++) { - Point<3,float> p = vd.getPos<0>(i); + Point<3,float> p = vd.getPos(i); // for each neighborhood particle for (size_t j = 0 ; j < verlet.get(i).size() ; j++) { auto & NN = verlet.get(i); - Point<3,float> q = vd.getPos<0>(NN.get(j)); + Point<3,float> q = vd.getPos(NN.get(j)); // some non-sense calculation as usage demo diff --git a/example/Vector/2_molecular_dynamic/Makefile b/example/Vector/2_molecular_dynamic/Makefile new file mode 100644 index 00000000..0944bdfc --- /dev/null +++ b/example/Vector/2_molecular_dynamic/Makefile @@ -0,0 +1,21 @@ +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: cell + +.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 new file mode 100644 index 00000000..1eecbac3 --- /dev/null +++ b/example/Vector/2_molecular_dynamic/config.cfg @@ -0,0 +1,2 @@ +[pack] +files = main.cpp Makefile diff --git a/example/Vector/2_molecular_dynamic/gc_plot2_out.html b/example/Vector/2_molecular_dynamic/gc_plot2_out.html new file mode 100644 index 00000000..a0d68196 --- /dev/null +++ b/example/Vector/2_molecular_dynamic/gc_plot2_out.html @@ -0,0 +1,129 @@ +<html> + <head> + <script type="text/javascript" src="https://www.gstatic.com/charts/loader.js"></script> + <script type="text/javascript"> + google.charts.load('current', {'packages':['corechart']}); + google.charts.setOnLoadCallback(drawVisualization); + + + function drawVisualization() { +var data0 = new google.visualization.DataTable(); +data0.addColumn('number','iteration'); +data0.addColumn('number','line0'); +data0.addRows([ +[0,18010.3], +[100,18010.5], +[200,16751.9], +[300,23684.3], +[400,23433.6], +[500,22997.1], +[600,24500.7], +[700,23194.2], +[800,23649.7], +[900,23663.8], +[1000,23838.5], +[1100,24004.5], +[1200,23994.6], +[1300,23702.6], +[1400,23636.1], +[1500,23877.8], +[1600,24171.7], +[1700,23985], +[1800,23507], +[1900,23308.7], +[2000,23820.8], +[2100,23417.1], +[2200,23607.6], +[2300,23670.5], +[2400,23722], +[2500,23757.4], +[2600,23914.1], +[2700,23541.1], +[2800,23037.1], +[2900,23520], +[3000,24058.2], +[3100,23087.6], +[3200,22901], +[3300,23182.9], +[3400,23352.6], +[3500,22819.6], +[3600,22970.7], +[3700,22831.1], +[3800,22506.7], +[3900,22614.4], +[4000,22716.8], +[4100,22969.7], +[4200,22976.2], +[4300,22690.9], +[4400,22382], +[4500,22731.3], +[4600,22547.5], +[4700,22194], +[4800,22350.9], +[4900,22686.1], +[5000,22328.7], +[5100,22629.4], +[5200,22597.8], +[5300,22272.9], +[5400,22649.2], +[5500,22057.9], +[5600,22137.9], +[5700,22347.8], +[5800,22465.9], +[5900,22530.8], +[6000,22393.6], +[6100,22117.9], +[6200,22008.9], +[6300,21868.5], +[6400,21330.7], +[6500,21390.1], +[6600,21782.7], +[6700,21393.1], +[6800,21718.1], +[6900,21874.9], +[7000,21867.5], +[7100,21453.3], +[7200,21734.5], +[7300,21439.8], +[7400,21454.8], +[7500,21700.6], +[7600,21638.9], +[7700,20965.1], +[7800,20846.4], +[7900,20689.3], +[8000,20758.1], +[8100,21154.9], +[8200,20726.4], +[8300,20594.4], +[8400,20460.2], +[8500,20694], +[8600,20509.1], +[8700,20492.3], +[8800,20342.2], +[8900,20397.8], +[9000,20673.6], +[9100,20500.1], +[9200,20203.3], +[9300,20221.9], +[9400,20219], +[9500,20021.4], +[9600,20164], +[9700,20237.4], +[9800,20089.2], +[9900,20061.7], +]); +var options0= { +title : 'Energy with time', +vAxis: {title: 'Energy'}, +hAxis: {title: 'iteration'}, +curveType: 'function', +lineWidth: 1, +intervals: { 'style':'area' }, +}; +var chart = new google.visualization.ComboChart(document.getElementById('chart_div0'));chart.draw(data0, options0); +}</script> +</head> +<body> +<div id="chart_div0" style="width: 900px; height: 500px;"></div> +</body> +</html> diff --git a/example/Vector/2_molecular_dynamic/main.cpp b/example/Vector/2_molecular_dynamic/main.cpp new file mode 100644 index 00000000..c5ae6136 --- /dev/null +++ b/example/Vector/2_molecular_dynamic/main.cpp @@ -0,0 +1,376 @@ + +#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<0>(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<0>(p))); + + while (Np.isNext()) + { + // Neighborhood particle q + auto q = Np.get(); + + if (q == p.getKey()) {++Np; continue;}; + + // repulsive + Point<3,double> xq = vd.getPos<0>(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<0>(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<0>(p))); + + while (Np.isNext()) + { + // Neighborhood particle q + auto q = Np.get(); + + if (q == p.getKey()) {++Np; continue;}; + + Point<3,double> xq = vd.getPos<0>(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.template getLastPos<0>()[0] = key.get(0) * it.getSpacing(0); + vd.template getLastPos<0>()[1] = key.get(1) * it.getSpacing(1); + vd.template getLastPos<0>()[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.template getPos<0>(p)[0] += vd.template getProp<velocity>(p)[0]*dt; + vd.template getPos<0>(p)[1] += vd.template getProp<velocity>(p)[1]*dt; + vd.template getPos<0>(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.write("particles",f); + double energy = calc_energy(vd,NN,L); + auto & vcl = create_vcluster(); + vcl.sum(energy); + vcl.execute(); + + 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/images/vector.cpp b/images/vector.cpp index 6c0e7a21..d824ef40 100644 --- a/images/vector.cpp +++ b/images/vector.cpp @@ -46,11 +46,11 @@ int main(int argc, char* argv[]) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); - vd.template getProp<1>(key)[0] = sin(10.0*vd.template getPos<s::x>(key)[0]); - vd.template getProp<1>(key)[1] = sin(10.0*vd.template getPos<s::x>(key)[1]); + vd.template getProp<1>(key)[0] = sin(10.0*vd.getPos(key)[0]); + vd.template getProp<1>(key)[1] = sin(10.0*vd.getPos(key)[1]); ++it; } @@ -77,8 +77,8 @@ int main(int argc, char* argv[]) { auto key = it.get(); - vd.template getPos<0>(key)[0] += 0.005; - vd.template getPos<0>(key)[1] += 0.005; + vd.getPos(key)[0] += 0.005; + vd.getPos(key)[1] += 0.005; vd.template getProp<1>(key)[0] = 0.005; vd.template getProp<1>(key)[1] = 0.005; diff --git a/src/Decomposition/nn_processor.hpp b/src/Decomposition/nn_processor.hpp index e7497c7c..bec37988 100755 --- a/src/Decomposition/nn_processor.hpp +++ b/src/Decomposition/nn_processor.hpp @@ -128,7 +128,7 @@ class nn_prcs break; case -1: bp.setLow(k,domain.getLow(k)); - bp.setHigh(k,ghost.getHigh(k)); + bp.setHigh(k,domain.getLow(k)+ghost.getHigh(k)); shift.get(k) = domain.getHigh(k)-domain.getLow(k); break; } diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp index c36a4746..6eb6d579 100644 --- a/src/Vector/vector_dist.hpp +++ b/src/Vector/vector_dist.hpp @@ -734,9 +734,9 @@ public: * \return the position of the element in space * */ - template<unsigned int id> inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<id>(vec_key.getKey())) + inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey())) { - return v_pos.template get<id>(vec_key.getKey()); + return v_pos.template get<0>(vec_key.getKey()); } /*! \brief Get the property of an element @@ -982,9 +982,9 @@ public: * \return the position of the element in space * */ - template<unsigned int id> inline auto getLastPos() -> decltype(v_pos.template get<id>(0)) + inline auto getLastPos() -> decltype(v_pos.template get<0>(0)) { - return v_pos.template get<id>(g_m - 1); + return v_pos.template get<0>(g_m - 1); } /*! \brief Get the property of the last element @@ -1020,6 +1020,32 @@ public: return getCellList(r_cut, g); } + /*! \brief Update a cell list using the stored particles + * + * \tparam CellL CellList type to construct + * + * \param cell_list Cell list to update + * + */ + template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > void updateCellList(CellL & cell_list) + { + // Clear the cell list from the previous particles + cell_list.clear(); + + // for each particle add the particle to the cell list + + auto it = getIterator(); + + while (it.isNext()) + { + auto key = it.get(); + + cell_list.add(this->template getPos(key), key.getKey()); + + ++it; + } + } + /*! \brief Construct a cell list starting from the stored particles * * It differ from the get getCellList for an additional parameter, in case the @@ -1055,18 +1081,7 @@ public: cell_list.Initialize(pbox, div); - // for each particle add the particle to the cell list - - auto it = getIterator(); - - while (it.isNext()) - { - auto key = it.get(); - - cell_list.add(this->template getPos<0>(key), key.getKey()); - - ++it; - } + updateCellList(cell_list); return cell_list; } @@ -1096,7 +1111,7 @@ public: vect_dist_key_dx key = it_p.get(); // Get the position of the particles - Point<dim, St> p = this->template getPos<0>(key); + Point<dim, St> p = this->template getPos(key); // Clear the neighborhood of the particle verlet.get(key.getKey()).clear(); @@ -1114,7 +1129,7 @@ public: continue; } - Point<dim, St> q = this->template getPos<0>(nnp); + Point<dim, St> q = this->template getPos(nnp); if (p.distance2(q) < r_cut2) verlet.get(key.getKey()).add(nnp); @@ -1285,7 +1300,7 @@ public: while (it.isNext()) { - size_t v = cdsm.getCell(this->template getPos<0>(it.get())); + size_t v = cdsm.getCell(this->template getPos(it.get())); dec.addComputationCost(v, 1); diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp index fc1ebc61..35259073 100644 --- a/src/Vector/vector_dist_unit_test.hpp +++ b/src/Vector/vector_dist_unit_test.hpp @@ -32,7 +32,7 @@ template<unsigned int dim> size_t total_n_part_lc(vector_dist<dim,float, Point_t { auto key = it2.get(); - noOut &= ct.isLocal(vd.template getPos<s::x>(key)); + noOut &= ct.isLocal(vd.getPos(key)); cnt++; @@ -69,10 +69,10 @@ template<unsigned int dim> inline void count_local_n_local(vector_dist<dim,float { auto key = it.get(); // Check if it is in the domain - if (box.isInsideNP(vd.template getPos<s::x>(key)) == true) + if (box.isInsideNP(vd.getPos(key)) == true) { // Check if local - if (ct.isLocalBC(vd.template getPos<s::x>(key),bc) == true) + if (ct.isLocalBC(vd.getPos(key),bc) == true) l_cnt++; else nl_cnt++; @@ -83,7 +83,7 @@ template<unsigned int dim> inline void count_local_n_local(vector_dist<dim,float } // Check that all particles are inside the Domain + Ghost part - if (dom_ext.isInside(vd.template getPos<s::x>(key)) == false) + if (dom_ext.isInside(vd.getPos(key)) == false) n_out++; ++it; @@ -164,8 +164,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost ) // set the particle position - vd.template getPos<s::x>(key_v)[0] = key.get(0) * spacing[0] + m_spacing[0]; - vd.template getPos<s::x>(key_v)[1] = key.get(1) * spacing[1] + m_spacing[1]; + vd.getPos(key_v)[0] = key.get(0) * spacing[0] + m_spacing[0]; + vd.getPos(key_v)[1] = key.get(1) * spacing[1] + m_spacing[1]; cobj++; @@ -193,7 +193,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost ) auto key = v_it2.get(); // fill with the processor ID where these particle live - vd.template getProp<p::s>(key) = vd.getPos<s::x>(key)[0] + vd.getPos<s::x>(key)[1] * 16; + vd.template getProp<p::s>(key) = vd.getPos(key)[0] + vd.getPos(key)[1] * 16; vd.template getProp<p::v>(key)[0] = v_cl.getProcessUnitID(); vd.template getProp<p::v>(key)[1] = v_cl.getProcessUnitID(); vd.template getProp<p::v>(key)[2] = v_cl.getProcessUnitID(); @@ -223,7 +223,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost ) auto key = g_it.get(); // Check the received data - BOOST_REQUIRE_EQUAL(vd.getPos<s::x>(key)[0] + vd.getPos<s::x>(key)[1] * 16,vd.template getProp<p::s>(key)); + BOOST_REQUIRE_EQUAL(vd.getPos(key)[0] + vd.getPos(key)[1] * 16,vd.template getProp<p::s>(key)); bool is_in = false; size_t b = 0; @@ -232,7 +232,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_ghost ) // check if the received data are in one of the ghost boxes for ( ; b < dec.getNEGhostBox() ; b++) { - if (dec.getEGhostBox(b).isInside(vd.getPos<s::x>(key)) == true ) + if (dec.getEGhostBox(b).isInside(vd.getPos(key)) == true ) { is_in = true; @@ -329,8 +329,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_2d ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); ++it; } @@ -349,7 +349,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_2d ) auto key = it2.get(); // Check if local - BOOST_REQUIRE_EQUAL(ct.isLocal(vd.template getPos<s::x>(key)),true); + BOOST_REQUIRE_EQUAL(ct.isLocal(vd.getPos(key)),true); cnt++; @@ -402,9 +402,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_3d ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); - vd.template getPos<s::x>(key)[2] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); + vd.getPos(key)[2] = ud(eg); ++it; } @@ -423,7 +423,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_3d ) auto key = it2.get(); // Check if local - BOOST_REQUIRE_EQUAL(ct.isLocal(vd.template getPos<s::x>(key)),true); + BOOST_REQUIRE_EQUAL(ct.isLocal(vd.getPos(key)),true); cnt++; @@ -486,8 +486,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_2d ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); ++it; } @@ -596,9 +596,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_3d ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); - vd.template getPos<s::x>(key)[2] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); + vd.getPos(key)[2] = ud(eg); ++it; } @@ -699,9 +699,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_random_walk ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); - vd.template getPos<s::x>(key)[2] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); + vd.getPos(key)[2] = ud(eg); ++it; } @@ -718,9 +718,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_random_walk ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] += 0.02 * ud(eg); - vd.template getPos<s::x>(key)[1] += 0.02 * ud(eg); - vd.template getPos<s::x>(key)[2] += 0.02 * ud(eg); + vd.getPos(key)[0] += 0.02 * ud(eg); + vd.getPos(key)[1] += 0.02 * ud(eg); + vd.getPos(key)[2] += 0.02 * ud(eg); ++it; } @@ -763,9 +763,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_map ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = 1.0; - vd.template getPos<s::x>(key)[1] = 1.0; - vd.template getPos<s::x>(key)[2] = 1.0; + vd.getPos(key)[0] = 1.0; + vd.getPos(key)[1] = 1.0; + vd.getPos(key)[2] = 1.0; ++it; } @@ -778,11 +778,11 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_map ) { auto key = it2.get(); - float f = vd.template getPos<s::x>(key)[0]; + float f = vd.getPos(key)[0]; BOOST_REQUIRE_EQUAL(f, 0.0); - f = vd.template getPos<s::x>(key)[1]; + f = vd.getPos(key)[1]; BOOST_REQUIRE_EQUAL(f, 0.0); - f = vd.template getPos<s::x>(key)[2]; + f = vd.getPos(key)[2]; BOOST_REQUIRE_EQUAL(f, 0.0); ++it2; @@ -815,9 +815,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_not_periodic_map ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = 1.0; - vd.template getPos<s::x>(key)[1] = 1.0; - vd.template getPos<s::x>(key)[2] = 1.0; + vd.getPos(key)[0] = 1.0; + vd.getPos(key)[1] = 1.0; + vd.getPos(key)[2] = 1.0; ++it; } @@ -830,11 +830,11 @@ BOOST_AUTO_TEST_CASE( vector_dist_not_periodic_map ) { auto key = it2.get(); - float f = vd.template getPos<s::x>(key)[0]; + float f = vd.getPos(key)[0]; BOOST_REQUIRE_EQUAL(f, 1.0); - f = vd.template getPos<s::x>(key)[1]; + f = vd.getPos(key)[1]; BOOST_REQUIRE_EQUAL(f, 1.0); - f = vd.template getPos<s::x>(key)[2]; + f = vd.getPos(key)[2]; BOOST_REQUIRE_EQUAL(f, 1.0); ++it2; @@ -876,15 +876,15 @@ BOOST_AUTO_TEST_CASE( vector_dist_out_of_bound_policy ) if (cnt < 1) { - vd.template getPos<s::x>(key)[0] = -0.06; - vd.template getPos<s::x>(key)[1] = -0.06; - vd.template getPos<s::x>(key)[2] = -0.06; + vd.getPos(key)[0] = -0.06; + vd.getPos(key)[1] = -0.06; + vd.getPos(key)[2] = -0.06; } else { - vd.template getPos<s::x>(key)[0] = 0.06; - vd.template getPos<s::x>(key)[1] = 0.06; - vd.template getPos<s::x>(key)[2] = 0.06; + vd.getPos(key)[0] = 0.06; + vd.getPos(key)[1] = 0.06; + vd.getPos(key)[2] = 0.06; } cnt++; @@ -954,9 +954,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] = ud(eg); - vd.template getPos<s::x>(key)[1] = ud(eg); - vd.template getPos<s::x>(key)[2] = ud(eg); + vd.getPos(key)[0] = ud(eg); + vd.getPos(key)[1] = ud(eg); + vd.getPos(key)[2] = ud(eg); ++it; } @@ -975,9 +975,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles ) { auto key = it.get(); - vd.template getPos<s::x>(key)[0] += 0.02 * ud(eg); - vd.template getPos<s::x>(key)[1] += 0.02 * ud(eg); - vd.template getPos<s::x>(key)[2] += 0.02 * ud(eg); + vd.getPos(key)[0] += 0.02 * ud(eg); + vd.getPos(key)[1] += 0.02 * ud(eg); + vd.getPos(key)[2] += 0.02 * ud(eg); ++it; } @@ -1000,9 +1000,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles ) { auto p = it2.get(); - Point<3,float> xp = vd.getPos<0>(p); + Point<3,float> xp = vd.getPos(p); - auto Np = NN.getIterator(NN.getCell(vd.getPos<0>(p))); + auto Np = NN.getIterator(NN.getCell(vd.getPos(p))); while (Np.isNext()) { @@ -1010,7 +1010,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles ) // repulsive - Point<3,float> xq = vd.getPos<0>(q); + Point<3,float> xq = vd.getPos(q); Point<3,float> f = (xp - xq); float distance = f.norm(); @@ -1081,9 +1081,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test ) auto key = it.get(); - vd.template getLastPos<s::x>()[0] = key.get(0) * it.getSpacing(0); - vd.template getLastPos<s::x>()[1] = key.get(1) * it.getSpacing(1); - vd.template getLastPos<s::x>()[2] = key.get(2) * it.getSpacing(2); + 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); ++it; } @@ -1130,14 +1130,14 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test ) size_t second_NN = 0; size_t third_NN = 0; - Point<3,float> p = vd.getPos<0>(i); + Point<3,float> p = vd.getPos(i); // for each neighborhood particle for (size_t j = 0 ; j < verlet.get(i).size() ; j++) { auto & NN = verlet.get(i); - Point<3,float> q = vd.getPos<0>(NN.get(j)); + Point<3,float> q = vd.getPos(NN.get(j)); float dist = p.distance(q); -- GitLab