merge

Author: dmarce1

CMakeLists.txt

@@ -1,9 +1,10 @@
-
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -O3 -march=native")
+set(CMAKE_CXX_COMPILER "g++")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -O3 -march=native -fopenmp")
 
 add_executable( fmm_test
 	src/main.cpp
 	src/multipole.cpp
 	src/expansion.cpp
+	src/lane_emden.cpp
 )
 

src/defs.hpp

@@ -23,15 +23,15 @@
 
 #include <array>
 
-#define CORRECTION_ON
+//#define CORRECTION_ON
 //#define CORRECTION_OPTIMIZE
 
 using real = double;
 using integer = long long;
 
 constexpr integer NDIM = 3;
 constexpr integer ncrit = 8;
-constexpr integer nparts = 100000;
+constexpr integer nparts = 20000;
 constexpr integer nchild = 8;
 constexpr integer LP = 20;
 

src/expansion.cpp

@@ -1,22 +1,19 @@
-
-
 /*  
-    Copyright (c) 2016 Dominic C. Marcello
-
-    This program is free software: you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
+ Copyright (c) 2016 Dominic C. Marcello
 
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
 
-    You should have received a copy of the GNU General Public License
-    along with this program.  If not, see <http://www.gnu.org/licenses/>.
-*/
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
 
+ You should have received a copy of the GNU General Public License
+ along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
 
 #include "expansion.hpp"
 #include <cmath>
