mirror/monofonIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00
Code Releases Activity

cleanup of Adrian's PANPHASIA_ho code to avoid compiler warnings in pedantic mode

Browse source
This commit is contained in:
Oliver Hahn 2021-09-01 21:30:34 +02:00
parent cb8c6bf209
commit ffcc03831f
6 changed files with 1466 additions and 1195 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

35
external/panphasia_ho/LICENSE vendored Normal file
View file

@ -0,0 +1,35 @@
The code in this subdirectory is part of Adrian Jenkins' PANPHASIA packet,
obtained from here http://icc.dur.ac.uk/Panphasia.php
PANPHASIA is not published under the GPL but has its own proprietary license,
make sure to visit the website before using the PANPHASIA functionality of
MUSIC2 and register your name.
We reproduce the licensing requirements for PANPHASIA from the above website
as retrieved on 2020/08/23:
We make our software available for free but with a licence that includes the
condition that users make sure the phases of any new simulation volumes set up
using Panphasia are published.
We are happy to collaborate with others on improving the software and providing
support for languages other than fortran. Contact: A.R.Jenkins@durham.ac.uk
LICENCE:
You are licensed to use this software free of charge on condition that:
- you will publish the phase descriptors and reference Jenkins (13) for any new
simulations that use Panphasia phases. You will pass on this condition to others
for any software or data you make available publically or privately that makes
use of Panphasia.
- that you will ensure any publications using results derived from Panphasia will
be submitted as a final version to arXiv prior to or coincident with publication
in a journal.
- that you report any bugs in this software as soon as confirmed to
A.R.Jenkins@durham.ac.uk
- that you understand that the software comes with no warranty and that is your
responsibility to ensure that it is suitable for the purpose that you intend.
- that you agree to having your name and email address stored for an indefinite
period in the future electronically in a database as a record that you agreed
the licence conditions.

230
external/panphasia_ho/high_order_panphasia_routines.c vendored
View file

