diff --git a/CMakeLists.txt b/CMakeLists.txt index a5088db..8d7f33e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,39 +16,36 @@ if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) enable_testing() add_executable(polesitter_tests tests/main.c) - target_link_libraries(polesitter_tests PRIVATE polesitter) - if(NOT MSVC) - target_link_libraries(polesitter_tests PRIVATE m) - endif() - - if(MSVC) - target_compile_options(polesitter_tests PRIVATE - /W4 /WX /O2 - ) - else() - target_compile_options(polesitter_tests PRIVATE - -Wall -Wextra -pedantic -O3 -ffast-math -march=native - ) - endif() - - add_test(NAME core_tests COMMAND polesitter_tests) + add_executable(polesitter_accuracy_tests tests/accuracy.c) + target_link_libraries(polesitter_accuracy_tests PRIVATE polesitter) add_executable(polesitter_benchmark tests/benchmark.c) target_link_libraries(polesitter_benchmark PRIVATE polesitter) if(NOT MSVC) + target_link_libraries(polesitter_tests PRIVATE m) + target_link_libraries(polesitter_accuracy_tests PRIVATE m) target_link_libraries(polesitter_benchmark PRIVATE m) endif() if(MSVC) - target_compile_options(polesitter_benchmark PRIVATE - /W4 /WX /O2 - ) + target_compile_options(polesitter_tests PRIVATE /W4 /WX /O2) + target_compile_options(polesitter_accuracy_tests PRIVATE /W4 /WX /O2) + target_compile_options(polesitter_benchmark PRIVATE /W4 /WX /O2) else() + target_compile_options(polesitter_tests PRIVATE + -Wall -Wextra -pedantic -O3 -ffast-math -march=native + ) + target_compile_options(polesitter_accuracy_tests PRIVATE + -Wall -Wextra -pedantic -O3 -ffast-math -march=native + ) target_compile_options(polesitter_benchmark PRIVATE -Wall -Wextra -pedantic -O3 -ffast-math -march=native ) endif() + + add_test(NAME core_tests COMMAND polesitter_tests) + add_test(NAME accuracy_tests COMMAND polesitter_accuracy_tests) endif() diff --git a/src/polesitter.h b/src/polesitter.h index a59fac4..7dca388 100644 --- a/src/polesitter.h +++ b/src/polesitter.h @@ -238,7 +238,6 @@ ps_result_t ps_prepare_particles(ps_particle_arrs_t* arrs, #endif // POLESITTER_H -#define POLESITTER_IMPLEMENTATION #ifdef POLESITTER_IMPLEMENTATION // ===================================================================== @@ -519,17 +518,19 @@ static void ps_impl_fmm_interaction_pass(ps_node_t* target, ps_node_t* src) { float dz = src->z - target->z; float dist_sq = (dx * dx) + (dy * dy) + (dz * dz); - // multipole acceptance criterion - // if dist > sum of half-widths, they are well-separated + // Conservative multipole acceptance criterion. The 12x squared + // separation factor was selected from direct-reference accuracy + // validation for this first-order expansion. // theta dictates accuracy vs speed float theta = 1.0F; float hw_sum = target->half_width + src->half_width; // well-separated - if (dist_sq > (theta * theta) * (hw_sum * hw_sum)) { - float dist = sqrtf(dist_sq); - float inv_r3 = 1.0F / (dist * dist_sq); - float inv_r5 = inv_r3 / dist_sq; + if (dist_sq > 12.0F * (theta * theta) * (hw_sum * hw_sum)) { + float soft_dist_sq = dist_sq + 0.1F; + float dist = sqrtf(soft_dist_sq); + float inv_r3 = 1.0F / (dist * soft_dist_sq); + float inv_r5 = inv_r3 / soft_dist_sq; float m0 = src->multipole[0]; float mx = src->multipole[1]; @@ -545,9 +546,9 @@ static void ps_impl_fmm_interaction_pass(ps_node_t* target, ps_node_t* src) { float m_dot_r = (mx * dx) + (my * dy) + (mz * dz); float dipole_coeff = 3.0F * m_dot_r * inv_r5; - float force_dip_x = (dx * dipole_coeff) - (mx * inv_r3); - float force_dip_y = (dy * dipole_coeff) - (my * inv_r3); - float force_dip_z = (dz * dipole_coeff) - (mz * inv_r3); + float force_dip_x = (mx * inv_r3) - (dx * dipole_coeff); + float force_dip_y = (my * inv_r3) - (dy * dipole_coeff); + float force_dip_z = (mz * inv_r3) - (dz * dipole_coeff); // accumulate into target's local expansion target->local[1] += force_m0_x + force_dip_x; @@ -724,7 +725,7 @@ static void ps_impl_fmm_p2p_pass(ps_node_t* target, ps_node_t* src, float hw_sum = target->half_width + src->half_width; // if well-separated M2L already handled it - if (dist_sq > (theta * theta) * (hw_sum * hw_sum)) { + if (dist_sq > 12.0F * (theta * theta) * (hw_sum * hw_sum)) { return; } @@ -814,7 +815,7 @@ static void ps_impl_fmm_p2p_pass(ps_node_t* target, ps_node_t* src, float32x4_t t_x_vec = vdupq_n_f32(t_x); float32x4_t t_y_vec = vdupq_n_f32(t_y); float32x4_t t_z_vec = vdupq_n_f32(t_z); - float32x4_t eps_vec = vdupq_n_f32(2.0F); + float32x4_t eps_vec = vdupq_n_f32(0.1F); // accumulators for target particles forces float32x4_t f_x_vec = vdupq_n_f32(0.0F); diff --git a/tests/accuracy.c b/tests/accuracy.c new file mode 100644 index 0000000..f4c7484 --- /dev/null +++ b/tests/accuracy.c @@ -0,0 +1,279 @@ +// clang-format off +#include +#include +#include +#include + +#define POLESITTER_IMPLEMENTATION +#include "../src/polesitter.h" + +#define ARENA_SIZE (1024ULL * 1024ULL * 32ULL) +#define SOFTENING_SQ 0.1F + +static int g_failures = 0; + +static void failf(const char* name, double actual, double expected, + double tolerance) { + (void)fprintf(stderr, + "[FAIL] %s: actual=%.9g expected=%.9g tolerance=%.9g\n", + name, actual, expected, tolerance); + ++g_failures; +} + +static void check_close(const char* name, double actual, double expected, + double tolerance) { + if (fabs(actual - expected) > tolerance) { + failf(name, actual, expected, tolerance); + } else { + (void)printf("[PASS] %s: %.9g\n", name, actual); + } +} + +static uint32_t lcg_next(uint32_t* state) { + *state = (*state * 1664525U) + 1013904223U; + return *state; +} + +static float rand_unit(uint32_t* state) { + return (float)(lcg_next(state) >> 8) / 16777215.0F; +} + +static void direct_reference(const float* x, const float* y, const float* z, + const float* mass, size_t count, double* fx, + double* fy, double* fz) { + for (size_t i = 0; i < count; ++i) { + double ax = 0.0; + double ay = 0.0; + double az = 0.0; + + for (size_t j = 0; j < count; ++j) { + if (i == j) { + continue; + } + + double dx = (double)x[j] - (double)x[i]; + double dy = (double)y[j] - (double)y[i]; + double dz = (double)z[j] - (double)z[i]; + double r2 = (dx * dx) + (dy * dy) + (dz * dz) + SOFTENING_SQ; + double inv_r = 1.0 / sqrt(r2); + double inv_r3 = inv_r * inv_r * inv_r; + double scale = (double)mass[i] * (double)mass[j] * inv_r3; + + ax += dx * scale; + ay += dy * scale; + az += dz * scale; + } + + fx[i] = ax; + fy[i] = ay; + fz[i] = az; + } +} + +static void run_dipole_sign_test(void) { + ps_node_t target; + ps_node_t source; + ps_impl_node_init(&target); + ps_impl_node_init(&source); + + target.x = 0.0F; + target.y = 0.0F; + target.z = 0.0F; + target.half_width = 1.0F; + + source.x = 100.0F; + source.y = 0.0F; + source.z = 0.0F; + source.half_width = 1.0F; + + /* One unit mass displaced +1 from the source-cell center. */ + source.multipole[0] = 1.0F; + source.multipole[1] = 1.0F; + + ps_impl_fmm_interaction_pass(&target, &source); + + /* + * First-order softened expansion for R=100 and delta=+1 is + * approximately 9.799856e-05. The large separation keeps this + * analytic check independent of the production MAC threshold. + */ + check_close("M2L dipole orientation", target.local[1], 0.00009799856, 5e-9); +} + +static void fill_uniform(float* x, float* y, float* z, float* mass, + uint32_t* id, size_t count) { + uint32_t state = 7331U; + for (size_t i = 0; i < count; ++i) { + x[i] = (rand_unit(&state) * 100.0F) - 50.0F; + y[i] = (rand_unit(&state) * 100.0F) - 50.0F; + z[i] = (rand_unit(&state) * 100.0F) - 50.0F; + mass[i] = 0.5F + (rand_unit(&state) * 1.5F); + id[i] = (uint32_t)i; + } +} + +static void fill_two_clusters(float* x, float* y, float* z, float* mass, + uint32_t* id, size_t count) { + uint32_t state = 20260803U; + size_t half = count / 2; + + for (size_t i = 0; i < count; ++i) { + float center = i < half ? -12.0F : 12.0F; + float skew = i < half ? -0.75F : 0.75F; + x[i] = center + skew + ((rand_unit(&state) - 0.5F) * 2.0F); + y[i] = (rand_unit(&state) - 0.5F) * 2.0F; + z[i] = (rand_unit(&state) - 0.5F) * 2.0F; + mass[i] = 0.5F + (rand_unit(&state) * 2.0F); + id[i] = (uint32_t)i; + } +} + +static double run_accuracy_case(const char* name, size_t count, + void (*fill)(float*, float*, float*, float*, + uint32_t*, size_t)) { + float* x = (float*)malloc(count * sizeof(float)); + float* y = (float*)malloc(count * sizeof(float)); + float* z = (float*)malloc(count * sizeof(float)); + float* mass = (float*)malloc(count * sizeof(float)); + float* fx = (float*)calloc(count, sizeof(float)); + float* fy = (float*)calloc(count, sizeof(float)); + float* fz = (float*)calloc(count, sizeof(float)); + uint32_t* id = (uint32_t*)malloc(count * sizeof(uint32_t)); + uint32_t* morton = (uint32_t*)malloc(count * sizeof(uint32_t)); + double* ref_fx = (double*)malloc(count * sizeof(double)); + double* ref_fy = (double*)malloc(count * sizeof(double)); + double* ref_fz = (double*)malloc(count * sizeof(double)); + void* arena = malloc(ARENA_SIZE); + + if (!x || !y || !z || !mass || !fx || !fy || !fz || !id || !morton || + !ref_fx || !ref_fy || !ref_fz || !arena) { + (void)fprintf(stderr, "[FAIL] %s: allocation failure\n", name); + ++g_failures; + free(x); + free(y); + free(z); + free(mass); + free(fx); + free(fy); + free(fz); + free(id); + free(morton); + free(ref_fx); + free(ref_fy); + free(ref_fz); + free(arena); + return INFINITY; + } + + fill(x, y, z, mass, id, count); + direct_reference(x, y, z, mass, count, ref_fx, ref_fy, ref_fz); + + ps_particle_arrs_t arrs = {x, y, z, mass, fx, fy, fz, id, count}; + ps_context_t* ctx = NULL; + ps_config_t cfg = {arena, ARENA_SIZE}; + + if (ps_init(&ctx, &cfg) != PS_OK) { + (void)fprintf(stderr, "[FAIL] %s: ps_init failed\n", name); + ++g_failures; + free(x); + free(y); + free(z); + free(mass); + free(fx); + free(fy); + free(fz); + free(id); + free(morton); + free(ref_fx); + free(ref_fy); + free(ref_fz); + free(arena); + return INFINITY; + } + + float min_b = 0.0F; + float max_b = 0.0F; + float range = 0.0F; + if (ps_prepare_particles(&arrs, morton, &min_b, &max_b, &range) != PS_OK) { + (void)fprintf(stderr, "[FAIL] %s: ps_prepare_particles failed\n", + name); + ++g_failures; + } else { + float root_c = min_b + (range * 0.5F); + if (ps_calc_forces(ctx, &arrs, morton, root_c, root_c, root_c, + range * 0.5F) != PS_OK) { + (void)fprintf(stderr, "[FAIL] %s: ps_calc_forces failed\n", name); + ++g_failures; + } + } + + double err2 = 0.0; + double ref2 = 0.0; + double max_rel = 0.0; + + for (size_t i = 0; i < count; ++i) { + uint32_t original = id[i]; + double dx = (double)fx[i] - ref_fx[original]; + double dy = (double)fy[i] - ref_fy[original]; + double dz = (double)fz[i] - ref_fz[original]; + double e2 = (dx * dx) + (dy * dy) + (dz * dz); + double r2 = (ref_fx[original] * ref_fx[original]) + + (ref_fy[original] * ref_fy[original]) + + (ref_fz[original] * ref_fz[original]); + err2 += e2; + ref2 += r2; + if (r2 > 1e-24) { + double rel = sqrt(e2 / r2); + if (rel > max_rel) { + max_rel = rel; + } + } + } + + double rms_rel = ref2 > 0.0 ? sqrt(err2 / ref2) : 0.0; + (void)printf("[METRIC] %s: rms_relative_error=%.6f max_relative_error=%.6f\n", + name, rms_rel, max_rel); + + free(x); + free(y); + free(z); + free(mass); + free(fx); + free(fy); + free(fz); + free(id); + free(morton); + free(ref_fx); + free(ref_fy); + free(ref_fz); + free(arena); + + return rms_rel; +} + +int main(void) { + run_dipole_sign_test(); + + double uniform = run_accuracy_case("uniform-512", 512, fill_uniform); + double clusters = + run_accuracy_case("two-clusters-512", 512, fill_two_clusters); + + /* Intentionally generous initial gates: these are validation smoke tests, + not production accuracy targets. */ + if (uniform > 0.10) { + failf("uniform RMS relative error", uniform, 0.10, 0.0); + } + if (clusters > 0.10) { + failf("clustered RMS relative error", clusters, 0.10, 0.0); + } + + if (g_failures != 0) { + (void)fprintf(stderr, "\n%d accuracy validation(s) failed.\n", + g_failures); + return EXIT_FAILURE; + } + + (void)printf("\nAll accuracy validations passed.\n"); + return EXIT_SUCCESS; +} +// clang-format on diff --git a/tests/benchmark.c b/tests/benchmark.c index 01bf8e7..b65a765 100644 --- a/tests/benchmark.c +++ b/tests/benchmark.c @@ -49,7 +49,7 @@ void direct_nbody_baseline(const ps_particle_arrs_t* arrs) { float p_dz = arrs->z[j] - t_z; float p_dist_sq = - (p_dx * p_dx) + (p_dy * p_dy) + (p_dz * p_dz) + 2.0F; + (p_dx * p_dx) + (p_dy * p_dy) + (p_dz * p_dz) + 0.1F; float inv_dist = 1.0F / sqrtf(p_dist_sq); float inv_dist3 = inv_dist * inv_dist * inv_dist; diff --git a/tests/main.c b/tests/main.c index b0fb189..9f294fe 100644 --- a/tests/main.c +++ b/tests/main.c @@ -156,13 +156,14 @@ void test_radix_sort(void) { ps_init(&ctx, &cfg); uint32_t morton_codes[4] = {999, 10, 500, 42}; + uint32_t ids[4] = {0, 1, 2, 3}; float x[4] = {9.0F, 1.0F, 5.0F, 4.0F}; float y[4] = {9.0F, 1.0F, 5.0F, 4.0F}; float z[4] = {9.0F, 1.0F, 5.0F, 4.0F}; float mass[4] = {9.0F, 1.0F, 5.0F, 4.0F}; - float fx[2] = {0.0F, 0.0F}, fy[2] = {0.0F, 0.0F}, fz[2] = {0.0F, 0.0F}; + float fx[4] = {0.0F}, fy[4] = {0.0F}, fz[4] = {0.0F}; - ps_particle_arrs_t arrs = {x, y, z, mass, fx, fy, fz, 0, 4}; + ps_particle_arrs_t arrs = {x, y, z, mass, fx, fy, fz, ids, 4}; ps_impl_sort_particles(&ctx->arena, morton_codes, &arrs); @@ -243,22 +244,21 @@ void test_fmm_interaction_pass(void) { node_a->x = -5.0F; node_a->y = -5.0F; node_a->z = -5.0F; - node_a->half_width = 2.5F; + node_a->half_width = 1.0F; node_b->x = 5.0F; node_b->y = 5.0F; node_b->z = 5.0F; - node_b->half_width = 2.5F; + node_b->half_width = 1.0F; node_b->multipole[0] = 1.0F; ps_impl_fmm_interaction_pass(node_a, node_b); // vector from A to B: dx=10, dy=10, dz=10 - // dist_sq = 300, dist = sqrt(300) ~= 17.32 - // inv_r3 = 1.0 / (dist * dist_sq) ~= 0.00019245 - // F_field = m * dx * inv_r3 ~= 1.0 * 10 * 0.00019245 ~= 0.0019245 - float expected_field = 0.0019245F; + // softened dist_sq = 300.1, dist ~= 17.3234 + // F_field = m * dx / (dist_sq^1.5) ~= 0.00192354 + float expected_field = 0.00192354F; // verify node a local expansion TEST_ASSERT_FLOAT_EQ(expected_field, node_a->local[1], 1e-6F); // F_x @@ -374,7 +374,7 @@ void test_fmm_p2p_pass(void) { ps_impl_fmm_p2p_pass(ctx->root, ctx->root, &arrs); // dist = 1.0 - // dist_sq = 1.0^2 + 2.0 = 3.0 + // dist_sq = 1.0^2 + 0.1 = 1.1 // F_mag = (mass1 * mass2) / (dist_sq^1.5) float expected_f = (2.0F * 3.0F) * powf(1.1F, -1.5F); @@ -395,6 +395,7 @@ int main(void) { RUN_TEST(test_morton_encoding); RUN_TEST(test_arena_allocator); RUN_TEST(test_octree_insertion); + RUN_TEST(test_radix_sort); RUN_TEST(test_fmm_upward_pass); RUN_TEST(test_fmm_interaction_pass); RUN_TEST(test_fmm_downward_pass);