diff --git a/tests/testNonSymInt.c b/tests/testNonSymInt.c new file mode 100644 index 0000000000000000000000000000000000000000..92d8ff24761afc669658daa7c0db6a040416c708 --- /dev/null +++ b/tests/testNonSymInt.c @@ -0,0 +1,257 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@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 <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ + +#include <fenv.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <unistd.h> +#include "swift.h" + +enum velocity_types { + velocity_zero, + velocity_random, + velocity_divergent, + velocity_rotating +}; + +char *serial_filename = "test_nonsym_serial.dat"; +char *vec_filename = "test_nonsym_vec.dat"; + +/** + * @brief Constructs an array of particles in a valid state prior to + * a IACT_NONSYM and IACT_NONSYM_VEC call. + * + * @param count No. of particles to create + * @param offset The position of the particle offset from (0,0,0). + * @param spacing Particle spacing. + * @param h The smoothing length of the particles in units of the inter-particle + *separation. + * @param partId The running counter of IDs. + * @param vel The type of velocity field (0, random, divergent, rotating) + */ +struct part *make_particles(int count, double *offset, double spacing, double h, + long long *partId, enum velocity_types vel) { + + struct part *particles; + if (posix_memalign((void **)&particles, part_align, + count * sizeof(struct part)) != 0) { + error("couldn't allocate particles, no. of particles: %d", (int)count); + } + + /* Construct the particles */ + struct part *p; + for (size_t i = 0; i < VEC_SIZE + 1; ++i) { + p = &particles[i]; + p->x[0] = offset[0] + spacing * i; + p->x[1] = offset[1] + spacing * i; + p->x[2] = offset[2] + spacing * i; + switch (vel) { + case velocity_zero: + p->v[0] = 0.f; + p->v[1] = 0.f; + p->v[2] = 0.f; + break; + case velocity_random: + p->v[0] = random_uniform(-0.05, 0.05); + p->v[1] = random_uniform(-0.05, 0.05); + p->v[2] = random_uniform(-0.05, 0.05); + break; + case velocity_divergent: + p->v[0] = p->x[0] - 1.5 * spacing; + p->v[1] = p->x[1] - 1.5 * spacing; + p->v[2] = p->x[2] - 1.5 * spacing; + break; + case velocity_rotating: + p->v[0] = p->x[1]; + p->v[1] = -p->x[0]; + p->v[2] = 0.f; + break; + } + p->h = h; + p->id = ++(*partId); + p->mass = 1.0f; + } + return particles; +} + +/** + * @brief Dumps all particle information to a file + */ +void dump_indv_particle_fields(char *fileName, struct part *p) { + + FILE *file = fopen(fileName, "a"); + + fprintf( + file, + "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e " + "%13e %13e %13e\n", + p->id, p->x[0], p->x[1], + p->x[2], p->v[0], p->v[1], + p->v[2], p->rho, p->rho_dh, + p->density.wcount, p->density.wcount_dh, +#ifdef GADGET2_SPH + p->div_v, p->density.rot_v[0], + p->density.rot_v[1], p->density.rot_v[2] +#else + 0., 0., 0., 0. +#endif + ); + fclose(file); +} + +/** + * @brief Creates a header for the output file + */ +void write_header(char *fileName) { + + FILE *file = fopen(fileName, "w"); + /* Write header */ + fprintf(file, + "# %4s %10s %10s %10s %10s %10s %10s %13s %13s %13s %13s %13s " + "%13s %13s %13s\n", + "ID", "pos_x", "pos_y", "pos_z", "v_x", "v_y", "v_z", "rho", "rho_dh", + "wcount", "wcount_dh", "div_v", "curl_vx", "curl_vy", "curl_vz"); + fprintf(file,"\nPARTICLES BEFORE INTERACTION:\n"); + fclose(file); +} + +/* + * @brief Calls the serial and vectorised version of the non-symmetrical density interaction. + * + * @param parts Particle array to be interacted + * @param count No. of particles to be interacted + * + */ +void test_nonsym_density_interaction(struct part *parts, int count) { + + /* Use the first particle in the array as the one that gets updated. */ + struct part pi = parts[0]; + //const float hig2 = hi * hi * kernel_gamma2; + + FILE *file; + write_header(serial_filename); + write_header(vec_filename); + + /* Dump state of particles before serial interaction. */ + dump_indv_particle_fields(serial_filename,&pi); + for(size_t i = 1; i < count; i++) + dump_indv_particle_fields(serial_filename,&parts[i]); + + /* Make copy of pi to be used in vectorised version. */ + struct part pi_vec = pi; + struct part pj_vec[VEC_SIZE]; + for(int i=0; i<VEC_SIZE; i++) + pj_vec[i] = parts[i + 1]; + + float r2q[VEC_SIZE] __attribute__((aligned(sizeof(float) * VEC_SIZE))); + float hiq[VEC_SIZE] __attribute__((aligned(sizeof(float) * VEC_SIZE))); + float hjq[VEC_SIZE] __attribute__((aligned(sizeof(float) * VEC_SIZE))); + float dxq[3 * VEC_SIZE] __attribute__((aligned(sizeof(float) * VEC_SIZE))); + struct part *piq[VEC_SIZE], *pjq[VEC_SIZE]; + + /* Perform serial interaction */ + for(int i=1; i<count; i++) { + /* Compute the pairwise distance. */ + float r2 = 0.0f; + float dx[3]; + for (int k = 0; k < 3; k++) { + dx[k] = pi.x[k] - parts[i].x[k]; + r2 += dx[k] * dx[k]; + } + + runner_iact_nonsym_density(r2, dx, pi.h, parts[i].h, &pi, &parts[i]); + } + + file = fopen(serial_filename, "a"); + fprintf(file,"\nPARTICLES AFTER INTERACTION:\n"); + fclose(file); + + /* Dump result of serial interaction. */ + dump_indv_particle_fields(serial_filename,&pi); + for(size_t i = 1; i < count; i++) + dump_indv_particle_fields(serial_filename,&parts[i]); + + /* Setup arrays for vector interaction. */ + for(int i=0; i<VEC_SIZE; i++) { + /* Compute the pairwise distance. */ + float r2 = 0.0f; + float dx[3]; + for (int k = 0; k < 3; k++) { + dx[k] = pi_vec.x[k] - pj_vec[i].x[k]; + r2 += dx[k] * dx[k]; + } + r2q[i] = r2; + dxq[3 * i + 0] = dx[0]; + dxq[3 * i + 1] = dx[1]; + dxq[3 * i + 2] = dx[2]; + hiq[i] = pi_vec.h; + hjq[i] = pj_vec[i].h; + piq[i] = &pi_vec; + pjq[i] = &pj_vec[i]; + } + + /* Dump state of particles before vector interaction. */ + dump_indv_particle_fields(vec_filename,piq[0]); + for(size_t i = 0; i < VEC_SIZE; i++) + dump_indv_particle_fields(vec_filename,pjq[i]); + + /* Perform vector interaction. */ + runner_iact_nonsym_vec_density(r2q, dxq, hiq, hjq, piq, pjq); + + file = fopen(vec_filename, "a"); + fprintf(file,"\nPARTICLES AFTER INTERACTION:\n"); + fclose(file); + + /* Dump result of serial interaction. */ + dump_indv_particle_fields(vec_filename,piq[0]); + for(size_t i = 0; i < VEC_SIZE; i++) + dump_indv_particle_fields(vec_filename,pjq[i]); +} + +/* And go... */ +int main(int argc, char *argv[]) { + double h = 1.2348, spacing = 0.5; + char outputFileNameExtension[200] = ""; + char vecOutputFileName[200] = ""; + double offset[3] = {0.0,0.0,0.0}; + int vel = velocity_zero; + int count = VEC_SIZE + 1; + + /* Get some randomness going */ + srand(0); + + /* Help users... */ + message("Smoothing length: h = %f", h); + message("Kernel: %s", kernel_name); + message("Neighbour target: N = %f", + h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0); + message("div_v target: div = %f", vel == 2 ? 3.f : 0.f); + message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f); + printf("\n"); + + /* Build the infrastructure */ + static long long partId = 0; + struct part *particles = make_particles(count,offset,spacing,h,&partId,vel); + + /* Call the test non-sym density test. */ + test_nonsym_density_interaction(particles,count); + + return 0; +}