diff --git a/src/grid_map_GPU.cu b/src/grid_map_GPU.cu index 0dd0bb6..face872 100644 --- a/src/grid_map_GPU.cu +++ b/src/grid_map_GPU.cu @@ -1,877 +1,1011 @@ /** * @Author Christoph Schaefer, EPFL (christophernstrerne.schaefer@epfl.ch), Gilles Fourestey (gilles.fourestey@epfl.ch) * @date July 2017 * @version 0,1 * */ #include #include "grid_map_GPU.cuh" #include "gradient2_GPU.cuh" #include #define BLOCK_SIZE_X 32 #define BLOCK_SIZE_Y 16 //#define ROT #define _SHARED_MEM #ifdef _SHARED_MEM #define SHARED __shared__ #warning "shared memory" extern __shared__ type_t shared[]; #else #define SHARED #endif #define Nx 1 #define Ny 0 #define cudasafe extern "C" { type_t myseconds(); } __device__ struct ellipse formeli_HPC(type_t a, type_t b, type_t c); //GPU mapping function declaration to change when figured out linkage problems __global__ void amplif_1_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); __global__ void amplif_2_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); __global__ void amplif_3_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); __global__ void amplif_4_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); __global__ void amplif_5_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); __global__ void amplif_6_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); +__global__ void shear_1_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); +__global__ void shear_2_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); +__global__ void shear_3_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); +__global__ void shear_4_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t z,int nbgridcells); __global__ void mass_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t mult, type_t ds, type_t dl, type_t h, type_t z,int nbgridcells); ////Map function selection map_gpu_function_t select_map_function(std::string mode, const struct runmode_param* runmode){ if (mode == "ampli"){ if(runmode->amplif == 1){ return &lif_1_grid_CPU_GPU; } else if(runmode->amplif == 2){ return &lif_2_grid_CPU_GPU; } else if(runmode->amplif == 3){ return &lif_3_grid_CPU_GPU; } else if(runmode->amplif == 4){ return &lif_4_grid_CPU_GPU; } else if(runmode->amplif == 5){ return &lif_5_grid_CPU_GPU; } else if(runmode->amplif == 6){ return &lif_6_grid_CPU_GPU; } else{ fprintf(stderr, "ERROR: Amplif mode %d not supported yet \n",runmode->amplif); exit(-1); } } else if(mode == "mass"){ if(runmode->mass == 1){ return &mass_1_grid_CPU_GPU; } else if(runmode->mass == 2){ return &mass_2_grid_CPU_GPU; } else if(runmode->mass == 3){ return &mass_3_grid_CPU_GPU; } else if(runmode->mass == 4){ return &mass_4_grid_CPU_GPU; } else{ fprintf(stderr, "ERROR: Mass mode %d not supported yet \n",runmode->mass); exit(-1); } } + else if(mode == "shear"){ + if(runmode->shear == 1){ + return &shear_1_grid_CPU_GPU; + } + else if(runmode->shear == 2){ + return &shear_2_grid_CPU_GPU; + } + else if(runmode->shear == 3){ + return &shear_3_grid_CPU_GPU; + } + else if(runmode->shear == 4){ + return &shear_4_grid_CPU_GPU; + } + else{ + fprintf(stderr, "ERROR: Mass mode %d not supported yet \n",runmode->mass); + exit(-1); + } + } + + + + else{ fprintf(stderr, "ERROR: No mode recognised \n"); exit(-1); } return 0; } ////Mass Map calculation, doesnt fit the bloody template... void map_mass_grid_GPU(map_gpu_function_t mapfunction, type_t *map, const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int nhalos ,int nbgridcells,int mode_amp, type_t zl, type_t zs ) { type_t dx = (frame->xmax - frame->xmin)/(nbgridcells - 1); type_t dy = (frame->ymax - frame->ymin)/(nbgridcells - 1); // map_mass_grid_GPU(mapfunction,map,cosmo, frame, lens, nhalos,mode_amp,zl, zs, dx, dy, nbgridcells, nbgridcells, 0, 0); } // void map_mass_grid_GPU(map_gpu_function_t mapfunction, type_t *map,const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int nhalos, int mode_amp, type_t zl, type_t zs, type_t dx, type_t dy, int nbgridcells_x, int nbgridcells_y, int istart, int jstart) { - int nBlocks_gpu = 0; - // Define the number of threads per block the GPU will use - cudaDeviceProp properties_gpu; - cudaGetDeviceProperties(&properties_gpu, 0); // Get properties of 0th GPU in use - - if (properties_gpu.maxThreadsDim[0]z[0]; } type_t dl0s = module_cosmodistances_objectObject(zl, zs, *cosmo); type_t dos = module_cosmodistances_observerObject(zs, *cosmo); type_t dol = module_cosmodistances_observerObject(zl, *cosmo); //select_ratio_function(std::string mode, const struct runmode_param* runmode, type_t dls, type_t ds); lens_gpu = (Potential_SOA *) malloc(sizeof(Potential_SOA)); lens_gpu->type = (int *) malloc(sizeof(int)); // Allocate variables on the GPU cudasafe(cudaMalloc( (void**)&(lens_kernel), sizeof(Potential_SOA)),"Gradient2gpu.cu : Alloc Potential_SOA: " ); cudasafe(cudaMalloc( (void**)&(type_gpu), nhalos*sizeof(int)),"Gradient2gpu.cu : Alloc type_gpu: " ); cudasafe(cudaMalloc( (void**)&(lens_x_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc x_gpu: " ); cudasafe(cudaMalloc( (void**)&(lens_y_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc y_gpu: " ); cudasafe(cudaMalloc( (void**)&(b0_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc b0_gpu: " ); cudasafe(cudaMalloc( (void**)&(angle_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc angle_gpu: " ); cudasafe(cudaMalloc( (void**)&(epot_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc epot_gpu: " ); cudasafe(cudaMalloc( (void**)&(rcore_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc rcore_gpu: " ); cudasafe(cudaMalloc( (void**)&(rcut_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc rcut_gpu: " ); cudasafe(cudaMalloc( (void**)&(anglecos_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc anglecos_gpu: " ); cudasafe(cudaMalloc( (void**)&(anglesin_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc anglesin_gpu: " ); cudasafe(cudaMalloc( (void**)&(dlsds_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc dlsds_gpu: " ); cudasafe(cudaMalloc( (void**)&(frame_gpu), sizeof(grid_param)),"Gradient2gpu.cu : Alloc frame_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_a_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_a_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_b_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_b_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_c_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_c_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_d_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_d_gpu: " ); cudasafe(cudaMalloc( (void**)&(map_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc map: " ); // Copy values to the GPU // cudasafe(cudaMemcpy(type_gpu,lens->type , nhalos*sizeof(int),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy type_gpu: " ); cudasafe(cudaMemcpy(lens_x_gpu,lens->position_x , nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy x_gpu: " ); cudasafe(cudaMemcpy(lens_y_gpu,lens->position_y , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy y_gpu: " ); cudasafe(cudaMemcpy(b0_gpu,lens->b0 , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2pu.cu : Copy b0_gpu: " ); cudasafe(cudaMemcpy(angle_gpu,lens->ellipticity_angle , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy angle_gpu: " ); cudasafe(cudaMemcpy(epot_gpu, lens->ellipticity_potential, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy epot_gpu: " ); cudasafe(cudaMemcpy(rcore_gpu, lens->rcore, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy rcore_gpu: " ); cudasafe(cudaMemcpy(rcut_gpu, lens->rcut, nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy rcut_gpu: " ); cudasafe(cudaMemcpy(anglecos_gpu, lens->anglecos, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy anglecos: " ); cudasafe(cudaMemcpy(anglesin_gpu, lens->anglesin, nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy anglesin: " ); cudasafe(cudaMemcpy(dlsds_gpu, lens->dlsds, nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy dlsds: " ); cudasafe(cudaMemcpy(frame_gpu, frame, sizeof(grid_param), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy fame_gpu: " ); // lens_gpu->type = type_gpu; lens_gpu->position_x = lens_x_gpu; lens_gpu->position_y = lens_y_gpu; lens_gpu->b0 = b0_gpu; lens_gpu->ellipticity_angle = angle_gpu; lens_gpu->ellipticity_potential = epot_gpu; lens_gpu->rcore = rcore_gpu; lens_gpu->rcut = rcut_gpu; lens_gpu->anglecos = anglecos_gpu; lens_gpu->anglesin = anglesin_gpu; lens_gpu->dlsds = dlsds_gpu; // cudaMemcpy(lens_kernel, lens_gpu, sizeof(Potential_SOA), cudaMemcpyHostToDevice); // type_t time = -myseconds(); // #if 0 int shalos = 0; while (shalos < nhalos) { int lens_type = lens->type[shalos]; type_t z = lens->z[shalos]; int count = 1; if (shalos + count < nhalos){ std::cerr << shalos << " " << count << " " << lens->type[shalos + count] << " " << lens->z[shalos + count] << " " << std::endl; while (lens->type[shalos + count] == lens_type and lens->z[shalos + count] == z ){ count++; if(shalos + count >= nhalos) break; //std::cerr << shalos << " " << count << " " << lens->type[shalos + count] << " " << lens->z[shalos + count] << " " << std::endl; } std::cerr << shalos << " " << count << " " << lens->type[shalos + count] << " " << lens->z[shalos + count] << " " << std::endl; } //if (shalos < nhalos) std::cerr << shalos << " " << count << " " << lens->type[shalos + count] << " " << lens->z[shalos + count] << " " << std::endl; shalos += count; } #endif //module_potentialDerivatives_totalGradient2_SOA_CPU_GPU(grid_grad2_a_gpu, grid_grad2_b_gpu, grid_grad2_c_gpu, grid_grad2_d_gpu, frame_gpu, lens_kernel, nhalos, dx, dy, nbgridcells_x, nbgridcells_y, istart, jstart); module_potentialDerivatives_Kmap_SOA_CPU_GPU(grid_grad2_a_gpu, grid_grad2_b_gpu, grid_grad2_c_gpu, grid_grad2_d_gpu, frame_gpu, lens_kernel, nhalos, dx, dy, nbgridcells_x, nbgridcells_y, istart, jstart); // cudasafe(cudaGetLastError(), "module_potentialDerivative_totalGradient_SOA_CPU_GPU"); //std::cerr << "ZS " << zs << " "<< dl0s << " "<< dos <h,zs,nbgridcells_x,nbgridcells_y,frame); cudasafe(cudaGetLastError(), "module_potentialDerivative_totalGradient_SOA_CPU_GPU"); cudaDeviceSynchronize(); // cudasafe(cudaMemcpy( map, map_gpu, (nbgridcells_x)*(nbgridcells_y)*sizeof(type_t), cudaMemcpyDeviceToHost )," --- Gradient2gpu.cu : Copy source_a_gpu: " ); // time += myseconds(); std::cout << " kernel time = " << time << " s." << std::endl; //printf("-----> %f %f \n",grid_grad_x[Nx], grid_grad_y[Ny]); // Free GPU memory cudaFree(lens_gpu); cudaFree(type_gpu); cudaFree(lens_x_gpu); cudaFree(lens_y_gpu); cudaFree(b0_gpu); cudaFree(angle_gpu); cudaFree(epot_gpu); cudaFree(rcore_gpu); cudaFree(rcut_gpu); cudaFree(anglecos_gpu); cudaFree(anglesin_gpu); cudaFree(dlsds_gpu); cudaFree(grid_grad2_a_gpu); cudaFree(grid_grad2_b_gpu); cudaFree(grid_grad2_c_gpu); cudaFree(grid_grad2_d_gpu); cudaFree(map_gpu); } ////General Map calculation void map_grid_GPU(map_gpu_function_t mapfunction, type_t *map, const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int nhalos ,int nbgridcells,int mode_amp, type_t z ) { type_t dx = (frame->xmax - frame->xmin)/(nbgridcells - 1); type_t dy = (frame->ymax - frame->ymin)/(nbgridcells - 1); // map_grid_GPU(mapfunction,map,cosmo, frame, lens, nhalos,mode_amp,z, dx, dy, nbgridcells, nbgridcells, 0, 0); } // void map_grid_GPU(map_gpu_function_t mapfunction, type_t *map,const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int nhalos, int mode_amp, type_t z, type_t dx, type_t dy, int nbgridcells_x, int nbgridcells_y, int istart, int jstart) { - int nBlocks_gpu = 0; - // Define the number of threads per block the GPU will use - cudaDeviceProp properties_gpu; - cudaGetDeviceProperties(&properties_gpu, 0); // Get properties of 0th GPU in use - - if (properties_gpu.maxThreadsDim[0]z[0], z, *cosmo); type_t dos = module_cosmodistances_observerObject(z, *cosmo); type_t dol = module_cosmodistances_observerObject(lens->z[0], *cosmo); //select_ratio_function(std::string mode, const struct runmode_param* runmode, type_t dls, type_t ds) lens_gpu = (Potential_SOA *) malloc(sizeof(Potential_SOA)); lens_gpu->type = (int *) malloc(sizeof(int)); // Allocate variables on the GPU cudasafe(cudaMalloc( (void**)&(lens_kernel), sizeof(Potential_SOA)),"Gradient2gpu.cu : Alloc Potential_SOA: " ); cudasafe(cudaMalloc( (void**)&(type_gpu), nhalos*sizeof(int)),"Gradient2gpu.cu : Alloc type_gpu: " ); cudasafe(cudaMalloc( (void**)&(lens_x_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc x_gpu: " ); cudasafe(cudaMalloc( (void**)&(lens_y_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc y_gpu: " ); cudasafe(cudaMalloc( (void**)&(b0_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc b0_gpu: " ); cudasafe(cudaMalloc( (void**)&(angle_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc angle_gpu: " ); cudasafe(cudaMalloc( (void**)&(epot_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc epot_gpu: " ); cudasafe(cudaMalloc( (void**)&(rcore_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc rcore_gpu: " ); cudasafe(cudaMalloc( (void**)&(rcut_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc rcut_gpu: " ); cudasafe(cudaMalloc( (void**)&(anglecos_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc anglecos_gpu: " ); cudasafe(cudaMalloc( (void**)&(anglesin_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc anglesin_gpu: " ); cudasafe(cudaMalloc( (void**)&(frame_gpu), sizeof(grid_param)),"Gradient2gpu.cu : Alloc frame_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_a_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_a_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_b_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_b_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_c_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_c_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_d_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_d_gpu: " ); cudasafe(cudaMalloc( (void**)&(map_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc map: " ); // Copy values to the GPU // cudasafe(cudaMemcpy(type_gpu,lens->type , nhalos*sizeof(int),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy type_gpu: " ); cudasafe(cudaMemcpy(lens_x_gpu,lens->position_x , nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy x_gpu: " ); cudasafe(cudaMemcpy(lens_y_gpu,lens->position_y , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy y_gpu: " ); cudasafe(cudaMemcpy(b0_gpu,lens->b0 , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2pu.cu : Copy b0_gpu: " ); cudasafe(cudaMemcpy(angle_gpu,lens->ellipticity_angle , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy angle_gpu: " ); cudasafe(cudaMemcpy(epot_gpu, lens->ellipticity_potential, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy epot_gpu: " ); cudasafe(cudaMemcpy(rcore_gpu, lens->rcore, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy rcore_gpu: " ); cudasafe(cudaMemcpy(rcut_gpu, lens->rcut, nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy rcut_gpu: " ); cudasafe(cudaMemcpy(anglecos_gpu, lens->anglecos, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy anglecos: " ); cudasafe(cudaMemcpy(anglesin_gpu, lens->anglesin, nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy anglesin: " ); cudasafe(cudaMemcpy(frame_gpu, frame, sizeof(grid_param), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy fame_gpu: " ); // lens_gpu->type = type_gpu; lens_gpu->position_x = lens_x_gpu; lens_gpu->position_y = lens_y_gpu; lens_gpu->b0 = b0_gpu; lens_gpu->ellipticity_angle = angle_gpu; lens_gpu->ellipticity_potential = epot_gpu; lens_gpu->rcore = rcore_gpu; lens_gpu->rcut = rcut_gpu; lens_gpu->anglecos = anglecos_gpu; lens_gpu->anglesin = anglesin_gpu; // cudaMemcpy(lens_kernel, lens_gpu, sizeof(Potential_SOA), cudaMemcpyHostToDevice); // type_t time = -myseconds(); - //std::cerr << "DLS " << dl0s << " " << dos << " " << dol << " " << lens->z[0] << " " << lens->z[1] << " " << z << std::endl; - //std::cerr <<"BLAAAAAAAAAAAA " << grid_grad2_a_gpu[0] << std::endl; // module_potentialDerivatives_totalGradient2_SOA_CPU_GPU(grid_grad2_a_gpu, grid_grad2_b_gpu, grid_grad2_c_gpu, grid_grad2_d_gpu, frame_gpu, lens_kernel, nhalos, dx, dy, nbgridcells_x, nbgridcells_y, istart, jstart); // mapfunction(map_gpu,grid_grad2_a_gpu, grid_grad2_b_gpu, grid_grad2_c_gpu, grid_grad2_d_gpu,dl0s,dos,dol,cosmo->h,z,nbgridcells_x,nbgridcells_y,frame); //cudasafe(cudaGetLastError(), "module_potentialDerivative_totalGradient_SOA_CPU_GPU"); cudaDeviceSynchronize(); // cudasafe(cudaMemcpy( map, map_gpu, (nbgridcells_x)*(nbgridcells_y)*sizeof(type_t), cudaMemcpyDeviceToHost )," --- Gradient2gpu.cu : Copy source_a_gpu: " ); // time += myseconds(); std::cout << " kernel time = " << time << " s." << std::endl; // Free GPU memory cudaFree(lens_gpu); cudaFree(type_gpu); cudaFree(lens_x_gpu); cudaFree(lens_y_gpu); cudaFree(b0_gpu); cudaFree(angle_gpu); cudaFree(epot_gpu); cudaFree(rcore_gpu); cudaFree(rcut_gpu); cudaFree(anglecos_gpu); cudaFree(anglesin_gpu); cudaFree(grid_grad2_a_gpu); cudaFree(grid_grad2_b_gpu); cudaFree(grid_grad2_c_gpu); cudaFree(grid_grad2_d_gpu); cudaFree(map_gpu); } //allows for resizing of the source void map_resizedgrid_GPU(map_gpu_function_t mapfunction, type_t *map,const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int nhalos, int mode_amp, type_t z, type_t dx, type_t dy, int nbgridcells_x, int nbgridcells_y, int istart, int jstart) { - int nBlocks_gpu = 0; - // Define the number of threads per block the GPU will use - cudaDeviceProp properties_gpu; - cudaGetDeviceProperties(&properties_gpu, 0); // Get properties of 0th GPU in use - - if (properties_gpu.maxThreadsDim[0]xmax - frame->xmin) / 6.; type_t resized_dy = (frame->ymax - frame->ymin) / 6.; resized_frame.xmax = frame->xmax - resized_dx; resized_frame.xmin = frame->xmin + resized_dx; resized_frame.ymax = frame->ymax - resized_dy; resized_frame.ymin = frame->xmin + resized_dy; grid_param *frame_gpu; Potential_SOA *lens_gpu,*lens_kernel; int *type_gpu; type_t *lens_x_gpu, *lens_y_gpu, *b0_gpu, *angle_gpu, *epot_gpu, *rcore_gpu, *rcut_gpu, *anglecos_gpu, *anglesin_gpu; type_t *grid_grad2_a_gpu, *grid_grad2_b_gpu , *grid_grad2_c_gpu, *grid_grad2_d_gpu, *map_gpu; type_t dl0s = module_cosmodistances_objectObject(lens->z[0], z, *cosmo); type_t dos = module_cosmodistances_observerObject(z, *cosmo); //select_ratio_function(std::string mode, const struct runmode_param* runmode, type_t dls, type_t ds) lens_gpu = (Potential_SOA *) malloc(sizeof(Potential_SOA)); lens_gpu->type = (int *) malloc(sizeof(int)); // Allocate variables on the GPU cudasafe(cudaMalloc( (void**)&(lens_kernel), sizeof(Potential_SOA)),"Gradient2gpu.cu : Alloc Potential_SOA: " ); cudasafe(cudaMalloc( (void**)&(type_gpu), nhalos*sizeof(int)),"Gradient2gpu.cu : Alloc type_gpu: " ); cudasafe(cudaMalloc( (void**)&(lens_x_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc x_gpu: " ); cudasafe(cudaMalloc( (void**)&(lens_y_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc y_gpu: " ); cudasafe(cudaMalloc( (void**)&(b0_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc b0_gpu: " ); cudasafe(cudaMalloc( (void**)&(angle_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc angle_gpu: " ); cudasafe(cudaMalloc( (void**)&(epot_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc epot_gpu: " ); cudasafe(cudaMalloc( (void**)&(rcore_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc rcore_gpu: " ); cudasafe(cudaMalloc( (void**)&(rcut_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc rcut_gpu: " ); cudasafe(cudaMalloc( (void**)&(anglecos_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc anglecos_gpu: " ); cudasafe(cudaMalloc( (void**)&(anglesin_gpu), nhalos*sizeof(type_t)),"Gradient2gpu.cu : Alloc anglesin_gpu: " ); cudasafe(cudaMalloc( (void**)&(frame_gpu), sizeof(grid_param)),"Gradient2gpu.cu : Alloc frame_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_a_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_a_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_b_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_b_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_c_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_c_gpu: " ); cudasafe(cudaMalloc( (void**)&(grid_grad2_d_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc source_d_gpu: " ); cudasafe(cudaMalloc( (void**)&(map_gpu), (nbgridcells_x) * (nbgridcells_y) *sizeof(type_t)),"Gradient2gpu.cu : Alloc map: " ); // Copy values to the GPU // cudasafe(cudaMemcpy(type_gpu,lens->type , nhalos*sizeof(int),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy type_gpu: " ); cudasafe(cudaMemcpy(lens_x_gpu,lens->position_x , nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy x_gpu: " ); cudasafe(cudaMemcpy(lens_y_gpu,lens->position_y , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy y_gpu: " ); cudasafe(cudaMemcpy(b0_gpu,lens->b0 , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2pu.cu : Copy b0_gpu: " ); cudasafe(cudaMemcpy(angle_gpu,lens->ellipticity_angle , nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy angle_gpu: " ); cudasafe(cudaMemcpy(epot_gpu, lens->ellipticity_potential, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy epot_gpu: " ); cudasafe(cudaMemcpy(rcore_gpu, lens->rcore, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy rcore_gpu: " ); cudasafe(cudaMemcpy(rcut_gpu, lens->rcut, nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy rcut_gpu: " ); cudasafe(cudaMemcpy(anglecos_gpu, lens->anglecos, nhalos*sizeof(type_t),cudaMemcpyHostToDevice ),"Gradient2gpu.cu : Copy anglecos: " ); cudasafe(cudaMemcpy(anglesin_gpu, lens->anglesin, nhalos*sizeof(type_t), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy anglesin: " ); cudasafe(cudaMemcpy(frame_gpu, frame, sizeof(grid_param), cudaMemcpyHostToDevice),"Gradient2gpu.cu : Copy fame_gpu: " ); // lens_gpu->type = type_gpu; lens_gpu->position_x = lens_x_gpu; lens_gpu->position_y = lens_y_gpu; lens_gpu->b0 = b0_gpu; lens_gpu->ellipticity_angle = angle_gpu; lens_gpu->ellipticity_potential = epot_gpu; lens_gpu->rcore = rcore_gpu; lens_gpu->rcut = rcut_gpu; lens_gpu->anglecos = anglecos_gpu; lens_gpu->anglesin = anglesin_gpu; // cudaMemcpy(lens_kernel, lens_gpu, sizeof(Potential_SOA), cudaMemcpyHostToDevice); // type_t time = -myseconds(); // module_potentialDerivatives_totalGradient2_SOA_CPU_GPU(grid_grad2_a_gpu, grid_grad2_b_gpu, grid_grad2_c_gpu, grid_grad2_d_gpu, frame_gpu, lens_kernel, nhalos, dx, dy, nbgridcells_x, nbgridcells_y, istart, jstart); // //mapfunction(map_gpu,grid_grad2_a_gpu, grid_grad2_b_gpu, grid_grad2_c_gpu, grid_grad2_d_gpu,dl0s,dos,z,mode_amp,nhalos,nbgridcells_x,nbgridcells_y); //amplif_grid_CPU_GPU(map_gpu,grid_grad2_a_gpu, grid_grad2_b_gpu, grid_grad2_c_gpu, grid_grad2_d_gpu,dl0s,z,mode_amp,nhalos,nbgridcells_x,nbgridcells_y); //cudasafe(cudaGetLastError(), "module_potentialDerivative_totalGradient_SOA_CPU_GPU"); cudaDeviceSynchronize(); // cudasafe(cudaMemcpy( map, map_gpu, (nbgridcells_x)*(nbgridcells_y)*sizeof(type_t), cudaMemcpyDeviceToHost )," --- Gradient2gpu.cu : Copy source_a_gpu: " ); // time += myseconds(); std::cout << " kernel time = " << time << " s." << std::endl; //printf("-----> %f %f \n",grid_grad_x[Nx], grid_grad_y[Ny]); // Free GPU memory cudaFree(lens_gpu); cudaFree(type_gpu); cudaFree(lens_x_gpu); cudaFree(lens_y_gpu); cudaFree(b0_gpu); cudaFree(angle_gpu); cudaFree(epot_gpu); cudaFree(rcore_gpu); cudaFree(rcut_gpu); cudaFree(anglecos_gpu); cudaFree(anglesin_gpu); cudaFree(grid_grad2_a_gpu); cudaFree(grid_grad2_b_gpu); cudaFree(grid_grad2_c_gpu); cudaFree(grid_grad2_d_gpu); cudaFree(map_gpu); } ////Map functions //Amplification NR 1 void amplif_1_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // //std::cerr << "Amplif 1 :" <>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } // __global__ void amplif_1_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) { type_t A,B,C; ellipse amp; type_t dlsds= dl0s/ds; //// int col = blockIdx.x*blockDim.x + threadIdx.x; int row = blockIdx.y*blockDim.y + threadIdx.y; // if ((row < nbgridcells) && (col < nbgridcells)) { int index = row*nbgridcells + col; A = 1. - grid_grad2_a[index]*dlsds; // 1 - DLS/DS * d2phixx B = - grid_grad2_b[index]*dlsds; // - DLS/DS * d2phixy C = 1. - grid_grad2_c[index]*dlsds; // 1 - DLS/DS * d2phiyy amp = formeli_HPC(A, B, C); ampli[index] = 1. / (amp.a * amp.b); //if(col == 0 and row == 0)printf(" GPUUUU: Grad %f %f %f ABC %f %f %f dlsds %f ab %f %f %f \n",grid_grad2_a[index]*dlsds,grid_grad2_b[index]*dlsds, grid_grad2_c[index]*dlsds,A,B,C,dlsds,amp.a,amp.b,(A - C)*(A - C) + 4*B*B); } } //Amplification NR 2 void amplif_2_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // amplif_2_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } // __global__ void amplif_2_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) { type_t A,B,C; ellipse amp; type_t dlsds= dl0s/ds; //// int col = blockIdx.x*blockDim.x + threadIdx.x; int row = blockIdx.y*blockDim.y + threadIdx.y; // if ((row < nbgridcells) && (col < nbgridcells)) { int index = row*nbgridcells + col; A = 1. - grid_grad2_a[index]*dlsds; // 1 - DLS/DS * d2phixx B = - grid_grad2_b[index]*dlsds; // - DLS/DS * d2phixy C = 1. - grid_grad2_c[index]*dlsds; // 1 - DLS/DS * d2phiyy amp = formeli_HPC(A, B, C); ampli[index] = 1. / fabs(amp.a * amp.b); } } //Amplification NR 3 void amplif_3_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // amplif_3_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } // __global__ void amplif_3_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) { type_t A,B,C; ellipse amp; type_t dlsds= dl0s/ds; //// int col = blockIdx.x*blockDim.x + threadIdx.x; int row = blockIdx.y*blockDim.y + threadIdx.y; // if ((row < nbgridcells) && (col < nbgridcells)) { int index = row*nbgridcells + col; A = 1. - grid_grad2_a[index]*dlsds; // 1 - DLS/DS * d2phixx B = - grid_grad2_b[index]*dlsds; // - DLS/DS * d2phixy C = 1. - grid_grad2_c[index]*dlsds; // 1 - DLS/DS * d2phiyy amp = formeli_HPC(A, B, C); ampli[index] = -2.5 * log10(fabs(amp.a * amp.b)); } } //Amplification NR 4 void amplif_4_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // amplif_4_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } // __global__ void amplif_4_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) { //type_t kappa,ga1,ga2,gam,gp; //// type_t dlsds = dl0s/ds; int col = blockIdx.x*blockDim.x + threadIdx.x; int row = blockIdx.y*blockDim.y + threadIdx.y; // if ((row < nbgridcells) && (col < nbgridcells)) { int index = row*nbgridcells + col; type_t kappa = (grid_grad2_a[index] + grid_grad2_c[index])*dlsds / 2.; type_t ga1 = (grid_grad2_a[index] - grid_grad2_c[index])*dlsds / 2.; type_t ga2 = grid_grad2_b[index]*dlsds; type_t gam = sqrt(ga1 * ga1 + ga2 * ga2); type_t gp = gam / (1 - kappa); ampli[index] = (1 - kappa) * (1 + gp * gp) / (1 - gp * gp); } } //Amplification NR 5 void amplif_5_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // amplif_5_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } // __global__ void amplif_5_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) { //// int col = blockIdx.x*blockDim.x + threadIdx.x; int row = blockIdx.y*blockDim.y + threadIdx.y; // if ((row < nbgridcells) && (col < nbgridcells)) { int index = row*nbgridcells + col; type_t kappa = (grid_grad2_a[index] + grid_grad2_c[index]) / 2.; ampli[index] = kappa; //if(col == 0 and row == 0)printf(" GPUUUU: Grad %f %f %f ",grid_grad2_a[index], grid_grad2_b[index],grid_grad2_c[index]); } } //Amplification NR 6 void amplif_6_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // amplif_6_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s, ds,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } // __global__ void amplif_6_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) { //// int col = blockIdx.x*blockDim.x + threadIdx.x; int row = blockIdx.y*blockDim.y + threadIdx.y; // if ((row < nbgridcells) && (col < nbgridcells)) { int index = row*nbgridcells + col; type_t ga1 = (grid_grad2_a[index] - grid_grad2_c[index]) / 2.; type_t ga2 = grid_grad2_b[index]; type_t gam = sqrt(ga1 * ga1 + ga2 * ga2); ampli[index] = gam; } } +//Shear NR 1 +void shear_1_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) +{ + int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks + int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; + // + //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); + // + dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); + dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); + // + cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); + // + shear_1_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); + cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); + // + cudaDeviceSynchronize(); + printf("GPU kernel done...\n"); +} +// +__global__ void shear_1_grid_GPU(type_t *shear,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) +{ + type_t A,B,C; + ellipse amp; + type_t dlsds= dl0s/ds; + //// + int col = blockIdx.x*blockDim.x + threadIdx.x; + int row = blockIdx.y*blockDim.y + threadIdx.y; + // + if ((row < nbgridcells) && (col < nbgridcells)) + { + + int index = row*nbgridcells + col; + A = 1. - grid_grad2_a[index]*dlsds; // 1 - DLS/DS * d2phixx + B = - grid_grad2_b[index]*dlsds; // - DLS/DS * d2phixy + C = 1. - grid_grad2_c[index]*dlsds; // 1 - DLS/DS * d2phiyy + amp = formeli_HPC(A, B, C); + shear[index] = (amp.a - amp.b) /2. ; + } +} +//Shear NR 2 +void shear_2_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) +{ + int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks + int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; + // + //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); + // + dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); + dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); + // + cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); + // + shear_2_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); + cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); + // + cudaDeviceSynchronize(); + printf("GPU kernel done...\n"); +} +// +__global__ void shear_2_grid_GPU(type_t *shear,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) +{ + type_t A,B,C; + ellipse amp; + type_t dlsds= dl0s/ds; + type_t q; + //// + int col = blockIdx.x*blockDim.x + threadIdx.x; + int row = blockIdx.y*blockDim.y + threadIdx.y; + // + if ((row < nbgridcells) && (col < nbgridcells)) + { + + int index = row*nbgridcells + col; + A = 1. - grid_grad2_a[index]*dlsds; // 1 - DLS/DS * d2phixx + B = - grid_grad2_b[index]*dlsds; // - DLS/DS * d2phixy + C = 1. - grid_grad2_c[index]*dlsds; // 1 - DLS/DS * d2phiyy + amp = formeli_HPC(A, B, C); + q = amp.a / amp.b; + shear[index] = fabs((q * q - 1.) / (q * q + 1.)); + } +} +//Shear NR 3 +void shear_3_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) +{ + int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks + int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; + // + //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); + // + dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); + dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); + // + cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); + // + shear_3_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); + cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); + // + cudaDeviceSynchronize(); + printf("GPU kernel done...\n"); +} +// +__global__ void shear_3_grid_GPU(type_t *shear,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) +{ + type_t dlsds= dl0s/ds; + //// + int col = blockIdx.x*blockDim.x + threadIdx.x; + int row = blockIdx.y*blockDim.y + threadIdx.y; + // + if ((row < nbgridcells) && (col < nbgridcells)) + { + + int index = row*nbgridcells + col; + shear[index] = (grid_grad2_a[index] - grid_grad2_c[index])*dlsds / 2. ; + //if(col == 0 and row == 0)printf(" GPUUUU: Grad %f %f %f ABC %f %f %f dlsds %f ab %f %f %f \n",grid_grad2_a[index]*dlsds,grid_grad2_b[index]*dlsds, grid_grad2_c[index]*dlsds,A,B,C,dlsds,amp.a,amp.b,(A - C)*(A - C) + 4*B*B); + } +} +//Shear NR 4 +void shear_4_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) +{ + int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks + int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; + // + //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); + // + dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); + dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); + // + //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); + // + //std::cerr << "Amplif 1 :" <>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dl0s,ds,z,nbgridcells_x); + cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); + // + cudaDeviceSynchronize(); + printf("GPU kernel done...\n"); +} +// +__global__ void shear_4_grid_GPU(type_t *shear,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t z,int nbgridcells) +{ + type_t dlsds= dl0s/ds; + //// + int col = blockIdx.x*blockDim.x + threadIdx.x; + int row = blockIdx.y*blockDim.y + threadIdx.y; + // + if ((row < nbgridcells) && (col < nbgridcells)) + { + + int index = row*nbgridcells + col; + shear[index] = grid_grad2_b[index] * dlsds ; + //if(col == 0 and row == 0)printf(" GPUUUU: Grad %f %f %f ABC %f %f %f dlsds %f ab %f %f %f \n",grid_grad2_a[index]*dlsds,grid_grad2_b[index]*dlsds, grid_grad2_c[index]*dlsds,A,B,C,dlsds,amp.a,amp.b,(A - C)*(A - C) + 4*B*B); + } +} + + //Mass NR 1 void mass_1_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // mass_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,1, dl0s,ds,h,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } //Mass NR 2 void mass_2_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); type_t dcrit = cH4piG * h / dl / dl0s * ds; // in g/cm^2 // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // mass_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dcrit, dl0s,ds,h,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } //Mass NR 3 void mass_3_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // type_t dx = (frame->xmax - frame->xmin)/(nbgridcells_x - 1); type_t dy = (frame->ymax - frame->ymin)/(nbgridcells_y - 1); //printf("grid_size_x = %d, grid_size_y = %d, nbgridcells_x = %d, nbgridcells_y = %d, istart = %d, jstart = %d (split)\n", GRID_SIZE_X, GRID_SIZE_Y, nbgridcells_x, nbgridcells_y, istart, jstart); type_t conv = MCRIT12 / h * dx * dy * dl / dl0s * ds; // in 10^12 M_sol/pixel // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // mass_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,conv, dl0s,ds,h,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } //Mass NR 4 void mass_4_grid_CPU_GPU(type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t ds, type_t dl, type_t h, type_t z, int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame) { int GRID_SIZE_X = (nbgridcells_x + BLOCK_SIZE_X - 1)/BLOCK_SIZE_X; // number of blocks int GRID_SIZE_Y = (nbgridcells_y + BLOCK_SIZE_Y - 1)/BLOCK_SIZE_Y; // type_t dcritA = cH0_4piG * h / dl / dl0s * ds; // in 10^12 M_sol/kpc^2 // dim3 threads(BLOCK_SIZE_X, BLOCK_SIZE_Y/1); dim3 grid (GRID_SIZE_X , GRID_SIZE_Y); // //printf("nhalos = %d, size of shared memory = %lf (split)\n", nhalos, (type_t) (8*nhalos + BLOCK_SIZE_X*BLOCK_SIZE_Y)*sizeof(type_t)); // cudaMemset(map, 0, nbgridcells_x*nbgridcells_y*sizeof(type_t)); // mass_grid_GPU<<>> (map,grid_grad2_a, grid_grad2_b,grid_grad2_c, grid_grad2_d,dcritA, dl0s,ds,h,z,nbgridcells_x); cudasafe(cudaGetLastError(), "amplif_grid_CPU_GPU"); // cudaDeviceSynchronize(); printf("GPU kernel done...\n"); } // __global__ void mass_grid_GPU(type_t *ampli,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t mult, type_t dls, type_t ds,type_t h, type_t z,int nbgridcells) { //// int col = blockIdx.x*blockDim.x + threadIdx.x; int row = blockIdx.y*blockDim.y + threadIdx.y; //type_t dlsds = dls / ds ; // if ((row < nbgridcells) && (col < nbgridcells)) { int index = row*nbgridcells + col; //if(col == 0 and row == 0)printf(" MASS GPU %f %f\n",grid_grad2_a[index], grid_grad2_c[index]); type_t ga1 = (grid_grad2_a[index] + grid_grad2_c[index]) * 0.5 ;//* dlsds ; ampli[index] = ga1*mult; } } __device__ struct ellipse formeli_HPC(type_t a, type_t b, type_t c) { struct ellipse eli; type_t e, delta, lambda, mu; // eq carateristique : det(M-xI) = 0 delta = (a - c)*(a - c) + 4*b*b; // 4*gamma^2 (cf phd_JPK eq 2.56) e = sqrt(delta); /*e is 2 * shear, ie 2*gamma*/ lambda = .5*(a + c + e); // 1 - k + gamma mu = .5*(a + c - e); // 1 - k - gamma eli.a = lambda; eli.b = mu; if (lambda != mu && fabs(b) > 1e-5) eli.theta = atan2(lambda - a, b); // cf phd_JPK eq 2.58, and // tan(theta)= ( -cos(2theta) +- 1 ) / sin(2theta) // ADDED by EJ 29/11/2007 else if ( a >= c ) // ellipse aligned along the major axis of magnification eli.theta = 0.; else eli.theta = acos(-1.) / 2.; // ellipse aligned along the minor axis of magnification return(eli); } diff --git a/src/grid_map_GPU.cuh b/src/grid_map_GPU.cuh index 9734fc4..a8d1b04 100644 --- a/src/grid_map_GPU.cuh +++ b/src/grid_map_GPU.cuh @@ -1,48 +1,53 @@ /** * @Author Christoph Schaefer, EPFL (christophernstrerne.schaefer@epfl.ch), Gilles Fourestey (gilles.fourestey@epfl.ch) * @date July 2017 * @version 0,1 * */ #ifndef GRID_MAP_GPU_CUH_ #define GRID_MAP_GPU_CUH_ //#include //#include "gradient2_GPU.cuh" #include "grid_gradient2_GPU.cuh" #include "module_cosmodistances.hpp" #include //#include ////type for mapping functions typedef void (*map_gpu_function_t) (type_t *map,type_t *grid_grad2_a,type_t *grid_grad2_b,type_t *grid_grad2_c,type_t *grid_grad2_d, type_t dl0s, type_t dos, type_t dl, type_t h, type_t z,int nbgridcells_x, int nbgridcells_y, const struct grid_param *frame); ////Map function selection map_gpu_function_t select_map_function(std::string mode,const struct runmode_param* runmode); ////Map functions //Amplification void amplif_1_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void amplif_2_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void amplif_3_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void amplif_4_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void amplif_5_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void amplif_6_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); //Mass void mass_1_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void mass_2_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void mass_3_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); void mass_4_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); +//Shear +void shear_1_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); +void shear_2_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); +void shear_3_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); +void shear_4_grid_CPU_GPU(type_t *,type_t *,type_t *,type_t *,type_t *, type_t , type_t, type_t , type_t , type_t ,int ,int ,const struct grid_param*); ////General mapping function void map_grid_GPU(map_gpu_function_t mapfunction, type_t *amplif,const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int Nlens,int nbgridcells, int mode_amp, type_t z ); void map_grid_GPU(map_gpu_function_t mapfunction, type_t *amplif,const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int Nlens, int mode_amp, type_t z, type_t dx, type_t dy, int nbgridcells_x, int nbgridcells_y, int istart = 0, int jstart = 0); void map_mass_grid_GPU(map_gpu_function_t mapfunction, type_t *map, const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int nhalos ,int nbgridcells,int mode_amp, type_t zm, type_t zs ); void map_mass_grid_GPU(map_gpu_function_t mapfunction, type_t *map,const struct cosmo_param *cosmo, const struct grid_param *frame, const struct Potential_SOA *lens, int nhalos, int mode_amp, type_t zm, type_t zs, type_t dx, type_t dy, int nbgridcells_x, int nbgridcells_y, int istart, int jstart); // #endif /* GRADIENTGPU_CUH_ */ diff --git a/src/module_readParameters.cpp b/src/module_readParameters.cpp index 26bbc13..67ab1eb 100644 --- a/src/module_readParameters.cpp +++ b/src/module_readParameters.cpp @@ -1,2925 +1,2932 @@ /** * @file module_readParameters.cpp * @Author Christoph Schaefer, EPFL (christophernstrerne.schaefer@epfl.ch) * @date July 2015 * @version 0,1 * @brief read parameters, images, sources and any other relevant information * * read parameter file * */ // Include //=========================================================================================================== #include #include #include #include #include #include #include #include #include "module_readParameters.hpp" #include "convert_coordinates.hpp" #include //=========================================================================================================== // Function definitions /** * * bayes.dat each # corresponds to one param available, it has to be in the same folder (can be changed) * * bayespot = array with bayes values * nparam = number of param available * nvalues = number of bayes potentials * */ void module_readParameters_preparebayes(int &nparam, int &nvalues){ nparam = 0; nvalues = 0; //Counting lines std::string line; std::ifstream IM("bayes.dat",std::ios::in); if ( IM ) { while( std::getline(IM,line) ) // Read every line { if ( strncmp(line.c_str(), "#", 1) == 0){ //Skipping commented lines nparam += 1; continue; } nvalues +=1; } } //nvalues -= 1; //exclude last empty line //nparam -=1; // Nparam is excluded //std::cerr << "nvalues :" << nvalues << "nparam :" << nparam << std::endl; IM.close(); } void module_readParameters_bayesmodels(double * bayespot, int nparam, int nvalues){ std::string line; std::ifstream IM("bayes.dat",std::ios::in); std::stringstream streamline; int j = 0; if ( IM ) { while( std::getline(IM,line) ) // Read every line { if ( strncmp(line.c_str(), "#", 1) == 0){ //Skipping commented lines continue; } streamline << line; for(int i = 0; i < nparam; i++){ streamline >> bayespot[j * nparam + i]; //IM >> bayespot[j * nparam + i]; if(j > nvalues){ fprintf(stderr, "ERROR: The bayes.dat file has an invalid line. please correct\nTypical problems: Whitespace on last bayes.dat line"); exit(-1); } //std::cerr << bayespot[j * nparam + i] << " " ; } streamline.clear(); //std::cerr << std::endl; j += 1; } } IM.close(); } /**setting potential using bayes[index] values. Does not support changing the bayes files config output defined by // Parameter constants in structure.h #define CX 0 #define CY 1 #define EPOT 2 #define EMASS 3 #define THETA 4 #define PHI 5 #define RC 6 #define B0 7 #define ALPHA 8 #define BETA 9 #define RCUT 10 #define MASSE 11 #define ZLENS 12 #define RCSLOPE 13 #define PMASS 14 */ void module_readParameters_setbayesmapmodels(const runmode_param* runmode, const cosmo_param* cosmology, const potentialoptimization* limit, potfile_param* potfile, Potential_SOA* lenses, double * bayespot, int nparam, int index){ int param_index = 2; int nhalo_index = 0; int SOA_index = 0; double DTR=acos(-1.)/180.; /* 1 deg in rad = pi/180 */ for(int i = 0; i < runmode->nhalos; i++){ //std::cerr << "Lenses SOA index: " << lenses->SOA_index[i] << std::endl; SOA_index = lenses->SOA_index[i]; if(limit[i].position.x.block >= 1){ lenses->position_x[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << " X : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; param_index++; } if(limit[i].position.y.block >= 1){ lenses->position_y[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << " Y : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; param_index++; } if(limit[i].ellipticity_potential.block >= 1){ lenses->ellipticity_potential[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << " epot : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; param_index++; } if(limit[i].ellipticity.block >= 1){ lenses->ellipticity[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << " emass : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; param_index++; } if(limit[i].ellipticity_angle.block >= 1){ lenses->ellipticity_angle[SOA_index] = bayespot[index*nparam+param_index]* DTR; //std::cerr << " X : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; lenses->anglecos[SOA_index] = cos(lenses->ellipticity_angle[SOA_index]); lenses->anglesin[SOA_index] = sin(lenses->ellipticity_angle[SOA_index]); param_index++; } if(limit[i].rcore.block >= 1){ lenses->rcore[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << " X : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; param_index++; } if(limit[i].vdisp.block >= 1){ lenses->vdisp[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << "VDISBLOC" << limit[i].vdisp.block <<" X : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << " " << nparam <= 1){ lenses->rcut[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << " X : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; param_index++; } if(limit[i].z.block >= 1){ lenses->z[SOA_index] = bayespot[index*nparam+param_index]; //std::cerr << " X : "<< index*nparam+param_index << " " << bayespot[index*nparam+param_index] << std::endl; param_index++; } module_readParameters_calculatePotentialparameter_SOA(lenses, SOA_index); } //std::cerr << "Potfile! " << runmode->potfile << std::endl; if(runmode->potfile != 0){ //Skip redshift image optimisation param_index += runmode->N_z_param; nhalo_index = runmode->nhalos; for(int ii = 0; ii < runmode->Nb_potfile; ii++){ //std::cerr << "ii " << ii << " " << param_index << " " << runmode->Nb_potfile << std::endl; //printf("DNDBSFB ircut %d\n",potfile->ircut); //Start potfile updating if(potfile[ii].ircut > 0){ potfile[ii].cut1 = bayespot[index*nparam+param_index]; //printf("cur %f ircut %d\n",potfile->cut,potfile->ircut); //std::cerr << index*nparam+param_index << " "<< bayespot[index*nparam+param_index] << std::endl; param_index++; } //std::cerr << "param " << param_index_pot << std::endl; if(potfile[ii].isigma > 0){ potfile[ii].sigma1 = bayespot[index*nparam+param_index]; //std::cerr << index*nparam+param_index << " "<< bayespot[index*nparam+param_index] << std::endl; param_index++; } setScalingRelations(runmode,cosmology,&potfile[ii],lenses,nhalo_index); nhalo_index += potfile[ii].npotfile; } /* for(int i = runmode->nhalos; i < runmode->nhalos + runmode->npotfile; i++){ param_index_pot = param_index; //std::cerr << "param " << param_index_pot << std::endl; //std::cerr << "param " << param_index_pot << std::endl; module_readParameters_calculatePotentialparameter_SOA(lenses, i); }*/ //update potential parameters } } //determine lens specific cosmological parameters needed for Kmap void module_readParameters_lens_dslds_calculation(const runmode_param* runmode, const cosmo_param* cosmo, Potential_SOA* lens){ type_t lens_z = lens->z[0]; type_t dl0s = module_cosmodistances_objectObject(lens_z, runmode->z_mass_s, *cosmo); type_t dos = module_cosmodistances_observerObject(runmode->z_mass_s, *cosmo); type_t dol = module_cosmodistances_observerObject(lens_z, *cosmo); lens->dlsds[0] = dl0s/dos; for(int ii = 0;ii n_tot_halos; ii++){ //std::cerr << ii << std::endl; //std::cerr << lens->z[ii] << std::endl; if(lens->z[ii] == lens_z){ lens->dlsds[ii] = dl0s/dos; } else{ lens_z = lens->z[ii]; dl0s = module_cosmodistances_objectObject(lens_z, runmode->z_mass_s, *cosmo); dos = module_cosmodistances_observerObject(runmode->z_mass_s, *cosmo); dol = module_cosmodistances_observerObject(lens_z, *cosmo); lens->dlsds[ii] = dl0s/dos; } //printf("dlsds %f\n",lens->dlsds[ii]); } } /** @brief This module function reads the cosmology parameters from the parameter file * * This module function reads the cosmology parameters from the parameter file. The given pointer is * updated with the results. * Read an infile: The structure of the file is as follows: | . | . | . |cosmology | model x <---- this is a value | H0 x | . . | . . | end | . | . | finish * * @param Cosmology cosmology * @param infile path to file */ void module_readParameters_readCosmology(std::string infile, cosmo_param &Cosmology) { std::string first, second, third, line1, line2; Cosmology.model = 1; Cosmology.H0 = 50; Cosmology.h = 1.; Cosmology.omegaM = 1.; Cosmology.omegaX = 0; Cosmology.wX = -1.; Cosmology.wa = 0.; Cosmology.curvature = 0.; std::ifstream IN(infile.c_str(),std::ios::in); // open file if ( IN ) { std::getline(IN,line1); std::istringstream read1(line1); // create a stream for the line read1 >> first; while(strncmp(first.c_str(), "fini",4) != 0 ) // read line by line until finish is reached { if ( strncmp(first.c_str(), "cosmolog", 8) == 0){ // if the line contains information about cosmology std::getline(IN,line2); // read the line word by word std::istringstream read2(line2); read2 >> second >> third; while(strncmp(second.c_str(), "end",3) != 0) // Go ahead until "end" is reached { if ( !strcmp(second.c_str(), "model") ) // set model of universe { Cosmology.model=atoi(third.c_str()); } else if ( !strcmp(second.c_str(), "H0") ) // set Hubble constant { Cosmology.H0=atof(third.c_str()); Cosmology.h = Cosmology.H0 / 50.; } else if ( !strcmp(second.c_str(), "omegaM") || !strcmp(second.c_str(), "omega") ) // set density of matter { Cosmology.omegaM=atof(third.c_str()); } else if ( !strcmp(second.c_str(), "omegaX") || !strcmp(second.c_str(), "lambda") ) // set cosmological constant { Cosmology.omegaX=atof(third.c_str()); } else if ( !strcmp(second.c_str(), "wX") || !strcmp(second.c_str(), "q") || !strcmp(second.c_str(), "w0") ) // set "q" for Model 2, "wX" for Model 3, "w0" for Model 4 { Cosmology.wX=atof(third.c_str()); } else if ( !strcmp(second.c_str(), "wa") || !strcmp(second.c_str(), "n") || !strcmp(second.c_str(), "delta") || !strcmp(second.c_str(), "w1") ) // set "n" for Model 2, "delta" for model 3, "w1" for model 4 { Cosmology.wa=atof(third.c_str()); } else if ( !strcmp(second.c_str(), "omegaK") ) // set universe curvature { Cosmology.curvature=atof(third.c_str()); } // read next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } // if a flat Universe if ( Cosmology.curvature == 0. ) { Cosmology.omegaX = 1 - Cosmology.omegaM; } else Cosmology.curvature = Cosmology.omegaM + Cosmology.omegaX - 1; } // read next line std::getline(IN,line1); std::istringstream read1(line1); read1 >> first; } IN.close(); } else { fprintf(stderr, "ERROR: file %s not found\n", infile.c_str()); // Exit if file was not found exit(-1); } } /** @brief This module function reads the number of sources, arclets, etc. in the parameter file. We need to know this to allocate the memory * * Function to read the number of multiple images and clumps * Check if there is 1 limit for each clump set * The program reads the number of potentials and limits defined, and checks whether there are the same number * Then it opens the imfile to read the numbers of images and set of images * Reads also the size of the multiple images area, the size of a pixel, and the runmode. * * @param infile path to file * @param runmode runmode parameter */ void read_runmode(std::istream &IN, struct runmode_param *runmode){ std::string first, second, third, fourth, fifth, line1, line2; //sscanf variables double in1, in2; //%lf in1, then (type_t)in1 ; std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third >> fourth; // Read in 4 words while (strncmp(second.c_str(), "end", 3)) // Read until we reach end { if ( !strcmp(second.c_str(), "nbgridcells") ) { sscanf(line2.c_str(), " %*s %d ", &runmode->nbgridcells); } if ( !strcmp(second.c_str(), "source") ) { char filename[FILENAME_SIZE]; sscanf(line2.c_str(), " %*s %d %s ", &runmode->source, &filename); runmode->sourfile = filename; } if ( !strcmp(second.c_str(), "image") ) { char filename[FILENAME_SIZE]; sscanf(line2.c_str(), " %*s %d %s ", &runmode->image, &filename); runmode->imagefile = filename; } if ( !strcmp(second.c_str(), "inverse") ) { sscanf(line2.c_str(), " %*s %d ", &runmode->inverse); } if ( !strcmp(second.c_str(), "mass") ) { char filename[FILENAME_SIZE]; sscanf(line2.c_str(), " %*s %d %d %lf %s", &runmode->mass, &runmode->mass_gridcells, &in1, &filename); runmode->z_mass = (type_t)in1; runmode->mass_name = filename; //runmode->z_mass_s =(type_t)in2; } if ( !strcmp(second.c_str(), "poten") ) { sscanf(line2.c_str(), " %*s %d %d %lf", &runmode->potential, &runmode->pot_gridcells, &in1); runmode->z_pot = (type_t)in1; } if ( !strcmp(second.c_str(), "dpl") ) { sscanf(line2.c_str(), " %*s %d %d %lf", &runmode->dpl, &runmode->dpl_gridcells, &in1); runmode->z_dpl = (type_t)in1; } if ( !strcmp(second.c_str(), "grid") ) { sscanf(line2.c_str(), " %*s %d %d %lf", &runmode->grid, &runmode->gridcells, &in1); runmode->zgrid = (type_t)in1; } if ( !strcmp(second.c_str(), "ampli") ) { char filename[FILENAME_SIZE]; sscanf(line2.c_str(), " %*s %d %d %lf %s", &runmode->amplif, &runmode->amplif_gridcells, &in1, &filename); runmode->z_amplif = (type_t)in1; runmode->amplif_name = filename; //std::cerr<ampli << << <shear, &runmode->shear_gridcells, &in1, &filename); + runmode->z_shear = (type_t)in1; + runmode->shear_name = filename; + } if ( !strcmp(second.c_str(), "arclets") ) { runmode->arclet = 1; // Not supported yet } if (!strcmp(second.c_str(), "reference")) { sscanf(line2.c_str(), "%*s %*d %lf %lf", &in1, &in2); runmode->ref_ra = (type_t)in1; runmode->ref_dec =(type_t)in2; //std::cerr << line2 << std::endl; //std::cout << "Reference: Ra " << runmode->ref_ra << " Dec:" << runmode->ref_dec <time); } if( !strncmp(second.c_str(), "Debug",5) ) // Read in if we are in debug mode { runmode->debug = 1; } // Read the next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } } void read_runmode_potential(std::istream &IN, int &numberPotentials){ std::string first, second, third, fourth, fifth, line1, line2; numberPotentials += 1; std::getline(IN,line2); std::istringstream read2(line2); double z(0); read2 >> second >> third; //std::cout << second << third << std::endl; while (strncmp(second.c_str(), "end", 3)) // Read until we reach end { if (!strcmp(second.c_str(), "z_lens")) // Get redshift { z=atof(third.c_str()); } // Read the next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } // Check if redshift from current halo was initialized if ( z == 0. ) { fprintf(stderr, "ERROR: A redshift is not defined for a potential \n"); exit(-1); } } void read_runmode_image(std::istream &IN, struct runmode_param *runmode){ std::string first, second, third, fourth, fifth, line1, line2; std::getline(IN,line2); std::istringstream read2(line2); double z(0); read2 >> second >> third; while (strncmp(second.c_str(), "end", 3)) // Read until we reach end { if ( !strcmp(second.c_str(), "multfile") ) { char filename[FILENAME_SIZE]; sscanf(line2.c_str(), " %*s %d %s ", &runmode->multi, &filename); runmode->imagefile = filename; } if ( !strcmp(second.c_str(), "z_m_limit") ) { runmode->N_z_param += 1 ; } //std::cout << runmode->multi << runmode->imagefile << std::endl; // Read the next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } } void read_runmode_source(std::istream &IN, struct runmode_param *runmode){ std::string first, second, third, fourth, fifth, line1, line2; std::getline(IN,line2); std::istringstream read2(line2); double z(0); read2 >> second >> third; while (strncmp(second.c_str(), "end", 3)) // Read until we reach end { if ( !strcmp(second.c_str(), "z_source") ) { sscanf(line2.c_str(), " %*s %lf ", &runmode->z_mass_s); } // Read the next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } } void read_runmode_potfile(std::istream &IN, struct runmode_param *runmode){ std::string first, second, third, fourth, fifth, line1, line2; runmode->potfile = 1; std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third >> fourth; // Read in 4 words while (strncmp(second.c_str(), "end", 3)) // Read until we reach end { if ( !strcmp(second.c_str(), "filein") ) { runmode->potfilename[runmode->Nb_potfile] = fourth; runmode->Nb_potfile = runmode->Nb_potfile + 1; break; } // Read the next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } } void read_runmode_countimages(struct runmode_param *runmode){ std::string line2; int old_j = 0; int j = 0; int imageIndex = 0; if (runmode->image == 1 or runmode->inverse == 1 or runmode->time >= 1 or runmode->multi >= 1){ std::string imageFile_name = runmode->imagefile; std::ifstream IM(imageFile_name.c_str(), std::ios::in); if ( IM ) { while( std::getline(IM,line2) ) // Read every line { if ( strncmp(line2.c_str(), "#", 1) == 0){ //Skipping commented lines continue; } else{ int j=atoi(line2.c_str()); if(j != old_j){ //If a new set has been reached, increase the nset iterator and update old j runmode->nsets +=1; old_j = j; } imageIndex++; } } } else{ fprintf(stderr, "ERROR: file %s not found\n", imageFile_name.c_str()); exit(-1); } runmode->nimagestot=imageIndex; IM.close(); } } void read_runmode_countsources(struct runmode_param *runmode){ std::string imageFile_name,line1, line2; if (runmode->source == 1 and runmode->image == 0 and runmode->multi == 0 ){ imageFile_name = runmode->sourfile; //printf("Source to image mode activated\n"); std::ifstream IM(imageFile_name.c_str(), std::ios::in); //printf(" Booo says hello again \n"); if ( IM ){ int i = 0; while( std::getline(IM,line1) ){ // Read until we reach the end i++; } runmode->nsets = i ; } else{ fprintf(stderr, "ERROR: file %s not found\n", imageFile_name.c_str()); exit(-1); } runmode->nimagestot= 0; // Important IM.close(); } } void read_runmode_countpotfile(struct runmode_param *runmode){ std::string imageFile_name,line1, line2; if (runmode->potfile == 1){ for(int ii = 0; ii < runmode->Nb_potfile; ii++){ imageFile_name = runmode->potfilename[ii]; std::ifstream IM(imageFile_name.c_str(), std::ios::in); if ( IM ){ int i = 0; while( std::getline(IM,line1) ){ // Read until we reach the end if ( line1[0] == '#' ){ continue; } i++; } runmode->npotfile[ii] = i ; } else{ fprintf(stderr, "ERROR: file %s not found\n", imageFile_name.c_str()); exit(-1); } IM.close(); runmode->n_tot_halos += runmode->npotfile[ii] ; } } } void module_readParameters_readRunmode(std::string infile, struct runmode_param *runmode) { std::string first, second, third, fourth, fifth, line1, line2; int Nsis(0), Npiemd(0); /// Set default values runmode->nbgridcells = 1000; runmode->source = 0; runmode->image = 0; runmode->N_z_param = 0; runmode->nsets = 0; runmode->nhalos = 0; runmode->multi = 0; runmode->amplif = 0; runmode->mass = 0; runmode->mass_gridcells = 1000; runmode->z_mass = 0.4; runmode->z_mass_s = 0.8; runmode->potential = 0; runmode->pot_gridcells = 1000; runmode->potfile = 0; //weird seg fault due to this line, haven't figured out why //runmode->npotfile = 0; runmode->z_pot = 0.8; runmode->dpl = 0; runmode->dpl_gridcells = 1000; runmode->z_dpl = 0.8; runmode->inverse = 0; runmode->arclet = 0; runmode->debug = 0; runmode->nimagestot = 0; runmode->nsets = 0; runmode->gridcells = 1000; //std::cerr << sizeof(*runmode) << std::endl; runmode->cline = 0; runmode->time = 0; runmode->inverse = 0; runmode->arclet = 0; runmode->ref_ra = 0; runmode->ref_dec = 0; runmode->Nb_potfile = 0; int j=0; std::string imageFile_name; int imageIndex=0; int numberPotentials=0, numberLimits=0; /*************************** read nhalos, imfile_name, pixel_size, multipleimagesarea_size and runmode from configuration file *********************/ std::ifstream IN(infile.c_str(), std::ios::in); if ( IN ) { std::getline(IN,line1); std::istringstream read1(line1); // create a stream for the line read1 >> first; // Read the first word //std::cout<cline = 1; } else if (!strncmp(first.c_str(), "source", 6)) { read_runmode_source(IN,runmode); } else if (!strncmp(first.c_str(), "potfile", 7)) { read_runmode_potfile(IN,runmode); } // read the next line std::getline(IN,line1); std::istringstream read1(line1); read1 >> first; //std::cout<nhalos=numberPotentials; runmode->n_tot_halos =numberPotentials; } else { fprintf(stderr, "ERROR: file %s not found\n", infile.c_str()); exit(-1); } // //getting nimage value (number of images), nset value (number of sources) and npotfile //if image or multi mode is activated get nimage and nset read_runmode_countimages(runmode); //if source mode is activated, get nset read_runmode_countsources(runmode); //if potfile mode, count number of potential in potfile read_runmode_countpotfile(runmode); //std::cerr <<"nsets: " <nsets <<" nhalos: " << runmode->nhalos << " nimagestot: " << runmode->nimagestot << " npotfile 1: " << runmode->npotfile[0] << " npotfile 2: " << runmode->npotfile[1] << " multi: " << runmode->multi<< std::endl; } /** @brief read the positions, redshifts and numbers of multiple images from the images file * * This module function reads in the strong lensing images * Output: image coordinates, image shapes (semi-minor, semi-major of ellipse and orientation angle), source redshifts, number of images per set * ** the images file must contain 1 x_center y_center major axis minor axis ellipticity angle redshift 1 . . . . . . 1 . . . . . . 2 . . . . . . 2 . . . . . . 2 . . . . . . 2 . . . . . . 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . nsetofimages . . . . . . At each line we store in j the index of the set of images to store the redshifts and the number of images of each set At each line we add an images in the position of images' array * * @param runmode runmode parameter * @param image array where images will be stored * @param nImagesSet array where the number of images per Set will be stored */ /// read the positions, redshifts and numbers of multiple images from the images file void module_readParameters_readImages(const struct runmode_param *runmode, struct galaxy image[], int nImagesSet[]) { std::string second, line1; int imageIndex=0; int j=0; //cast variables double cast_x,cast_y,cast_a,cast_b,cast_theta,cast_z,cast_mag; for(int i=0; insets; i++){ nImagesSet[i]=0; } /*********initialisation of nimages array*********/ for(int i=0; inimagestot; i++){ /* Initialise here the variables of the images*/ image[i].center.x = image[i].center.y = 0; image[i].shape.a = image[i].shape.b = image[i].shape.theta = 0; image[i].redshift = 0; } //printf("imagefile :%d \n", runmode->imagefile.c_str // Read in images std::ifstream IM(runmode->imagefile.c_str(),std::ios::in); if ( IM ) { int i = 0; int old_j = 1; int source_index = 0; while( std::getline(IM,line1) ) // Read until we reach the end { // Read in the parameters, * means we read in a parameter but do not store it sscanf(line1.c_str(), "%d %lf %lf %lf %lf %lf %lf %lf", &j, &cast_x, &cast_y, &cast_a, &cast_b, &cast_theta, &cast_z, &cast_mag); //Casting image[i].center.x =(type_t)cast_x; image[i].center.y =(type_t)cast_y; image[i].shape.a =(type_t)cast_a; image[i].shape.b =(type_t)cast_b; image[i].shape.theta =(type_t)cast_theta; image[i].redshift =(type_t)cast_z; image[i].mag =(type_t)cast_mag; //Variables int j=atoi(line1.c_str()); if(j != old_j){ //If a new set has been reached, increase the nset iterator and update old j source_index +=1; old_j = j; } nImagesSet[source_index]++; // Increase the counter of the number of system for system with number j-1 by 1 imageIndex++; i++; } } else{ std::cout << "Error : file " << runmode->imagefile << " not found" << std::endl; exit(-1); } IM.close(); } /** @brief read the positions, redshifts and numbers of multiple sources from the sources file * * This module function reads in the strong lensing sources * Output: sources coordinates, sources shapes (semi-minor, semi-major of ellipse and orientation angle), source redshifts * * the images file must contain x_center y_center major axis minor axis ellipticity angle redshift 1 . . . . . . 2 . . . . . . 3 . . . . . . * * @param runmode runmode parameter * @param source array where sources will be stored */ /// read the positions, redshifts and numbers of multiple images from the images file void module_readParameters_readSources(struct runmode_param *runmode, struct galaxy source[]) { std::string second, line1; int j=0; //cast variables double cast_x,cast_y,cast_a,cast_b,cast_theta,cast_z,cast_mag; //Calculating nset std::ifstream IM(runmode->sourfile.c_str(),std::ios::in); if ( IM ){ int i = 0; while( std::getline(IM,line1) ){ // Read until we reach the end i++; } runmode->nsets = i ; } else{ std::cout << "Error : file " << runmode->sourfile << " not found" << std::endl; exit(-1); } IM.close(); /*********initialisation of nimages array*********/ for(int i=0; insets; i++){ /* Initialise here the variables of the images*/ source[i].center.x = source[i].center.y = 0; source[i].shape.a = source[i].shape.b = source[i].shape.theta = 0; source[i].redshift = 0; source[i].mag = 0; } std::ifstream IM2(runmode->sourfile.c_str(),std::ios::in); // Read in images if ( IM2 ) { int i = 0; while( std::getline(IM2,line1) ) // Read until we reach the end { // Read in the parameters, * means we read in a parameter but do not store it sscanf(line1.c_str(), "%d %lf %lf %lf %lf %lf %lf %lf", &j, &cast_x, &cast_y, &cast_a, &cast_b, &cast_theta, &cast_z, &cast_mag); //Casting source[i].center.x =(type_t)cast_x; source[i].center.y =(type_t)cast_y; source[i].shape.a =(type_t)cast_a; source[i].shape.b =(type_t)cast_b; source[i].shape.theta =(type_t)cast_theta; source[i].redshift =(type_t)cast_z; source[i].mag =(type_t)cast_mag; i++; } } else{ std::cout << "Error : file " << runmode->sourfile << " not found" << std::endl; exit(-1); } IM2.close(); } /** @brief This module function reads the critical and caustic line information *@param infile path to file * @param cline cline parameter variable */ #if 0 void module_readParameters_CriticCaustic(std::string infile, cline_param *cline){ std::string first, second, third, line1, line2; //cast variables double cast_1; //Default value initialiasation cline->nplan = 0; for(int i =0; i < NPZMAX; ++i){ cline->cz[i] = 0; cline->dos[i] = 0; cline->dls[i] = 0; cline->dlsds[i] = 0; } cline->limitLow = 1; cline->dmax = 1; cline->limitHigh = 10; cline->nbgridcells = 1000; std::ifstream IN(infile.c_str(), std::ios::in); if ( IN ){ while(std::getline(IN,line1)){ std::istringstream read1(line1); // create a stream for the line read1 >> first; if ( strncmp(first.c_str(), "cline", 5) == 0){ // Read in runmode information std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; // Read in 4 words while (strncmp(second.c_str(), "end", 3)) // Read until we reach end { if ( !strcmp(second.c_str(), "nplan") ){ sscanf(third.c_str(), "%d", &cline->nplan); if ( cline->nplan > NPZMAX ){ cline->nplan = NPZMAX; } int j = 0; while ( read2 >> third ) { sscanf(third.c_str(), "%lf", &cast_1); cline->cz[j] =(type_t)cast_1; //printf(" zf %f \n",cline->cz[j]); j++; } } if ( !strcmp(second.c_str(), "dmax") ) { sscanf(third.c_str(), "%lf", &cast_1); cline->dmax =(type_t)cast_1; cline->xmax = cline->dmax; cline->xmin = -cline->dmax; cline->ymax = cline->dmax; cline->ymin = -cline->dmax; } if ( !strcmp(second.c_str(), "pas") || !strcmp(second.c_str(), "step") || !strcmp(second.c_str(), "limitLow") ) { sscanf(third.c_str(), "%lf", &cast_1); cline->limitLow =(type_t)cast_1; } if ( !strcmp(second.c_str(), "limitHigh") ) { sscanf(third.c_str(), "%lf", &cast_1); cline->limitHigh =(type_t)cast_1; } if ( !strcmp(second.c_str(), "nbgridcells") ) { sscanf(third.c_str(), "%d", &cline->nbgridcells); } // Read the next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } } } // closes while loop } // closes if(IN) IN.close(); } #endif /** @brief This module function reads the potfile information *@param infile path to file * @param potfile_param potfile parameter variable */ void module_readParameters_readpotfiles_param(std::string infile, potfile_param potfile[], cosmo_param cosmology){ std::string first, second, third, line1, line2; //cast variables double cast_1, cast_2; //Default value potfile initialiasation int Int_potfile = 0; std::ifstream IN(infile.c_str(), std::ios::in); if ( IN ){ while(std::getline(IN,line1)){ //std::getline(IN,line1); std::istringstream read1(line1); // create a stream for the line read1 >> first; if ( strncmp(first.c_str(), "potfile", 7) == 0){ // Read in potfile information potfile[Int_potfile].potid = 0; potfile[Int_potfile].ftype = 0; potfile[Int_potfile].type = 0; potfile[Int_potfile].zlens = 0; potfile[Int_potfile].mag0 = 0; potfile[Int_potfile].isigma = 0; potfile[Int_potfile].sigma = -1; potfile[Int_potfile].sigma1 = 0; potfile[Int_potfile].sigma2 = 0; potfile[Int_potfile].core = -1.0; potfile[Int_potfile].corekpc = -1; potfile[Int_potfile].ircut = 0; potfile[Int_potfile].cut1 = DBL_MAX; potfile[Int_potfile].cut2 = 0; potfile[Int_potfile].cutkpc1 = DBL_MAX; potfile[Int_potfile].cutkpc2 = 0; potfile[Int_potfile].islope = 0; potfile[Int_potfile].slope1 = 0; potfile[Int_potfile].slope2 = 0; potfile[Int_potfile].npotfile = 0; potfile[Int_potfile].ivdscat = 0; potfile[Int_potfile].vdscat1 =0; potfile[Int_potfile].vdscat2 = 0; potfile[Int_potfile].ircutscat = 0; std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; // Read in 4 words while (strncmp(second.c_str(), "end", 3)) // Read until we reach end { if ( !strcmp(second.c_str(), "filein") ) { sscanf(line2.c_str(), " %*s %d %s ", &potfile[Int_potfile].ftype, potfile[Int_potfile].potfile); //runmode->potfilename[Int_potfile] = fourth; } else if ( !strcmp(second.c_str(), "type") ) { sscanf(third.c_str(), "%d", &potfile[Int_potfile].type); } else if ( !strcmp(second.c_str(), "zlens") || !strcmp(second.c_str(), "z_lens")) { sscanf(third.c_str(), "%lf", &cast_1); potfile[Int_potfile].zlens =(type_t)cast_1; } else if (!strcmp(second.c_str(), "mag0") || !strcmp(second.c_str(), "r200")) { sscanf(third.c_str(), "%lf", &cast_1); potfile[Int_potfile].mag0 =(type_t)cast_1; } else if (!strcmp(second.c_str(), "select")) { sscanf(third.c_str(), "%d", &potfile[Int_potfile].select); } else if (!strcmp(second.c_str(), "core")) { sscanf(third.c_str(), "%lf", &cast_1); potfile[Int_potfile].core =(type_t)cast_1; } else if (!strcmp(second.c_str(), "corekpc")) { sscanf(third.c_str(), "%lf", &cast_1); potfile[Int_potfile].corekpc =(type_t)cast_1; } else if (!strcmp(second.c_str(), "cut")) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].ircut, &cast_1, &cast_2); potfile[Int_potfile].cut1 =(type_t)cast_1; potfile[Int_potfile].cut2 =(type_t)cast_2; } else if (!strcmp(second.c_str(), "cutkpc")) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].ircut, &cast_1, &cast_2); potfile[Int_potfile].cutkpc1 =(type_t)cast_1; potfile[Int_potfile].cutkpc2 =(type_t)cast_2; //std::cerr << potfile[Int_potfile].cutkpc1 << potfile[Int_potfile].cutkpc2 << std::endl; } else if (!strcmp(second.c_str(), "slope") || !strcmp(second.c_str(), "m200slope")) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].islope, &cast_1, &cast_2); potfile[Int_potfile].slope1 =(type_t)cast_1; potfile[Int_potfile].slope2 =(type_t)cast_2; } else if (!strcmp(second.c_str(), "sigma")) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].isigma, &cast_1, &cast_2); potfile[Int_potfile].sigma1 =(type_t)cast_1; potfile[Int_potfile].sigma2 =(type_t)cast_2; } else if (!strcmp(second.c_str(), "vdslope") || !strcmp(second.c_str(), "c200slope")) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].ivdslope, &cast_1, &cast_2); potfile[Int_potfile].vdslope1 =(type_t)cast_1; potfile[Int_potfile].vdslope2 =(type_t)cast_2; } else if (!strncmp(second.c_str(), "vdscat", 6)) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].ivdscat, &cast_1, &cast_2); potfile[Int_potfile].vdscat1 =(type_t)cast_1; potfile[Int_potfile].vdscat2 =(type_t)cast_2; } else if (!strncmp(second.c_str(), "rcutscat", 8)) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].ircutscat, &cast_1, &cast_2); potfile[Int_potfile].rcutscat1 =(type_t)cast_1; potfile[Int_potfile].rcutscat2 =(type_t)cast_2; } else if (!strcmp(second.c_str(), "a") || !strcmp(second.c_str(), "m200")) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].ia, &cast_1, &cast_2); potfile[Int_potfile].a1 =(type_t)cast_1; potfile[Int_potfile].a2 =(type_t)cast_2; } else if (!strcmp(second.c_str(), "b") || !strcmp(second.c_str(), "c200")) { sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile[Int_potfile].ib, &cast_1, &cast_2); potfile[Int_potfile].b1 =(type_t)cast_1; potfile[Int_potfile].b2 =(type_t)cast_2; } // Read the next line std::getline(IN,line2); std::istringstream read2(line2); read2 >> second >> third; } //Calculate vdisp, rcut and rcore //********************************************************************* // Set the Potfile current and limiting values //********************************************************************* if ( potfile[Int_potfile].ftype <= 4 ) { // Scale potfile SIGMA potfile[Int_potfile].sigma = potfile[Int_potfile].sigma1; // ... and potfile RCUT if ( potfile[Int_potfile].cut1 == DBL_MAX && potfile[Int_potfile].cutkpc1 != DBL_MAX ) { potfile[Int_potfile].cut1 = potfile[Int_potfile].cutkpc1 / (d0 / cosmology.h * module_cosmodistances_observerObject(potfile[Int_potfile].zlens,cosmology)); potfile[Int_potfile].cut2 = potfile[Int_potfile].cutkpc2 / (d0 / cosmology.h * module_cosmodistances_observerObject(potfile[Int_potfile].zlens,cosmology)); } potfile[Int_potfile].cut = potfile[Int_potfile].cut1; // ... and potfile RCORE if ( potfile[Int_potfile].core == -1.0 && potfile[Int_potfile].corekpc != -1 ) potfile[Int_potfile].core = potfile[Int_potfile].corekpc / (d0 / cosmology.h * module_cosmodistances_observerObject(potfile[Int_potfile].zlens,cosmology)); // ... and potfile RCUT SLOPE potfile[Int_potfile].slope = potfile[Int_potfile].slope1; // ... and potfile VDSLOPE potfile[Int_potfile].vdslope = potfile[Int_potfile].vdslope1; } //std::cerr << "Int potfile : " << Int_potfile << std::endl; Int_potfile += 1; } } // closes while loop } // closes if(IN) IN.close(); } /** @brief This module function loads the potential from the potfile into a lens *@param infile path to file * @param cline cline parameter variable */ /* void module_readParameters_readpotfiles(const runmode_param *runmode, potfile_param *potfile, Potential *lens){ std::string first, second, line1; double cast_1, cast_2; double aa,bb; double DTR=acos(-1.)/180.; //cast variables double cast_x,cast_y,cast_theta,cast_lum,cast_mag; potfile->reference_ra = potfile->reference_dec = 0; std::cerr << "Ref pot: "<< potfile->reference_ra << " " << potfile->reference_dec << std::endl; // Read in potentials std::ifstream IM(potfile->potfile,std::ios::in); if ( IM ) { int i = runmode->nhalos; while( std::getline(IM,line1) ) // Read until we reach the end { std::istringstream read1(line1); // create a stream for the line read1 >> first; // Skip commented lines with # //std::cerr << first << std::endl; if (!strncmp(first.c_str(), "#REFERENCE", 10) ){ sscanf(line1.c_str(), " %*s %*d%lf%lf", &cast_1, &cast_2); potfile->reference_ra = (type_t) cast_1; potfile->reference_dec = (type_t) cast_2; std::cerr << "Ref potfiles: "<< potfile->reference_ra << " " << potfile->reference_dec << std::endl; continue; } else if ( line1[0] == '#' ){ continue; } // Default initialisation of clump lens[i].type = potfile->type; lens[i].z = potfile->zlens; lens[i].ellipticity_potential = lens[i].ellipticity = 0.; lens[i].alpha = 0.; lens[i].rcut = 0.; lens[i].rcore = 0.; lens[i].rscale = 0.; lens[i].mag = 0.; lens[i].lum = 0.; lens[i].vdisp = 0.; lens[i].position.x = lens[i].position.y = 0.; lens[i].ellipticity_angle = 0.; lens[i].weight = 0; lens[i].exponent = 0; lens[i].einasto_kappacritic = 0; // Read a line of the catalog if ( potfile->ftype == 1 || potfile->ftype == 3 ) { sscanf( line1.c_str(), "%s%lf%lf%lf%lf%lf%lf%lf", &lens[i].name, &cast_x, &cast_y, &aa, &bb, &cast_theta, &cast_mag, &cast_lum); lens[i].ellipticity = (type_t) (aa*aa-bb*bb)/(aa*aa+bb*bb); if ( lens[i].ellipticity < 0 ) { fprintf( stderr, "ERROR: The potfile clump %s has a negative ellipticity.\n", lens[i].name ); exit(-1); } //goto NEXT; } //Casting lens[i].position.x =(type_t)cast_x; lens[i].position.y =(type_t)cast_y; lens[i].theta =(type_t)cast_theta; lens[i].lum =(type_t)cast_lum; lens[i].mag =(type_t)cast_mag; //general parameters lens[i].vdisp = potfile->sigma; lens[i].rcore = potfile->core; lens[i].rcut = potfile->cut; lens[i].ellipticity_angle = lens[i].theta* DTR; //Calculate parameters like b0, potential ellipticity and anyother parameter depending on the profile module_readParameters_calculatePotentialparameter(&lens[i]); //Variables potfile->npotfile++; i++; } } else{ std::cout << "Error : file " << potfile->potfile << " not found" << std::endl; exit(-1); } IM.close(); } */ void module_readParameters_readpotfiles_SOA(const runmode_param *runmode, const cosmo_param *cosmology, potfile_param potfile[], Potential_SOA *lens){ std::string first, second, line1; double cast_1, cast_2; double aa,bb; double DTR=acos(-1.)/180.; /* 1 deg in rad = pi/180 */ //cast variables double cast_x,cast_y,cast_theta,cast_lum,cast_mag, cast_name; int index = runmode->nhalos; // Read in potentials for(int jj = 0; jj < runmode->Nb_potfile ; jj++){ std::ifstream IM(potfile[jj].potfile,std::ios::in); potfile[jj].npotfile = 0; if ( IM ) { int i = index; while( std::getline(IM,line1) ) // Read until we reach the end { std::istringstream read1(line1); // create a stream for the line read1 >> first; // Skip commented lines with # //std::cerr << first << std::endl; if (!strncmp(first.c_str(), "#REFERENCE", 10)){ sscanf(line1.c_str(), " %*s %d%lf%lf", &potfile[jj].reference_mode, &cast_1, &cast_2); potfile[jj].reference_ra = (type_t) cast_1; potfile[jj].reference_dec = (type_t) cast_2; //std::cerr << "Ref pot: "<< potfile[jj].reference_ra << " " << potfile[jj].reference_dec << std::endl; continue; } // Skip commented lines with # else if ( line1[0] == '#' ) continue; cast_x = cast_y = cast_theta = cast_lum = cast_mag = cast_name = DBL_MAX; //std::cerr << "Turn " << i << std::endl; // Default initialisation of clump lens->type[i] = potfile[jj].type; lens->z[i] = potfile[jj].zlens; lens->ellipticity_potential[i] = lens->ellipticity[i] = 0.; lens->alpha[i] = 0.; lens->rcut[i] = 0.; lens->rcore[i] = 0.; lens->rscale[i] = 0.; lens->mag[i] = 0.; lens->lum[i] = 0.; lens->vdisp[i] = 0.; lens->position_x[i] = lens->position_y[i] = 0.; lens->ellipticity_angle[i] = 0.; lens->weight[i] = 0; lens->exponent[i] = 0; lens->einasto_kappacritic[i] = 0; //std::cerr << "Init finished "<< std::endl; //std::cerr << line1 << std::endl; // Read a line of the catalog if ( potfile[jj].ftype == 1 || potfile[jj].ftype == 3 ) { sscanf( line1.c_str(), "%*s%lf%lf%lf%lf%lf%lf%lf", &cast_x, &cast_y, &aa, &bb, &cast_theta, &cast_mag, &cast_lum); if (aa == bb){ lens->ellipticity[i] = 0; } else{ lens->ellipticity[i] = (type_t) (aa*aa-bb*bb)/(aa*aa+bb*bb); } //std::cerr << aa << bb <ellipticity[i] <ellipticity[i] < 0 ) { fprintf( stderr, "ERROR: The potfile clump %d has a negative ellipticity.\n", i ); exit(-1); } //goto NEXT; } else{ fprintf( stderr, "ERROR: Unknown ftype %d.\n", potfile[jj].ftype ); exit(-1); } //Casting lens->name[i] =(type_t)cast_name; lens->position_x[i] =(type_t)cast_x; lens->position_y[i] =(type_t)cast_y; //Cette partie me fait maaaaaaal au yeux .... buhu. Lenstool-HPC introduit une erreur sur le x a cause de la multiplication par un cosinus pour rester compatible lenstool convertXY_to_abs(&lens->position_x[i], &lens->position_y[i], potfile[jj].reference_mode, potfile[jj].reference_ra, potfile[jj].reference_dec ); convertXY_to_rela(&lens->position_x[i], &lens->position_y[i], potfile[jj].reference_mode, potfile[jj].reference_ra, potfile[jj].reference_dec ); //module_cosmodistances_relativecoordinates_XY( &lens->position_x[i], &lens->position_y[i], potfile[jj].reference_mode, potfile[jj].reference_ra, potfile[jj].reference_dec ); lens->theta[i] =(type_t)cast_theta; lens->lum[i] =(type_t)cast_lum; lens->mag[i] =(type_t)cast_mag; //general parameters //lens->vdisp[i] = potfile[jj].sigma; //lens->rcore[i] = potfile[jj].core; //lens->rcut[i] = potfile[jj].cut; lens->ellipticity_angle[i] = lens->theta[i]* DTR; lens->anglecos[i] = cos(lens->ellipticity_angle[i]); lens->anglesin[i] = sin(lens->ellipticity_angle[i]); if(cast_x == DBL_MAX or cast_y == DBL_MAX or cast_theta == DBL_MAX or cast_lum == DBL_MAX or cast_mag == DBL_MAX){ std::cerr << "Corrupted potfile: Potential " << i << " of potfile " << potfile[jj].potfile << " has missing values" << std::endl; exit(-1); } //Variables potfile[jj].npotfile++; i++; } setScalingRelations(runmode,cosmology,&potfile[jj],lens,index); index += potfile[jj].npotfile; //std::cerr << } else{ std::cout << "Error : file " << potfile->potfile << " not found" << std::endl; exit(-1); } IM.close(); } } void setScalingRelations(const runmode_param *runmode, const cosmo_param *cosmology, potfile_param *pot, Potential_SOA* lenses, int index){ //********************************************************************* // Check if the scaling relations are defined //********************************************************************* if ( pot->ftype == 1 ) { if ( pot->sigma1 == -1 ) { fprintf(stderr, "ERROR: potfile: sigma not defined\n"); exit(-1); } if ( pot->cutkpc1 == DBL_MAX && pot->cut1 == DBL_MAX ) { fprintf(stderr, "ERROR: potfile: cut length not defined\n"); exit(-1); } if ( pot->corekpc == -1 && pot->core == -1 ) { fprintf(stderr, "ERROR: potfile: core length not defined\n"); exit(-1); } } else{ fprintf(stderr, "ERROR: potfile: potfile type %d not supported\n", pot->ftype); exit(-1); } //********************************************************************* // Set the Potfile current and limiting values //********************************************************************* if ( pot->ftype <= 4 ) { // Scale potfile SIGMA pot->sigma = pot->sigma1; // ... and potfile RCUT if ( pot->cut1 == DBL_MAX && pot->cutkpc1 != DBL_MAX ) { pot->cut1 = pot->cutkpc1 / (d0 / cosmology->h * module_cosmodistances_observerObject(pot->zlens,*cosmology)); pot->cut2 = pot->cutkpc2 / (d0 / cosmology->h * module_cosmodistances_observerObject(pot->zlens,*cosmology)); } pot->cut = pot->cut1; // ... and potfile RCORE if ( pot->core == -1.0 && pot->corekpc != -1 ) pot->core = pot->corekpc / (d0 / cosmology->h * module_cosmodistances_observerObject(pot->zlens,*cosmology)); // ... and potfile RCUT SLOPE pot->slope = pot->slope1; // ... and potfile VDSLOPE pot->vdslope = pot->vdslope1; } // set potfile VDSCAT for all potfile scaling relations pot->vdscat = pot->vdscat1; // ... and potfile RCUTSCAT pot->rcutscat = pot->rcutscat1; for ( int i = index ; i < index + pot->npotfile ; i++ ){ if ( lenses->mag[i] != 0 ){ lenses->rcore[i] = pot->core * pow(10., 0.4 * (pot->mag0 - lenses->mag[i]) / 2.);} /* * Scale the sigma and rcut of a potfile clump according to the potfile parameters */ // loop over the potfile clumps to scale if ( pot->ftype <= 4 ) { //std::cerr << lenses->mag[i] << std::endl; if ( lenses->mag[i] != 0 ) { lenses->vdisp[i] = pot->sigma * pow(10., 0.4 * (pot->mag0 - lenses->mag[i]) / pot->vdslope); //std::cerr << " "<< pot->sigma1 <<" "<< pot->vdslope << " "<< lenses->mag[i] << " "<< pot->mag0 << " "<< lenses->vdisp[i] << " " << i << std::endl; /* The factor of 2 so that with slope1 = 4, we have * 2/slope1=1/2, then Brainerd, Blandford, Smail, 1996 */ lenses->rcut[i] = pot->cut * pow(10., 0.4 * (pot->mag0 - lenses->mag[i]) * 2. / pot->slope); } if ( pot->ivdscat != 0 ){ fprintf(stderr, "ERROR: potfile: ivdscat not supported yet\n"); exit(-1); } // Convert sigma to b0 //set_dynamics(i); if ( pot->ircutscat != 0 ){ fprintf(stderr, "ERROR: potfile: ircutscat not supported yet\n"); exit(-1); } } //Calculate parameters like b0, potential ellipticity and anyother parameter depending on the profile module_readParameters_calculatePotentialparameter_SOA(lenses,i); } } /** @brief read the information about arclets * !Not used! Will be reworked * This module function reads in the arclet images for weak lensing * @param Arclet file */ void module_readParameters_arclets(std::string arclets_filename, point arclets_position[], ellipse arclets_shape[], double arclets_redshift[]) { std::string second, line1; int j=0; //cast variables double cast_x,cast_y,cast_a,cast_b,cast_theta,cast_z; /******************** read the arclets file *******************/ /* the arclets file must contain id x_center y_center a b theta redshift line : 1 . . . . . . . 2 3 . . narclets */ std::ifstream IM(arclets_filename.c_str(),std::ios::in); if ( IM ) { while( std::getline(IM,line1)) // Read every line { // Read in parameters, * means we read the parameter but don't store it (*s) sscanf(line1.c_str(), "%*s %lf %lf %lf %lf %lf %lf", &cast_x, &cast_y, &cast_a, &cast_b, &cast_theta, &cast_z); //Casting arclets_position[j].x =(type_t)cast_x; arclets_position[j].y =(type_t)cast_y; arclets_shape[j].a =(type_t)cast_a; arclets_shape[j].b =(type_t)cast_b; arclets_shape[j].theta =(type_t)cast_theta; arclets_redshift[j] =(type_t)cast_z; j++; } IM.close(); } else { printf("Error : file %s not found\n",arclets_filename.c_str()); exit(-1); } } /** @brief This module function reads in if a parameter will be optimized by the MCMC or stay fixed. * * This module function reads in if a parameter will be optimized by the MCMC or stay fixed. * If it will be optimized, it specifies its minimum and maximum allowed values. Unless declared otherwise by the user, input values are fixed and won't be optimized. * * read an infile : | . | . |limit | x_center x x x x <--- these values contains : block min max accuracy | y_center x x x x if block=1 this is a free parameter, otherwise not | . . . . . min and max define the extremal value of this parameter | . . . . . accuracy is a criterium of convergence for the MCMC | end | . | . |limit | x_center x x x x | y_center x x x x | . . . . . | . . . . . | end | . | . |finish and fills the variables with these values * @param infile path to input file * @param host_potentialoptimization array where limits will be stored * @param nhalos number of mass distributions */ void module_readParameters_limit(std::string infile, struct potentialoptimization host_potentialoptimization[], int nhalos ) { std::string first, second, line1, line2; int i=0; //cast variables double cast_min,cast_max,cast_sigma; //double d1 = d0 / cosmology.h * module_cosmodistances_observerObject(lens_temp.z,cosmology); type_t DTR=acos(-1.)/180.; /* 1 deg in rad = pi/180 */ /*** initialize the block variables to zero (= not to be optimized) ***/ for(int index=0; index> first; //td::cerr << " 1: "<< first << std::endl; if (!strncmp(first.c_str(), "limit", 5)) // Read the limits { while(std::getline(IN,line2)) { std::istringstream read2(line2); read2 >> second; // Read in 1 word //std::cerr << " 2: "<< second << std::endl; if (strcmp(second.c_str(), "end") == 0) break; // stop reading at "end" and move to next potential limit section if (!strcmp(second.c_str(), "x_centre") || // Read in for x center !strcmp(second.c_str(), "x_center") ) { sscanf(line2.c_str(), "%*s %d %lf %lf %lf", &host_potentialoptimization[i].position.x.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].position.x.min =(type_t)cast_min; host_potentialoptimization[i].position.x.max =(type_t)cast_max; host_potentialoptimization[i].position.x.sigma =(type_t)cast_sigma; } else if (!strcmp(second.c_str(), "y_centre") || // Read in for y center !strcmp(second.c_str(), "y_center") ) { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].position.y.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].position.y.min =(type_t)cast_min; host_potentialoptimization[i].position.y.max =(type_t)cast_max; host_potentialoptimization[i].position.y.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "v_disp") ) // Read in for ellipticity { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].vdisp.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].vdisp.min =(type_t)cast_min; host_potentialoptimization[i].vdisp.max =(type_t)cast_max; host_potentialoptimization[i].vdisp.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "ellip_pot")) // Read in for ellipticity { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].ellipticity_potential.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].ellipticity_potential.min =(type_t)cast_min; host_potentialoptimization[i].ellipticity_potential.max =(type_t)cast_max; host_potentialoptimization[i].ellipticity_potential.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "ellipticitymass") || !strcmp(second.c_str(), "ellipticity") || !strcmp(second.c_str(), "ellipticite") || !strcmp(second.c_str(), "gamma") ) // Read in for ellipticity { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].ellipticity.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].ellipticity.min =(type_t)cast_min; host_potentialoptimization[i].ellipticity.max =(type_t)cast_max; host_potentialoptimization[i].ellipticity.sigma =(type_t)cast_sigma; } else if (!strcmp(second.c_str(), "ellipticity_angle") || !strcmp(second.c_str(), "angle_pos")) // Read in for ellipticity angle { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].ellipticity_angle.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].ellipticity_angle.min =(type_t)cast_min; host_potentialoptimization[i].ellipticity_angle.max =(type_t)cast_max; host_potentialoptimization[i].ellipticity_angle.sigma =(type_t)cast_sigma; host_potentialoptimization[i].ellipticity_angle.min *= DTR; host_potentialoptimization[i].ellipticity_angle.max *= DTR; host_potentialoptimization[i].ellipticity_angle.sigma *= DTR; } else if ( !strcmp(second.c_str(), "rcut") || !strcmp(second.c_str(), "cut_radius")) // Read in for r cut { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rcut.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].rcut.min =(type_t)cast_min; host_potentialoptimization[i].rcut.max =(type_t)cast_max; host_potentialoptimization[i].rcut.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "cut_radius_kpc")) // Read in for r cut { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rcut.block, &cast_min, &cast_max, &cast_sigma); std::cerr << "LIMIT: Rcut with kpc units not supported for the moment" << std::endl; } else if ( !strcmp(second.c_str(), "rcore") || !strcmp(second.c_str(), "core_radius")) // Read in for core radius { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rcore.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].rcore.min =(type_t)cast_min; host_potentialoptimization[i].rcore.max =(type_t)cast_max; host_potentialoptimization[i].rcore.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "core_radius_kpc")) // Read in for core radius { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rcore.block, &cast_min, &cast_max, &cast_sigma); std::cerr << "LIMIT: Rcore with kpc units not supported for the moment" << std::endl; } else if ( !strcmp(second.c_str(), "NFW_rs") || // Read in for NFW scale radius !strcmp(second.c_str(), "rscale") ) { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rscale.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].rscale.min =(type_t)cast_min; host_potentialoptimization[i].rscale.max =(type_t)cast_max; host_potentialoptimization[i].rscale.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "exponent") ) // Read in for exponent { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].exponent.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].exponent.min =(type_t)cast_min; host_potentialoptimization[i].exponent.max =(type_t)cast_max; host_potentialoptimization[i].exponent.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "alpha") ) // Read in for alpha { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].alpha.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].alpha.min =(type_t)cast_min; host_potentialoptimization[i].alpha.max =(type_t)cast_max; host_potentialoptimization[i].alpha.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "einasto_kappacritic") ) // Read in for critical kappa { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].einasto_kappacritic.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].einasto_kappacritic.min =(type_t)cast_min; host_potentialoptimization[i].einasto_kappacritic.max =(type_t)cast_max; host_potentialoptimization[i].einasto_kappacritic.sigma =(type_t)cast_sigma; } else if (!strcmp(second.c_str(), "virial_mass") || // Read in for virial mass !strcmp(second.c_str(), "masse") || !strcmp(second.c_str(), "m200") || !strcmp(second.c_str(), "mass")) { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].weight.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].weight.min =(type_t)cast_min; host_potentialoptimization[i].weight.max =(type_t)cast_max; host_potentialoptimization[i].weight.sigma =(type_t)cast_sigma; } else if ( !strcmp(second.c_str(), "z_lens") ) // Read in for redshift { sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].z.block, &cast_min, &cast_max, &cast_sigma); host_potentialoptimization[i].z.min =(type_t)cast_min; host_potentialoptimization[i].z.max =(type_t)cast_max; host_potentialoptimization[i].z.sigma =(type_t)cast_sigma; } } // end of inner while loop i++; // Move to next potential } } } IN.close(); } /** @brief This module function reads in the potential form and its parameters (e.g. NFW). * * @param infile path to input file * @param lens array where mass distributions will be stored * @param nhalos number of mass distributions */ void module_readParameters_Potential(std::string infile, Potential lens[], int nhalos) { double DTR=acos(-1.)/180.; /* 1 deg in rad = pi/180 */ Potential *ilens; std::string first, second, third, line1, line2; int i=0; std::ifstream IN(infile.c_str(), std::ios::in); if ( IN ) { while(std::getline(IN,line1)) { std::istringstream read1(line1); // create a stream for the line read1 >> first; if (!strncmp(first.c_str(), "potent", 6)) // Read in potential { /***********************************************/ ilens = &lens[i]; ilens->position.x = ilens->position.y = 0.; ilens->ellipticity = 0; ilens->ellipticity_potential = 0.; ilens->ellipticity_angle = 0.; ilens->rcut = 0.; ilens->rcore = 0; ilens->weight = 0; ilens->rscale = 0; ilens->exponent = 0; ilens->alpha = 0.; ilens->einasto_kappacritic = 0; ilens->z = 0; while(std::getline(IN,line2)) { std::istringstream read2(line2); read2 >> second >> third; if (strcmp(second.c_str(), "end") == 0) // Move to next potential and initialize it { if ( ilens->z == 0. ) // Check if redshift from current halo was initialized { fprintf(stderr, "ERROR: No redshift defined for potential %d\n", i); exit(-1); } break; // Break while loop and move to next potential } // Read in values if ( !strcmp(second.c_str(), "profil") || // Get profile !strcmp(second.c_str(), "profile") ) { if(!strcmp(third.c_str(), "PIEMD") || !strcmp(third.c_str(), "1") ) { ilens->type=1; strcpy(ilens->type_name,"PIEMD");//ilens->type_name="PIEMD"; } if(!strcmp(third.c_str(), "NFW") || !strcmp(third.c_str(), "2") ) { ilens->type=2; strcpy(ilens->type_name,"NFW");//ilens->type_name="NFW"; } if(!strcmp(third.c_str(), "SIES") || !strcmp(third.c_str(), "3") ) { ilens->type=3; strcpy(ilens->type_name,"SIES");//ilens->type_name="SIES"; } if(!strncmp(third.c_str(), "point", 5) || !strcmp(third.c_str(), "4") ) { ilens->type=4; strcpy(ilens->type_name,"point");//ilens->type_name="point"; } if(!strcmp(third.c_str(), "SIE") || !strcmp(third.c_str(), "5") ) { ilens->type=5; strcpy(ilens->type_name,"SIE");//ilens->type_name="point"; } if(!strcmp(third.c_str(), "8") ) { ilens->type=8; strcpy(ilens->type_name,"PIEMD1");//ilens->type_name="point"; } if(!strcmp(third.c_str(), "81") ) { ilens->type=81; strcpy(ilens->type_name,"PIEMD81");//ilens->type_name="point"; } } else if (!strcmp(second.c_str(), "name")) // Get name of lens { sscanf(third.c_str(),"%s",ilens->name); } else if (!strcmp(second.c_str(), "x_centre") || // Get x center !strcmp(second.c_str(), "x_center") ) { ilens->position.x=atof(third.c_str()); //std::cout << "PositionX : " << std::setprecision(15) << ilens->position.x << std::endl; } else if (!strcmp(second.c_str(), "y_centre") || // Get y center !strcmp(second.c_str(), "y_center") ) { ilens->position.y=atof(third.c_str()); } else if ( !strcmp(second.c_str(), "ellipticitymass") || !strcmp(second.c_str(), "ellipticity") ) // Get ellipticity { ilens->ellipticity=atof(third.c_str()); //ilens->ellipticity=ilens->ellipticity/3.; } else if (!strcmp(second.c_str(), "ellipticity_angle") || !strcmp(second.c_str(), "angle_pos")) // Get ellipticity angle { ilens->ellipticity_angle=atof(third.c_str()); ilens->ellipticity_angle *= DTR; } else if ( !strcmp(second.c_str(), "rcore") || !strcmp(second.c_str(), "core_radius")) // Get core radius { ilens->rcore=atof(third.c_str()); } else if (!strcmp(second.c_str(), "rcut") || !strcmp(second.c_str(), "cut_radius")) // Get cut radius { ilens->rcut=atof(third.c_str()); } else if (!strcmp(second.c_str(), "NFW_rs") || // Get scale radius of NFW !strcmp(second.c_str(), "rscale")) { ilens->rscale=atof(third.c_str()); } else if (!strcmp(second.c_str(), "exponent") ) // Get exponent { ilens->exponent=atof(third.c_str()); } else if (!strcmp(second.c_str(), "alpha") ) // Get alpha { ilens->alpha=atof(third.c_str()); } else if (!strcmp(second.c_str(), "einasto_kappacritic") || // Get critical kappa !strcmp(second.c_str(), "kappacritic")) { ilens->einasto_kappacritic=atof(third.c_str()); } else if (!strcmp(second.c_str(), "z_lens")) // Get redshift { ilens->z=atof(third.c_str()); } else if (!strcmp(second.c_str(), "v_disp")) // Get Dispersion velocity { ilens->vdisp=atof(third.c_str()); } else if ( !strncmp(second.c_str(), "virial_mass", 6) || // Get virial mass !strcmp(second.c_str(), "masse") || !strcmp(second.c_str(), "m200") || !strcmp(second.c_str(), "mass") ) { ilens->weight=atof(third.c_str()); } } // closes inner while loop //Calculate parameters like b0, potential ellipticity and anyother parameter depending on the profile module_readParameters_calculatePotentialparameter(ilens); i++; // Set counter to next potential } // closes if loop } // closes while loop if(i==0){ printf("Parameter potential not found in the file %s",infile.c_str()); exit(-1); } if(i>1){ for(int j=1; j> first; if (!strncmp(first.c_str(), "potent", 6)) // Read in potential { while(std::getline(IN,line2)) { std::istringstream read2(line2); read2 >> second >> third; if (strcmp(second.c_str(), "end") == 0) // Move to next potential and initialize it { break; // Break while loop and move to next potential } if ( !strcmp(second.c_str(), "profil") || // Get profile !strcmp(second.c_str(), "profile") ) { if(!strcmp(third.c_str(), "PIEMD") || !strcmp(third.c_str(), "1") ) { ind=atoi(third.c_str()); N_type[ind] += 1; } else if(!strcmp(third.c_str(), "NFW") || !strcmp(third.c_str(), "2") ) { ind=atoi(third.c_str()); N_type[ind] += 1; } else if(!strcmp(third.c_str(), "SIES") || !strcmp(third.c_str(), "3") ) { ind=atoi(third.c_str()); N_type[ind] += 1; } else if(!strncmp(third.c_str(), "point", 5) || !strcmp(third.c_str(), "4") ) { ind=atoi(third.c_str()); N_type[ind] += 1; } else if(!strcmp(third.c_str(), "SIE") || !strcmp(third.c_str(), "5") ) { ind=atoi(third.c_str()); N_type[ind] += 1; } else if(!strcmp(third.c_str(), "8") ) //PIEMD { ind=atoi(third.c_str()); N_type[ind] += 1; } else if(!strcmp(third.c_str(), "81") ) //PIEMD81 { ind=atoi(third.c_str()); N_type[ind] += 1; //std::cerr << "Type First: " << ind << std::endl; } else if(!strcmp(third.c_str(), "14") ) //PIEMD81 { ind=atoi(third.c_str()); N_type[ind] += 1; //std::cerr << "Type First: " << ind << std::endl; } else{ printf( "ERROR: Unknown Lensprofile, Emergency stop\n"); exit (EXIT_FAILURE); } } } } } } IN.close(); IN.clear(); IN.open(infile.c_str(), std::ios::in); } void module_readParameters_PotentialSOA_direct(std::string infile, Potential_SOA *lens_SOA, int nhalos, int n_tot_halos, cosmo_param cosmology){ double DTR=acos(-1.)/180.; /* 1 deg in rad = pi/180 */ double core_radius_kpc = 0.; double cut_radius_kpc = 0.; int N_type[100]; int Indice_type[100]; int ind, initial_index; Potential lens_temp; //Init of lens_SOA lens_SOA->name = new type_t[n_tot_halos]; lens_SOA->type = new int[n_tot_halos]; lens_SOA->position_x = new type_t[n_tot_halos]; lens_SOA->position_y = new type_t[n_tot_halos]; lens_SOA->b0 = new type_t[n_tot_halos]; lens_SOA->vdisp = new type_t[n_tot_halos]; lens_SOA->ellipticity_angle = new type_t[n_tot_halos]; lens_SOA->ellipticity = new type_t[n_tot_halos]; lens_SOA->ellipticity_potential = new type_t[n_tot_halos]; lens_SOA->rcore = new type_t[n_tot_halos]; lens_SOA->rcut = new type_t[n_tot_halos]; lens_SOA->z = new type_t[n_tot_halos]; lens_SOA->anglecos = new type_t[n_tot_halos]; lens_SOA->anglesin = new type_t[n_tot_halos]; lens_SOA->mag = new type_t[n_tot_halos]; lens_SOA->lum = new type_t[n_tot_halos]; lens_SOA->weight = new type_t[n_tot_halos]; lens_SOA->exponent = new type_t[n_tot_halos]; lens_SOA->einasto_kappacritic = new type_t[n_tot_halos]; lens_SOA->rscale = new type_t[n_tot_halos]; lens_SOA->alpha = new type_t[n_tot_halos]; lens_SOA->theta = new type_t[n_tot_halos]; lens_SOA->dlsds = new type_t[n_tot_halos]; lens_SOA->SOA_index = new int[n_tot_halos]; //Used to store the initial index of lenses initial_index = 0; //Init of N_types and Indice_type (Number of lenses of a certain type) for(int i=0;i < 100; ++i){ N_type[i] = 0; Indice_type[i] = 0; } //First sweep through the runmode file to find N_type (number of types) read_potentialSOA_ntypes(infile,N_type); //Calcuting starting points for each type in lens array for(int i=1;i < 100; ++i){ Indice_type[i] = N_type[i]+Indice_type[i-1]; //printf("%d %d \n ",N_type[i], Indice_type[i]); } std::string first, second, third, line1, line2; std::ifstream IN(infile.c_str(), std::ios::in); if(IN){ while(std::getline(IN,line1)) { first = ""; std::istringstream read1(line1); // create a stream for the line read1 >> first; //std::cerr << " 1: " << first << std::endl; if (!strncmp(first.c_str(), "potent", 6)) // Read in potential { lens_temp.position.x = lens_temp.position.y = 0.; lens_temp.ellipticity = 0; lens_temp.ellipticity_potential = 0.; lens_temp.ellipticity_angle = 0.; lens_temp.vdisp = 0.; lens_temp.rcut = 0.; lens_temp.rcore = 0; lens_temp.b0 = 0; core_radius_kpc = 0.; cut_radius_kpc = 0; lens_temp.weight = 0; lens_temp.rscale = 0; lens_temp.exponent = 0; lens_temp.alpha = 0.; lens_temp.einasto_kappacritic = 0; lens_temp.z = 0; while(std::getline(IN,line2)) { //Init temp potential std::istringstream read2(line2); read2 >> second >> third; //std::cerr << line2 << std::endl; //std::cerr << " 2: " << second << std::endl; if (strcmp(second.c_str(), "end") == 0) // Move to next potential and initialize it { if ( lens_temp.z == 0. ) // Check if redshift from current halo was initialized { fprintf(stderr, "ERROR: No redshift defined for potential at position x: %f and y: %f \n", lens_temp.position.x , lens_temp.position.y); exit(-1); } break; // Break while loop and move to next potential } //Find profile if ( !strcmp(second.c_str(), "profil") || // Get profile !strcmp(second.c_str(), "profile") ) { lens_temp.type=atoi(third.c_str()); //std::cerr << lens_temp.type << std::endl; } else if (!strcmp(second.c_str(), "name")) // Get name of lens { sscanf(third.c_str(),"%s",lens_temp.name); } else if (!strcmp(second.c_str(), "x_centre") || // Get x center !strcmp(second.c_str(), "x_center") ) { lens_temp.position.x=atof(third.c_str()); //std::cout << "PositionX : " << std::setprecision(15) << lens_temp.position.x << std::endl; } else if (!strcmp(second.c_str(), "y_centre") || // Get y center !strcmp(second.c_str(), "y_center") ) { lens_temp.position.y=atof(third.c_str()); } else if ( !strcmp(second.c_str(), "ellipticitymass") || !strcmp(second.c_str(), "ellipticity") || !strcmp(second.c_str(), "ellipticite") ) // Get ellipticity { lens_temp.ellipticity=atof(third.c_str()); //lens_temp.ellipticity=lens_temp.ellipticity/3.; } else if (!strcmp(second.c_str(), "ellipticity_angle") || !strcmp(second.c_str(), "angle_pos")) // Get ellipticity angle { lens_temp.ellipticity_angle=atof(third.c_str()); lens_temp.ellipticity_angle *= DTR; } else if ( !strcmp(second.c_str(), "rcore") || !strcmp(second.c_str(), "core_radius")) // Get core radius { lens_temp.rcore=atof(third.c_str()); } else if (!strcmp(second.c_str(), "rcut") || !strcmp(second.c_str(), "cut_radius")) // Get cut radius { lens_temp.rcut=atof(third.c_str()); } else if ( !strcmp(second.c_str(), "core_radius_kpc")) // Get core radius { core_radius_kpc=atof(third.c_str()); } else if (!strcmp(second.c_str(), "cut_radius_kpc")) // Get cut radius { cut_radius_kpc=atof(third.c_str()); } else if (!strcmp(second.c_str(), "NFW_rs") || // Get scale radius of NFW !strcmp(second.c_str(), "rscale")) { lens_temp.rscale=atof(third.c_str()); } else if (!strcmp(second.c_str(), "exponent") ) // Get exponent { lens_temp.exponent=atof(third.c_str()); } else if (!strcmp(second.c_str(), "alpha") ) // Get alpha { lens_temp.alpha=atof(third.c_str()); } else if (!strcmp(second.c_str(), "einasto_kappacritic") || // Get critical kappa !strcmp(second.c_str(), "kappacritic")) { lens_temp.einasto_kappacritic=atof(third.c_str()); } else if (!strcmp(second.c_str(), "z_lens")) // Get redshift { lens_temp.z=atof(third.c_str()); //std::cerr << lens_temp.z << std::endl; } else if (!strcmp(second.c_str(), "v_disp")) // Get Dispersion velocity { lens_temp.vdisp=atof(third.c_str()); //std::cerr << "vdisp : "<< third << " " << lens_temp.vdisp << std::endl; } else if ( !strncmp(second.c_str(), "virial_mass", 6) || // Get virial mass !strcmp(second.c_str(), "masse") || !strcmp(second.c_str(), "m200") || !strcmp(second.c_str(), "mass") ) { lens_temp.weight=atof(third.c_str()); } } // closes inner while loop // Converting distance in kpc to arcsec. double d1 = d0 / cosmology.h * module_cosmodistances_observerObject(lens_temp.z,cosmology); //printf(" D1 HPC : %f %f %f %f\n",d1, d0,cosmology.h,lens_temp.z ); // Set rcore value in kpc or in arcsec. if ( core_radius_kpc != 0. ) lens_temp.rcore = core_radius_kpc / d1; else core_radius_kpc = lens_temp.rcore * d1; // Set rcut value in kpc or in arcsec. if ( cut_radius_kpc != 0. ) { //std::cerr << "d1 " << d1 << std::endl; lens_temp.rcut = cut_radius_kpc / d1;} else cut_radius_kpc = lens_temp.rcut * d1; //Calculate parameters like b0, potential ellipticity and anyother parameter depending on the profile module_readParameters_calculatePotentialparameter(&lens_temp); //assign value to SOA //std::cerr << "Type + indice :" << lens_temp.type << Indice_type[lens_temp.type-1] << std::endl; if(Indice_type[lens_temp.type-1] type[ind] = lens_temp.type; lens_SOA->position_x[ind] = lens_temp.position.x; lens_SOA->position_y[ind] = lens_temp.position.y; lens_SOA->b0[ind] = lens_temp.b0; lens_SOA->vdisp[ind] = lens_temp.vdisp; lens_SOA->ellipticity_angle[ind] = lens_temp.ellipticity_angle; lens_SOA->ellipticity[ind] = lens_temp.ellipticity; lens_SOA->ellipticity_potential[ind] = lens_temp.ellipticity_potential; lens_SOA->rcore[ind] = lens_temp.rcore; lens_SOA->rcut[ind] = lens_temp.rcut; lens_SOA->z[ind] = lens_temp.z; lens_SOA->anglecos[ind] = cos(lens_temp.ellipticity_angle); lens_SOA->anglesin[ind] = sin(lens_temp.ellipticity_angle); //Store new index for bayes map purposes lens_SOA->SOA_index[initial_index] = ind; initial_index += 1; Indice_type[lens_temp.type-1] += 1; } } // closes if loop } // closes while loop } IN.close(); } /** @brief This module function converts potentials to potentieal_SOA (AOS to SOA conversion) * */ void module_readParameters_PotentialSOA(std::string infile, Potential *lens, Potential_SOA *lens_SOA, int nhalos){ lens_SOA->type = new int[nhalos]; lens_SOA->position_x = new type_t[nhalos]; lens_SOA->position_y = new type_t[nhalos]; lens_SOA->b0 = new type_t[nhalos]; lens_SOA->ellipticity_angle = new type_t[nhalos]; lens_SOA->ellipticity = new type_t[nhalos]; lens_SOA->ellipticity_potential = new type_t[nhalos]; lens_SOA->rcore = new type_t[nhalos]; lens_SOA->rcut = new type_t[nhalos]; lens_SOA->z = new type_t[nhalos]; lens_SOA->anglecos = new type_t[nhalos]; lens_SOA->anglesin = new type_t[nhalos]; int N_type[100]; int Indice_type[100]; int ind; for(int i=0;i < 100; ++i){ N_type[i] = 0; Indice_type[i] = 0; } for (int i = 0; i < nhalos; ++i){ N_type[lens[i].type] += 1; } for(int i=1;i < 100; ++i){ Indice_type[i] = N_type[i]+Indice_type[i-1]; //printf("%d %d \n ",N_type[i], Indice_type[i]); } for (int i = 0; i < nhalos; ++i){ if(Indice_type[lens[i].type-1] type[ind] = lens[i].type; lens_SOA->position_x[ind] = lens[i].position.x; //std::cerr << lens_SOA[1].position_x[*i_point] << " " << lens[i].position.x << std::endl; lens_SOA->position_y[ind] = lens[i].position.y; lens_SOA->b0[ind] = lens[i].b0; lens_SOA->ellipticity_angle[ind] = lens[i].ellipticity_angle; lens_SOA->ellipticity[ind] = lens[i].ellipticity; lens_SOA->ellipticity_potential[ind] = lens[i].ellipticity_potential; lens_SOA->rcore[ind] = lens[i].rcore; lens_SOA->rcut[ind] = lens[i].rcut; lens_SOA->z[ind] = lens[i].z; lens_SOA->anglecos[ind] = cos(lens[i].ellipticity_angle); lens_SOA->anglesin[ind] = sin(lens[i].ellipticity_angle); Indice_type[lens[i].type-1] += 1; } } //printf("Bla anglecos = %f\n", lens_SOA->anglecos[0]); } /** @brief This module function calculates profile depended information like the impactparameter b0 and the potential ellipticity epot * * @param lens: mass distribution for which to calculate parameters */ void module_readParameters_calculatePotentialparameter(Potential *lens){ switch (lens->type) { case(5): /*Elliptical Isothermal Sphere*/ //impact parameter b0 lens->b0 = 4* pi_c2 * lens->vdisp * lens->vdisp ; //ellipticity_potential lens->ellipticity_potential = lens->ellipticity/3 ; break; case(8): /* PIEMD */ //impact parameter b0 lens->b0 = 6.*pi_c2 * lens->vdisp * lens->vdisp; //ellipticity_parameter if ( lens->ellipticity == 0. && lens->ellipticity_potential != 0. ){ // emass is (a2-b2)/(a2+b2) lens->ellipticity = 2.*lens->ellipticity_potential / (1. + lens->ellipticity_potential * lens->ellipticity_potential); //printf("1 : %f %f \n",lens->ellipticity,lens->ellipticity_potential); } else if ( lens->ellipticity == 0. && lens->ellipticity_potential == 0. ){ lens->ellipticity_potential = 0.00001; //printf("2 : %f %f \n",lens->ellipticity,lens->ellipticity_potential); } else{ // epot is (a-b)/(a+b) lens->ellipticity_potential = (1. - sqrt(1 - lens->ellipticity * lens->ellipticity)) / lens->ellipticity; //printf("3 : %f %f \n",lens->ellipticity,lens->ellipticity_potential); } break; case(81): /* PIEMD */ //impact parameter b0 lens->b0 = 6.*pi_c2 * lens->vdisp * lens->vdisp; //ellipticity_parameter if ( lens->ellipticity == 0. && lens->ellipticity_potential != 0. ){ // emass is (a2-b2)/(a2+b2) lens->ellipticity = 2.*lens->ellipticity_potential / (1. + lens->ellipticity_potential * lens->ellipticity_potential); //printf("1 : %f %f \n",lens->ellipticity,lens->ellipticity_potential); } else if ( lens->ellipticity == 0. && lens->ellipticity_potential == 0. ){ lens->ellipticity_potential = 0.00001; //printf("2 : %f %f \n",lens->ellipticity,lens->ellipticity_potential); } else{ // epot is (a-b)/(a+b) lens->ellipticity_potential = (1. - sqrt(1 - lens->ellipticity * lens->ellipticity)) / lens->ellipticity; //printf("3 : %f %f \n",lens->ellipticity,lens->ellipticity_potential); } break; case(14): lens->ellipticity_potential = lens->ellipticity / 3; break; default: printf( "ERROR: LENSPARA profil type of clump %s unknown : %d\n",lens->name, lens->type); exit (EXIT_FAILURE); break; }; } void module_readParameters_calculatePotentialparameter_SOA(Potential_SOA *lens, int ind){ switch (lens->type[ind]) { case(5): /*Elliptical Isothermal Sphere*/ //impact parameter b0 lens->b0[ind] = 4* pi_c2 * lens->vdisp[ind] * lens->vdisp[ind] ; //ellipticity_potential if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] != 0. ){ lens->ellipticity_potential[ind] = lens->ellipticity[ind]/3 ; } else{ lens->ellipticity[ind] = lens->ellipticity_potential[ind]*3; } break; case(8): /* PIEMD */ //impact parameter b0 lens->b0[ind] = 6.*pi_c2 * lens->vdisp[ind] * lens->vdisp[ind]; //ellipticity_parameter if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] != 0. ){ // emass is (a2-b2)/(a2+b2) lens->ellipticity[ind] = 2.*lens->ellipticity_potential[ind] / (1. + lens->ellipticity_potential[ind] * lens->ellipticity_potential[ind]); //printf("1 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]); } else if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] == 0. ){ lens->ellipticity_potential[ind] = 0.00001; //printf("2 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]); } else{ // epot is (a-b)/(a+b) lens->ellipticity_potential[ind] = (1. - sqrt(1 - lens->ellipticity[ind] * lens->ellipticity[ind])) / lens->ellipticity[ind]; //printf("3 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]); } break; case(81): /* PIEMD */ //impact parameter b0 lens->b0[ind] = 6.*pi_c2 * lens->vdisp[ind] * lens->vdisp[ind]; //ellipticity_parameter if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] != 0. ){ // emass is (a2-b2)/(a2+b2) lens->ellipticity[ind] = 2.*lens->ellipticity_potential[ind] / (1. + lens->ellipticity_potential[ind] * lens->ellipticity_potential[ind]); //printf("1 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]); } else if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] == 0. ){ lens->ellipticity_potential[ind] = 0.00001; //printf("2 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]); } else{ // epot is (a-b)/(a+b) lens->ellipticity_potential[ind] = (1. - sqrt(1 - lens->ellipticity[ind] * lens->ellipticity[ind])) / lens->ellipticity[ind]; //printf("3 : %f %f \n",lens->ellipticity[ind],lens->ellipticity_potential[ind]); } break; default: if ( lens->ellipticity[ind] == 0. && lens->ellipticity_potential[ind] != 0. ){ lens->ellipticity[ind] = lens->ellipticity_potential[ind]*3; } else{ lens->ellipticity_potential[ind] = lens->ellipticity[ind]/3 ; } break; }; } /** @brief This module function reads the grid information * * @param infile path to input file * @param grid array where grid information will be stored */ void module_readParameters_Grid(std::string infile, grid_param *grid) { std::string first, second, third, line1, line2; double cast_1; double dmax ; std::ifstream IN(infile.c_str(), std::ios::in); if ( IN ) { while(std::getline(IN,line1)) { std::istringstream read1(line1); // create a stream for the line read1 >> first; if (!strncmp(first.c_str(), "grille", 7) || !strncmp(first.c_str(), "frame", 5) || !strncmp(first.c_str(), "champ", 5)) // Read in potential { while(std::getline(IN,line2)) { std::istringstream read2(line2); read2 >> second >> third; if (strcmp(second.c_str(), "end") == 0) // Move to next potential and initialize it { break; // Break while loop and move to next potential } // Read in values if (!strcmp(second.c_str(), "dmax")) { sscanf(third.c_str(), "%lf", &dmax); //printf( "\t%s\t%lf\n", second.c_str(), dmax); //grid->dmax = (type_t) cast_1; grid->xmin = -dmax; grid->xmax = dmax; grid->ymin = -dmax; grid->ymax = dmax; grid->rmax = dmax * sqrt(2.); } else if (!strcmp(second.c_str(), "xmin")) { sscanf(third.c_str(), "%lf", &cast_1); grid->xmin = (type_t) cast_1; //printf("\t%s\t%lf\n", second.c_str(), grid->xmin); } else if (!strcmp(second.c_str(), "xmax")) { sscanf(third.c_str(), "%lf", &cast_1); grid->xmax = (type_t) cast_1; //printf("\t%s\t%lf\n", second.c_str(), grid->xmax); } else if (!strcmp(second.c_str(), "ymin")) { sscanf(third.c_str(), "%lf", &cast_1); grid->ymin = (type_t) cast_1; //printf( "\t%s\t%lf\n", second.c_str(), grid->ymin); } else if (!strcmp(second.c_str(), "ymax")) { sscanf(third.c_str(), "%lf", &cast_1); grid->ymax = (type_t) cast_1; //printf( "\t%s\t%lf\n", second.c_str(), grid->ymax); } } } // closes if loop } // closes while loop } // closes if(IN) IN.close(); } /* @brief Get number of sets for cleanlens mode * !!Not used. Will be reworked!! * * This module function reads in the "clean lens" mode sources: These sources are read in and lensed only once assuming a fixed mass model // Then we can see the predicted location of the images and compare with their real positions * * @param clean lens file, number of clean lens images * @return coordinates for each image, shape of each image, redshifts, number of sets, number of images per set */ // Get number of sets for cleanlens mode void module_readParameters_SingleLensingSourcesNumberSets(std::string infile, int &nsetofimages_cleanlens ) { std::string line1; int setNumber=0; nsetofimages_cleanlens=0; // Get number of sets of images std::ifstream IM(infile.c_str(),std::ios::in); if ( IM ) { while( std::getline(IM,line1)) { // Read in parameters sscanf(line1.c_str(), "%d", &setNumber); if(setNumber > nsetofimages_cleanlens) nsetofimages_cleanlens = setNumber; } } IM.close(); } /** @brief Prints out cosmology */ void module_readParameters_debug_cosmology(int DEBUG, cosmo_param cosmology) { if (DEBUG == 1) // If we are in debug mode { printf("\n\n####################### Summary of Input Parameters #######################\n"); printf("Cosmology: Cosmology model = %d, omegaM = %lf, omegaX = %lf, curvature = %lf, wX = %lf, wa = %lf, H0 = %lf, h = %lf\n", cosmology.model, cosmology.omegaM, cosmology.omegaX, cosmology.curvature, cosmology.wX, cosmology.wa, cosmology.H0, cosmology.h); } } /** @brief Prints out runmode */ void module_readParameters_debug_runmode(int DEBUG, runmode_param runmode) { if (DEBUG == 1) // If we are in debug mode { printf("Runmode: nhalos = %d, nsets = %d, nimagestot = %d, source = %d, image = %d, arclet = %d, cline = %d, inverse= %d, multi= %d ampli= %d DEBUG = %d\n", runmode.nhalos, runmode.nsets, runmode.nimagestot, runmode.source, runmode.image, runmode.arclet, runmode.cline, runmode.inverse, runmode.multi, runmode.amplif, runmode.debug); } } /** @brief Prints out images */ void module_readParameters_debug_image(int DEBUG, galaxy image[], int nImagesSet[], int nsets) { if (DEBUG == 1) // If we are in debug mode { int index = 0; for ( int i = 0; i < nsets; ++i){ for( int j = 0; j < nImagesSet[i]; ++j){ printf("Image [%d]: x = %lf, y = %lf, shape: a = %f, b = %f, theta = %lf, redshift = %lf, nImagesSet = %d,\n",index , image[index].center.x, image[index].center.y, image[index].shape.a, image[index].shape.b, image[index].shape.theta, image[index].redshift, nImagesSet[i]); index +=1; //printf( "%d \n", index ); } } } } /** @brief Prints out source */ void module_readParameters_debug_source(int DEBUG, galaxy source[], int nsets) { if (DEBUG == 1) // If we are in debug mode { for ( int i = 0; i < nsets; ++i){ printf("Source[%d]: x = %lf, y = %lf, shape: a = %f, b = %f, theta = %lf, redshift = %lf. \n\n", i, source[i].center.x, source[i].center.y, source[i].shape.a, source[i].shape.b, source[i].shape.theta, source[i].redshift); } } } /** @brief Prints out massdistributions */ void module_readParameters_debug_potential(int DEBUG, Potential lenses[], int nhalos) { if (DEBUG == 1) // If we are in debug mode { for ( int i = 0; i < nhalos; ++i){ printf("Potential[%d]: x = %f, y = %f, vdisp = %f, type = %d, type_name = %s, name = %s,\n \t ellipticity = %f, ellipticity_pot = %f, ellipticity angle (radians) = %f, rcore = %f, rcut = %f,\n \t rscale = %f, exponent = %f, alpha = %f, einasto kappa critical = %f, z = %f\n", i,lenses[i].position.x, lenses[i].position.y, lenses[i].vdisp, lenses[i].type, lenses[i].type_name, lenses[i].name, lenses[i].ellipticity, lenses[i].ellipticity_potential, lenses[i].ellipticity_angle, lenses[i].rcore, lenses[i].rcut, lenses[i].rscale, lenses[i].exponent, lenses[i].alpha, lenses[i].einasto_kappacritic, lenses[i].z); } } } void module_readParameters_debug_potential_SOA(int DEBUG, Potential_SOA lenses, int nhalos) { if (DEBUG == 1) // If we are in debug mode { for ( int i = 0; i < nhalos; ++i){ printf("Potential[%d]: x = %f, y = %f, b0 = %f, vdisp = %f, type = %d, \n \t ellipticity = %f, ellipticity_pot = %f, ellipticity angle (radians) = %f, rcore = %f, rcut = %f,\n \t z = %f\n, angle cos %f, sin %f \n", i,lenses.position_x[i], lenses.position_y[i], lenses.b0[i],lenses.vdisp[i], lenses.type[i], lenses.ellipticity[i], lenses.ellipticity_potential[i], lenses.ellipticity_angle[i], lenses.rcore[i], lenses.rcut[i], lenses.z[i], lenses.anglecos[i], lenses.anglesin[i]); } } } /** @brief Prints out potfile_param */ void module_readParameters_debug_potfileparam(int DEBUG, potfile_param *potfile) { if (DEBUG == 1) // If we are in debug mode { printf("Potfile: potid = %d, filename = %s, ftype = %d, type = %d, zlens = %f, mag0 = %f, sigma = %f,\n \t core = %f, corekpc = %f, ircut = %d, cut1 = %f, cut2 = %f,\n \t cutkpc1 = %f, cutkpc2 = %f, islope = %d, slope1 = %f, slope2 = %f\n", potfile->potid,potfile->potfile,potfile->ftype,potfile->type,potfile->zlens,potfile->mag0,potfile->sigma,potfile->core,potfile->corekpc,potfile->ircut,potfile->cut1,potfile->cut2,potfile->cutkpc1,potfile->cutkpc2,potfile->islope,potfile->slope1,potfile->slope2); } } /** @brief Prints out cline */ void module_readParameters_debug_criticcaustic(int DEBUG, cline_param cline) { if (DEBUG == 1) // If we are in debug mode { printf("Cline: Number of planes = %d ", cline.nplan); for( int i=0; i < cline.nplan; ++i){ printf(" z[%d] = %f, ",i ,cline.cz[i]); } printf(" dmax = %f, Low = %f, High= %f, Nb of gridcells %d \n", cline.dmax , cline.limitLow ,cline.limitHigh, cline.nbgridcells); } } /** @brief Prints out limits */ void module_readParameters_debug_limit(int DEBUG, struct potentialoptimization host_potentialoptimization) { if (DEBUG == 1) // If we are in debug mode { printf("DEBUG: Center.x B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.position.x.block, host_potentialoptimization.position.x.min, host_potentialoptimization.position.x.max, host_potentialoptimization.position.x.sigma); printf("DEBUG: Center.y B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.position.y.block, host_potentialoptimization.position.y.min, host_potentialoptimization.position.y.max, host_potentialoptimization.position.y.sigma); printf("DEBUG: weight.y B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.weight.block, host_potentialoptimization.weight.min, host_potentialoptimization.weight.max, host_potentialoptimization.weight.sigma); printf("DEBUG: ellipticity_angle B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.ellipticity_angle.block, host_potentialoptimization.ellipticity_angle.min, host_potentialoptimization.ellipticity_angle.max, host_potentialoptimization.ellipticity_angle.sigma); printf("DEBUG: ellipticity B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.ellipticity.block, host_potentialoptimization.ellipticity.min, host_potentialoptimization.ellipticity.max, host_potentialoptimization.ellipticity.sigma); printf("DEBUG: ellipticity_potential B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.ellipticity_potential.block, host_potentialoptimization.ellipticity_potential.min, host_potentialoptimization.ellipticity_potential.max, host_potentialoptimization.ellipticity_potential.sigma); printf("DEBUG: rcore B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.rcore.block, host_potentialoptimization.rcore.min, host_potentialoptimization.rcore.max, host_potentialoptimization.rcore.sigma); printf("DEBUG: rcut B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.rcut.block, host_potentialoptimization.rcut.min, host_potentialoptimization.rcut.max, host_potentialoptimization.rcut.sigma); printf("DEBUG: rscale B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.rscale.block, host_potentialoptimization.rscale.min, host_potentialoptimization.rscale.max, host_potentialoptimization.rscale.sigma); printf("DEBUG: exponent B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.exponent.block, host_potentialoptimization.exponent.min, host_potentialoptimization.exponent.max, host_potentialoptimization.exponent.sigma); printf("DEBUG: vdisp B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.vdisp.block, host_potentialoptimization.vdisp.min, host_potentialoptimization.vdisp.max, host_potentialoptimization.vdisp.sigma); printf("DEBUG: alpha B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.alpha.block, host_potentialoptimization.alpha.min, host_potentialoptimization.alpha.max, host_potentialoptimization.alpha.sigma); printf("DEBUG: z B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.z.block, host_potentialoptimization.z.min, host_potentialoptimization.z.max, host_potentialoptimization.z.sigma); } } diff --git a/src/structure_hpc.hpp b/src/structure_hpc.hpp index de54923..f81a74b 100644 --- a/src/structure_hpc.hpp +++ b/src/structure_hpc.hpp @@ -1,562 +1,567 @@ /** * @file structure.h * @Author Thomas Jalabert, EPFL (me@example.com) , Christoph Schaefer, EPFL (christophernstrerne.schaefer@epfl.ch) * @date July 2015 * @version 0,1 * @brief Header file to define the used structures (e.g. defined structs) * * @param configuration file (parameters.par) * @return Depends on choice in the configuration file, e.g. least chi2 model */ // Header guard #ifndef STRUCTURE_HPC_H #define STRUCTURE_HPC_H #include #include "type.h" /*****************************************************************/ /* */ /* Constants: Will be migrated to constants.h when there are too many of them*/ /* */ /*****************************************************************/ // GPU definitions #define threadsPerBlock 512 // threads for each set of images #define MAXIMPERSOURCE 20 // maximum number of multiple images for one source #define MAXIM 200 // maximum total number of images treated // Dimensions definitions #define NPZMAX 9 /* maximum number of critical lines in g_cline struct*/ //#define CLINESIZE 500 /* maximum number of critical and caustic points for Cline mode. Adjust depending on RAM*/ #define NPOTFILESIZE 2000 //maximum number of potential in potfiles #define N_MAX_POTFILE 10 //maximum number of potfiles #define DMIN 1e-4 // distance minimale de convergence dans le plan image (in arcsec) #define NITMAX 100 #define NTYPES 2 #define DBL_MAX 1.7976931348623157e+308 // gNFW definitions #define gNFW_ARRAY_SIZE 1800 // Set the dimension of the gnfw table gNFW_ARRAY_SIZE, must be 1800 for the current table file // Filename definition #define FILENAME_SIZE 50 // size of a filename in .par file //constants definition #define pi_c2 7.209970e-06 /* pi en arcsecond/ c^2 = 648000/vol/vol */ #define cH2piG 0.11585881 /* cH0/2piG en g/cm^2 (H0=50) */ #define cH4piG 0.057929405 /* cH0/4piG en g/cm^2 (H0=50) */ #define cH0_4piG 2.7730112e-4 /* cH0/4piG en 10^12 M_Sol/kpc^2 (H0=50) */ #define d0 29.068701 /* vol/h0*1000/rta -- c/H0 en h-1.kpc/arcsecond (h0=50)*/ #define MCRIT12 .2343165 /* c^3/4piGh0/RTA/RTA in 1e12 M_sol/arcsec^2 (h0=50) */ /*****************************************************************/ /* */ /* Types definition */ /* */ /*****************************************************************/ /** @brief Point: Structure of 2 coordinates * * @param x: X coordinate * @param y: Y coordinate */ #ifdef __WITH_LENSTOOL #include "structure.h" #else struct point { type_t x; type_t y; }; /** @brief Complex: Structure of 2 doubles * @param re: Real Part * @param im: Imaginary Part */ struct complex { type_t re; type_t im; }; /** @brief Segment: Structure of two points */ struct segment { point a; point b; }; /** @brief triplet: Structure of three points defining a triangle */ struct triplet { struct point a; struct point b; struct point c; }; /** @brief bitriplet: Defines two linked triangles (one in source plane and one in image plane) * @param i: Triangle on image plane * @param s: Triangle on source plane */ struct bitriplet { struct triplet i; struct triplet s; }; /** @brief contains the table needed to compute the potential derivatives of general NFW profile */ typedef struct { double alpha_now, x_now, kappa, dpl; } gNFWdata; /** @brief Matrix: 2by2 doubles */ struct matrix { type_t a; type_t b; type_t c; type_t d; }; /** @brief ellipse: for shape computation * @param a: semimajor axis * @param b: semiminor axis * @param theta: shape ellipticity */ struct ellipse { type_t a; type_t b; type_t theta; }; #endif /** @brief Storage type for sources, lenses and arclets * @param center: position of the center of galaxy * @param shape: shape of galaxy * @param mag: magnitude * @param redshift: redshift * @param dls: Distance lens-source * @param dos: Distance observer-source * @param dr: dls/dos */ struct galaxy { //char name[IDSIZE]; struct point center; struct ellipse shape; type_t mag; type_t redshift; type_t dls; type_t dos; type_t dr; }; /** @brief Contains the information for optimising a parameter in the inverse mode * @param block: blockorfree variable (whether a parameter is blocked or free for the mcmc algorithm) * @param min: lower optimisation value * @param max: upper optimisation value * @param sigma: optimisation step (MIGHT NOT BE TAKEN INTO ACCOUNT) */ struct optimize_block { int block; type_t min; type_t max; type_t sigma; }; /** @brief two optimize_block to simulate a point */ struct optimize_point // blockorfree for the point structure { struct optimize_block x; struct optimize_block y; }; /** @brief Contains the information for optimising the potential in the inverse mode * @param position: position of the center of the halo * @param weight: weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW) * @param b0: Impact parameter * @param ellipticity_angle: orientation of the clump * @param ellipticity: Mass ellipticity * @param ellipticity_potential: Potential ellipticity * @param rcore: PIEMD specific value * @param rcut: PIEMD specific value * @param rscale: scale radius for NFW, Einasto * @param exponent: exponent for Einasto * @param vdisp: Dispersion velocity * @param alpha: exponent for general NFW * @param einasto_kappacritic: critical kappa for Einasto profile * @param z: redshift */ struct potentialoptimization // block or free variable for the MCMC for the potential { struct optimize_point position; struct optimize_block weight; struct optimize_block b0; struct optimize_block ellipticity_angle; struct optimize_block ellipticity; struct optimize_block ellipticity_potential; struct optimize_block rcore; struct optimize_block rcut; struct optimize_block rscale; struct optimize_block exponent; struct optimize_block vdisp; struct optimize_block alpha; struct optimize_block einasto_kappacritic; struct optimize_block z; }; /** @brief Contains the information of the lens potential * @param type: 1=PIEMD , 2=NFW, 3=SIES, 4=point mass, 5=SIS, 8=PIEMD * @param type_name: IEMD, NFW, SIES, point * @param name: name of the clump (e.g. name of the galaxy) : not compulsory * @param position: position of the center of the halo * @param weight: weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW) * @param b0: Impact parameter * @param ellipticity_angle: * @param ellipticity: Mass ellipticity * @param ellipticity_potential: Potential ellipticity * @param rcore: PIEMD specific value * @param rcut: PIEMD specific value * @param rscale: scale radius for NFW, Einasto * @param exponent: exponent for Einasto * @param vdisp: Dispersion velocity * @param alpha: exponent for general NFW * @param einasto_kappacritic: critical kappa for Einasto profile * @param z: redshift */ struct Potential_SOA { int* type; // 1=PIEMD ; 2=NFW; 3=SIES, 4=point mass char type_name[10]; // PIEMD, NFW, SIES, point type_t* name; // name of the clump (e.g. name of the galaxy) : not compulsory //struct point position; // position of the center of the halo type_t* position_x; // position of the center of the halo type_t* position_y; // position of the center of the halo type_t* weight; // weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW) type_t* b0; // Impact parameter type_t* vdisp; //Dispersion velocity type_t* ellipticity_angle; // orientation of the clump type_t* ellipticity; // ellipticity of the mass distribition type_t* ellipticity_potential; //ellipticity of the potential type_t* rcore; // core radius type_t* rcut; // cut radius type_t* rscale; // scale radius for NFW, Einasto type_t* exponent; // exponent for Einasto type_t* alpha; // exponent for general NFW type_t* einasto_kappacritic; // critical kappa for Einasto profile type_t* z; // redshift of the clump type_t* mag; //magnitude type_t* lum; //luminosity type_t* theta; //theta type_t* anglecos; //theta precomputation of cosinus and sinus values type_t* anglesin; //theta type_t* sigma; // sigma type_t* dlsds; // dls/dos (for Kmapping purposes int* SOA_index; //for bayes map purposes, it contains at the ith position the position in the SOA list of potential[i] }; struct Potential { int type; // 1=PIEMD ; 2=NFW; 3=SIES, 4=point mass char type_name[10]; // PIEMD, NFW, SIES, point char name[20]; // name of the clump (e.g. name of the galaxy) : not compulsory struct point position; // position of the center of the halo type_t weight; // weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW) type_t b0; // Impact parameter type_t vdisp; //Dispersion velocity type_t ellipticity_angle; // orientation of the clump type_t ellipticity; // ellipticity of the mass distribition type_t ellipticity_potential; //ellipticity of the potential type_t rcore; // core radius type_t rcut; // cut radius type_t rscale; // scale radius for NFW, Einasto type_t exponent; // exponent for Einasto type_t alpha; // exponent for general NFW type_t einasto_kappacritic; // critical kappa for Einasto profile type_t z; // redshift of the clump type_t mag; //magnitude type_t lum; //luminosity type_t theta; //theta type_t sigma; // sigma type_t dlsds; // dls/dos (for Kmapping purposes }; /*****************************************************************/ /* */ /* Control structure definition */ /* */ /*****************************************************************/ /** @brief Control structure for runmode parameters * * Default values are set in module_readParameters_readRunmode * * @param nbgridcells: Number of grid cells * @param source: flag for simple source to image conversion * @param sourfile: file name for source information * @param image: flag for simple image to source conversion * @param imafile: file name for image information * @param mass: flag for mass fitsfile * @param massgridcells: Nb of cell for fitsfile * @param z_mass: redshift for which to be computed * @param z_mass_s: redshift of source for which effect of mass will be computed * @param potential: flag for potential fitsfile * @param potgridcells: Nb of cell for fitsfile * @param z_pot: redshift for which to be computed * @param dpl: flag for displacement fitsfile * @param dplgridcells: Nb of cell for fitsfile * @param z_dpl: redshift for which to be computed * @param inverse: flag for inversion mode (MCMC etc,) * @param arclet: flag for arclet mode * @param debug: flag for debug mode * @param nimagestotal: total number of lensed images in file * @param nsets: number of sources attributed to the lensed images * @param nhalo: Number of halos * @param grid: 0 for automatic grid (not implemented), 1 for grid infor applying on source plane, 2 for grid information applying on image plane * @param nbgridcells: Number of grid cells * @param zgrid: redshift of grid * @param cline: flag for critical and caustic line calculation */ struct runmode_param { int nbgridcells; //Source Mode int source; int nsets; //Image Mode int image; int N_z_param; int nimagestot; //Mult Mode int multi; //reference type_t ref_ra ; type_t ref_dec; //Mass Mode int mass; int mass_gridcells; type_t z_mass; type_t z_mass_s; std::string mass_name; //Potential Mode int potential; int pot_gridcells; type_t z_pot; int nhalos; int n_tot_halos; //Potfile Mode int potfile; int npotfile[N_MAX_POTFILE]; //number of halos inside 1 file int Nb_potfile; //number of potfiles //displacement Mode int dpl; int dpl_gridcells; type_t z_dpl; //Inverse Mode int inverse; //Arclet Mode int arclet; //Debug Mode int debug; //Grid Mode int grid; int gridcells; type_t zgrid; //Critic and caustic mode int cline; //Amplification Mode int amplif; int amplif_gridcells; type_t z_amplif; std::string amplif_name; + //Shear Mode + int shear; + int shear_gridcells; + type_t z_shear; + std::string shear_name; //Time/Benchmark mode int time; //SOA variables int Nlens[NTYPES]; std::string imagefile; std::string potfilename[N_MAX_POTFILE]; std::string sourfile; }; /** @brief Not used yet * */ struct image_param { }; /** @brief Not used yet * */ struct source_param { }; /** @brief Contains Grid information */ struct grid_param { type_t xmin; type_t xmax; type_t ymin; type_t ymax; type_t lmin; type_t lmax; type_t rmax; }; /** @brief Control structure for cosmological parameters * * @param model: Cosmological model * @param omegaM: * @param omegaX: * @param curvature: curvature parameter * @param wX: * @param wa: * @param H0: Hubble constant * @param h: H0/50 */ struct cosmo_param { int model; type_t omegaM; type_t omegaX; type_t curvature; type_t wX; type_t wa; type_t H0; type_t h; }; /** @brief Control structure for potfile parameters * * @param potid: 1: pot P, 2: pot Q @param ftype: @param potfile[FILENAME_SIZE]; @param type; @param zlens; @param core;CCC @param corekpc; @param mag0; @param select; @param ircut; @param cut, cut1, cut2; @param cutkpc1, cutkpc2; @param isigma; @param sigma, sigma1, sigma2; @param islope; @param slope, slope1, slope2; @param ivdslope; @param vdslope, vdslope1, vdslope2; @param ivdscat; @param vdscat, vdscat1, vdscat2; @param ircutscat; @param rcutscat, rcutscat1, rcutscat2; @param ia; // scaling relation of msm200 @param a, a1, a2; @param ib; // scaling relation of msm200 @param b, b1, b2; */ struct potfile_param { int potid; // 1: pot P, 2: pot Q int ftype; char potfile[FILENAME_SIZE]; int reference_mode; type_t reference_ra; type_t reference_dec; int type; type_t zlens; type_t core; type_t corekpc; type_t mag0; int select; int ircut; type_t cut, cut1, cut2; type_t cutkpc1, cutkpc2; int isigma; type_t sigma, sigma1, sigma2; int islope; type_t slope, slope1, slope2; int ivdslope; type_t vdslope, vdslope1, vdslope2; int ivdscat; type_t vdscat, vdscat1, vdscat2; int ircutscat; type_t rcutscat, rcutscat1, rcutscat2; int ia; // scaling relation of msm200 type_t a, a1, a2; int ib; // scaling relation of msm200 type_t b, b1, b2; int npotfile; }; /** @brief Control structure for caustic and critical lines * * @param nplan: number of sourceplanes for which the caustic lines will be computed * @param cz: redshift values array for respective sourceplanes * @param dos: distcosmo1 to redshift z * @param dls: distcosmo2 between lens[0] and z * @param dlsds: ratio of dl0s/dos * @param limitLow: minimum size of the squares in MARCHINGSQUARES * @param dmax: Size of search area * @param xmin: * @param xmax: * @param ymin: * @param ymax: * @param limithigh: maximum size of the squares in MARCHINGSQUARES algorithm * @param nbgridcells: nbgridcells * nbgridcells = number of pixels for critical line computation */ struct cline_param { int nplan; type_t cz[NPZMAX]; type_t dos[NPZMAX]; // distcosmo1 to redshift z type_t dls[NPZMAX]; // distcosmo2 between lens[0] and z type_t dlsds[NPZMAX]; // ratio of dl0s/dos type_t limitLow; // minimum size of the squares in MARCHINGSQUARES or initial step size in SNAKE type_t dmax; type_t xmin; type_t xmax; type_t ymin; type_t ymax; type_t limitHigh; // maximum size of the squares in MARCHINGSQUARES algorithm int nbgridcells; // nbgridcells * nbgridcells = number of pixels for critical line computation }; #endif // if STRUCTURE_H diff --git a/utils/exec/m1931.par b/utils/exec/m1931.par index 1db03a5..17c4c27 100644 --- a/utils/exec/m1931.par +++ b/utils/exec/m1931.par @@ -1,128 +1,129 @@ runmode inverse 0 0.2 350 reference 3 292.95682 -26.575718 verbose 0 - mass 4 200 0.35 mass_4 - ampli 5 200 1000 ampli_5 + mass 0 200 0.35 mass_4 + ampli 0 200 1000 ampli_5 + shear 4 200 5 shear_1 end image multfile 1 immul_ongoing.cat mult_wcs 1 forme -1 z_m_limit 1 2.0 1 1.00 6.0 0.10 z_m_limit 2 4.0 1 1.00 6.0 0.10 z_m_limit 3 6.0 1 1.00 8.0 0.10 z_m_limit 4 5.0 1 1.00 8.0 0.10 sigposArcsec 0.5 end grille nombre 128 polaire 0 nlentille 201 nlens_opt 3 nlens_crit 3 end source z_source 1.25 end potentiel 1 # halo A profil 81 x_centre 0. y_centre 0. ellipticite 0.3 angle_pos 55.0 v_disp 1000. #rcore 50 core_radius_kpc 50. cut_radius_kpc 1000. z_lens 0.35 end limit 1 x_centre 1 -5. 5. 0.1 y_centre 1 -5. 5. 0.1 ellipticite 1 0.0 0.7 0.1 angle_pos 1 0.0 180.0 0.1 core_radius 1 1. 35. 0.1 v_disp 1 800. 1300. 0.1 end potentiel 2 # halo B profil 81 x_centre -4.03 y_centre -40.76 ellipticite 0.3 angle_pos 55.0 v_disp 1000. core_radius_kpc 50. cut_radius_kpc 1000. z_lens 0.35 end limit 2 x_centre 1 -10. 0. 0.1 y_centre 1 -45. -35. 0.1 ellipticite 1 0.0 0.9 0.1 angle_pos 1 0.0 180.0 0.1 core_radius 1 1. 35. 0.1 v_disp 1 100. 800. 0.1 end potentiel 3 # BCG profil 81 x_centre 0. y_centre 0. ellipticite 0.273 angle_pos 58.82 v_disp 100. core_radius_kpc 0.115 cut_radius_kpc 100. z_lens 0.35 end limit 3 cut_radius_kpc 1 0.000 200.000 0.100 v_disp 1 0.000 300.000 0.100 end potfile galaxies filein 1 m1931_CM_zCLASH-MUSE.cat zlens 0.35 profil 81 type 81 corekpc 0.15 mag0 19.65 sigma 3 158. 27. 250. cutkpc 1 10.0 100. 0.1 slope 0 4. vdslope 0 4. end potfile filein 1 potfile_cl_N-NW-S3.txt zlens 0.35 type 81 corekpc 0 150.0 mag0 27.0 sigma 1 100.0 300.0 cutkpc 0 10.0 100. 0.1 slope 0 4.000000 4.200000 vdslope 0 4.000000 0.000000 vdscatter 0 0.000000 0.000000 rcutscatter 0 0.000000 0.000000 end cline nplan 0 1.49 1.894 2.497 dmax 150.0 pas .1 algorithm marchingsquare end grande iso 0 0 0.0 0.0 0.0 name best profil 0 0 contour 1 0 large_dist 0.3 end cosmologie H0 70. omega 0.3 lambda 0.7 end champ dmax 136.5 end fini diff --git a/utils/exec/main.cpp b/utils/exec/main.cpp index 9baea43..34225a4 100644 --- a/utils/exec/main.cpp +++ b/utils/exec/main.cpp @@ -1,421 +1,469 @@ /** * @file main.cpp * @Author Christoph Schaaefer, EPFL (christophernstrerne.schaefer@epfl.ch) * @date October 2016 * @brief Benchmark for gradhalo function */ #include #include #include #include #include #include #include #include // //#include #include // //#include #include //#include #include "timer.h" #include "gradient.hpp" #include "gradient2.hpp" #include "chi_CPU.hpp" #include "module_cosmodistances.hpp" #include "module_readParameters.hpp" #include "grid_gradient2_CPU.hpp" #include "grid_amplif_CPU.hpp" #include "module_writeFits.hpp" #ifdef __WITH_GPU #include "grid_gradient_GPU.cuh" #include "grid_map_GPU.cuh" #include "grid_gradient2_GPU.cuh" //#include "gradient_GPU.cuh" #endif #ifdef __WITH_LENSTOOL //#include "setup.hpp" #warning "linking with libtool..." #include #include #include #include #include #include // // struct g_mode M; struct g_pot P[NPOTFILE]; struct g_pixel imFrame, wFrame, ps, PSF; struct g_cube cubeFrame; struct g_dyn Dy; // //TV // struct g_source S; struct g_image I; struct g_grille G; struct g_msgrid H; // multi-scale grid struct g_frame F; struct g_large L; struct g_cosmo C; struct g_cline CL; struct g_observ O; struct pot lens[NLMAX]; struct pot lmin[NLMAX], lmax[NLMAX], prec[NLMAX]; struct g_cosmo clmin, clmax; /*cosmological limits*/ struct galaxie smin[NFMAX], smax[NFMAX]; // limits on source parameters struct ipot ip; struct MCarlo mc; struct vfield vf; struct vfield vfmin,vfmax; // limits on velocity field parameters struct cline cl[NIMAX]; lensdata *lens_table; // int block[NLMAX][NPAMAX]; /*switch for the lens optimisation*/ int cblock[NPAMAX]; /*switch for the cosmological optimisation*/ int sblock[NFMAX][NPAMAX]; /*switch for the source parameters*/ int vfblock[NPAMAX]; /*switch for the velocity field parameters*/ double excu[NLMAX][NPAMAX]; double excd[NLMAX][NPAMAX]; /* supplments tableaux de valeurs pour fonctions g pour Einasto * * Ce sont trois variables globales qu'on pourra utiliser dans toutes les fonctions du projet * */ #define CMAX 20 #define LMAX 80 float Tab1[LMAX][CMAX]; float Tab2[LMAX][CMAX]; float Tab3[LMAX][CMAX]; int nrline, ntline, flagr, flagt; long int narclet; struct point gimage[NGGMAX][NGGMAX], gsource_global[NGGMAX][NGGMAX]; struct biline radial[NMAX], tangent[NMAX]; struct galaxie arclet[NAMAX], source[NFMAX], image[NFMAX][NIMAX]; struct galaxie cimage[NFMAX]; struct pointgal gianti[NPMAX][NIMAX]; struct point SC; double elix; double alpha_e; double *v_xx; double *v_yy; double **map_p; double **tmp_p; double **map_axx; double **map_ayy; #endif void gradient_grid_GPU_sorted(type_t *grid_grad_x, type_t *grid_grad_y, const struct grid_param *frame, const struct Potential_SOA *lens, int Nlens, int nbgridcells); // // int module_readCheckInput_readInput(int argc, char *argv[], std::string *outdir) { /// check if there is a correct number of arguments, and store the name of the input file in infile char* infile; struct stat file_stat; // If we do not have 3 arguments, stop if ( argc != 3 ) { fprintf(stderr, "\nUnexpected number of arguments\n"); fprintf(stderr, "\nUSAGE:\n"); fprintf(stderr, "lenstool input_file output_directorypath [-n]\n\n"); exit(-1); } else if ( argc == 3 ) infile=argv[1]; std::ifstream ifile(infile,std::ifstream::in); // Open the file int ts = (int) time (NULL); char buffer[10]; std::stringstream ss; ss << ts; std::string trimstamp = ss.str(); // //std::string outdir = argv[2]; *outdir = argv[2]; //*outdir += "-"; //*outdir += trimstamp; std::cout << *outdir << std::endl; // check whether the output directory already exists if (stat(outdir->c_str(), &file_stat) < 0){ mkdir(outdir->c_str(), S_IRUSR | S_IWUSR | S_IXUSR | S_IRGRP | S_IWGRP | S_IXGRP | S_IROTH ); } else { printf("Error : Directory %s already exists. Specify a non existing directory.\n",argv[2]); exit(-1); } // check whether the input file exists. If it could not be opened (ifile = 0), it does not exist if(ifile){ ifile.close(); } else{ printf("The file %s does not exist, please specify a valid file name\n",infile); exit(-1); } return 0; } // // // int main(int argc, char *argv[]) { // // Setting Up the problem // // This module function reads the terminal input when calling LENSTOOL and checks that it is correct // Otherwise it exits LENSTOOL // char cwd[1024]; if (getcwd(cwd, sizeof(cwd)) != NULL) fprintf(stdout, "Current working dir: %s\n", cwd); // std::string path; module_readCheckInput_readInput(argc, argv, &path); // // This module function reads the cosmology parameters from the parameter file // Input: struct cosmologicalparameters cosmology, parameter file // Output: Initialized cosmology struct cosmo_param cosmology; // Cosmology struct to store the cosmology data from the file std::string inputFile = argv[1]; // Input file module_readParameters_readCosmology(inputFile, cosmology); // // This module function reads the runmode paragraph and the number of sources, arclets, etc. in the parameter file. // The runmode_param stores the information of what exactly the user wants to do with lenstool. struct runmode_param runmode; module_readParameters_readRunmode(inputFile, &runmode); module_readParameters_debug_cosmology(runmode.debug, cosmology); module_readParameters_debug_runmode(1, runmode); // //=== Declaring variables // struct grid_param frame; struct galaxy images[runmode.nimagestot]; struct galaxy sources[runmode.nsets]; //struct Potential lenses[runmode.n_tot_halos]; struct Potential_SOA lenses_SOA_table[NTYPES]; struct Potential_SOA lenses_SOA; struct cline_param cline; struct potfile_param potfile[runmode.Nb_potfile]; //struct Potential potfilepotentials[runmode.npotfile]; struct potentialoptimization host_potentialoptimization[runmode.nhalos]; int nImagesSet[runmode.nsets]; // Contains the number of images in each set of images //Bayesmap specific variables type_t* bayespot; int nparam, nvalues; // This module function reads in the potential form and its parameters (e.g. NFW) // Input: input file // Output: Potentials and its parameters module_readParameters_PotentialSOA_direct(inputFile, &lenses_SOA, runmode.nhalos, runmode.n_tot_halos, cosmology); module_readParameters_debug_potential_SOA(0, lenses_SOA, runmode.nhalos); module_readParameters_limit(inputFile, host_potentialoptimization, runmode.nhalos ); #if 0 for(int ii = 0; ii .fits char fname2[50]; // .fits FILE *pFile; int i; double *index; // list of bayes.dat lines int seed; // random seed int tmp; printf("Setting up lenstool using %d lenses...", runmode.n_tot_halos); fflush(stdout); //convert_to_LT(&lenses_SOA, runmode.nhalos+runmode.npotfile); // Read the .par file init_grille(argv[1], 1); // remove the .fits extension tcpo filename if( M.imass ) M.massfile[strlen(M.massfile)-5]=0; if( M.ishear ) M.shearfile[strlen(M.shearfile)-5]=0; if( M.iampli ) M.amplifile[strlen(M.amplifile)-5]=0; if( M.idpl ) { M.dplxfile[strlen(M.dplxfile)-5]=0; M.dplyfile[strlen(M.dplyfile)-5]=0; } if( M.pixel ) M.pixelfile[strlen(M.pixelfile)-5]=0; if( M.iclean ) ps.pixfile[strlen(ps.pixfile)-5]=0; // Read catalog of multiple images readConstraints(); // Initialise the grid if( G.pol != 0 ) gridp(); else grid(); // Switch to silent mode M.verbose = 0; printf("ok\n"); #if 1 for(int i = 0; i < G.nlens; i++){ printf("Lenstool Potential[%d]: x = %f, y = %f, vdisp = %f, type = %d \n \t ellipticity = %f, ellipticity_pot = %f, ellipticity angle (radians) = %f, rcore = %f, rcut = %f,\n z = %f\n", i,lens[i].C.x, lens[i].C.y, lens[i].sigma, lens[i].type, lens[i].emass, lens[i].epot, lens[i].theta, lens[i].rc, lens[i].rcut, lens[i].z); } #endif // Create the ./tmp directory i = system("mkdir -p tmp"); // Prepare the index list index = (double *) malloc((unsigned) nVal*sizeof(double)); for( i = 0 ; i < nVal ; i++ ) index[i]=i; seed = -2; std::cerr << " Finished setting up" << std::endl; /* antecedant d'une grille source ou d'une grille image */ if ( M.grille != 0 ){ t_lt = -myseconds(); g_grid(M.grille, M.ngrille, M.zgrille); t_lt += myseconds();} /* grille du amplification (keyword ampli)*/ if ( M.iampli != 0 ){ t_lt = -myseconds(); g_ampli(M.iampli, M.nampli, M.zampli, M.amplifile); t_lt += myseconds();} /* grille du potential */ if ( M.ipoten != 0 ){ t_lt = -myseconds(); g_poten(M.ipoten, M.npoten, M.zpoten, M.potenfile); t_lt += myseconds();} /* grille du mass */ if ( M.imass != 0 ){ t_lt = -myseconds(); g_mass(M.imass, M.nmass, M.zmass, S.zs, M.massfile); t_lt += myseconds();} /* grille du dpl */ if ( M.idpl != 0 ){ t_lt = -myseconds(); g_dpl(M.idpl, M.ndpl, M.zdpl, M.dplxfile, M.dplyfile); t_lt += myseconds();} + /* grille du curv */ + if ( M.icurv != 0 ) + g_curv(M.icurv, M.ncurv, M.zcurv, M.cxxfile, M.cxyfile, M.cyyfile); + + /* grille du shear */ + if ( M.ishear != 0 ) + g_shear(M.ishear, M.nshear, M.zshear, M.shearfile); + + /* grille du time-delay */ + if ( M.itime != 0 ) + g_time(M.itime, M.ntime, M.ztime, M.timefile); + + /* grille du shear_field */ + if ( M.ishearf != 0 ) + g_shearf(M.ishearf, M.zshearf, M.shearffile, M.nshearf); + + /* amplification field grid */ + if ( M.iamplif != 0 ) + g_amplif(M.iamplif, M.zamplif, M.ampliffile); + #endif #ifdef __WITH_GPU double t_1, t_2; t_1 = -myseconds(); std::cerr << runmode.mass << " " << runmode.amplif << " " << std::endl; if (runmode.mass > 0){ //Allocation type_t* mass_GPU = (type_t *) malloc((int) (runmode.mass_gridcells) * (runmode.mass_gridcells) * sizeof(type_t)); ////calculate maps std::cout << " GPU launching for map mass " << std::endl; t_2 = -myseconds(); module_readParameters_debug_potential_SOA(0, lenses_SOA, runmode.n_tot_halos); //Init memset(mass_GPU, 0, (runmode.mass_gridcells) * (runmode.mass_gridcells) * sizeof(type_t)); //Choosing Function definition map_gpu_function_t map_gpu_func; map_gpu_func = select_map_function("mass",&runmode); //calculating map using defined function map_mass_grid_GPU(map_gpu_func,mass_GPU,&cosmology, &frame, &lenses_SOA, runmode.n_tot_halos, runmode.mass_gridcells ,runmode.mass, runmode.z_mass, runmode.z_mass_s); //writing module_writeFits(path,runmode.mass_name,mass_GPU,&runmode,runmode.mass_gridcells,&frame, runmode.ref_ra, runmode.ref_dec ); t_2 += myseconds(); std::cout << " Time " << std::setprecision(15) << t_2 << std::endl; std::cerr << "**" << mass_GPU[0] << std::endl; free(mass_GPU); } if (runmode.amplif > 0){ //Allocation type_t* ampli_GPU = (type_t *) malloc((int) (runmode.amplif_gridcells) * (runmode.amplif_gridcells) * sizeof(type_t)); ////calculate maps std::cout << " GPU launching for map amplif " << std::endl; t_2 = -myseconds(); module_readParameters_debug_potential_SOA(0, lenses_SOA, runmode.n_tot_halos); //Init memset(ampli_GPU, 0, (runmode.amplif_gridcells) * (runmode.amplif_gridcells) * sizeof(type_t)); //Choosing Function definition map_gpu_function_t map_gpu_func; map_gpu_func = select_map_function("ampli",&runmode); //calculating map using defined function map_grid_GPU(map_gpu_func,ampli_GPU,&cosmology, &frame, &lenses_SOA, runmode.n_tot_halos, runmode.amplif_gridcells ,runmode.amplif, runmode.z_amplif); //writing //std::cerr << runmode.amplif_name << std::endl; module_writeFits(path,runmode.amplif_name,ampli_GPU,&runmode,runmode.amplif_gridcells,&frame, runmode.ref_ra, runmode.ref_dec ); t_2 += myseconds(); std::cout << " Time " << std::setprecision(15) << t_2 << std::endl; std::cerr << "**" << ampli_GPU[0] << std::endl; //std::cerr << "**" << ampli_GPU[0] << std::endl; free(ampli_GPU); } + if (runmode.shear > 0){ + //Allocation + type_t* shear_GPU = (type_t *) malloc((int) (runmode.shear_gridcells) * (runmode.shear_gridcells) * sizeof(type_t)); + + ////calculate maps + std::cout << " GPU launching for map shear " << std::endl; + t_2 = -myseconds(); + module_readParameters_debug_potential_SOA(0, lenses_SOA, runmode.n_tot_halos); + //Init + memset(shear_GPU, 0, (runmode.shear_gridcells) * (runmode.shear_gridcells) * sizeof(type_t)); + + //Choosing Function definition + map_gpu_function_t map_gpu_func; + map_gpu_func = select_map_function("shear",&runmode); + + //calculating map using defined function + map_grid_GPU(map_gpu_func,shear_GPU,&cosmology, &frame, &lenses_SOA, runmode.n_tot_halos, runmode.shear_gridcells ,runmode.shear, runmode.z_shear); + + //writing + //std::cerr << runmode.amplif_name << std::endl; + module_writeFits(path,runmode.shear_name,shear_GPU,&runmode,runmode.shear_gridcells,&frame, runmode.ref_ra, runmode.ref_dec ); + t_2 += myseconds(); + std::cout << " Time " << std::setprecision(15) << t_2 << std::endl; + std::cerr << "**" << shear_GPU[0] << std::endl; + + //std::cerr << "**" << ampli_GPU[0] << std::endl; + free(shear_GPU); + } // t_1 += myseconds(); std::cout << "Lenstool 1 Map Time " << std::setprecision(15) << t_lt << std::endl; //std::cout << "HPC Total Time " << std::setprecision(15) << t_1 << std::endl; #endif } diff --git a/utils/maps/m1931.par b/utils/maps/m1931.par index 1db03a5..26284bc 100644 --- a/utils/maps/m1931.par +++ b/utils/maps/m1931.par @@ -1,128 +1,129 @@ runmode inverse 0 0.2 350 reference 3 292.95682 -26.575718 verbose 0 mass 4 200 0.35 mass_4 ampli 5 200 1000 ampli_5 + shear 1 200 0.3 shear_1 end image multfile 1 immul_ongoing.cat mult_wcs 1 forme -1 z_m_limit 1 2.0 1 1.00 6.0 0.10 z_m_limit 2 4.0 1 1.00 6.0 0.10 z_m_limit 3 6.0 1 1.00 8.0 0.10 z_m_limit 4 5.0 1 1.00 8.0 0.10 sigposArcsec 0.5 end grille nombre 128 polaire 0 nlentille 201 nlens_opt 3 nlens_crit 3 end source z_source 1.25 end potentiel 1 # halo A profil 81 x_centre 0. y_centre 0. ellipticite 0.3 angle_pos 55.0 v_disp 1000. #rcore 50 core_radius_kpc 50. cut_radius_kpc 1000. z_lens 0.35 end limit 1 x_centre 1 -5. 5. 0.1 y_centre 1 -5. 5. 0.1 ellipticite 1 0.0 0.7 0.1 angle_pos 1 0.0 180.0 0.1 core_radius 1 1. 35. 0.1 v_disp 1 800. 1300. 0.1 end potentiel 2 # halo B profil 81 x_centre -4.03 y_centre -40.76 ellipticite 0.3 angle_pos 55.0 v_disp 1000. core_radius_kpc 50. cut_radius_kpc 1000. z_lens 0.35 end limit 2 x_centre 1 -10. 0. 0.1 y_centre 1 -45. -35. 0.1 ellipticite 1 0.0 0.9 0.1 angle_pos 1 0.0 180.0 0.1 core_radius 1 1. 35. 0.1 v_disp 1 100. 800. 0.1 end potentiel 3 # BCG profil 81 x_centre 0. y_centre 0. ellipticite 0.273 angle_pos 58.82 v_disp 100. core_radius_kpc 0.115 cut_radius_kpc 100. z_lens 0.35 end limit 3 cut_radius_kpc 1 0.000 200.000 0.100 v_disp 1 0.000 300.000 0.100 end potfile galaxies filein 1 m1931_CM_zCLASH-MUSE.cat zlens 0.35 profil 81 type 81 corekpc 0.15 mag0 19.65 sigma 3 158. 27. 250. cutkpc 1 10.0 100. 0.1 slope 0 4. vdslope 0 4. end potfile filein 1 potfile_cl_N-NW-S3.txt zlens 0.35 type 81 corekpc 0 150.0 mag0 27.0 sigma 1 100.0 300.0 cutkpc 0 10.0 100. 0.1 slope 0 4.000000 4.200000 vdslope 0 4.000000 0.000000 vdscatter 0 0.000000 0.000000 rcutscatter 0 0.000000 0.000000 end cline nplan 0 1.49 1.894 2.497 dmax 150.0 pas .1 algorithm marchingsquare end grande iso 0 0 0.0 0.0 0.0 name best profil 0 0 contour 1 0 large_dist 0.3 end cosmologie H0 70. omega 0.3 lambda 0.7 end champ dmax 136.5 end fini diff --git a/utils/maps/main.cpp b/utils/maps/main.cpp index a35b552..a7b0b2e 100644 --- a/utils/maps/main.cpp +++ b/utils/maps/main.cpp @@ -1,482 +1,510 @@ /** * @file main.cpp * @Author Christoph Schaaefer, EPFL (christophernstrerne.schaefer@epfl.ch) * @date October 2016 * @brief Benchmark for gradhalo function */ #include #include #include #include #include #include #include #include // //#include #include // //#include #include //#include #include "timer.h" #include "gradient.hpp" #include "gradient2.hpp" #include "chi_CPU.hpp" #include "module_cosmodistances.hpp" #include "module_readParameters.hpp" #include "grid_gradient2_CPU.hpp" #include "grid_amplif_CPU.hpp" #include "module_writeFits.hpp" #ifdef __WITH_GPU #include "grid_gradient_GPU.cuh" #include "grid_map_GPU.cuh" #include "grid_gradient2_GPU.cuh" //#include "gradient_GPU.cuh" #endif #ifdef __WITH_LENSTOOL //#include "setup.hpp" #warning "linking with libtool..." #include #include #include #include #include #include // // struct g_mode M; struct g_pot P[NPOTFILE]; struct g_pixel imFrame, wFrame, ps, PSF; struct g_cube cubeFrame; struct g_dyn Dy; // //TV // struct g_source S; struct g_image I; struct g_grille G; struct g_msgrid H; // multi-scale grid struct g_frame F; struct g_large L; struct g_cosmo C; struct g_cline CL; struct g_observ O; struct pot lens[NLMAX]; struct pot lmin[NLMAX], lmax[NLMAX], prec[NLMAX]; struct g_cosmo clmin, clmax; /*cosmological limits*/ struct galaxie smin[NFMAX], smax[NFMAX]; // limits on source parameters struct ipot ip; struct MCarlo mc; struct vfield vf; struct vfield vfmin,vfmax; // limits on velocity field parameters struct cline cl[NIMAX]; lensdata *lens_table; // int block[NLMAX][NPAMAX]; /*switch for the lens optimisation*/ int cblock[NPAMAX]; /*switch for the cosmological optimisation*/ int sblock[NFMAX][NPAMAX]; /*switch for the source parameters*/ int vfblock[NPAMAX]; /*switch for the velocity field parameters*/ double excu[NLMAX][NPAMAX]; double excd[NLMAX][NPAMAX]; /* supplments tableaux de valeurs pour fonctions g pour Einasto * * Ce sont trois variables globales qu'on pourra utiliser dans toutes les fonctions du projet * */ #define CMAX 20 #define LMAX 80 float Tab1[LMAX][CMAX]; float Tab2[LMAX][CMAX]; float Tab3[LMAX][CMAX]; int nrline, ntline, flagr, flagt; long int narclet; struct point gimage[NGGMAX][NGGMAX], gsource_global[NGGMAX][NGGMAX]; struct biline radial[NMAX], tangent[NMAX]; struct galaxie arclet[NAMAX], source[NFMAX], image[NFMAX][NIMAX]; struct galaxie cimage[NFMAX]; struct pointgal gianti[NPMAX][NIMAX]; struct point SC; double elix; double alpha_e; double *v_xx; double *v_yy; double **map_p; double **tmp_p; double **map_axx; double **map_ayy; #endif void gradient_grid_GPU_sorted(type_t *grid_grad_x, type_t *grid_grad_y, const struct grid_param *frame, const struct Potential_SOA *lens, int Nlens, int nbgridcells); // // int module_readCheckInput_readInput(int argc, char *argv[], std::string *outdir) { /// check if there is a correct number of arguments, and store the name of the input file in infile char* infile; struct stat file_stat; // If we do not have 3 arguments, stop if ( argc != 3 ) { fprintf(stderr, "\nUnexpected number of arguments\n"); fprintf(stderr, "\nUSAGE:\n"); fprintf(stderr, "lenstool input_file output_directorypath [-n]\n\n"); exit(-1); } else if ( argc == 3 ) infile=argv[1]; std::ifstream ifile(infile,std::ifstream::in); // Open the file int ts = (int) time (NULL); char buffer[10]; std::stringstream ss; ss << ts; std::string trimstamp = ss.str(); // //std::string outdir = argv[2]; *outdir = argv[2]; //*outdir += "-"; //*outdir += trimstamp; std::cout << *outdir << std::endl; // check whether the output directory already exists if (stat(outdir->c_str(), &file_stat) < 0){ mkdir(outdir->c_str(), S_IRUSR | S_IWUSR | S_IXUSR | S_IRGRP | S_IWGRP | S_IXGRP | S_IROTH ); } else { printf("Error : Directory %s already exists. Specify a non existing directory.\n",argv[2]); exit(-1); } // check whether the input file exists. If it could not be opened (ifile = 0), it does not exist if(ifile){ ifile.close(); } else{ printf("The file %s does not exist, please specify a valid file name\n",infile); exit(-1); } return 0; } // // // int main(int argc, char *argv[]) { // // Setting Up the problem // // This module function reads the terminal input when calling LENSTOOL and checks that it is correct // Otherwise it exits LENSTOOL // char cwd[1024]; if (getcwd(cwd, sizeof(cwd)) != NULL) fprintf(stdout, "Current working dir: %s\n", cwd); // std::string path; module_readCheckInput_readInput(argc, argv, &path); // // This module function reads the cosmology parameters from the parameter file // Input: struct cosmologicalparameters cosmology, parameter file // Output: Initialized cosmology struct cosmo_param cosmology; // Cosmology struct to store the cosmology data from the file std::string inputFile = argv[1]; // Input file module_readParameters_readCosmology(inputFile, cosmology); // // This module function reads the runmode paragraph and the number of sources, arclets, etc. in the parameter file. // The runmode_param stores the information of what exactly the user wants to do with lenstool. struct runmode_param runmode; module_readParameters_readRunmode(inputFile, &runmode); module_readParameters_debug_cosmology(runmode.debug, cosmology); module_readParameters_debug_runmode(1, runmode); // //=== Declaring variables // struct grid_param frame; struct galaxy images[runmode.nimagestot]; struct galaxy sources[runmode.nsets]; //struct Potential lenses[runmode.n_tot_halos]; struct Potential_SOA lenses_SOA_table[NTYPES]; struct Potential_SOA lenses_SOA; struct cline_param cline; struct potfile_param potfile[runmode.Nb_potfile]; //struct Potential potfilepotentials[runmode.npotfile]; struct potentialoptimization host_potentialoptimization[runmode.nhalos]; int nImagesSet[runmode.nsets]; // Contains the number of images in each set of images //Bayesmap specific variables type_t* bayespot; int nparam, nvalues; // This module function reads in the potential form and its parameters (e.g. NFW) // Input: input file // Output: Potentials and its parameters module_readParameters_PotentialSOA_direct(inputFile, &lenses_SOA, runmode.nhalos, runmode.n_tot_halos, cosmology); module_readParameters_debug_potential_SOA(1, lenses_SOA, runmode.nhalos); module_readParameters_limit(inputFile, host_potentialoptimization, runmode.nhalos ); #if 0 for(int ii = 0; ii .fits char fname2[50]; // .fits FILE *pFile; int i; double *index; // list of bayes.dat lines int seed; // random seed int tmp; printf("Setting up lenstool using %d lenses...", runmode.n_tot_halos); fflush(stdout); //convert_to_LT(&lenses_SOA, runmode.nhalos+runmode.npotfile); // Read the .par file init_grille(argv[1], 1); // remove the .fits extension tcpo filename if( M.imass ) M.massfile[strlen(M.massfile)-5]=0; if( M.ishear ) M.shearfile[strlen(M.shearfile)-5]=0; if( M.iampli ) M.amplifile[strlen(M.amplifile)-5]=0; if( M.idpl ) { M.dplxfile[strlen(M.dplxfile)-5]=0; M.dplyfile[strlen(M.dplyfile)-5]=0; } if( M.pixel ) M.pixelfile[strlen(M.pixelfile)-5]=0; if( M.iclean ) ps.pixfile[strlen(ps.pixfile)-5]=0; // Read catalog of multiple images readConstraints(); // Initialise the grid if( G.pol != 0 ) gridp(); else grid(); // Switch to silent mode M.verbose = 0; printf("ok\n"); std::cerr << " Read Bayes models" << std::endl; // Read the bayes.dat file array = readBayesModels(&nParam, &nVal); if( array == NULL ) { fprintf(stderr, "ERROR: bayes.dat file not found\n"); return -1; } #if 0 for(int i = 0; i < G.nlens; i++){ printf("Lenstool Potential[%d]: x = %f, y = %f, vdisp = %f, type = %d \n \t ellipticity = %f, ellipticity_pot = %f, ellipticity angle (radians) = %f, rcore = %f, rcut = %f,\n z = %f\n", i,lens[i].C.x, lens[i].C.y, lens[i].sigma, lens[i].type, lens[i].emass, lens[i].epot, lens[i].theta, lens[i].rc, lens[i].rcut, lens[i].z); } #endif // Create the ./tmp directory i = system("mkdir -p tmp"); // Prepare the index list index = (double *) malloc((unsigned) nVal*sizeof(double)); for( i = 0 ; i < nVal ; i++ ) index[i]=i; seed = -2; std::cerr << " Finished setting up" << std::endl; //Defining maps int ampli = 1; t_lt_total = -myseconds(); // Loop over each line for( i = 0; i < nVal && i < 2000; i++ ) { // Randomly draw a line from index array tmp = (int) floor(d_random(&seed) * (nVal - i)); iVal = index[i+tmp]; // and swap the indexes in the index list index[i+tmp] = index[i]; // Set the lens parameters from setBayesModel( iVal, nVal, array ); #if 0 for(int i = 0; i < G.nlens; i++){ printf("Lenstool Potential[%d]: x = %f, y = %f, vdisp = %f, type = %d \n \t ellipticity = %f, ellipticity_pot = %f, ellipticity angle (radians) = %f, rcore = %f, rcut = %f,\n z = %f\n", i,lens[i].C.x, lens[i].C.y, lens[i].sigma, lens[i].type, lens[i].emass, lens[i].epot, lens[i].theta, lens[i].rc, lens[i].rcut, lens[i].z); } #endif if( M.imass != 0 ) { sprintf( fname, "tmp/%s%ld.fits",M.massfile, iVal ); printf("INFO: Compute file %d/%ld : %s [CTRL-C to interrupt]\n",i+1, nVal,fname); fflush(stdout); pFile = fopen( fname, "r" ); if( pFile == NULL ) { pFile = fopen( fname, "w"); fprintf( pFile, "busy\n" ); fclose(pFile); t_lt = -myseconds(); g_mass( M.imass, M.nmass, M.zmass, S.zs, fname ); t_lt += myseconds(); std::cout << " Time = " << std::setprecision(15) << t_lt << " " << turn < 0){ //Allocation type_t* mass_GPU = (type_t *) malloc((int) (runmode.mass_gridcells) * (runmode.mass_gridcells) * sizeof(type_t)); for(int ii = 0; ii < nvalues; ii++){ ////calculate maps std::cout << " GPU launching for map mass " << ii << std::endl; t_2 = -myseconds(); ////set bayes potential module_readParameters_setbayesmapmodels(&runmode, &cosmology, host_potentialoptimization, potfile, &lenses_SOA,bayespot,nparam, ii); //std::cerr << runmode.n_tot_halos << std::endl; module_readParameters_debug_potential_SOA(0, lenses_SOA, runmode.n_tot_halos); //Init memset(mass_GPU, 0, (runmode.mass_gridcells) * (runmode.mass_gridcells) * sizeof(type_t)); //Choosing Function definition map_gpu_function_t map_gpu_func; map_gpu_func = select_map_function("mass",&runmode); //calculating map using defined function map_mass_grid_GPU(map_gpu_func,mass_GPU,&cosmology, &frame, &lenses_SOA, runmode.n_tot_halos, runmode.mass_gridcells ,runmode.mass, runmode.z_mass, runmode.z_mass_s); std::cerr <<" para : " << runmode.z_mass << " " << runmode.z_mass_s << std::endl; //writing //std::cerr << runmode.amplif_name << std::endl; module_writeFits(path,runmode.mass_name,ii,mass_GPU,&runmode,runmode.mass_gridcells,&frame, runmode.ref_ra, runmode.ref_dec ); t_2 += myseconds(); std::cout << " Time " << std::setprecision(15) << t_2 << std::endl; std::cerr << "**" << mass_GPU[0] << std::endl; } //std::cerr << "**" << ampli_GPU[0] << std::endl; free(mass_GPU); } if (runmode.amplif > 0){ //Allocation type_t* ampli_GPU = (type_t *) malloc((int) (runmode.amplif_gridcells) * (runmode.amplif_gridcells) * sizeof(type_t)); for(int ii = 0; ii < nvalues; ii++){ ////calculate maps std::cout << " GPU launching for map amplif " << ii << std::endl; t_2 = -myseconds(); ////set bayes potential module_readParameters_setbayesmapmodels(&runmode, &cosmology, host_potentialoptimization, potfile, &lenses_SOA,bayespot,nparam, ii); module_readParameters_debug_potential_SOA(0, lenses_SOA, runmode.n_tot_halos); //Init memset(ampli_GPU, 0, (runmode.amplif_gridcells) * (runmode.amplif_gridcells) * sizeof(type_t)); //Choosing Function definition map_gpu_function_t map_gpu_func; map_gpu_func = select_map_function("ampli",&runmode); //calculating map using defined function map_grid_GPU(map_gpu_func,ampli_GPU,&cosmology, &frame, &lenses_SOA, runmode.n_tot_halos, runmode.amplif_gridcells ,runmode.amplif, runmode.z_amplif); //writing //std::cerr << runmode.amplif_name << std::endl; module_writeFits(path,runmode.amplif_name,ii,ampli_GPU,&runmode,runmode.amplif_gridcells,&frame, runmode.ref_ra, runmode.ref_dec ); t_2 += myseconds(); std::cout << " Time " << std::setprecision(15) << t_2 << std::endl; std::cerr << "**" << ampli_GPU[0] << std::endl; } //std::cerr << "**" << ampli_GPU[0] << std::endl; free(ampli_GPU); } + if (runmode.shear > 0){ + //Allocation + type_t* shear_GPU = (type_t *) malloc((int) (runmode.shear_gridcells) * (runmode.shear_gridcells) * sizeof(type_t)); + for(int ii = 0; ii < nvalues; ii++){ + ////calculate maps + std::cout << " GPU launching for map amplif " << ii << std::endl; + t_2 = -myseconds(); + ////set bayes potential + module_readParameters_setbayesmapmodels(&runmode, &cosmology, host_potentialoptimization, potfile, &lenses_SOA,bayespot,nparam, ii); + module_readParameters_debug_potential_SOA(0, lenses_SOA, runmode.n_tot_halos); + //Init + memset(shear_GPU, 0, (runmode.shear_gridcells) * (runmode.shear_gridcells) * sizeof(type_t)); + + //Choosing Function definition + map_gpu_function_t map_gpu_func; + map_gpu_func = select_map_function("shear",&runmode); + + //calculating map using defined function + map_grid_GPU(map_gpu_func,shear_GPU,&cosmology, &frame, &lenses_SOA, runmode.n_tot_halos, runmode.shear_gridcells ,runmode.shear, runmode.z_shear); + + //writing + module_writeFits(path,runmode.shear_name,ii,shear_GPU,&runmode,runmode.shear_gridcells,&frame, runmode.ref_ra, runmode.ref_dec ); + t_2 += myseconds(); + std::cout << " Time " << std::setprecision(15) << t_2 << std::endl; + std::cerr << "**" << shear_GPU[0] << std::endl; + } + free(shear_GPU); + } // t_1 += myseconds(); std::cout << "Lenstool Total Time " << std::setprecision(15) << t_lt_total << std::endl; std::cout << "HPC Total Time " << std::setprecision(15) << t_1 << std::endl; #endif }