diff --git a/Build/Makefile.M100 b/Build/Makefile.M100 index 1fa008e8cf29ab9f112d3ad399a6fe4dd5faa901..eaf92e21428d06d38ce46146f1cee1677a972341 100644 --- a/Build/Makefile.M100 +++ b/Build/Makefile.M100 @@ -1,3 +1,4 @@ + CC = gcc CXX = g++ @@ -5,8 +6,11 @@ MPICC = mpicc MPIC++ = mpiCC -#FFTW_INCL= -I/home/taffoni/sw/include -#FFTW_LIB= -L/home/taffoni/sw/lib +FFTW_INCL= -I/cineca/prod/opt/libraries/fftw/3.3.8/spectrum_mpi--10.4.0--binary/include/ +FFTW_LIB= /cineca/prod/opt/libraries/fftw/3.3.8/spectrum_mpi--10.4.0--binary/lib + +FITSIO_INCL= -I/m100_work/IscrC_CD-DLS/cfitsio-3.49 +FITSIO_LIB= /m100_work/IscrC_CD-DLS/cfitsio-3.49 NVCC = nvcc @@ -15,6 +19,6 @@ NVLIB = -L/cineca/prod/opt/compilers/cuda/10.1/none/lib64/ -lcudart -lcuda OMP= -fopenmp -CFLAGS += -I. $(FFTW_INCL) $(GSL_INCL) $(MPI_INCL) +CFLAGS += -I. $(FFTW_INCL) $(GSL_INCL) $(MPI_INCL) $(FITSIO_INCL) OPTIMIZE = $(OMP) -O3 -mtune=native diff --git a/Makefile b/Makefile index 23d161d0f6b432626329c9ffa52f94d2e31fab2a..f632e6ed188e180fe79f0ccdb42d037fe29fb751 100644 --- a/Makefile +++ b/Makefile @@ -27,14 +27,21 @@ endif #OPT += -DNVIDIA # perform one-side communication (suggested) instead of reduce (only if MPI is active) -OPT += -DONE_SIDE +#OPT += -DONE_SIDE # write the full 3D cube of gridded visibilities and its FFT transform #OPT += -DWRITE_DATA # write the final image OPT += -DWRITE_IMAGE # perform w-stacking phase correction OPT += -DPHASE_ON +# Support CFITSIO +OPT += -DFITSIO +# Perform true parallel images writing +#OPT += -DPARALLELIO +ifeq (FITSIO,$(findstring FITSIO,$(OPT))) + LIBS += -L$(FITSIO_LIB) -lcfitsio +endif DEPS = w-stacking.h w-stacking-fftw.c w-stacking.cu phase_correction.cu COBJ = w-stacking.o w-stacking-fftw.o phase_correction.o diff --git a/scripts/plotgrid_fits.py b/scripts/plotgrid_fits.py new file mode 100644 index 0000000000000000000000000000000000000000..0b4583ea6d8c48febb643e62dcbcce294a1b0f46 --- /dev/null +++ b/scripts/plotgrid_fits.py @@ -0,0 +1,32 @@ +#!/usr/bin/python3 +import numpy as np +import matplotlib.pyplot as plt +from astropy.io import fits +real_fits = "/home/emanuele/hpc_imaging/test_fits_real.fits" +imag_fits = "/home/emanuele/hpc_imaging/test_fits_img.fits" +nplanes = 1 + + +with fits.open(real_fits) as hdu_real: + img_hdu_real = np.array(hdu_real[0].data) + +with fits.open(imag_fits) as hdu_imag: + img_hdu_imag = np.array(hdu_imag[0].data) + + +xaxis = int(np.sqrt(img_hdu_real.size)) +yaxes = xaxis +residual = np.vectorize(complex)(img_hdu_real, img_hdu_imag) + +cumul2d = residual.reshape((xaxis,yaxes,nplanes), order='F') + +for i in range(nplanes): + gridded = np.squeeze(cumul2d[:,:,i]) + ax = plt.subplot() + img = ax.imshow(np.abs(np.fft.fftshift(gridded)), aspect='auto', interpolation='none', origin='lower') + ax.set_xlabel('cell') + ax.set_ylabel('cell') + cbar = plt.colorbar(img) + cbar.set_label('norm(FFT)',size=18) + figname='fits_image_' + str(i) + '.png' + plt.savefig(figname) diff --git a/w-stacking-fftw.c b/w-stacking-fftw.c index 5469014518e9e03a3d65ee9461bccb9c58cebfc1..55244baa73b6378fde3c111568963353bc33ee24 100644 --- a/w-stacking-fftw.c +++ b/w-stacking-fftw.c @@ -1,6 +1,9 @@ #include #include #include +#ifdef FITSIO +#include "fitsio.h" +#endif #ifdef USE_MPI #include #ifdef USE_FFTW @@ -44,6 +47,7 @@ int main(int argc, char * argv[]) int rank; int size; + // Define main filenames FILE * pFile; FILE * pFile1; @@ -51,7 +55,7 @@ int main(int argc, char * argv[]) FILE * pFileimg; // Global filename to be composed char filename[1000]; - + // MS paths char datapath[900]; char datapath_multi[NFILES][900]; @@ -123,6 +127,19 @@ int main(int argc, char * argv[]) // Half support size double w_supporth = (double)((w_support-1)/2)*dx; + + // Initialize FITS image parameters + fitsfile *fptreal; + fitsfile *fptrimg; + int status; + char testfitsreal[FILENAMELENGTH] = "parallel_np2_real.fits"; + char testfitsimag[FILENAMELENGTH] = "parallel_np2_img.fits"; + + long naxis = 2; + long naxes[2] = { grid_size_x, grid_size_y }; + + + // Internal profiling parameters clock_t start, end, start0, startk, endk; double setup_time, process_time, mpi_time, fftw_time, tot_time, kernel_time, reduce_time, compose_time, phase_time; @@ -195,15 +212,14 @@ if(rank == 0){ // LOCAL grid size xaxis = local_grid_size_x; yaxis = local_grid_size_y; - clock_gettime(CLOCK_MONOTONIC, &begin); start = clock(); // INPUT FILES (only the first ndatasets entries are used) int ndatasets = 1; - //strcpy(datapath_multi[0],"data/newgauss2noconj_t201806301100_SBL180.binMS/"); + strcpy(datapath_multi[0],"/m100_scratch/userexternal/ederubei/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.binMS/"); //strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/newgauss4_t201806301100_SBL180.binMS/"); - strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.binMS/"); + //strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.binMS/"); //strcpy(datapath_multi[1],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_134MHz.pre-cal.binMS/"); strcpy(datapath,datapath_multi[0]); @@ -917,8 +933,8 @@ if(rank == 0){ clock_gettime(CLOCK_MONOTONIC, &begin); start = clock(); if(rank == 0)printf("PHASE CORRECTION\n"); - double* image_real = (double*) calloc(xaxis*yaxis,sizeof(double)); - double* image_imag = (double*) calloc(xaxis*yaxis,sizeof(double)); + double* image_real = (double*) calloc(xaxis*yaxis,sizeof(double)); + double* image_imag = (double*) calloc(xaxis*yaxis,sizeof(double)); phase_correction(gridss,image_real,image_imag,xaxis,yaxis,num_w_planes,grid_size_x,grid_size_y,resolution,wmin,wmax,num_threads); @@ -931,23 +947,106 @@ if(rank == 0){ if(rank == 0) { - pFilereal = fopen (fftfile2,"wb"); - pFileimg = fopen (fftfile3,"wb"); - fclose(pFilereal); - fclose(pFileimg); + #ifdef FITSIO + printf("REMOVING RESIDUAL FITS FILE\n"); + remove(testfitsreal); + remove(testfitsimag); + + + printf("FITS CREATION\n"); + status = 0; + + fits_create_file(&fptrimg, testfitsimag, &status); + fits_create_img(fptrimg, DOUBLE_IMG, naxis, naxes, &status); + fits_close_file(fptrimg, &status); + + status = 0; + + fits_create_file(&fptreal, testfitsreal, &status); + fits_create_img(fptreal, DOUBLE_IMG, naxis, naxes, &status); + fits_close_file(fptreal, &status); + #endif + + pFilereal = fopen (fftfile2, "wb"); + pFileimg = fopen (fftfile3, "wb"); + fclose(pFilereal); + fclose(pFileimg); } #ifdef USE_MPI MPI_Barrier(MPI_COMM_WORLD); #endif if(rank == 0)printf("WRITING IMAGE\n"); + long * fpixel = (long *) malloc(sizeof(long)*naxis); + long * lpixel = (long *) malloc(sizeof(long)*naxis); + + #ifdef USE_MPI + MPI_Barrier(MPI_COMM_WORLD); + #endif + #ifdef PARALLELIO + #ifdef FITSIO + + fpixel[0] = 1; + fpixel[1] = rank*yaxis+1; + lpixel[0] = xaxis; + lpixel[1] = (rank+1)*yaxis; + + status = 0; + fits_open_image(&fptreal, testfitsreal, READWRITE, &status); + fits_write_subset(fptreal, TDOUBLE, fpixel, lpixel, image_real, &status); + fits_close_file(fptreal, &status); + + status = 0; + fits_open_image(&fptrimg, testfitsimag, READWRITE, &status); + fits_write_subset(fptrimg, TDOUBLE, fpixel, lpixel, image_imag, &status); + fits_close_file(fptrimg, &status); + + #endif //FITSIO + + pFilereal = fopen (fftfile2,"ab"); + pFileimg = fopen (fftfile3,"ab"); + + long global_index = rank*(xaxis*yaxis)*sizeof(double); + + fseek(pFilereal, global_index, SEEK_SET); + fwrite(image_real, xaxis*yaxis, sizeof(double), pFilereal); + fseek(pFileimg, global_index, SEEK_SET); + fwrite(image_imag, xaxis*yaxis, sizeof(double), pFileimg); + + fclose(pFilereal); + fclose(pFileimg); + + #ifdef USE_MPI + MPI_Barrier(MPI_COMM_WORLD); + #endif + #else for (int isector=0; isector