Commit 1688d72a authored by Claudio Gheller's avatar Claudio Gheller
Browse files

Merge branch 'test_derubeis' into 'main' 

CFITSIO and parallel images writing implementation

See merge request !3
parents 0da6cc99 b48befbc

Build/Makefile.M100

+7 −3
Original line number Diff line number Diff line 


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

Makefile

+8 −1
Original line number Diff line number Diff line
@@ -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

scripts/plotgrid_fits.py

0 → 100644
+32 −0
Original line number Diff line number Diff line 
#!/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)

w-stacking-fftw.c

+112 −11
Original line number Diff line number Diff line 
#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#ifdef FITSIO

#include "fitsio.h"

#endif

#ifdef USE_MPI

#include <mpi.h>

#ifdef USE_FFTW
@@ -44,6 +47,7 @@ int main(int argc, char * argv[]) 
	int rank;

	int size;





        // Define main filenames

	FILE * pFile;

	FILE * pFile1;
@@ -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]);
@@ -931,6 +947,26 @@ if(rank == 0){ 


        if(rank == 0)

        {

            #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);
@@ -940,6 +976,49 @@ if(rank == 0){ 
	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<size; isector++)

	{

	    #ifdef USE_MPI
@@ -947,7 +1026,27 @@ if(rank == 0){ 
            #endif

	    if(isector == rank)

	    {

               #ifdef FITSIO



	       printf("%d writing\n",isector);

               

               fpixel[0] = 1;

               fpixel[1] = isector*yaxis+1;

               lpixel[0] = xaxis;

               lpixel[1] = (isector+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");
@@ -962,10 +1061,12 @@ if(rank == 0){ 
               fclose(pFileimg);

	    }

	}

	#endif //PARALLELIO

	#ifdef USE_MPI

	MPI_Barrier(MPI_COMM_WORLD);

        #endif





#endif //WRITE_IMAGE