Add tests/gpu_sim_check

2 files changed, 148 insertions, 0 deletions

diff --git a/tests/.gitignore b/tests/.gitignore
index 3c726461..b855d2bd 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -1,4 +1,5 @@
 *.o
 .deps
 list_check
+gpu_sim_check
 .dirstamp
diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c
new file mode 100644
index 00000000..28e85bcc
--- /dev/null
+++ b/tests/gpu_sim_check.c
@@ -0,0 +1,147 @@
+/*
+ * gpu_sim_check.c
+ *
+ * Check that GPU simulation agrees with CPU version
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "../src/diffraction.h"
+#include "../src/diffraction-gpu.h"
+#include "../src/detector.h"
+#include "../src/beam-parameters.h"
+#include "../src/utils.h"
+
+
+int main(int argc, char *argv[])
+{
+	struct gpu_context *gctx;
+	struct image gpu_image;
+	struct image cpu_image;
+	UnitCell *cell;
+	UnitCell *cell_raw;
+	struct detector *det;
+	struct beam_params *beam;
+	int i;
+	double gpu_min, gpu_max, gpu_tot;
+	double cpu_min, cpu_max, cpu_tot;
+	double dev, perc;
+	const double sep = 20.0;
+
+	gctx = setup_gpu(1, NULL, NULL, NULL, 0);
+	if ( gctx == NULL ) {
+		ERROR("Couldn't set up GPU.\n");
+		return 1;
+	}
+
+	cell_raw = cell_new_from_parameters(28.1e-9, 28.1e-9, 16.5e-9,
+	                          deg2rad(90.0), deg2rad(90.0), deg2rad(120.0));
+
+	cell = cell_rotate(cell_raw, random_quaternion());
+
+	gpu_image.width = 1024;
+	gpu_image.height = 1024;
+	cpu_image.width = 1024;
+	cpu_image.height = 1024;
+	det = calloc(1, sizeof(struct detector));
+	det->n_panels = 2;
+	det->panels = calloc(2, sizeof(struct panel));
+
+	det->panels[0].min_fs = 0;
+	det->panels[0].max_fs = 1023;
+	det->panels[0].min_ss = 0;
+	det->panels[0].max_ss = 511;
+	det->panels[0].fsx = 1;
+	det->panels[0].fsy = 0;
+	det->panels[0].ssx = 0;
+	det->panels[0].ssy = 1;
+	det->panels[0].xfs = 1;
+	det->panels[0].yfs = 0;
+	det->panels[0].xss = 0;
+	det->panels[0].yss = 1;
+	det->panels[0].cnx = -512.0;
+	det->panels[0].cny = -512.0-sep;
+	det->panels[0].clen = 100.0e-3;
+	det->panels[0].res = 9090.91;
+
+	det->panels[1].min_fs = 0;
+	det->panels[1].max_fs = 1023;
+	det->panels[1].min_ss = 512;
+	det->panels[1].max_ss = 1023;
+	det->panels[1].fsx = 1;
+	det->panels[1].fsy = 0;
+	det->panels[1].ssx = 0;
+	det->panels[1].ssy = 1;
+	det->panels[1].xfs = 1;
+	det->panels[1].yfs = 0;
+	det->panels[1].xss = 0;
+	det->panels[1].yss = 1;
+	det->panels[1].cnx = -512.0;
+	det->panels[1].cny = sep;
+	det->panels[1].clen = 100.0e-3;
+	det->panels[1].res = 9090.91;
+
+	cpu_image.det = det;
+	gpu_image.det = det;
+
+	beam = calloc(1, sizeof(struct beam_params));
+	beam->fluence = 1.0e15;
+	beam->beam_radius = 1.0e-6;
+	beam->photon_energy = 9000.0;
+	beam->bandwidth = 0.1 / 100.0;
+	beam->divergence = 0.0;
+	beam->dqe = 1.0;
+	beam->adu_per_photon = 1.0;
+	cpu_image.beam = beam;
+	gpu_image.beam = beam;
+
+	cpu_image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
+	gpu_image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
+
+	get_diffraction_gpu(gctx, &gpu_image, 8, 8, 8, cell);
+	get_diffraction(&cpu_image, 8, 8, 8, NULL, NULL, NULL, cell,
+	                GRADIENT_MOSAIC, "1");
+
+	gpu_min = +INFINITY;  gpu_max = -INFINITY;  gpu_tot = 0.0;
+	cpu_min = +INFINITY;  cpu_max = -INFINITY;  cpu_tot = 0.0;
+	dev = 0.0;
+	for ( i=0; i<1024*1024; i++ ) {
+
+		const double cpu = cpu_image.data[i];
+		const double gpu = gpu_image.data[i];
+
+		if ( cpu > cpu_max ) cpu_max = cpu;
+		if ( cpu < cpu_min ) cpu_min = cpu;
+		if ( gpu > gpu_max ) gpu_max = gpu;
+		if ( gpu < gpu_min ) gpu_min = gpu;
+		gpu_tot += gpu;
+		cpu_tot += cpu;
+		dev += fabs(gpu - cpu);
+
+	}
+	perc = 100.0*dev/cpu_tot;
+
+	STATUS("GPU: min=%8e, max=%8e, total=%8e\n", gpu_min, gpu_max, gpu_tot);
+	STATUS("CPU: min=%8e, max=%8e, total=%8e\n", cpu_min, cpu_max, cpu_tot);
+	STATUS("dev = %8e (%5.2f%% of CPU total)\n", dev, perc);
+
+	cell_free(cell);
+	free_detector_geometry(det);
+	free(beam);
+
+	if ( perc > 1.0 ) return 1;
+
+	return 0;
+}