diff --git a/tests/testKernelGrav.c b/tests/testKernelGrav.c index c05d15bcb4930f8430c4d4371620623e07605f34..42ce6066b605c8cc9ac78a6e565dd59a500112f4 100644 --- a/tests/testKernelGrav.c +++ b/tests/testKernelGrav.c @@ -17,14 +17,16 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ +#include "const.h" #include "kernel_gravity.h" +#include "kernel_long_gravity.h" #include <math.h> #include <stdio.h> #include <stdlib.h> #include <strings.h> -#define numPoints (1 << 6) +#define numPoints (1 << 7) /** * @brief The Gadget-2 gravity kernel function @@ -54,7 +56,7 @@ float gadget(float r, float h) { int main() { const float h = 3.f; - const float r_max = 5.f; + const float r_max = 6.f; for (int k = 1; k < numPoints; ++k) { @@ -82,5 +84,29 @@ int main() { } printf("\nAll values are consistent\n"); + + /* Now test the long range function */ + const float a_smooth = const_gravity_a_smooth; + + 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 = erfc(u / 2) + u * exp(-u * u / 4) / sqrt(M_PI); + + 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); + + if (fabsf(gadget_w - swift_w) > 2e-7) { + printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w); + return 1; + } + } + return 0; }