bubble-dynamics

commit d8fbc4a5754fbdae871da60c87a734eb662b7d02
parent 0ad2649d61478ea0c045b2df3888c402a51705b7
Author: Ursula Rasthofer <urasthofer@ethz.ch>
Date:   Mon, 20 Feb 2017 15:06:18 +0100

added comments

Diffstat:
A
rebuild.sh
 | 
1
+
M
src/kernels/KMClusterPositions_D.h
 | 
23
++++++++++++++++++-----
M
src/kernels/KMCluster_FC.h
 | 
12
+++++++-----
M
src/kernels/KMCluster_TY.h
 | 
8
++++----
M
src/kernels/RPCluster.h
 | 
8
++++----
M
src/kernels/RPClusterPositions_D.h
 | 
8
++++----

6 files changed, 38 insertions(+), 22 deletions(-)
diff --git a/rebuild.sh b/rebuild.sh
@@ -0,0 +1 @@
+make clean; make eigen=1 pos=1
diff --git a/src/kernels/KMClusterPositions_D.h b/src/kernels/KMClusterPositions_D.h
@@ -9,14 +9,15 @@
 
 #include <cassert>
 #include <cstdlib>
+#include <iostream>
 #include <fstream>
 #include <cmath>
 #include <vector>
 
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
-//#include <Eigen/Core>
-//#include <Eigen/Dense>
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
+#include <Eigen/Core>
+#include <Eigen/Dense>
 
 #include "KernelBase.h"
 #include "BubbleBase.h"
@@ -134,7 +135,7 @@ public:
 
                     if (Dij < (Ui[0]+Uj[0]))
                     {
-                      std::cerr << "Colliding bubbles " << Dij << " " << Ui[0] << " " << Uj[0] << std::endl;
+                      std::cerr << "Colliding bubbles " << i << " " << j << " distance " << Dij << " radius i " << Ui[0] << " radius j " << Uj[0] << std::endl;
                       abort();
                     }
 
@@ -266,6 +267,18 @@ public:
         //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.fullPivHouseholderQr().solve(b);
         // this is notably faster and sufficiently accurate (rasthofer Feb 3, 2017)
         Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.colPivHouseholderQr().solve(b);
+        // don't use these solvers although they are fast (err e-2)
+        // to increase speed use less expensive time integration scheme
+        //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.llt().solve(b);
+        //Eigen::Matrix<Real, Eigen::Dynamic, 1> Rddot = A.ldlt().solve(b);
+#ifndef NDEBUG
+        Real err = (A*Rddot - b).norm() / b.norm();
+
+        ofstream errfile;
+        errfile.open("error.dat", ios::out | ios::app);
+        errfile << scientific << setprecision(6) << err << endl;
+        errfile.close();
+#endif /* NDEBUG */
 
         for (size_t i = 0; i < _N; ++i)
         {
diff --git a/src/kernels/KMCluster_FC.h b/src/kernels/KMCluster_FC.h
@@ -12,10 +12,10 @@
 #include <cmath>
 #include <vector>
 
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
-//#include <Eigen/Core>
-//#include <Eigen/Dense>
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
+#include <Eigen/Core>
+#include <Eigen/Dense>
 
 #include "KernelBase.h"
 #include "BubbleBase.h"
@@ -109,7 +109,9 @@ public:
                 else
                 {
                     const Real UDij = Uj[0]/bd.Dij[j + i*_N];
-                    A(i,j) = (static_cast<Real>(1) + Ui[1]*clInv)*Uj[0]*UDij + static_cast<Real>(6.0)*Uj[1]*UDij*Ui[0]*clInv;
+                    A(i,j) = (static_cast<Real>(1) + Ui[1]*clInv)*Uj[0]*UDij + static_cast<Real>(1.0)*Uj[1]*UDij*Ui[0]*clInv;
+                    // extension of Matrix: keeping all terms would yield a factor of 6 for the second term, 
+                    // following simplifications of Fuster and Colonius results in a factor of 1
                     bnbr += Uj[1]*Uj[1]*UDij;
                     bnbr2 += Uj[1]*Uj[1]*Uj[1]/bd.Dij[j + i*_N];
                 }
diff --git a/src/kernels/KMCluster_TY.h b/src/kernels/KMCluster_TY.h
@@ -12,10 +12,10 @@
 #include <cmath>
 #include <vector>
 
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
-//#include <Eigen/Core>
-//#include <Eigen/Dense>
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
+#include <Eigen/Core>
+#include <Eigen/Dense>
 
 #include "KernelBase.h"
 #include "BubbleBase.h"
diff --git a/src/kernels/RPCluster.h b/src/kernels/RPCluster.h
@@ -12,10 +12,10 @@
 #include <cmath>
 #include <vector>
 
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
-//#include <Eigen/Core>
-//#include <Eigen/Dense>
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
+#include <Eigen/Core>
+#include <Eigen/Dense>
 
 #include "KernelBase.h"
 #include "BubbleBase.h"
diff --git a/src/kernels/RPClusterPositions_D.h b/src/kernels/RPClusterPositions_D.h
@@ -13,10 +13,10 @@
 #include <cmath>
 #include <vector>
 
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
-#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
-//#include <Eigen/Core>
-//#include <Eigen/Dense>
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Core"
+//#include "/cluster/home/ursular/research/codes/eigen/Eigen/Dense"
+#include <Eigen/Core>
+#include <Eigen/Dense>
 
 #include "KernelBase.h"
 #include "BubbleBase.h"