/******************************************************************************* * This file is part of SWIFT. * Copyright (C) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl) * James Willis (james.s.willis@durham.ac.uk) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * 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 Lesser General Public License * along with this program. If not, see . * ******************************************************************************/ #include "const.h" #include "kernel_gravity.h" #include "kernel_long_gravity.h" #include #include #include #include #define numPoints (1 << 7) /** * @brief The Gadget-2 gravity kernel function * * Taken from Gadget-2.0.7's forcetree.c lines 2755-2800 * * @param r The distance between particles * @param epsilon The cut-off distance of the kernel */ float gadget(float r, float epsilon) { const float h = epsilon; const float h_inv = 1.f / h; const float u = r * h_inv; if (u >= 1) { const float r_inv = 1. / r; return r_inv * r_inv * r_inv; } else { if (u < 0.5) return h_inv * h_inv * h_inv * (10.666666666667 + u * u * (32.0 * u - 38.4)); else return h_inv * h_inv * h_inv * (21.333333333333 - 48.0 * u + 38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u)); } } int main(int argc, char *argv[]) { const float h = 3.f; const float r_max = 6.f; for (int k = 1; k < numPoints; ++k) { const float r = (r_max * k) / numPoints; const float gadget_w = gadget(r, h); const float h_inv = 1.f / h; const float h_inv3 = h_inv * h_inv * h_inv; const float u = r * h_inv; float swift_w; if (r >= h) { swift_w = 1 / (r * r * r); } else { kernel_grav_eval(u, &swift_w); swift_w *= h_inv3; } if (fabsf(gadget_w - swift_w) > 1e-5 * fabsf(gadget_w)) { printf("%2d: r= %f h= %f u= %f Wg(r,h)= %f Ws(r,h)= %f\n", k, r, h, u, gadget_w, swift_w); printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w); return 1; } } printf("\nAll values are consistent\n"); /* Now test the long range function */ /* const float a_smooth = 4.5f; */ /* for (int k = 1; k < numPoints; ++k) { */ /* const float r = (r_max * k) / numPoints; */ /* const float u = r / a_smooth; */ /* float swift_w; */ /* kernel_long_grav_eval(u, &swift_w); */ /* float gadget_w = erfcf(u / 2) + u * expf(-u * u / 4) / sqrtf(M_PI); */ /* if (fabsf(gadget_w - swift_w) > 1e-4 * fabsf(gadget_w)) { */ /* printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r, * a_smooth, */ /* u, swift_w, gadget_w); */ /* printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w); */ /* return 1; */ /* } */ /* } */ return 0; }