@ -195,7 +195,7 @@ for (i=0; i<8*Nbasis; i++) output_vec_children[i] = work_vec2[i];
void box_muller_(PAN_REAL *unif_rand, PAN_REAL *gvar)
{
int i, j, k, count;
int i, count;
const PAN_REAL pi = 4.0 * atan(1.0);
const PAN_REAL two_pi = 2.0 * pi;
@ -264,7 +264,7 @@ void speed_test2_()
size_t N_cells = 1e6;
PAN_REAL parent[Nbasis];
PAN_REAL child[8 * Nbasis];
// PAN_REAL child[8 * Nbasis];
PAN_REAL output[8 * Nbasis];
//ticks tick_start = getticks();
@ -285,7 +285,7 @@ void speed_test2_()
//tic_total = getticks() - tick_start;
//printf("Computed %ld cells in time %.3f %s\n",N_cells,clocks_from_ticks(tic_total),clocks_getunit());
};
}
//===================================================================================
void compute_all_properties_of_a_panphasia_cell_(size_t *level, size_t *j1, size_t *j2, size_t *j3,
@ -326,7 +326,7 @@ void compute_all_properties_of_a_panphasia_cell_(size_t *level, size_t *j1, size
//
//
// };
};
}
//==================================================================================
void test_random_dist_(size_t ishift)
@ -369,8 +369,8 @@ void test_random_dist_(size_t ishift)
char str1[100], str2[100];
sprintf(str1, "Gaussian_random_distribution_%llu.dat", ishift);
sprintf(str2, "Log_uniform_random_distribution_%llu.dat", ishift);
sprintf(str1, "Gaussian_random_distribution_%lu.dat", ishift);
sprintf(str2, "Log_uniform_random_distribution_%lu.dat", ishift);
FILE *file = fopen(str1, "w");
FILE *file2 = fopen(str2, "w");
@ -378,7 +378,7 @@ void test_random_dist_(size_t ishift)
for (size_t j3 = 0; j3 < NC; j3++)
{
if (j3 % 10000000 == 0)
printf("Looped over %lld\n", j3);
printf("Looped over %ld\n", j3);
return_uniform_pseudo_rands_threefry4x64_(l, j1, j2, j3, unif_randoms, seed_value, allow_non_zero_seed_saftey_catch);
@ -402,9 +402,9 @@ void test_random_dist_(size_t ishift)
rms_value = sqrt(sum_squares / (double)nrand);
printf("Number of rands %ld RMS = %12.10lg Deviation %lg \n",
printf("Number of rands %lld RMS = %12.10lg Deviation %lg \n",
nrand, rms_value, (rms_value - 1.0) * sqrt((double)nrand));
fprintf(file, "Number of rands %ld RMS = %12.10lg Deviation %lg \n",
fprintf(file, "Number of rands %lld RMS = %12.10lg Deviation %lg \n",
nrand, rms_value, (rms_value - 1.0) * sqrt((double)nrand));
for (int i = 0; i < 100; i++)
@ -414,9 +414,9 @@ void test_random_dist_(size_t ishift)
{
g_expected = 0.5 * (erf(0.2 * sqrt(0.5) * ((PAN_REAL)(i - array_offset) + 0.5)) - erf(0.2 * sqrt(0.5) * ((PAN_REAL)(i - array_offset) - 0.5))) * (PAN_REAL)nrand;
printf("%d %ld %f %f \n", i - array_offset, gauss_dist[i], g_expected, (gauss_dist[i] - g_expected) / sqrt(gauss_dist[i]));
printf("%d %lld %f %f \n", i - array_offset, gauss_dist[i], g_expected, (gauss_dist[i] - g_expected) / sqrt(gauss_dist[i]));
fprintf(file, "%d %ld %f %f \n", i - array_offset, gauss_dist[i], g_expected, (gauss_dist[i] - g_expected) / sqrt(gauss_dist[i]));
fprintf(file, "%d %lld %f %f \n", i - array_offset, gauss_dist[i], g_expected, (gauss_dist[i] - g_expected) / sqrt(gauss_dist[i]));
};
if (log_uniform_dist[i] != 0)
{
@ -612,20 +612,20 @@ int demo_descriptor_()
// char str[200] = "[Panph6,L21,(1136930,890765,1847934),S3,CH2414478110,Auriga_volume2]";
// char str[200] = "[Panph6,L21,(1136930,890765,1847934),S3,CH-999,Auriga_volume2]";
char copy[200];
const char s[20] = "[,L,(),S,CH,]";
char *token;
// char copy[200];
// const char s[20] = "[,L,(),S,CH,]";
// char *token;
size_t desc_level, desc_x, desc_y, desc_z, desc_size;
long long int desc_ch;
char desc_name[100];
char desc_iden[8];
// size_t desc_level, desc_x, desc_y, desc_z, desc_size;
// long long int desc_ch;
// char desc_name[100];
// char desc_iden[8];
int error_code;
int pan_mode;
descriptor_read_in = 0;
if (error_code = parse_and_validate_descriptor_(str,&pan_mode))
if ((error_code = parse_and_validate_descriptor_(str, &pan_mode)))
{
printf("Invalid descriptor %s\n", str);
@ -639,12 +639,12 @@ int demo_descriptor_()
if (descriptor_read_in)
{
printf("-----------------------------------------\n");
printf("Descriptor order: %llu\n", descriptor_order);
printf("Descriptor base level: %llu\n", descriptor_base_level);
printf("Descriptor x-origin: %llu\n", descriptor_xorigin);
printf("Descriptor y-origin: %llu\n", descriptor_yorigin);
printf("Descriptor z-origin: %llu\n", descriptor_zorigin);
printf("Descriptor base size: %llu\n", descriptor_base_size);
printf("Descriptor order: %lu\n", descriptor_order);
printf("Descriptor base level: %lu\n", descriptor_base_level);
printf("Descriptor x-origin: %lu\n", descriptor_xorigin);
printf("Descriptor y-origin: %lu\n", descriptor_yorigin);
printf("Descriptor z-origin: %lu\n", descriptor_zorigin);
printf("Descriptor base size: %lu\n", descriptor_base_size);
printf("Descriptor check digit:%lld\n", descriptor_check_digit);
printf("Descriptor name %s\n", descriptor_name);
printf("-----------------------------------------\n");
@ -674,8 +674,8 @@ int demo_descriptor_()
verbose = 0;
if (error_code = PANPHASIA_init_level_(&rel_lev,
&rel_orig_x, &rel_orig_y, &rel_orig_z, &verbose))
if ((error_code = PANPHASIA_init_level_(&rel_lev,
&rel_orig_x, &rel_orig_y, &rel_orig_z, &verbose)))
{
printf("Error %d in initialing PANPHASIA_init_level_\n",
error_code);
@ -704,9 +704,9 @@ int demo_descriptor_()
for (int i = 0; i < Nbasis / 3; i++)
copy_list[i] = 3 * i;
if (error_code = PANPHASIA_compute_coefficients_(&xstart, &ystart, &zstart,
&xextent, &yextent, &zextent, copy_list, &ncopy,
output_values, &flag_output_mode, &verbose))
if ((error_code = PANPHASIA_compute_coefficients_(&xstart, &ystart, &zstart,
&xextent, &yextent, &zextent, copy_list, &ncopy,
output_values, &flag_output_mode, &verbose)))
{
printf("Error %d in PANPHASIA_compute_coefficients \n", error_code);
@ -721,13 +721,13 @@ int demo_descriptor_()
for (size_t xco = 0; xco < xextent; xco++)
for (size_t yco = 0; yco < yextent; yco++)
for (size_t zco = 0; zco < zextent; zco++)
fprintf(file, "%llu %llu %llu %f\n", xco, yco, zco, output_values[ncopy * (xco * yextent * zextent + yco * zextent + zco)]);
fprintf(file, "%lu %lu %lu %f\n", xco, yco, zco, output_values[ncopy * (xco * yextent * zextent + yco * zextent + zco)]);
fclose(file);
};
return (0);
};
}
int PANPHASIA_init_descriptor_(const char *descriptor, int *verbose)
{
@ -758,7 +758,7 @@ int PANPHASIA_init_descriptor_(const char *descriptor, int *verbose)
check_panphasia_key_(verb);
int pan_mode;
if (error = parse_and_validate_descriptor_(descriptor,&pan_mode))
if ((error = parse_and_validate_descriptor_(descriptor, &pan_mode)))
{
printf("-----------------------------------------\n");
printf("Error initating start-up Panphasia routines \n");
@ -771,7 +771,7 @@ int PANPHASIA_init_descriptor_(const char *descriptor, int *verbose)
if (*verbose)
printf("Sucessfully started Panphasia with the descriptor:\n%s\n", descriptor);
return (0);
};
}
/////////////////////////////////////////////////////////////////////////////////
@ -838,7 +838,7 @@ void PANPHASIA_init_descriptor_checks()
printf("===================================================\n");
panphasia_rel_origin_set = 0; // Force user to set rel origin themselves.
};
}
int PANPHASIA_init_level_(size_t *rel_lev,
size_t *rel_orig_x, size_t *rel_orig_y,
@ -869,16 +869,16 @@ int PANPHASIA_init_level_(size_t *rel_lev,
{
printf("-----------------------------------------------------------------\n");
printf("Initialising a Panphasia subgrid\n");
printf("Relative level %llu\n", rel_level);
printf("Relative origin (%llu,%llu,%llu)\n", rel_origin_x, rel_origin_y, rel_origin_z);
printf("The maximum possible extent of this subgrid is %llu cells\n", rel_coord_max);
printf("Relative level %lu\n", rel_level);
printf("Relative origin (%lu,%lu,%lu)\n", rel_origin_x, rel_origin_y, rel_origin_z);
printf("The maximum possible extent of this subgrid is %lu cells\n", rel_coord_max);
printf("-----------------------------------------------------------------\n");
};
panphasia_rel_origin_set = 1;
return (0);
};
}
//======================================================================================
//======================================================================================
@ -910,9 +910,12 @@ int PANPHASIA_compute_coefficients_(size_t *xstart, size_t *ystart, size_t *zsta
if (*zstart >= rel_coord_max)
return (203);
if (*xextent > rel_coord_max) return (204);
if (*yextent > rel_coord_max) return (205);
if (*zextent > rel_coord_max) return (206);
if (*xextent > rel_coord_max)
return (204);
if (*yextent > rel_coord_max)
return (205);
if (*zextent > rel_coord_max)
return (206);
if ((*ncopy < 0) || (*ncopy > Nbasis))
return (207);
@ -920,7 +923,8 @@ int PANPHASIA_compute_coefficients_(size_t *xstart, size_t *ystart, size_t *zsta
if ((copy_list[0] < 0) || (copy_list[*ncopy - 1] >= Nbasis))
return (208);
if ((*xextent==0)||(*yextent==0)||(*zextent==0)) return(0);
if ((*xextent == 0) || (*yextent == 0) || (*zextent == 0))
return (0);
for (int i = 1; i < *ncopy; i++)
if (copy_list[i] <= copy_list[i - 1])
@ -1009,17 +1013,17 @@ int PANPHASIA_compute_coefficients_(size_t *xstart, size_t *ystart, size_t *zsta
{
int error_code;
if (error_code = return_binary_tree_cell_lists(level_max, list_cell_x_coord,
*xextent, ret_x_list_coords, nreturn_x, child_pointer_x,
level_count_x, level_begin_x, index_perm_x))
if ((error_code = return_binary_tree_cell_lists(level_max, list_cell_x_coord,
*xextent, ret_x_list_coords, nreturn_x, child_pointer_x,
level_count_x, level_begin_x, index_perm_x)))
return (error_code);
if (error_code = return_binary_tree_cell_lists(level_max, list_cell_y_coord,
*yextent, ret_y_list_coords, nreturn_y, child_pointer_y,
level_count_y, level_begin_y, index_perm_y))
if ((error_code = return_binary_tree_cell_lists(level_max, list_cell_y_coord,
*yextent, ret_y_list_coords, nreturn_y, child_pointer_y,
level_count_y, level_begin_y, index_perm_y)))
return (error_code);
if (error_code = return_binary_tree_cell_lists(level_max, list_cell_z_coord,
*zextent, ret_z_list_coords, nreturn_z, child_pointer_z,
level_count_z, level_begin_z, index_perm_z))
if ((error_code = return_binary_tree_cell_lists(level_max, list_cell_z_coord,
*zextent, ret_z_list_coords, nreturn_z, child_pointer_z,
level_count_z, level_begin_z, index_perm_z)))
return (error_code);
};
//===================================================================
@ -1035,7 +1039,7 @@ int PANPHASIA_compute_coefficients_(size_t *xstart, size_t *ystart, size_t *zsta
};
if (*verbose)
printf("Total number cells: %llu \n", number_of_cells);
printf("Total number cells: %lu \n", number_of_cells);
cell_memory_to_allocate = sizeof(PAN_REAL) * number_of_cells * Nbasis;
};
@ -1047,13 +1051,13 @@ int PANPHASIA_compute_coefficients_(size_t *xstart, size_t *ystart, size_t *zsta
//========================================================================================
// Loop over octree starting at the root, for all relevant cells at each level
//========================================================================================
size_t total_number_cells = 0;
size_t num_cell_compute = 0;
size_t num_level_max_cells = 0;
size_t total_num_children = 0;
// size_t total_number_cells = 0;
// size_t num_cell_compute = 0;
// size_t num_level_max_cells = 0;
// size_t total_num_children = 0;
{
size_t cell_index, j1, j2, j3;
size_t child_cells[8];
// size_t child_cells[8];
size_t xoffset, yoffset, zoffset;
size_t ix, iy, iz;
size_t xco, yco, zco;
@ -1135,7 +1139,7 @@ int PANPHASIA_compute_coefficients_(size_t *xstart, size_t *ystart, size_t *zsta
}; // end loop over possible children
if (*verbose > 1)
printf("Cell: L%llu %llu %llu %llu\n", level, j1, j2, j3);
printf("Cell: L%lu %lu %lu %lu\n", level, j1, j2, j3);
}; // z/y/x-coordinate/level
@ -1262,31 +1266,31 @@ int parse_and_validate_descriptor_(const char *descriptor, int *pan_mode)
switch (nelement)
{
case 1:
if (sscanf(token, "Panph%llu", &desc_order) != 1)
if (sscanf(token, "Panph%lu", &desc_order) != 1)
return (440001);
break;
case 2:
if (sscanf(token, "L%llu", &desc_level) != 1)
if (sscanf(token, "L%lu", &desc_level) != 1)
return 440002;
break;
case 3:
if (sscanf(token, "%llu", &desc_x) != 1)
if (sscanf(token, "%lu", &desc_x) != 1)
return 440003;
break;
case 4:
if (sscanf(token, "%llu", &desc_y) != 1)
if (sscanf(token, "%lu", &desc_y) != 1)
return 440004;
break;
case 5:
if (sscanf(token, "%llu", &desc_z) != 1)
if (sscanf(token, "%lu", &desc_z) != 1)
return 440005;
break;
case 6:
if (sscanf(token, "S%llu", &desc_size) != 1)
if (sscanf(token, "S%lu", &desc_size) != 1)
return 440005;
break;
case 7:
if (sscanf(token, "KK%lld", &desc_kk_limit) == 1)
if (sscanf(token, "KK%lu", &desc_kk_limit) == 1)
{
kk_limit_set = 1;
token = strtok(NULL, split);
@ -1295,7 +1299,7 @@ int parse_and_validate_descriptor_(const char *descriptor, int *pan_mode)
return 440006;
break;
case 8:
if (sscanf(token, "%s", &desc_name) != 1)
if (sscanf(token, "%199s", desc_name) != 1)
return 440007;
break;
}
@ -1304,12 +1308,12 @@ int parse_and_validate_descriptor_(const char *descriptor, int *pan_mode)
if (kk_limit_set == 0)
{
sprintf(descriptor_as_read, "[Panph%llu,L%llu,(%llu,%llu,%llu),S%llu,CH%lld,%s]",
sprintf(descriptor_as_read, "[Panph%lu,L%lu,(%lu,%lu,%lu),S%lu,CH%lld,%s]",
desc_order, desc_level, desc_x, desc_y, desc_z, desc_size, desc_ch, desc_name);
}
else
{
sprintf(descriptor_as_read, "[Panph%llu,L%llu,(%llu,%llu,%llu),S%llu,KK%lld,CH%lld,%s]",
sprintf(descriptor_as_read, "[Panph%lu,L%lu,(%lu,%lu,%lu),S%lu,KK%ld,CH%lld,%s]",
desc_order, desc_level, desc_x, desc_y, desc_z, desc_size, desc_kk_limit, desc_ch, desc_name);
}
@ -1334,7 +1338,7 @@ int parse_and_validate_descriptor_(const char *descriptor, int *pan_mode)
strcpy(full_descriptor, descriptor);
descriptor_read_in = 1;
*pan_mode = (desc_order==1)? 0:1; // 0 - Old descriptor: 1 HO descriptor
*pan_mode = (desc_order == 1) ? 0 : 1; // 0 - Old descriptor: 1 HO descriptor
comp_ch = compute_check_digit_(); // check the check digit
@ -1359,7 +1363,7 @@ void calc_absolute_coordinates(size_t xrel, size_t yrel, size_t zrel, size_t *xa
// printf("descriptor_zorigin %llu rel_level %llu zrel %llu rel_origin_z %llu rel_coord_max %llu \n descriptor_base_level %llu, zabs %llu\n",
//y descriptor_zorigin,rel_level,zrel,rel_origin_z,rel_coord_max,descriptor_base_level,*zabs);
};
}
int cell_information(size_t cell_id, size_t *cumulative_cell_index, size_t *cuboid_x_dimen,
size_t *cuboid_y_dimen, size_t *cuboid_z_dimen, size_t *cell_lev,
@ -1452,44 +1456,44 @@ int return_binary_tree_cell_lists(size_t level_max, size_t *list_cell_coordinate
return (0);
}
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
//
// Test code for checking the appropriate moments are preserved
// between levels in Panphasia
//
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
//
// Test code for checking the appropriate moments are preserved
// between levels in Panphasia
//
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
#include <gsl/gsl_sf_legendre.h>
void integrate_cell(int, int, int, size_t, size_t, size_t, FFTW_REAL *, double *);
int compute_panphasia_(double, double, double, size_t, size_t, size_t, FFTW_REAL *, double *);
void test_cell_moments(char *, size_t, size_t, size_t, size_t, size_t, double *);
void test_cell_moments(const char *, size_t, size_t, size_t, size_t, size_t, double *);
////////////////////////////////////////////////////////////////////////////////
void test_moments_()
{
int lev = 10;
// int lev = 10;
char descriptor_demo[300] = "Hello!";
printf("Demo string %s\n", descriptor_demo);
// descriptor_pair_generate_();//, descriptor_demo);
printf("Parameters: %s\n", descriptor_demo);
size_t nlevel = 1;
// size_t nlevel = 1;
double coefficients1[Nbasis];
double coefficients2[Nbasis];
@ -1515,7 +1519,7 @@ void test_moments_()
yco = (yco_full) >> (63 - level);
zco = (zco_full) >> (63 - level);
sprintf(descriptor, "[Panph6,L%ld,(%llu,%llu,%llu),S1,CH-999,test]", level, xco, yco, zco);
sprintf(descriptor, "[Panph6,L%ld,(%lu,%lu,%lu),S1,CH-999,test]", level, xco, yco, zco);
// printf("%s\n",descriptor);
test_cell_moments(descriptor, 0, 0, 0, 0, 1, coefficients1);
@ -1560,7 +1564,7 @@ void test_moments_()
printf("Completed moment test successfully.\n");
}
void test_cell_moments(char root_descriptor[200], size_t rel_lev, size_t rel_orig_x,
void test_cell_moments(const char *root_descriptor, size_t rel_lev, size_t rel_orig_x,
size_t rel_orig_y, size_t rel_orig_z, size_t extent, double *coeff)
{
@ -1572,8 +1576,8 @@ void test_cell_moments(char root_descriptor[200], size_t rel_lev, size_t rel_ori
verbose = 0;
if (error_code = PANPHASIA_init_level_(&rel_lev,
&rel_orig_x, &rel_orig_y, &rel_orig_z, &verbose))
if ((error_code = PANPHASIA_init_level_(&rel_lev,
&rel_orig_x, &rel_orig_y, &rel_orig_z, &verbose)))
{
printf("Error %d in initialing PANPHASIA_init_level_\n",
error_code);
@ -1598,9 +1602,9 @@ void test_cell_moments(char root_descriptor[200], size_t rel_lev, size_t rel_ori
abort();
}
if (error_code = PANPHASIA_compute_coefficients_(&xstart, &ystart, &zstart,
&xextent, &yextent, &zextent, copy_list, &ncopy,
output_values, &flag_output_mode, &verbose))
if ((error_code = PANPHASIA_compute_coefficients_(&xstart, &ystart, &zstart,
&xextent, &yextent, &zextent, copy_list, &ncopy,
output_values, &flag_output_mode, &verbose)))
{
printf("Error %d in PANPHASIA_compute_coefficients_ \n", error_code);
@ -1622,7 +1626,8 @@ void test_cell_moments(char root_descriptor[200], size_t rel_lev, size_t rel_ori
*/
double sum_coefficients[Nbasis];
for( size_t i=0; i<Nbasis; ++i ) sum_coefficients[i] = 0.0;
for (size_t i = 0; i < Nbasis; ++i)
sum_coefficients[i] = 0.0;
double results[Nbasis];
size_t num_cells = 0;
@ -1647,15 +1652,15 @@ void test_cell_moments(char root_descriptor[200], size_t rel_lev, size_t rel_ori
void integrate_cell(int ix, int iy, int iz, size_t xextent, size_t yextent, size_t zextent, FFTW_REAL *output_values, double *results)
{
/////////////////////////////////////////////////////////////////////////////
//
// This function computes the integral over a cell of the product of the
// Panphasia field with an 'analysing' Legendre polynomial. As the
// integrand is a polynomial, Gaussian quadrature can be used for
// integration as it is exact up to rounding error provide p_order
// is less than 10.
//
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
//
// This function computes the integral over a cell of the product of the
// Panphasia field with an 'analysing' Legendre polynomial. As the
// integrand is a polynomial, Gaussian quadrature can be used for
// integration as it is exact up to rounding error provide p_order
// is less than 10.
//
/////////////////////////////////////////////////////////////////////////////
const double GQ_weights[5] = {0.2955242247147529, 0.2692667193099963,
0.2190863625159820, 0.1494513491505806,
@ -1695,9 +1700,10 @@ void integrate_cell(int ix, int iy, int iz, size_t xextent, size_t yextent, size
double range = 0.5 * (b - a);
double sum[Nbasis];
for( size_t i=0; i<Nbasis; ++i ) sum[i] = 0.0;
for (size_t i = 0; i < Nbasis; ++i)
sum[i] = 0.0;
double test_sum = 0.0;
// double test_sum = 0.0;
for (int i = 0; i < 10; i++)
{
@ -1840,7 +1846,7 @@ void compute_sph_bessel_coeffs(int nfft, int pmax, int n4dimen, int fdim, double
int j = (i <= nfft / 2) ? i : i - nfft;
int k = abs(j);
double sign = (j < 0) ? pow(-1.0, l) : 1.0;
double x = pi*(double)fdim*(double)k/(double)nfft;
double x = pi * (double)fdim * (double)k / (double)nfft;
double result;
spherical_bessel_(&l, &x, &result);

806
external/panphasia_ho/pan_mpi_routines.c vendored
View file

@ -13,7 +13,6 @@
#include <omp.h>
#endif
extern const int Nbasis;
extern const int irank_p[3][84];
@ -22,556 +21,513 @@ extern size_t descriptor_kk_limit;
extern size_t descriptor_base_size;
int PANPHASIA_compute_kspace_field_(size_t relative_level, ptrdiff_t N0_fourier_grid,
ptrdiff_t local_n0_fourier_return, ptrdiff_t local_0_start_fourier_return,
FFTW_COMPLEX *return_field)
ptrdiff_t local_n0_fourier_return, ptrdiff_t local_0_start_fourier_return,
FFTW_COMPLEX *return_field)
{
size_t copy_list[Nbasis];
size_t copy_list[Nbasis];
int pmax = 6;
int nsubdivide = 21; //(pmax%2==0)?pmax+1:pmax+2;
size_t ncopy = (pmax + 1) * (pmax + 2) * (pmax + 3) / 6;
int pmax = 6;
if (ncopy % nsubdivide != 0)
return (100010);
int nchunk = ncopy / nsubdivide;
int nsubdivide = 21; //(pmax%2==0)?pmax+1:pmax+2;
int verbose = 1;
int flag_output_mode = 2;
int error;
ptrdiff_t size_to_alloc_fourier;
ptrdiff_t size_to_alloc_pan;
ptrdiff_t local_fourier_x_start, local_fourier_x_end;
FFTW_PLAN output_coeff_forward_plan;
size_t ncopy = (pmax+1)*(pmax+2)*(pmax+3)/6;
ptrdiff_t N0_pan_grid = descriptor_base_size << relative_level;
if (ncopy%nsubdivide!=0) return(100010);
int nchunk = ncopy/nsubdivide;
if (N0_fourier_grid % N0_pan_grid != 0)
return (100015);
int verbose = 1;
int flag_output_mode=2;
int error;
ptrdiff_t size_to_alloc_fourier;
ptrdiff_t size_to_alloc_pan;
ptrdiff_t local_fourier_x_start, local_fourier_x_end;
FFTW_PLAN output_coeff_forward_plan;
int fdim = N0_fourier_grid / N0_pan_grid;
size_t nfft_dim = N0_fourier_grid;
size_t npan_dim = N0_pan_grid;
ptrdiff_t N0_pan_grid = descriptor_base_size<<relative_level;
int SHARED_FOUR_PAN_SPACE = (nsubdivide == 1) && (fdim == 1) && (sizeof(PAN_REAL) == sizeof(FFTW_REAL));
if (N0_fourier_grid%N0_pan_grid!=0) return (100015);
////////////////////////////////////////////////////////////////////////////////////
int fdim = N0_fourier_grid/N0_pan_grid;
size_t nfft_dim = N0_fourier_grid;
size_t npan_dim = N0_pan_grid;
if (pmax > descriptor_order)
return (100020);
int SHARED_FOUR_PAN_SPACE = (nsubdivide==1)&&(fdim==1)&&(sizeof(PAN_REAL)==sizeof(FFTW_REAL));
for (size_t i = 0; i < Nbasis; i++)
copy_list[i] = i;
//printf("Dimensions of FT (%td,%td,%td)\n",N0_fourier_grid,N0_fourier_grid,N0_fourier_grid);
//printf("Dimensions of PG (%td,%td,%td)\n",N0_pan_grid,N0_pan_grid,N0_pan_grid);
//printf("local_no %td local_0_start_fourier %td\n",local_n0_fourier_return, local_0_start_fourier_return);
// Compute 1-D Spherical Bessel coefficients for each order //////////////////
// These are needed for the convolutions below //////////////////
size_t n4dimen;
////////////////////////////////////////////////////////////////////////////////////
n4dimen = (nfft_dim % 4 == 0) ? 4 * (nfft_dim / 4) + 4 : 4 * (nfft_dim / 4) + 5;
if (pmax>descriptor_order) return(100020);
double complex *sph_bessel_coeff = FFTW_MALLOC(sizeof(double complex) * n4dimen * (pmax + 1));
for (size_t i=0; i<Nbasis; i++) copy_list[i]=i;
if (sph_bessel_coeff == NULL)
return (100030);
//printf("Dimensions of FT (%td,%td,%td)\n",N0_fourier_grid,N0_fourier_grid,N0_fourier_grid);
//printf("Dimensions of PG (%td,%td,%td)\n",N0_pan_grid,N0_pan_grid,N0_pan_grid);
//printf("local_no %td local_0_start_fourier %td\n",local_n0_fourier_return, local_0_start_fourier_return);
compute_sph_bessel_coeffs(nfft_dim, pmax, n4dimen, fdim, sph_bessel_coeff);
// Compute 1-D Spherical Bessel coefficients for each order //////////////////
// These are needed for the convolutions below //////////////////
size_t n4dimen;
//printf("Reached here! ndimen4 %ld\n",n4dimen);
///////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
// Determine sizes of Fourier and Panphasia coefficient 3-D arrays //
ptrdiff_t local_n0_pan;
ptrdiff_t local_0_start_pan;
{
n4dimen=(nfft_dim%4==0) ? 4*(nfft_dim/4)+4 : 4*(nfft_dim/4)+5;
int rank = 3;
double complex *sph_bessel_coeff = FFTW_MALLOC(sizeof(double complex)*n4dimen*(pmax+1));
ptrdiff_t local_n0_fourier;
ptrdiff_t local_0_start_fourier;
const ptrdiff_t ndimens_alloc_fourier[] = {N0_fourier_grid, N0_fourier_grid, N0_fourier_grid + 2}; // Allocated for r2c
ptrdiff_t howmany = nchunk;
size_to_alloc_fourier = FFTW_MPI_LOCAL_SIZE_MANY(rank, ndimens_alloc_fourier, howmany,
FFTW_MPI_DEFAULT_BLOCK, MPI_COMM_WORLD,
&local_n0_fourier, &local_0_start_fourier);
if (sph_bessel_coeff==NULL) return(100030);
if (local_0_start_fourier != local_0_start_fourier_return)
{
printf("Error local_0_start_fourier!=local_0_start_fourier_return\n");
return (100032);
};
compute_sph_bessel_coeffs(nfft_dim, pmax, n4dimen, fdim, sph_bessel_coeff);
if (local_n0_fourier != local_n0_fourier_return)
{
printf("Error local_n0_fourier!=local_n0_fourier_return\n");
return (100033);
};
//printf("Reached here! ndimen4 %ld\n",n4dimen);
///////////////////////////////////////////////////////////////////////////////
ptrdiff_t abs_fourier_x_start, abs_fourier_x_end;
if (local_0_start_fourier_return % fdim == 0)
{
abs_fourier_x_start = local_0_start_fourier_return;
}
else
{
abs_fourier_x_start = local_0_start_fourier_return + fdim - (local_0_start_fourier_return % fdim);
};
if ((local_0_start_fourier_return + local_n0_fourier_return - 1) % fdim == 0)
{
abs_fourier_x_end = local_0_start_fourier_return + local_n0_fourier_return - 1;
}
else
{
abs_fourier_x_end = (local_0_start_fourier_return + local_n0_fourier_return - 1) -
((local_0_start_fourier_return + local_n0_fourier_return - 1) % fdim);
};
if ((abs_fourier_x_end - abs_fourier_x_start) % fdim != 0)
return (100036);
//////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////
// Determine sizes of Fourier and Panphasia coefficient 3-D arrays //
local_0_start_pan = abs_fourier_x_start / fdim;
local_n0_pan = 1 + (abs_fourier_x_end - abs_fourier_x_start) / fdim;
local_fourier_x_start = abs_fourier_x_start - local_0_start_fourier_return;
local_fourier_x_end = abs_fourier_x_end - local_0_start_fourier_return;
ptrdiff_t local_n0_pan;
ptrdiff_t local_0_start_pan;
{
int rank =3;
// const ptrdiff_t ndimens_alloc_pan[] = {N0_pan_grid, N0_pan_grid, N0_pan_grid+2};
howmany = ncopy;
ptrdiff_t local_n0_fourier;
ptrdiff_t local_0_start_fourier;
const ptrdiff_t ndimens_alloc_fourier[] = {N0_fourier_grid, N0_fourier_grid, N0_fourier_grid+2}; // Allocated for r2c
ptrdiff_t howmany = nchunk;
size_to_alloc_fourier = FFTW_MPI_LOCAL_SIZE_MANY(rank, ndimens_alloc_fourier, howmany,
FFTW_MPI_DEFAULT_BLOCK,MPI_COMM_WORLD,
&local_n0_fourier,&local_0_start_fourier);
if (local_0_start_fourier!=local_0_start_fourier_return){
printf("Error local_0_start_fourier!=local_0_start_fourier_return\n");
return(100032);
};
size_to_alloc_pan = howmany * local_n0_pan * N0_pan_grid * (N0_pan_grid + 2);
if (local_n0_fourier!=local_n0_fourier_return){
printf("Error local_n0_fourier!=local_n0_fourier_return\n");
return(100033);
};
//printf("size_to_alloc_fourier = %td\n",size_to_alloc_fourier);
//printf("size_to_alloc_pan = %td\n",size_to_alloc_pan);
};
/////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
// Allocate arrays to store Panphasia coefficients and Fourier information
// If nsubdivide==1 then use the same structure to store both.
ptrdiff_t abs_fourier_x_start, abs_fourier_x_end;
void *panphasia_coefficients = FFTW_MALLOC(sizeof(PAN_REAL) * size_to_alloc_pan);
void *fourier_grids;
if (local_0_start_fourier_return%fdim==0){
abs_fourier_x_start = local_0_start_fourier_return;
}else{
abs_fourier_x_start = local_0_start_fourier_return + fdim - (local_0_start_fourier_return%fdim);
};
if (panphasia_coefficients == NULL)
return (100040);
void *mode_weightings = FFTW_MALLOC(sizeof(FFTW_REAL) * size_to_alloc_fourier / nchunk);
if ((local_0_start_fourier_return + local_n0_fourier_return - 1)%fdim==0){
abs_fourier_x_end = local_0_start_fourier_return + local_n0_fourier_return - 1;
}else{
abs_fourier_x_end = (local_0_start_fourier_return + local_n0_fourier_return - 1)-
((local_0_start_fourier_return + local_n0_fourier_return - 1)%fdim);
};
FFTW_REAL *ptr_mode_weightings;
ptr_mode_weightings = mode_weightings;
memset(ptr_mode_weightings, 0, sizeof(FFTW_REAL) * size_to_alloc_fourier / nchunk);
FFTW_REAL *ptr_real_fourier_grid;
FFTW_REAL *ptr_panphasia_coefficients = panphasia_coefficients;
FFTW_COMPLEX *ptr_cmplx_fourier_grid;
if (SHARED_FOUR_PAN_SPACE)
{
ptr_real_fourier_grid = panphasia_coefficients;
ptr_cmplx_fourier_grid = panphasia_coefficients;
}
else
{
fourier_grids = FFTW_MALLOC(sizeof(FFTW_REAL) * size_to_alloc_fourier);
if (fourier_grids == NULL)
return (100041);
ptr_real_fourier_grid = fourier_grids;
ptr_cmplx_fourier_grid = fourier_grids;
};
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
if ((abs_fourier_x_end-abs_fourier_x_start)%fdim!=0) return (100036);
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
//////////// Compute the Panphasia coefficients ////////////////////////////////
local_0_start_pan = abs_fourier_x_start/fdim;
local_n0_pan = 1 + (abs_fourier_x_end-abs_fourier_x_start)/fdim;
{
local_fourier_x_start = abs_fourier_x_start - local_0_start_fourier_return;
local_fourier_x_end = abs_fourier_x_end - local_0_start_fourier_return;
size_t xorigin = local_0_start_pan, yorigin = 0, zorigin = 0;
const ptrdiff_t ndimens_alloc_pan[] = {N0_pan_grid, N0_pan_grid, N0_pan_grid+2};
howmany = ncopy;
size_t xextent = local_n0_pan, yextent = N0_pan_grid, zextent = N0_pan_grid;
size_to_alloc_pan = howmany * local_n0_pan * N0_pan_grid * ( N0_pan_grid+2);
if ((error = PANPHASIA_compute_coefficients_(&xorigin, &yorigin, &zorigin, &xextent, &yextent,
&zextent, copy_list, &ncopy, panphasia_coefficients,
&flag_output_mode, &verbose)))
return (100100 + error);
};
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
////////// Output diagnostics for small grids only
//printf("size_to_alloc_fourier = %td\n",size_to_alloc_fourier);
//printf("size_to_alloc_pan = %td\n",size_to_alloc_pan);
{
};
/////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////
if (N0_pan_grid < -1)
{
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
char filename[100];
sprintf(filename, "output_real_space_field.%d", rank);
FILE *fp;
printf("local_n0_pan %ld npan_dim %ld N0_pan_grid %ld\n", local_n0_pan,
npan_dim, N0_pan_grid);
fp = fopen(filename, "w");
for (int ix = 0; ix < local_n0_pan; ix++)
for (int iy = 0; iy < npan_dim; iy++)
for (int iz = 0; iz < npan_dim; iz++)
{
int index = ix * N0_pan_grid * (N0_pan_grid + 2) + iy * (N0_pan_grid + 2) + iz;
fprintf(fp, "%6lu%6d%6d %14.8lf %d\n", ix + local_0_start_pan, iy, iz,
ptr_panphasia_coefficients[index], index);
};
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
// Allocate arrays to store Panphasia coefficients and Fourier information
// If nsubdivide==1 then use the same structure to store both.
fclose(fp);
};
};
void *panphasia_coefficients= FFTW_MALLOC(sizeof(PAN_REAL)*size_to_alloc_pan);
void *fourier_grids;
//----------------------------------------------------------------------------------
////////////// Set up FTTW plan //////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
if (panphasia_coefficients==NULL) return(100040);
// printf("Making the plan ... \n");
void *mode_weightings = FFTW_MALLOC(sizeof(FFTW_REAL)*size_to_alloc_fourier/nchunk);
//////////////////// Make plan for ncopy interleaved FTs ///////////////////////////
//////////////////////////////////////////////////////////////////////////
FFTW_REAL *ptr_mode_weightings;
ptr_mode_weightings = mode_weightings;
memset(ptr_mode_weightings, 0, sizeof(FFTW_REAL)*size_to_alloc_fourier/nchunk);
{
int rank = 3;
const ptrdiff_t ndimens[3] = {N0_fourier_grid, N0_fourier_grid, N0_fourier_grid};
ptrdiff_t howmany = nchunk;
ptrdiff_t block = FFTW_MPI_DEFAULT_BLOCK;
ptrdiff_t tblock = FFTW_MPI_DEFAULT_BLOCK;
unsigned flags = FFTW_ESTIMATE;
output_coeff_forward_plan = FFTW_MPI_PLAN_MANY_DTF_R2C(rank, ndimens,
howmany, block, tblock,
ptr_real_fourier_grid, ptr_cmplx_fourier_grid,
MPI_COMM_WORLD, flags);
if (output_coeff_forward_plan == NULL)
{
printf("Null plan\n");
return (100051);
};
};
FFTW_REAL *ptr_real_fourier_grid;
FFTW_REAL *ptr_panphasia_coefficients = panphasia_coefficients;
FFTW_COMPLEX *ptr_cmplx_fourier_grid;
//printf("Plan completed ... \n");
if (SHARED_FOUR_PAN_SPACE){
ptr_real_fourier_grid = panphasia_coefficients;
ptr_cmplx_fourier_grid = panphasia_coefficients;
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//----------------------------------------------------------------------------------
}else{
memset(return_field, 0, local_n0_fourier_return * N0_fourier_grid * (N0_fourier_grid + 2) * sizeof(FFTW_REAL));
fourier_grids= FFTW_MALLOC(sizeof(FFTW_REAL)*size_to_alloc_fourier);
if (fourier_grids==NULL) return(100041);
for (int iter = 0; iter < nsubdivide; iter++)
{
ptr_real_fourier_grid = fourier_grids;
ptr_cmplx_fourier_grid = fourier_grids;
};
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
int moffset = iter * nchunk;
if (!SHARED_FOUR_PAN_SPACE)
{
memset(ptr_real_fourier_grid, 0, sizeof(FFTW_REAL) * size_to_alloc_fourier);
// Copy Panphasia coefficients to Fourier grid with appropriate stride - fdim
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
//////////// Compute the Panphasia coefficients ////////////////////////////////
size_t index_p, index_f;
int m;
{
size_t xorigin= local_0_start_pan, yorigin=0, zorigin=0;
size_t xextent =local_n0_pan, yextent = N0_pan_grid, zextent = N0_pan_grid;
if (error = PANPHASIA_compute_coefficients_(&xorigin,&yorigin,&zorigin,&xextent,&yextent,
&zextent, copy_list, &ncopy, panphasia_coefficients,
&flag_output_mode,&verbose)) return(100100+error);
};
///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
////////// Output diagnostics for small grids only
{
if (N0_pan_grid<-1){
int rank;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
char filename[100];
sprintf(filename,"output_real_space_field.%d",rank);
FILE *fp;
printf("local_n0_pan %ld npan_dim %ld N0_pan_grid %ld\n",local_n0_pan,
npan_dim,N0_pan_grid);
fp = fopen(filename,"w");
for (int ix=0; ix<local_n0_pan; ix++)
for (int iy=0; iy < npan_dim; iy++)
for (int iz=0; iz < npan_dim; iz++){
int index = ix*N0_pan_grid*(N0_pan_grid+2) + iy*(N0_pan_grid+2) + iz;
fprintf(fp,"%6llu%6d%6d %14.8lf %d\n",ix+local_0_start_pan,iy,iz,
ptr_panphasia_coefficients[index],index);
};
fclose(fp);
};
};
//----------------------------------------------------------------------------------
////////////// Set up FTTW plan //////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
// printf("Making the plan ... \n");
//////////////////// Make plan for ncopy interleaved FTs ///////////////////////////
//////////////////////////////////////////////////////////////////////////
{
int rank = 3;
const ptrdiff_t ndimens[3] = {N0_fourier_grid, N0_fourier_grid, N0_fourier_grid};
ptrdiff_t howmany = nchunk;
ptrdiff_t block = FFTW_MPI_DEFAULT_BLOCK;
ptrdiff_t tblock = FFTW_MPI_DEFAULT_BLOCK;
unsigned flags = FFTW_ESTIMATE;
output_coeff_forward_plan = FFTW_MPI_PLAN_MANY_DTF_R2C(rank, ndimens,
howmany, block, tblock,
ptr_real_fourier_grid, ptr_cmplx_fourier_grid,
MPI_COMM_WORLD, flags);
if (output_coeff_forward_plan==NULL) {
printf("Null plan\n");
return(100051);
};
};
//printf("Plan completed ... \n");
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
//----------------------------------------------------------------------------------
memset(return_field, 0, local_n0_fourier_return*N0_fourier_grid *(N0_fourier_grid +2) * sizeof(FFTW_REAL));
for (int iter = 0; iter < nsubdivide; iter++){
int moffset = iter*nchunk;
if (!SHARED_FOUR_PAN_SPACE){
memset(ptr_real_fourier_grid, 0, sizeof(FFTW_REAL)*size_to_alloc_fourier);
// Copy Panphasia coefficients to Fourier grid with appropriate stride - fdim
size_t index_p,index_f;
int m;
for (int ix_f = local_fourier_x_start, ix_p=0; ix_f <=local_fourier_x_end; ix_f+=fdim,ix_p++)
for (int ix_f = local_fourier_x_start, ix_p = 0; ix_f <= local_fourier_x_end; ix_f += fdim, ix_p++)
#ifdef USE_OPENMP
#pragma omp parallel for collapse(2) private (index_p,index_f,m)
#pragma omp parallel for collapse(2) private(index_p, index_f, m)
#endif
for(int iy_f=0, iy_p=0 ; iy_f<nfft_dim; iy_f+=fdim, iy_p++)
for(int iz_f=0, iz_p=0; iz_f<nfft_dim; iz_f+=fdim, iz_p++){
index_p = ix_p*N0_pan_grid*(N0_pan_grid + 2) + iy_p*(N0_pan_grid + 2) + iz_p;
index_f = ix_f*N0_fourier_grid*(N0_fourier_grid + 2) + iy_f*(N0_fourier_grid + 2) + iz_f;
for (int iy_f = 0, iy_p = 0; iy_f < nfft_dim; iy_f += fdim, iy_p++)
for (int iz_f = 0, iz_p = 0; iz_f < nfft_dim; iz_f += fdim, iz_p++)
{
index_p = ix_p * N0_pan_grid * (N0_pan_grid + 2) + iy_p * (N0_pan_grid + 2) + iz_p;
index_f = ix_f * N0_fourier_grid * (N0_fourier_grid + 2) + iy_f * (N0_fourier_grid + 2) + iz_f;
for (m=0; m<nchunk; m++){
ptr_real_fourier_grid[nchunk*index_f + m] =
ptr_panphasia_coefficients[ncopy*index_p + moffset + m];
for (m = 0; m < nchunk; m++)
{
ptr_real_fourier_grid[nchunk * index_f + m] =
ptr_panphasia_coefficients[ncopy * index_p + moffset + m];
};
};
};
};
};
};
FFTW_MPI_EXECUTE_DFT_R2C(output_coeff_forward_plan,
ptr_real_fourier_grid, ptr_cmplx_fourier_grid);
FFTW_MPI_EXECUTE_DFT_R2C(output_coeff_forward_plan,
ptr_real_fourier_grid,ptr_cmplx_fourier_grid);
{
// Convolve and combine the FT of the Panphasia coefficient field
size_t index1,index2;
complex weight;
int m;
{
// Convolve and combine the FT of the Panphasia coefficient field
size_t index1, index2;
complex weight;
//int m;
#ifdef USE_OPENMP
#pragma omp parallel for collapse(3) \
private (index1,index2,weight,m)
#pragma omp parallel for collapse(3) private(index1, index2, weight, m)
#endif
for(int ix=0;ix<local_n0_fourier_return;ix++)
for(int iy=0;iy<nfft_dim;iy++)
for(int iz=0;iz<=nfft_dim/2;iz++){
index1 = ix*N0_fourier_grid*(N0_fourier_grid/2+1) + iy*(N0_fourier_grid/2+1) + iz;
for (int m=0; m<nchunk; m++){
index2 = nchunk*index1 + m;
for (int ix = 0; ix < local_n0_fourier_return; ix++)
for (int iy = 0; iy < nfft_dim; iy++)
for (int iz = 0; iz <= nfft_dim / 2; iz++)
{
index1 = ix * N0_fourier_grid * (N0_fourier_grid / 2 + 1) + iy * (N0_fourier_grid / 2 + 1) + iz;
for (int m = 0; m < nchunk; m++)
{
index2 = nchunk * index1 + m;
weight = sph_bessel_coeff[n4dimen*irank_p[0][m+moffset]+ix+local_0_start_fourier_return]*
sph_bessel_coeff[n4dimen*irank_p[1][m+moffset]+iy]*
sph_bessel_coeff[n4dimen*irank_p[2][m+moffset]+iz];
return_field[index1] += weight * ptr_cmplx_fourier_grid[index2];
ptr_mode_weightings[index1]+=cabs(weight)*cabs(weight);
};
};
};
weight = sph_bessel_coeff[n4dimen * irank_p[0][m + moffset] + ix + local_0_start_fourier_return] *
sph_bessel_coeff[n4dimen * irank_p[1][m + moffset] + iy] *
sph_bessel_coeff[n4dimen * irank_p[2][m + moffset] + iz];
return_field[index1] += weight * ptr_cmplx_fourier_grid[index2];
ptr_mode_weightings[index1] += cabs(weight) * cabs(weight);
};
};
};
}; // End loop over iter
}; // End loop over iter
// printf("Reached here 10!\n");
// printf("Reached here 10!\n");
//Add phase shift and normalise field
{
//Add phase shift and normalise field
{
double complex phase_shift_and_scale;
int kx, ky, kz;
const double pi = 4.0 * atan(1.0);
size_t index1;
double complex phase_shift_and_scale;
int kx,ky,kz;
const double pi = 4.0 * atan(1.0);
size_t index1;
double min_weight = 100.0;
double min_weight = 100.0;
#ifdef USE_OPENMP
#pragma omp parallel for collapse(3) \
private (index1,kx,ky,kz,phase_shift_and_scale)
#pragma omp parallel for collapse(3) private(index1, kx, ky, kz, phase_shift_and_scale)
#endif
for(int ix=0;ix<local_n0_fourier_return;ix++)
for(int iy=0;iy<nfft_dim;iy++)
for(int iz=0;iz<=nfft_dim/2;iz++){
index1 = ix*N0_fourier_grid*(N0_fourier_grid/2+1) + iy*(N0_fourier_grid/2+1) + iz;
kx = (ix+local_0_start_fourier_return>nfft_dim/2) ?
ix + local_0_start_fourier_return - nfft_dim : ix + local_0_start_fourier_return;
ky = (iy > nfft_dim/2) ? iy-nfft_dim : iy;
kz = (iz > nfft_dim/2) ? iz-nfft_dim : iz;
for (int ix = 0; ix < local_n0_fourier_return; ix++)
for (int iy = 0; iy < nfft_dim; iy++)
for (int iz = 0; iz <= nfft_dim / 2; iz++)
{
index1 = ix * N0_fourier_grid * (N0_fourier_grid / 2 + 1) + iy * (N0_fourier_grid / 2 + 1) + iz;
kx = (ix + local_0_start_fourier_return > nfft_dim / 2) ? ix + local_0_start_fourier_return - nfft_dim : ix + local_0_start_fourier_return;
ky = (iy > nfft_dim / 2) ? iy - nfft_dim : iy;
kz = (iz > nfft_dim / 2) ? iz - nfft_dim : iz;
if ( (kx==nfft_dim/2)||(ky==nfft_dim/2)||(kz==nfft_dim/2)){
// Set Nyquist modes to zero - not used by IC_Gen anyway.
phase_shift_and_scale = 0.0; //1.0/pow((double)nfft_dim,1.5); // No phase shift
ptr_mode_weightings[index1] = 0.0; // Set squared weight to zero
}else{
phase_shift_and_scale = sqrt( (double)(fdim*fdim*fdim))*
cexp( (double)fdim * (-I)*pi*(double)(kx + ky + kz)/sqrt(ptr_mode_weightings[index1])/
(double)nfft_dim)/pow((double)nfft_dim,1.5);
if ((kx == nfft_dim / 2) || (ky == nfft_dim / 2) || (kz == nfft_dim / 2))
{
// Set Nyquist modes to zero - not used by IC_Gen anyway.
phase_shift_and_scale = 0.0; //1.0/pow((double)nfft_dim,1.5); // No phase shift
ptr_mode_weightings[index1] = 0.0; // Set squared weight to zero
}
else
{
phase_shift_and_scale = sqrt((double)(fdim * fdim * fdim)) *
cexp((double)fdim * (-I) * pi * (double)(kx + ky + kz) / sqrt(ptr_mode_weightings[index1]) /
(double)nfft_dim) /
pow((double)nfft_dim, 1.5);
};
return_field[index1] *= phase_shift_and_scale;
if (ptr_mode_weightings[index1] < min_weight)
min_weight = ptr_mode_weightings[index1];
};
return_field[index1] *= phase_shift_and_scale;
if (ptr_mode_weightings[index1]<min_weight)
min_weight=ptr_mode_weightings[index1];
};
// printf("Minimum weight %lf\n",sqrt(min_weight));
};
// printf("Minimum weight %lf\n",sqrt(min_weight));
// printf("Reached here 11!\n");
// Rescale selected Fourier modes to unit amplitude.
// By default this part is not executed.
};
// printf("Reached here 11!\n");
// Rescale selected Fourier modes to unit amplitude.
// By default this part is not executed.
if (descriptor_kk_limit>0){
size_t index1;
complex weight;
size_t ksquared;
int kx,ky,kz;
if (descriptor_kk_limit > 0)
{
size_t index1;
complex weight;
size_t ksquared;
int kx, ky, kz;
#ifdef USE_OPENMP
#pragma omp parallel for collapse(3) \
private (index1,kx,ky,kz,ksquared,weight)
#pragma omp parallel for collapse(3) private(index1, kx, ky, kz, ksquared, weight)
#endif
for(int ix=0;ix<local_n0_fourier_return;ix++)
for(int iy=0;iy<nfft_dim;iy++)
for(int iz=0;iz<=nfft_dim/2;iz++){
kx = (ix+local_0_start_fourier_return>nfft_dim/2) ?
ix + local_0_start_fourier_return - nfft_dim : ix + local_0_start_fourier_return;
ky = (iy > nfft_dim/2) ? iy-nfft_dim : iy;
kz = (iz > nfft_dim/2) ? iz-nfft_dim : iz;
ksquared = kx*kx + ky*ky + kz*kz;
if ( (kx!=nfft_dim/2)&&(ky!=nfft_dim/2)&&(kz!=nfft_dim/2)){ //Omit Nyquist modes
for (int ix = 0; ix < local_n0_fourier_return; ix++)
for (int iy = 0; iy < nfft_dim; iy++)
for (int iz = 0; iz <= nfft_dim / 2; iz++)
{
kx = (ix + local_0_start_fourier_return > nfft_dim / 2) ? ix + local_0_start_fourier_return - nfft_dim : ix + local_0_start_fourier_return;
ky = (iy > nfft_dim / 2) ? iy - nfft_dim : iy;
kz = (iz > nfft_dim / 2) ? iz - nfft_dim : iz;
ksquared = kx * kx + ky * ky + kz * kz;
if ((kx != nfft_dim / 2) && (ky != nfft_dim / 2) && (kz != nfft_dim / 2))
{ //Omit Nyquist modes
if ((ksquared<=descriptor_kk_limit)&&(ksquared!=0)){
index1 = ix*N0_fourier_grid*(N0_fourier_grid/2+1) + iy*(N0_fourier_grid/2+1) + iz;
weight = cabs(return_field[index1]);
return_field[index1] /= weight;
};
};
};
};
//printf("Reached here 12!\n");
if ((ksquared <= descriptor_kk_limit) && (ksquared != 0))
{
index1 = ix * N0_fourier_grid * (N0_fourier_grid / 2 + 1) + iy * (N0_fourier_grid / 2 + 1) + iz;
weight = cabs(return_field[index1]);
return_field[index1] /= weight;
};
};
};
};
//printf("Reached here 12!\n");
int rank;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
char filename[100];
sprintf(filename,"output_k_space_alt.%d",rank);
sprintf(filename, "output_k_space_alt.%d", rank);
FILE *fp;
if (nfft_dim < -1)
{
FILE *fp;
if (nfft_dim <-1){
fp = fopen(filename, "w");
for (int ix = 0; ix < local_n0_fourier_return; ix++)
for (int iy = 0; iy < nfft_dim; iy++)
for (int iz = 0; iz <= nfft_dim / 2; iz++)
{
FILE *fp;
int index = ix * N0_fourier_grid * (N0_fourier_grid / 2 + 1) + iy * (N0_fourier_grid / 2 + 1) + iz;
fprintf(fp, "%6lu%6d%6d %14.8lf %14.8lf %14.8lf \n", ix + local_0_start_fourier_return, iy, iz,
creal(return_field[index]), cimag(return_field[index]), sqrt(ptr_mode_weightings[index]));
// ptr_mode_weightings[index]);
};
fclose(fp);
}
else
{
fp = fopen(filename,"w");
for (int ix=0; ix<local_n0_fourier_return; ix++)
for (int iy=0; iy < nfft_dim; iy++)
for (int iz=0; iz <= nfft_dim/2; iz++){
fp = fopen(filename, "w");
int index = ix*N0_fourier_grid*(N0_fourier_grid/2+1) + iy*(N0_fourier_grid/2+1) + iz;
fprintf(fp,"%6llu%6d%6d %14.8lf %14.8lf %14.8lf \n",ix+local_0_start_fourier_return,iy,iz,
creal(return_field[index]),cimag(return_field[index]),sqrt(ptr_mode_weightings[index]));
// ptr_mode_weightings[index]);
};
fclose(fp);
for (int ix = 0; ix < local_n0_fourier_return; ix++)
for (int iy = 0; iy < nfft_dim; iy++)
for (int iz = 0; iz <= nfft_dim / 2; iz++)
{
if (ix + iy + iz + local_0_start_fourier_return < 100)
{
int index = ix * N0_fourier_grid * (N0_fourier_grid / 2 + 1) + iy * (N0_fourier_grid / 2 + 1) + iz;
fprintf(fp, "%6lu%6d%6d %14.8lf %14.8lf %14.8lf \n", ix + local_0_start_fourier_return, iy, iz,
creal(return_field[index]), cimag(return_field[index]), sqrt(ptr_mode_weightings[index]));
// ptr_mode_weightings[index]);
};
};
fclose(fp);
};
// Transpose output field
{
}else{
FFTW_PLAN transpose_output_plan;
unsigned flags = FFTW_ESTIMATE;
fp = fopen(filename,"w");
for (int ix=0; ix<local_n0_fourier_return; ix++)
for (int iy=0; iy < nfft_dim; iy++)
for (int iz=0; iz <= nfft_dim/2; iz++){
if (ix+iy+iz+local_0_start_fourier_return<100){
int index = ix*N0_fourier_grid*(N0_fourier_grid/2+1) + iy*(N0_fourier_grid/2+1) + iz;
fprintf(fp,"%6llu%6d%6d %14.8lf %14.8lf %14.8lf \n",ix+local_0_start_fourier_return,iy,iz,
creal(return_field[index]),cimag(return_field[index]),sqrt(ptr_mode_weightings[index]));
// ptr_mode_weightings[index]);
};
};
fclose(fp);
void *ptr_inter = return_field;
FFTW_REAL *ptr_return_as_real_field;
ptr_return_as_real_field = ptr_inter;
};
// int rank = 2;
// ptrdiff_t size_to_transpose;
ptrdiff_t howmany = N0_fourier_grid + 2;
// ptrdiff_t local_n0, local_0_start;
// ptrdiff_t local_n1, local_1_start;
// Transpose output field
// const ptrdiff_t ndimens[] = {N0_fourier_grid, N0_fourier_grid};
{
// size_to_transpose = FFTW_MPI_LOCAL_SIZE_MANY_TRANSPOSED(rank, ndimens, howmany,
// FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
// MPI_COMM_WORLD,&local_n0, &local_0_start,&local_n1, &local_1_start);
FFTW_PLAN transpose_output_plan;
unsigned flags = FFTW_ESTIMATE;
// printf("size_to_transpose = %td\n",size_to_transpose);
void *ptr_inter = return_field;
transpose_output_plan = FFTW_MPI_PLAN_MANY_TRANSPOSE(N0_fourier_grid, N0_fourier_grid,
howmany, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
ptr_return_as_real_field, ptr_return_as_real_field,
MPI_COMM_WORLD, flags);
FFTW_REAL *ptr_return_as_real_field;
ptr_return_as_real_field = ptr_inter;
//printf("Transpose plan completed.\n");
int rank = 2;
FFTW_MPI_EXECUTE_R2R(transpose_output_plan, ptr_return_as_real_field, ptr_return_as_real_field);
ptrdiff_t size_to_transpose;
ptrdiff_t howmany = N0_fourier_grid + 2;
ptrdiff_t local_n0, local_0_start;
ptrdiff_t local_n1, local_1_start;
//printf("Transpose completed.\n");
};
const ptrdiff_t ndimens[] = {N0_fourier_grid, N0_fourier_grid};
// Free all memory assigned by FFTW_MALLOC
FFTW_FREE(mode_weightings);
FFTW_FREE(sph_bessel_coeff);
FFTW_FREE(panphasia_coefficients);
if (!SHARED_FOUR_PAN_SPACE)
FFTW_FREE(fourier_grids);
size_to_transpose = FFTW_MPI_LOCAL_SIZE_MANY_TRANSPOSED(rank, ndimens, howmany,
FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
MPI_COMM_WORLD,&local_n0, &local_0_start,&local_n1, &local_1_start);
// printf("size_to_transpose = %td\n",size_to_transpose);
transpose_output_plan = FFTW_MPI_PLAN_MANY_TRANSPOSE(N0_fourier_grid, N0_fourier_grid,
howmany, FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK,
ptr_return_as_real_field, ptr_return_as_real_field,
MPI_COMM_WORLD, flags);
//printf("Transpose plan completed.\n");
FFTW_MPI_EXECUTE_R2R(transpose_output_plan,ptr_return_as_real_field,ptr_return_as_real_field);
//printf("Transpose completed.\n");
};
// Free all memory assigned by FFTW_MALLOC
FFTW_FREE(mode_weightings);
FFTW_FREE(sph_bessel_coeff);
FFTW_FREE(panphasia_coefficients);
if (!SHARED_FOUR_PAN_SPACE) FFTW_FREE(fourier_grids);
FFTW_DESTROY_PLAN(output_coeff_forward_plan);
//printf("Reached end of PANPHASIA_compute_kspace_field_\n");
return(0);
};
FFTW_DESTROY_PLAN(output_coeff_forward_plan);
//printf("Reached end of PANPHASIA_compute_kspace_field_\n");
return (0);
}
//==========================================================================================
//==========================================================================================

2
external/panphasia_ho/panphasia_functions.h vendored
View file

@ -66,7 +66,7 @@ int PANPHASIA_compute_coefficients_(size_t *xstart, size_t *ystart, size_t*zstar
void test_moments_();
void test_propogation_of_moments_(int iterations);
void test_cell_moments(char *,size_t, size_t, size_t, size_t, size_t, double *);
void test_cell_moments(const char*,size_t, size_t, size_t, size_t, size_t, double *);
void spherical_bessel_(int *, double *, double *);

1586
external/panphasia_ho/uniform_rand_threefry4x64.c vendored
View file

File diff suppressed because it is too large Load diff

2
src/plugins/random_panphasia.cc
View file

@ -162,6 +162,8 @@ void RNG_panphasia::Run_Panphasia_Highorder(Grid_FFT<real_t> &g)
Grid_FFT<real_t> pan_grid({{N0, N0, N0}}, {{boxlength_, boxlength_, boxlength_}});
_unused(alloc_local);
assert(pan_grid.n_[0] == N0);
assert(pan_grid.n_[1] == N0);
assert(pan_grid.n_[2] == N0);