@@ -166,10 +163,11 @@ void expansion::invert() {
 expansion::~expansion() {
 }
 
-void multipole_interaction(expansion& L1, expansion& L2, const multipole& M1, const multipole& M2, space_vector dX) {
+void multipole_interaction(expansion& L1, expansion& L2, const multipole& M1,
+		const multipole& M2, space_vector dX) {
 
-	const real delta[3][3] = { { real(1.0), real(0.0), real(0.0) }, { real(0.0), real(1.0), real(0.0) }, { real(0.0),
-			real(0.0), real(1.0) } };
+	const real delta[3][3] = { { real(1.0), real(0.0), real(0.0) }, { real(0.0),
+			real(1.0), real(0.0) }, { real(0.0), real(0.0), real(1.0) } };
 
 	real y0 = 0.0;
 	for (integer d = 0; d != NDIM; ++d) {
@@ -187,16 +185,28 @@ void multipole_interaction(expansion& L1, expansion& L2, const multipole& M1, co
 #endif
 #endif
 	D() = d0;
+
+	for (integer a = 0; a < 3; a++) {
+		for (integer b = a; b < 3; b++) {
+			for (integer c = b; c < 3; c++) {
+				D(a, b, c) = 0.0;
+			}
+		}
+	}
+
 	for (integer a = 0; a < 3; a++) {
 		D(a) = dX[a] * d1;
+		D(a, a, a) += dX[a] * d2;
 		for (integer b = a; b < 3; b++) {
+			D(a, a, b) += dX[b] * d2;
+			D(b, b, a) += dX[a] * d2;
 			D(a, b) = dX[a] * dX[b] * d2 + delta[a][b] * d1;
 			for (integer c = b; c < 3; c++) {
-				D(a, b, c) = dX[a] * dX[b] * dX[c] * d3;
-				D(a, b, c) += (delta[a][b] * dX[c] + delta[b][c] * dX[a] + delta[c][a] * dX[b]) * d2;
+				D(a, b, c) += dX[a] * dX[b] * dX[c] * d3;
 			}
 		}
 	}
+
 #ifdef CORRECTION_ON
 	std::array<expansion, NDIM> D4;
 	for (integer i = 0; i != NDIM; ++i) {
@@ -226,44 +236,27 @@ void multipole_interaction(expansion& L1, expansion& L2, const multipole& M1, co
 
 	L1() += M2() * D();
 	L2() += M1() * D();
-	for (integer a = 0; a < 3; a++) {
-		for (integer b = 0; b < 3; b++) {
-			L1() += M2(a, b) * D(a, b) * 0.5;
-			L2() += M1(a, b) * D(a, b) * 0.5;
-		}
-	}
-
 	for (integer a = 0; a < 3; a++) {
 		L1(a) += M2() * D(a);
 		L2(a) -= M1() * D(a);
-	}
-	for (integer a = 0; a < 3; a++) {
 		for (integer b = 0; b < 3; b++) {
+			L1() += M2(a, b) * D(a, b) * 0.5;
+			L2() += M1(a, b) * D(a, b) * 0.5;
 			for (integer c = 0; c < 3; c++) {
 				L1(a) += M2(c, b) * D(a, b, c) * 0.5;
 				L2(a) -= M1(c, b) * D(a, b, c) * 0.5;
 			}
 		}
-	}
-
-	for (integer a = 0; a < 3; a++) {
 		for (integer b = a; b < 3; b++) {
 			L1(a, b) += M2() * D(a, b);
 			L2(a, b) += M1() * D(a, b);
-		}
-	}
-
-	for (integer a = 0; a < 3; a++) {
-		for (integer b = a; b < 3; b++) {
 			for (integer c = b; c < 3; c++) {
 				L1(a, b, c) += M2() * D(a, b, c);
 				L2(a, b, c) -= M1() * D(a, b, c);
 			}
 		}
 	}
 
-
-
 #ifdef CORRECTION_ON
 	for (integer i = 0; i != NDIM; ++i) {
 		for (integer j = 0; j != NDIM; ++j) {
@@ -277,64 +270,3 @@ void multipole_interaction(expansion& L1, expansion& L2, const multipole& M1, co
 	}
 #endif
 }
-
-/*
- std::array<expansion, NDIM> expansion::get_derivatives() const {
- std::array<expansion, NDIM> D;
- for (integer i = 0; i != NDIM; ++i) {
- D[i]() = (*this)(i);
- for (integer a = 0; a < NDIM; ++a) {
- D[i](a) = (*this)(i, a);
- for (integer b = a; b < NDIM; ++b) {
- D[i](a, b) = (*this)(i, a, b);
- for (integer c = b; c != NDIM; ++c) {
- D[i](a, b, c) = real(0.0);
- }
- }
- }
- }
- #ifdef CORRECTION_ON
- const real rinv = real(1) / r;
- const real D2 = -real(3.0) * std::pow(rinv, real(5));
- const real D3 = +real(15.) * std::pow(rinv, real(7));
- const real D4 = -real(105) * std::pow(rinv, real(9));
- std::array<expansion, NDIM> E1;
- for (integer i = 0; i != NDIM; ++i) {
- E1[i] = 0.0;
- for (integer j = 0; j != NDIM; ++j) {
- for (integer k = j; k != NDIM; ++k) {
- for (integer l = k; l != NDIM; ++l) {
- E1[i](j, k, l) += delta[i][j] * delta[k][l] * D2;
- E1[i](j, k, l) += delta[i][l] * delta[j][k] * D2;
- E1[i](j, k, l) += delta[i][k] * delta[l][j] * D2;
-
- E1[i](j, k, l) += delta[i][j] * R[k] * R[l] * D3;
- E1[i](j, k, l) += delta[i][l] * R[j] * R[k] * D3;
- E1[i](j, k, l) += delta[i][k] * R[l] * R[j] * D3;
-
- E1[i](j, k, l) += R[i] * R[j] * delta[k][l] * D3;
- E1[i](j, k, l) += R[i] * R[l] * delta[j][k] * D3;
- E1[i](j, k, l) += R[i] * R[k] * delta[l][j] * D3;
-
- E1[i](j, k, l) += R[i] * R[j] * R[k] * R[l] * D4;
-
-
- }
- }
- }
- }
- for (integer i = 0; i != NDIM; ++i) {
- for (integer j = 0; j != NDIM; ++j) {
- for (integer k = 0; k != NDIM; ++k) {
- for (integer l = 0; l != NDIM; ++l) {
- D[i]() -= E1[i](j, k, l) * M(j, k, l) * (real(1) / real(6));
- D[i]() += E1[i](j, k, l) * N(j, k, l) * M() / N() * (real(1) / real(6));
- }
- }
- }
- }
- #endif
- return D;
- }
- */
-