clean up

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;

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

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;

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

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

unit_tests/Perfectly_Elastic_Collision_nD_correctness_check.hpp

-//
-//  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

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 */