diff --git a/Particle_Method_Algorithms/Diffusion_PSE_Euler_1D.hpp b/Particle_Method_Algorithms/Diffusion_PSE_Euler_1D.hpp
index 8256bf0f9409b7c36f2d03ac136a310be028ecb0..b706d551d1e33f6518c701bacf8374cc721a7396 100644
--- a/Particle_Method_Algorithms/Diffusion_PSE_Euler_1D.hpp
+++ b/Particle_Method_Algorithms/Diffusion_PSE_Euler_1D.hpp
@@ -61,7 +61,7 @@ namespace Diffusion_PSE_Euler_1D
protected:
bool stop(const GlobalVariable& g)override{
- return (g.endT<=g.t);
+ return (g.endT<g.t);
}
void neighborhood(const GlobalVariable& g)override{
this->cutoffRadius=g.rc;
diff --git a/Particle_Method_Algorithms/PSE_1D.hpp b/Particle_Method_Algorithms/PSE_1D.hpp
index 6de8b636b35123ecbc79d05a4e8cca3623a8bca5..9bc9f4293989d20b1f5df12449f7a8fb5dbe546b 100644
--- a/Particle_Method_Algorithms/PSE_1D.hpp
+++ b/Particle_Method_Algorithms/PSE_1D.hpp
@@ -22,6 +22,7 @@ namespace PSE_1D
public:
double concentration;
double accumulator;
+ //int interactionCOunt;
//overritten helper functions for the print out
std::string getLables() const override {
@@ -67,11 +68,15 @@ namespace PSE_1D
}
void interact(const GlobalVariable& g, Particle& p, Particle& q)override{
double d = q.position[0]-p.position[0];//calulates the euclidian distance squared
- p.accumulator+= (q.concentration-p.concentration)*exp(-d*d/(4*g.epsilon*g.epsilon ));//calculate the excange of strangth (concentration) between the neighbors
+ //p.accumulator+= (q.concentration-p.concentration)*exp(-d*d/(4*g.epsilon*g.epsilon ));//calculate the excange of strangth (concentration) between the neighbors
+ p.accumulator+= (q.concentration-p.concentration)*exp(-(d/g.epsilon)*(d/(4.0*g.epsilon)));
+ //p.interactionCOunt++;
}
void evolve(GlobalVariable& g,Particle& p)override{
- p.concentration= g.meshSpacing/(g.epsilon*g.epsilon*g.epsilon*2*sqrt(M_PI)) * p.accumulator;
- p.accumulator=0;//set to 0 to have a fresh accumulator for the next time step
+ //p.concentration= g.meshSpacing/(g.epsilon*g.epsilon*g.epsilon*2*sqrt(M_PI)) * p.accumulator;
+ p.concentration= (g.meshSpacing/g.epsilon) * (p.accumulator/g.epsilon) /(g.epsilon*2.0*sqrt(M_PI));
+ p.accumulator=0.0;//set to 0 to have a fresh accumulator for the next time step
+ //cout<<p.interactionCOunt<<endl;
}
};
diff --git a/prototype_interface/main b/prototype_interface/main
index a5f1fef86f8d968680a7a6164cf99cca2e2229c3..40af43f622050f3c6afe21d12fe78522f6bf1d3d 100755
Binary files a/prototype_interface/main and b/prototype_interface/main differ
diff --git a/prototype_interface/main.cpp b/prototype_interface/main.cpp
index e3c6c0eb069c5dfcd586ef8399a6bf641bc9251b..d5a8fc6e7d5e0e37d2e63b2e12b0de9b48c819b7 100644
--- a/prototype_interface/main.cpp
+++ b/prototype_interface/main.cpp
@@ -10,7 +10,6 @@
#include "PM_Point.hpp"
#include "Perfectly_Elastic_Collision_nD_unit_test.hpp"
-#include "Perfectly_Elastic_Collision_nD_correctness_check.hpp"
#include "functions_unit_tests.hpp"
#include "Diffusion_PSE_Euler_unit_test.hpp"
#include "Diffusion_PSE_Euler_correctness_check.hpp"
@@ -24,27 +23,49 @@ int main(int argc, const char * argv[]) {
auto start = std::chrono::high_resolution_clock::now();//time mesurment start time
//test cases, they should be used one by one
-// test_PM_Point();
-// Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_unit_test();
-// Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_big_unit_test();
+ // test_PM_Point();
+
+ // Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_unit_test();
+
+ // Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_big_unit_test();
+
// SPH::sph_unit_test();
- //Diffusion_PSE_Euler_3D::diffusion_PSE_Euler_unit_test();
-// functions_unit_test();
- //Diffusion_PSE_Euler_3D::diffusion_PSE_Euler_correctness_check(32);
- //Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_correctness_check();
-
- // //Conovergence Study
+ // Diffusion_PSE_Euler_3D::diffusion_PSE_Euler_unit_test();
+
+ // functions_unit_test();
+
+ // Diffusion_PSE_Euler_3D::diffusion_PSE_Euler_correctness_check(32);
+
+ // Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_correctness_check();
+
+ // //Conovergence Study with sin() (paper)
// for(int i=1;i<100;i++){
// PSE_1D::pse_convergence_unit_test(power(2,i));
// }
- // Time Complexity of increasing domain size with const particle density
- for(int i=0;i<100;i++){
- Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_big_time_test(0.5*power(2,i));
- }
+ // //Conovergence Study with gaus-curve
+ // for(int i=1;i<100;i++){
+ // PSE_1D::pse_convergence_unit_test_gauss(power(2,i));
+ // }
+
+ // //Time Complexity of increasing domain size constant particle density (paper)
+ // for(int i=1;i<100;i++){
+ // Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_big_time_test_02(0.5*power(i,2));
+ // }
+
+ // //Time Complexity of increasing number of particles (paper)
+ // for(int i=0;i<100;i++){
+ // Diffusion_PSE_Euler_3D::diffusion_PSE_Euler_runtime_test(power(2,i));
+ // }
+
+//random other tests
+ // // Time Complexity of increasing domain size with const particle density
+ // for(int i=0;i<100;i++){
+ // Perfectly_Elastic_Collision::perfectly_elastic_collision_2D_big_time_test(0.5*power(2,i));
+ // }
// //Time Complexity of increasing particle number constant domain size
// for(int i=0;i<100;i++){
@@ -57,6 +78,7 @@ int main(int argc, const char * argv[]) {
// }
+
//time mesurment
auto finish = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed = finish - start;
diff --git a/unit_tests/Diffusion_PSE_Euler_unit_test.hpp b/unit_tests/Diffusion_PSE_Euler_unit_test.hpp
index c9a71340c6ad5983ad48fd74299ea7e820792e6f..be384e288f28bcb45f9a4b0437f36ddbe089409d 100644
--- a/unit_tests/Diffusion_PSE_Euler_unit_test.hpp
+++ b/unit_tests/Diffusion_PSE_Euler_unit_test.hpp
@@ -18,7 +18,7 @@
//Particle Method Instances for Diffusion
namespace Diffusion_PSE_Euler_3D {
- double exactSolution(PM_Point<double,3> pos, double t, double diffConst){
+ double exactSolution(PM_Point<double,3> pos, double t, double diffConst){
double c=4*t*diffConst;
return exp(-pos.abs2()/c)/sqrt(power(M_PI*c,3));
}
@@ -51,9 +51,117 @@ namespace Diffusion_PSE_Euler_3D {
}
Particle_Method diff(g, mesh);
- diff.run(1000);
+ diff.run(10);
}
+
+ void diffusion_PSE_Euler_runtime_test(double scaling=64){
+ //Domain width is 1.0 and the center is 0.0 in all dimentions
+
+ GlobalVariable g;
+ g.diffusionConst=0.01;
+ g.meshSpacing=1.0/scaling;
+ g.epsilon=g.meshSpacing;
+ g.endT=0.5;
+ g.t=0.0;
+ g.dt=0.002;
+ g.rc=g.meshSpacing*3;
+
+ int halfsize=scaling/2.0;
+ g.meshDim = PM_Point<int,3> (2*halfsize+1);
+ g.domain_min= -halfsize*g.meshSpacing*PM_Point<double,3>(1);
+ g.domain_max= halfsize*g.meshSpacing*PM_Point<double,3>(1);
+
+ vector<Particle> mesh;
+ int length=g.meshDim[0]*g.meshDim[1]*g.meshDim[2];
+ mesh.resize(length);
+ mesh[(length-1)/2].concentration=1.0/power(g.meshSpacing,3);//initial peak at the center
+
+ //just for storing the data
+ for (int i=0; i<mesh.size(); i++){
+ mesh[i].position=(PM_Point<double, 3>)toVectorialIndex(i, g.meshDim)*g.meshSpacing+g.domain_min;
+ }
+
+ Particle_Method diff(g, mesh);
+
+ auto start = std::chrono::high_resolution_clock::now();//time mesurment start time
+ diff.run(-1);
+ auto finish = std::chrono::high_resolution_clock::now();
+ std::chrono::duration<double> elapsed = finish - start;
+
+
+
+ //init native C++
+ double t=0;
+ double dt=0.002; //std=0.002
+ double tend=0.5;
+
+ double meshSpacing=1.0/scaling;
+ double rc=meshSpacing*3;//cutof radius std=3
+ double epsilon=meshSpacing;//kernal width std=meshSpacing
+ double diffusionConst=0.01;
+
+ int halfNo=scaling/2.0;;
+ int no1D=halfNo*2+1;
+ int noParticles=no1D*no1D*no1D;
+ double x[noParticles];
+ double y[noParticles];
+ double z[noParticles];
+ double concentration[noParticles];
+ double accumulator[noParticles];
+
+ for(int i=0;i<no1D;i++){
+ for(int j=0;j<no1D;j++){
+ for(int k=0;k<no1D;k++){
+ x[k+j*no1D+i*no1D*no1D]=k*meshSpacing-halfNo*meshSpacing;
+ y[k+j*no1D+i*no1D*no1D]=j*meshSpacing-halfNo*meshSpacing;
+ z[k+j*no1D+i*no1D*no1D]=i*meshSpacing-halfNo*meshSpacing;
+ concentration[k+j*no1D+i*no1D*no1D]=0;
+ accumulator[k+j*no1D+i*no1D*no1D]=0;
+ }
+ }
+ }
+ concentration[(noParticles-1)/2]=1.0/meshSpacing/meshSpacing/meshSpacing;//initial peak at the center
+
+ cout << noParticles<<"," << elapsed.count() <<",";
+ start = std::chrono::high_resolution_clock::now();//time mesurment start time
+ //caclulation
+ while(t<tend){
+ for(int i=0;i<noParticles;i++){
+ for(int j=0;j<noParticles;j++){
+ double d2= (x[i]-x[j])*(x[i]-x[j])
+ +(y[i]-y[j])*(y[i]-y[j])
+ +(z[i]-z[j])*(z[i]-z[j]);
+ if(d2<=rc*rc){
+ double factor=d2/epsilon/epsilon;
+ double factor5=factor*factor*factor*factor*factor;
+ accumulator[i]+= (concentration[j]-concentration[i]) / (factor5+1.0);
+ }
+ }
+ }
+ for(int i=0;i<noParticles;i++){
+ double fa=meshSpacing/epsilon;
+ double fa3=fa*fa*fa;
+ double fb=epsilon*M_PI*epsilon*M_PI;
+ concentration[i]+= dt* diffusionConst * 15.0 * fa3/fb * accumulator[i];//Euler time step
+ accumulator[i]=0;//set to 0 to have a fresh accumulator for the next time step
+ }
+ t+=dt;
+ }
+
+ finish = std::chrono::high_resolution_clock::now();
+ elapsed = finish - start;
+
+ cout << elapsed.count() << ",";
+
+
+ //for correctness studies
+ double difference=0;
+ for (int i=0; i < noParticles; i++){
+ difference+=abs(mesh[i].concentration-concentration[i]);
+ }
+ cout << difference << endl;
+ }
}
#endif /* Diffusion_PSE_Euler_unit_test_h */
\ No newline at end of file
diff --git a/unit_tests/PSE_1D_convergence_unit_test.hpp b/unit_tests/PSE_1D_convergence_unit_test.hpp
index 1e936af9147959062f35bed9dc823d1fc27317ee..381c8852219030a02478bf0b07f96325a95c2f30 100644
--- a/unit_tests/PSE_1D_convergence_unit_test.hpp
+++ b/unit_tests/PSE_1D_convergence_unit_test.hpp
@@ -16,13 +16,17 @@ namespace PSE_1D {
void pse_convergence_unit_test(int scaling=100){
GlobalVariable g;
-
+ //Settings for the convergence study in the aper
double domainSize=4.0;
int no_neigbors=12;
g.meshSpacing=domainSize/scaling;
+
+
g.epsilon=g.meshSpacing;
g.rc=g.meshSpacing*no_neigbors;
+
+
int length=scaling+2*no_neigbors+1;
g.meshDim = PM_Point<int,1> (length);
@@ -65,6 +69,88 @@ namespace PSE_1D {
// cout << s.str() << endl;
save_to_file <Particle, 1>(errorMesh,s.str());
}
+
+
+
+ double gauss(double x){
+ double sigma=0.05;
+ return exp(-(x*x)/(sigma*sigma))/(sqrt(M_PI)*sigma);
+ }
+ double gaussdx(double x){
+ double sigma=0.05;
+ return -2.0*x*exp(-x*x/sigma/sigma)/sqrt(M_PI)/power(sigma,3);
+ }
+ double gaussdx2(double x){
+ double sigma=0.05;
+ //return 2.0/sqrt(M_PI)/power(sigma,3)*(2*x*x/sigma/sigma-1)*exp(-x*x/sigma/sigma);
+ return exp(-(x*x)/(sigma*sigma))*(4.0*x*x-2.0*sigma*sigma)/(sqrt(M_PI)*power(sigma,5));
+ }
+ void pse_convergence_unit_test_gauss(int scaling=100){
+ GlobalVariable g;
+ //Settings for the convergence study in the aper
+ double domainSize=1.0;
+ double no_neigbors=9;
+
+ g.meshSpacing=domainSize/scaling;
+
+ double c=1.0;
+ g.epsilon=g.meshSpacing/c;
+ g.rc=g.epsilon*no_neigbors;
+
+
+
+ int length=scaling+2*no_neigbors+1;
+ g.meshDim = PM_Point<int,1> (length);
+ g.domain_min= PM_Point<double,1>(-domainSize/2.0-g.rc);
+ g.domain_max= PM_Point<double,1>(domainSize/2.0+g.rc);
+
+ vector<Particle> mesh;
+ mesh.resize(length);
+
+ //init
+ for (int i=0; i<length; i++){
+ mesh[i].position[0]=i*g.meshSpacing+g.domain_min[0];
+ mesh[i].concentration= gauss(mesh[i].position[0]);
+ //mesh[i].interactionCOunt=0;
+ //cout <<mesh[i].position[0]<<","<< mesh[i].concentration <<","<< gauss(mesh[i].position[0]) <<endl;
+ }
+ //cout <<endl <<endl;
+
+
+ Particle_Method diff(g, mesh);
+ auto start = std::chrono::high_resolution_clock::now();//time mesurment start time
+ diff.run(1000);
+ auto finish = std::chrono::high_resolution_clock::now();
+ auto runtime = finish - start;
+
+
+ //for convergence studies
+ double maxError=0.0;
+ double l2error=0.0;
+ double maxValue=0.0;
+ vector<Particle> errorMesh;
+ errorMesh.resize(length);
+ for (int i=no_neigbors; i < length-no_neigbors; i++){
+ errorMesh[i].position=mesh[i].position;
+ double error=mesh[i].concentration-gaussdx2(mesh[i].position[0]);
+ errorMesh[i].concentration=abs(error);
+ l2error+=error*error;
+ maxValue=max(maxValue,abs(gaussdx2(mesh[i].position[0])));
+ maxError=max(maxError,errorMesh[i].concentration);
+ //cout <<mesh[i].position[0]<<","<< mesh[i].concentration <<","<< gaussdx2(mesh[i].position[0]) <<endl;
+ }
+ //l2error=sqrt(l2error)/maxValue/(scaling+1);//sqrt(l2error)/maxValue/(scaling+1);
+ l2error=l2error/maxValue/maxValue/(scaling+1);
+ maxError=maxError/maxValue;
+
+ cout << g.meshSpacing<<","<< maxError <<","<< runtime.count() << endl;
+ std::stringstream s;
+ s.str("");
+ s << "PSE_error_"<< ".csv";
+ // cout << s.str() << endl;
+ save_to_file <Particle, 1>(errorMesh,s.str());
+ }
+
}
#endif /* Diffusion_PSE_Euler_unit_test_h */
\ No newline at end of file
diff --git a/unit_tests/Perfectly_Elastic_Collision_nD_correctness_check.hpp b/unit_tests/Perfectly_Elastic_Collision_nD_correctness_check.hpp
deleted file mode 100644
index e4705d610f587b1508249d4c3dc2d2864f9201b8..0000000000000000000000000000000000000000
--- a/unit_tests/Perfectly_Elastic_Collision_nD_correctness_check.hpp
+++ /dev/null
@@ -1,124 +0,0 @@
-//
-// Perfectly_Elastic_Collision_nD_correctness_check.hpp
-// prototype_interface
-//
-// Created by Johannes Pahlke on 15.06.2022.
-//
-
-#ifndef Perfectly_Elastic_Collision_nD_correctness_check_h
-#define Perfectly_Elastic_Collision_nD_correctness_check_h
-
-#include "Perfectly_Elastic_Collision_nD.hpp"
-#include <type_traits>
-#include <random>
-
-//Particle Method Instances for perfectly elastic collision
-namespace Perfectly_Elastic_Collision{
-
- void perfectly_elastic_collision_2D_correctness_check(){
- std::random_device rd;
- std::mt19937 e2(0);
- std::uniform_real_distribution<> rand(-0.5, 0.5);
-
- //invalid initial values (the particle are too close)
- // int size=10000;
- // vector<Particle<double,2>> balls;
- // balls.resize(size);
- // for (int i=0; i<size; i++){
- // balls[i].position[0]=rand(e2);
- // balls[i].position[1]=rand(e2);
- // balls[i].velocity[0]=rand(e2)/10;
- // balls[i].velocity[1]=rand(e2)/10;
- // }
- // GlobalVariable<double> g;
- // g.dt=0.1;
- // g.endT=10;
- // g.t=0.0;
- // g.rc=0.01;
- int size=3;
- vector<Particle<double,2>> balls;
- balls.resize(3);
- balls[0].position[0]=0.0;
- balls[0].position[1]=0.0;
- balls[1].position[0]=1.0;
- balls[1].position[1]=0.2;
- balls[2].position[0]=-0.1;
- balls[2].position[1]=0.8;
-
- balls[0].velocity[0]=0.25;
- balls[0].velocity[1]=0;
- balls[1].velocity[0]=-0.50;
- balls[1].velocity[1]=0;
- balls[2].velocity[0]=0.20;
- balls[2].velocity[1]=-0.70;
-
- GlobalVariable<double> g;
- g.dt=0.1;
- g.endT=10;
- g.t=0.0;
- g.rc=0.5;
-
- Particle_Method<double, 2> pec(g,balls);
-
-
-// Init native C++ variante
- double t=0;
- double dt=0.1;
- double tend=10;
- double rc=0.01;
-
- double x[size];
- double y[size];
- double vx[size];
- double vy[size];
- for(int i=0;i<size;i++){
- x[i]=balls[i].position[0];
- y[i]=balls[i].position[1];
- vx[i]=balls[i].velocity[0];
- vy[i]=balls[i].velocity[1];
- }
-
-// run framework variant
- pec.run();
-
-//run native C++ variante
- while(t<tend){
- for(int i=0;i<size;i++){
- for(int j=0;j<size;j++){
- // double d2= (x[i]-x[j])*(x[i]-x[j])
- // +(y[i]-y[j])*(y[i]-y[j]);
- double diffx=x[j]-x[i];
- double diffy=y[j]-y[i];
- double d2=diffx*diffx+diffy*diffy;
- if(d2<=rc*rc && (diffx>0.0 || (diffx==0.0 && diffy>0.0)) ){
- double diffx=x[j]-x[i];
- double diffy=y[j]-y[i];
- double diffvx=vx[j]-vx[i];
- double diffvy=vy[j]-vy[i];
- double temp=diffx*diffvx+diffy*diffvy;
-
- vx[i]+=diffx/d2*temp;
- vy[i]+=diffy/d2*temp;
- vx[j]-=diffx/d2*temp;
- vy[j]-=diffy/d2*temp;
- }
- }
- }
- for(int i=0;i<size;i++){
- x[i]+=dt*vx[i];
- y[i]+=dt*vy[i];
- }
- t+=dt;
- }
-
-//calculate the difference of the two versions
- double difference=0;
- for (int i=0; i < size; i++){
- difference+=abs(balls[i].position[0]-x[i])+abs(balls[i].position[1]-y[i]);
- }
- cout << "difference: " << difference << endl;
-
- }
-}
-
-#endif /* Perfectly_Elastic_Collision_nD_correctness_check_h */
\ No newline at end of file
diff --git a/unit_tests/Perfectly_Elastic_Collision_nD_unit_test.hpp b/unit_tests/Perfectly_Elastic_Collision_nD_unit_test.hpp
index c7d3a5736f7261b892cc1eb3450f1aa3538cecd1..bb8165e6a6a4333b8b52ff3e2b8858886a39b7f6 100644
--- a/unit_tests/Perfectly_Elastic_Collision_nD_unit_test.hpp
+++ b/unit_tests/Perfectly_Elastic_Collision_nD_unit_test.hpp
@@ -177,7 +177,204 @@ namespace Perfectly_Elastic_Collision{
int cellNo=power(domainSize/g.rc,dimension);
cout << cellNo << "," << particleNo << "," << elapsed.count() << endl;
}
+
+
+ //Instance with not valid initial values (the particle are too close)
+ void perfectly_elastic_collision_2D_big_time_test_02(double domainSize){
+ const int dimension=2;
+
+ int paNoPerDimAndCell=3;
+
+ GlobalVariable<double> g;
+ g.dt=0.1;
+ g.endT=10;
+ g.t=0.0;
+ g.rc=0.5;
+
+ double particleSPacing=g.rc/paNoPerDimAndCell;
+
+ //particleNo=10*power(domainSize,dimension);//this keeps the particle density on average the same low density
+ int particleNoPerDim=paNoPerDimAndCell*domainSize/g.rc;
+ int particleNo=power(particleNoPerDim,dimension);
+
+
+ std::random_device rd;
+ std::mt19937 e2(0);
+ std::uniform_real_distribution<> rand(-0.5, 0.5);
+
+
+ vector<Particle<double,dimension>> balls;
+ balls.resize(particleNo);
+ for (int i=0; i<particleNo; i++){
+ balls[i].position=(PM_Point<double, dimension>)toVectorialIndex(i, PM_Point<int,dimension> (particleNoPerDim))*particleSPacing;
+ for(int j=0; j<dimension; j++){
+ // balls[i].position[j]=rand(e2)/;
+ // balls[i].position[j]=
+ balls[i].velocity[j]=rand(e2)/10.0;
+ }
+ }
+
+
+
+
+
+ // Init native C++ variante
+ double t=g.t;
+ double dt=g.dt;
+ double tend=g.endT;
+ double rc=g.rc;
+
+ double x[particleNo];
+ double y[particleNo];
+ double vx[particleNo];
+ double vy[particleNo];
+ for(int i=0;i<particleNo;i++){
+ x[i]=balls[i].position[0];
+ y[i]=balls[i].position[1];
+ vx[i]=balls[i].velocity[0];
+ vy[i]=balls[i].velocity[1];
+ }
+
+ // Time mesurment // run framework variant
+ Particle_Method<double, dimension> pec(g,balls);
+ auto start = std::chrono::high_resolution_clock::now();//time mesurment start time
+ pec.run();
+ auto finish = std::chrono::high_resolution_clock::now();
+ std::chrono::duration<double> elapsed = finish - start;
+
+ int cellNo=power(domainSize/g.rc,dimension);
+ cout << cellNo << "," << particleNo << "," << elapsed.count() << ",";
+
+
+ start = std::chrono::high_resolution_clock::now();//time mesurment start time
+//run native C++ variante
+ while(t<tend){
+ for(int i=0;i<particleNo;i++){
+ for(int j=0;j<particleNo;j++){
+ // double d2= (x[i]-x[j])*(x[i]-x[j])
+ // +(y[i]-y[j])*(y[i]-y[j]);
+ double diffx=x[j]-x[i];
+ double diffy=y[j]-y[i];
+ double d2=diffx*diffx+diffy*diffy;
+ if(d2<=rc*rc && (diffx>0.0 || (diffx==0.0 && diffy>0.0)) ){
+ double diffx=x[j]-x[i];
+ double diffy=y[j]-y[i];
+ double diffvx=vx[j]-vx[i];
+ double diffvy=vy[j]-vy[i];
+ double temp=diffx*diffvx+diffy*diffvy;
+
+ vx[i]+=diffx/d2*temp;
+ vy[i]+=diffy/d2*temp;
+ vx[j]-=diffx/d2*temp;
+ vy[j]-=diffy/d2*temp;
+ }
+ }
+ }
+ for(int i=0;i<particleNo;i++){
+ x[i]+=dt*vx[i];
+ y[i]+=dt*vy[i];
+ }
+ t+=dt;
+ }
+ finish = std::chrono::high_resolution_clock::now();
+ elapsed = finish - start;
+
+
+ cout << elapsed.count() << endl;
+ }
+
+
+ void perfectly_elastic_collision_2D_correctness_check(){
+ std::random_device rd;
+ std::mt19937 e2(0);
+ std::uniform_real_distribution<> rand(-0.5, 0.5);
+
+ int size=3;
+ vector<Particle<double,2>> balls;
+ balls.resize(3);
+ balls[0].position[0]=0.0;
+ balls[0].position[1]=0.0;
+ balls[1].position[0]=1.0;
+ balls[1].position[1]=0.2;
+ balls[2].position[0]=-0.1;
+ balls[2].position[1]=0.8;
+
+ balls[0].velocity[0]=0.25;
+ balls[0].velocity[1]=0;
+ balls[1].velocity[0]=-0.50;
+ balls[1].velocity[1]=0;
+ balls[2].velocity[0]=0.20;
+ balls[2].velocity[1]=-0.70;
+
+ GlobalVariable<double> g;
+ g.dt=0.1;
+ g.endT=10;
+ g.t=0.0;
+ g.rc=0.5;
+
+ Particle_Method<double, 2> pec(g,balls);
+
+
+// Init native C++ variante
+ double t=g.t;
+ double dt=g.dt;
+ double tend=g.endT;
+ double rc=g.rc;
+
+ double x[size];
+ double y[size];
+ double vx[size];
+ double vy[size];
+ for(int i=0;i<size;i++){
+ x[i]=balls[i].position[0];
+ y[i]=balls[i].position[1];
+ vx[i]=balls[i].velocity[0];
+ vy[i]=balls[i].velocity[1];
+ }
+
+// run framework variant
+ pec.run();
+
+//run native C++ variante
+ while(t<tend){
+ for(int i=0;i<size;i++){
+ for(int j=0;j<size;j++){
+ double diffx=x[j]-x[i];
+ double diffy=y[j]-y[i];
+
+ double d2=diffx*diffx+diffy*diffy;
+ if(d2<=rc*rc && (diffx>0.0 || (diffx==0.0 && diffy>0.0)) ){
+ double diffvx=vx[j]-vx[i];
+ double diffvy=vy[j]-vy[i];
+ double temp=diffx*diffvx+diffy*diffvy;
+
+ vx[i]+=diffx/d2*temp;
+ vy[i]+=diffy/d2*temp;
+ vx[j]-=diffx/d2*temp;
+ vy[j]-=diffy/d2*temp;
+ }
+ }
+ }
+ for(int i=0;i<size;i++){
+ x[i]+=dt*vx[i];
+ y[i]+=dt*vy[i];
+ }
+ t+=dt;
+ }
+
+//calculate the difference of the two versions
+ double difference=0;
+ for (int i=0; i < size; i++){
+ difference+=abs(balls[i].position[0]-x[i])+abs(balls[i].position[1]-y[i]);
+ }
+ cout << "difference: " << difference << endl;
+
+ }
+
+
}
+
+
#endif /* Perfectly_Elastic_Collision_nD_unit_test_h */