Commit 1865c4ec authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Extended the gravityDerivatives test to also verify the M2P terms.

parent fe44f590
......@@ -935,13 +935,14 @@ int main(int argc, char* argv[]) {
message("Seed = %d", seed);
srand(seed);
/* Start by testing M2L */
for (int i = 0; i < 100; ++i) {
const double dx = 100. * ((double)rand() / (RAND_MAX));
const double dy = 100. * ((double)rand() / (RAND_MAX));
const double dz = 100. * ((double)rand() / (RAND_MAX));
message("Testing gravity for r=(%e %e %e)", dx, dy, dz);
message("Testing M2L gravity for r=(%e %e %e)", dx, dy, dz);
const double r_s = 100. * ((double)rand() / (RAND_MAX));
const double r_s_inv = 1. / r_s;
......@@ -967,89 +968,209 @@ int main(int argc, char* argv[]) {
/* Now check everything... */
/* 0th order terms */
test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "D_000");
test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "M2L D_000");
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 1st order terms */
test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "D_100");
test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "D_010");
test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "D_001");
test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "M2L D_100");
test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "M2L D_010");
test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "M2L D_001");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order terms */
test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "D_200");
test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "D_020");
test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "D_002");
test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "D_110");
test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "D_101");
test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "D_011");
test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "M2L D_200");
test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "M2L D_020");
test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "M2L D_002");
test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "M2L D_110");
test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "M2L D_101");
test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "M2L D_011");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
tol *= 2.5;
/* 3rd order terms */
test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "D_300");
test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "D_030");
test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "D_003");
test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "D_210");
test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "D_201");
test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "D_120");
test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "D_021");
test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "D_102");
test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "D_012");
test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "D_111");
test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "M2L D_300");
test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "M2L D_030");
test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "M2L D_003");
test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "M2L D_210");
test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "M2L D_201");
test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "M2L D_120");
test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "M2L D_021");
test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "M2L D_102");
test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "M2L D_012");
test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "M2L D_111");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order terms */
test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "D_400");
test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "D_040");
test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "D_004");
test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "D_310");
test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "D_301");
test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "D_130");
test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "D_031");
test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "D_103");
test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "D_013");
test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "D_220");
test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "D_202");
test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "D_022");
test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "D_211");
test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "D_121");
test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "D_112");
test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "M2L D_400");
test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "M2L D_040");
test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "M2L D_004");
test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "M2L D_310");
test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "M2L D_301");
test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "M2L D_130");
test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "M2L D_031");
test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "M2L D_103");
test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "M2L D_013");
test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "M2L D_220");
test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "M2L D_202");
test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "M2L D_022");
test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "M2L D_211");
test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "M2L D_121");
test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "M2L D_112");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
tol *= 2.5;
/* 5th order terms */
test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "D_500");
test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "D_050");
test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "D_005");
test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "D_410");
test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "D_401");
test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "D_140");
test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "D_041");
test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "D_104");
test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "D_014");
test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "D_320");
test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "D_302");
test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "D_230");
test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "D_032");
test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "D_203");
test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "D_023");
test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "D_311");
test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "D_131");
test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "D_113");
test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "D_122");
test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "D_212");
test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "D_221");
test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "M2L D_500");
test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "M2L D_050");
test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "M2L D_005");
test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "M2L D_410");
test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "M2L D_401");
test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "M2L D_140");
test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "M2L D_041");
test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "M2L D_104");
test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "M2L D_014");
test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "M2L D_320");
test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "M2L D_302");
test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "M2L D_230");
test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "M2L D_032");
test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "M2L D_203");
test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "M2L D_023");
test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "M2L D_311");
test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "M2L D_131");
test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "M2L D_113");
test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "M2L D_122");
test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "M2L D_212");
test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "M2L D_221");
#endif
message("All good!");
}
/* And now the M2P terms */
for (int i = 0; i < 100; ++i) {
const double dx = 100. * ((double)rand() / (RAND_MAX));
const double dy = 100. * ((double)rand() / (RAND_MAX));
const double dz = 100. * ((double)rand() / (RAND_MAX));
message("Testing M2P gravity for r=(%e %e %e)", dx, dy, dz);
const double r_s = 100. * ((double)rand() / (RAND_MAX));
const double r_s_inv = 1. / r_s;
const int periodic = 0;
message("Mesh scale r_s=%e periodic=%d", r_s, periodic);
/* Compute distance */
const double r2 = dx * dx + dy * dy + dz * dz;
const double r_inv = 1. / sqrt(r2);
const double r = r2 * r_inv;
const double eps = r / 10.;
/* Compute all derivatives */
struct potential_derivatives_M2P pot;
potential_derivatives_compute_M2P(dx, dy, dz, r2, r_inv, eps, periodic,
r_s_inv, &pot);
/* Minimal value we care about */
const double min = 1e-9;
/* Now check everything... */
/* 0th order terms */
test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "M2P D_000");
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 1st order terms */
test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "M2P D_100");
test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "M2P D_010");
test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "M2P D_001");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order terms */
test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "M2P D_200");
test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "M2P D_020");
test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "M2P D_002");
test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "M2P D_110");
test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "M2P D_101");
test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "M2P D_011");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
tol *= 2.5;
/* 3rd order terms */
test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "M2P D_300");
test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "M2P D_030");
test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "M2P D_003");
test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "M2P D_210");
test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "M2P D_201");
test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "M2P D_120");
test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "M2P D_021");
test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "M2P D_102");
test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "M2P D_012");
test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "M2P D_111");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order terms */
test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "M2P D_400");
test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "M2P D_040");
test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "M2P D_004");
test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "M2P D_310");
test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "M2P D_301");
test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "M2P D_130");
test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "M2P D_031");
test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "M2P D_103");
test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "M2P D_013");
test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "M2P D_220");
test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "M2P D_202");
test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "M2P D_022");
test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "M2P D_211");
test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "M2P D_121");
test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "M2P D_112");
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
tol *= 2.5;
/* 5th order terms */
test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "M2P D_500");
test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "M2P D_050");
test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "M2P D_005");
test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "M2P D_410");
test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "M2P D_401");
test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "M2P D_140");
test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "M2P D_041");
test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "M2P D_104");
test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "M2P D_014");
test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "M2P D_320");
test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "M2P D_302");
test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "M2P D_230");
test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "M2P D_032");
test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "M2P D_203");
test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "M2P D_023");
test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "M2P D_311");
test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "M2P D_131");
test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "M2P D_113");
test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "M2P D_122");
test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "M2P D_212");
test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "M2P D_221");
#endif
message("All good!");
}
/* All happy */
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment