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

WIP: some more cleanup

Browse source
This commit is contained in:
Oliver Hahn 2023-02-12 22:38:28 -08:00
parent bd12e40e22
commit 244f11a821
5 changed files with 69 additions and 126 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
ics_example.conf
View file

@ -12,8 +12,7 @@ align_top = no
baryons = no
use_2LPT = no
use_LLA = no
periodic_TF = yes
zero_zoom_velocity = no
[cosmology]
Omega_m = 0.276

113
src/convolution_kernel.cc
View file

@ -1,10 +1,9 @@
/*
convolution_kernel.cc - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
a code to generate multi-scale initial conditions for cosmological simulations
Copyright (C) 2010-19 Oliver Hahn
Copyright (C) 2010-23 Oliver Hahn
*/
@ -96,16 +95,23 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
std::complex<double> dcmode(RE(cdata[0]), IM(cdata[0]));
if (!pk->is_ksampled())
#pragma omp parallel
{
#pragma omp parallel for
const size_t veclen = cparam_.nz / 2 + 1;
double *kvec = new double[veclen];
double *Tkvec = new double[veclen];
double *argvec = new double[veclen];
#pragma omp for
for (int i = 0; i < cparam_.nx; ++i)
for (int j = 0; j < cparam_.ny; ++j)
{
for (int k = 0; k < cparam_.nz / 2 + 1; ++k)
{
size_t ii = (size_t)(i * cparam_.ny + j) * (size_t)(cparam_.nz / 2 + 1) + (size_t)k;
double kx, ky, kz;
kx = (double)i;
@ -117,92 +123,42 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
if (ky > cparam_.ny / 2)
ky -= cparam_.ny;
double arg = (kx + ky + kz) * dstag;
std::complex<double> carg(cos(arg), sin(arg));
kvec[k] = sqrt(kx * kx + ky * ky + kz * kz);
argvec[k] = (kx + ky + kz) * dstag;
}
std::complex<double>
ccdata(RE(cdata[ii]), IM(cdata[ii])),
cckernel(RE(ckernel[ii]), IM(ckernel[ii]));
pk->at_k(veclen, kvec, Tkvec);
for (int k = 0; k < cparam_.nz / 2 + 1; ++k)
{
size_t ii = (size_t)(i * cparam_.ny + j) * (size_t)(cparam_.nz / 2 + 1) + (size_t)k;
std::complex<double> carg(cos(argvec[k]), sin(argvec[k]));
std::complex<double> ccdata(RE(cdata[ii]), IM(cdata[ii]));
if( fix ){
ccdata = ccdata / std::abs(ccdata);
ccdata = ccdata / std::abs(ccdata) / fftnormp;
}
if( flip ){
ccdata = -ccdata;
}
ccdata = ccdata * cckernel * fftnorm * carg;
ccdata = ccdata * Tkvec[k] * fftnorm * carg;
RE(cdata[ii]) = ccdata.real();
IM(cdata[ii]) = ccdata.imag();
}
}
delete[] kvec;
delete[] Tkvec;
delete[] argvec;
}
else
{
#pragma omp parallel
{
// we now set the correct DC mode below...
RE(cdata[0]) = 0.0;
IM(cdata[0]) = 0.0;
const size_t veclen = cparam_.nz / 2 + 1;
double *kvec = new double[veclen];
double *Tkvec = new double[veclen];
double *argvec = new double[veclen];
#pragma omp for
for (int i = 0; i < cparam_.nx; ++i)
for (int j = 0; j < cparam_.ny; ++j)
{
for (int k = 0; k < cparam_.nz / 2 + 1; ++k)
{
double kx, ky, kz;
kx = (double)i;
ky = (double)j;
kz = (double)k;
if (kx > cparam_.nx / 2)
kx -= cparam_.nx;
if (ky > cparam_.ny / 2)
ky -= cparam_.ny;
kvec[k] = sqrt(kx * kx + ky * ky + kz * kz);
argvec[k] = (kx + ky + kz) * dstag;
}
pk->at_k(veclen, kvec, Tkvec);
for (int k = 0; k < cparam_.nz / 2 + 1; ++k)
{
size_t ii = (size_t)(i * cparam_.ny + j) * (size_t)(cparam_.nz / 2 + 1) + (size_t)k;
std::complex<double> carg(cos(argvec[k]), sin(argvec[k]));
std::complex<double> ccdata(RE(cdata[ii]), IM(cdata[ii]));
if( fix ){
ccdata = ccdata / std::abs(ccdata) / fftnormp;
}
if( flip ){
ccdata = -ccdata;
}
ccdata = ccdata * Tkvec[k] * fftnorm * carg;
RE(cdata[ii]) = ccdata.real();
IM(cdata[ii]) = ccdata.imag();
}
}
delete[] kvec;
delete[] Tkvec;
delete[] argvec;
}
// we now set the correct DC mode below...
RE(cdata[0]) = 0.0;
IM(cdata[0]) = 0.0;
}
LOGUSER("Performing backward FFT...");
@ -229,7 +185,6 @@ void perform(kernel *pk, void *pd, bool shift, bool fix, bool flip)
#endif
// set the DC mode here to avoid a possible truncation error in single precision
if (pk->is_ksampled())
{
size_t nelem = (size_t)cparam_.nx * (size_t)cparam_.ny * (size_t)cparam_.nz;
real_t mean = dcmode.real() * fftnorm / (real_t)nelem;

48
src/densities.cc
View file

@ -114,8 +114,13 @@ void fft_coarsen(m1 &v, m2 &V)
val_fine *= val_phas * fftnorm / 8.0;
RE(ccoarse[qc]) = val_fine.real();
IM(ccoarse[qc]) = val_fine.imag();
if( i!=(int)nxF/2 && j!=(int)nyF/2 && k!=(int)nzF/2 ){
RE(ccoarse[qc]) = val_fine.real();
IM(ccoarse[qc]) = val_fine.imag();
}else{
RE(ccoarse[qc]) = 0.0;//val_fine.real();
IM(ccoarse[qc]) = 0.0;//val_fine.imag();
}
}
delete[] rfine;
@ -449,27 +454,14 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
convolution::kernel_creator *the_kernel_creator;
if (kspaceTF)
{
std::cout << " - Using k-space transfer function kernel.\n";
LOGUSER("Using k-space transfer function kernel.");
std::cout << " - Using k-space transfer function kernel.\n";
LOGUSER("Using k-space transfer function kernel.");
#ifdef SINGLE_PRECISION
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_k_float"];
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_k_float"];
#else
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_k_double"];
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_k_double"];
#endif
}
else
{
std::cout << " - Using real-space transfer function kernel.\n";
LOGUSER("Using real-space transfer function kernel.");
#ifdef SINGLE_PRECISION
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_real_float"];
#else
the_kernel_creator = convolution::get_kernel_map()["tf_kernel_real_double"];
#endif
}
convolution::kernel *the_tf_kernel = the_kernel_creator->create(cf, ptf, refh, type);
@ -502,21 +494,13 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
refh.size(levelmin + i, 1), refh.size(levelmin + i, 2));
if( refh.get_margin() > 0 ){
fine = new PaddedDensitySubGrid<real_t>(refh.offset(levelmin + i, 0),
refh.offset(levelmin + i, 1),
refh.offset(levelmin + i, 2),
refh.size(levelmin + i, 0),
refh.size(levelmin + i, 1),
refh.size(levelmin + i, 2),
refh.get_margin(), refh.get_margin(), refh.get_margin() );
fine = new PaddedDensitySubGrid<real_t>( refh.offset(levelmin + i, 0), refh.offset(levelmin + i, 1), refh.offset(levelmin + i, 2),
refh.size(levelmin + i, 0), refh.size(levelmin + i, 1), refh.size(levelmin + i, 2),
refh.get_margin(), refh.get_margin(), refh.get_margin() );
LOGUSER(" margin = %d",refh.get_margin());
}else{
fine = new PaddedDensitySubGrid<real_t>(refh.offset(levelmin + i, 0),
refh.offset(levelmin + i, 1),
refh.offset(levelmin + i, 2),
refh.size(levelmin + i, 0),
refh.size(levelmin + i, 1),
refh.size(levelmin + i, 2));
fine = new PaddedDensitySubGrid<real_t>( refh.offset(levelmin + i, 0), refh.offset(levelmin + i, 1), refh.offset(levelmin + i, 2),
refh.size(levelmin + i, 0), refh.size(levelmin + i, 1), refh.size(levelmin + i, 2));
LOGUSER(" margin = %d",refh.size(levelmin + i, 0)/2);
}
/////////////////////////////////////////////////////////////////////////

21
src/main.cc
View file

@ -421,7 +421,8 @@ int main(int argc, const char *argv[])
bool
do_baryons = cf.getValue<bool>("setup", "baryons"),
do_2LPT = cf.getValueSafe<bool>("setup", "use_2LPT", false),
do_LLA = cf.getValueSafe<bool>("setup", "use_LLA", false);
do_LLA = cf.getValueSafe<bool>("setup", "use_LLA", false),
do_counter_mode = cf.getValueSafe<bool>("setup", "zero_zoom_velocity", false);
transfer_function_plugin *the_transfer_function_plugin = select_transfer_function_plugin(cf);
@ -460,7 +461,7 @@ int main(int argc, const char *argv[])
}
the_region_generator = select_region_generator_plugin(cf);
//------------------------------------------------------------------------------
//... determine run parameters
//------------------------------------------------------------------------------
@ -619,7 +620,7 @@ int main(int argc, const char *argv[])
//... compute counter-mode to minimize advection errors
counter_mode_amp[icoord] = compute_finest_mean(data_forIO);
add_constant_value( data_forIO, -counter_mode_amp[icoord] );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord] );
LOGUSER("Writing CDM displacements");
the_output_plugin->write_dm_position(icoord, data_forIO);
@ -756,7 +757,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity-mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
LOGUSER("Writing CDM velocities");
the_output_plugin->write_dm_velocity(icoord, data_forIO);
@ -826,7 +827,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
LOGUSER("Writing CDM velocities");
the_output_plugin->write_dm_velocity(icoord, data_forIO);
@ -885,7 +886,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]*cosmo.vfact );
LOGUSER("Writing baryon velocities");
the_output_plugin->write_gas_velocity(icoord, data_forIO);
@ -1008,7 +1009,7 @@ int main(int argc, const char *argv[])
//... compute counter-mode to minimize advection errors
counter_mode_amp[icoord] = compute_finest_mean(data_forIO);
add_constant_value( data_forIO, -counter_mode_amp[icoord] );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord] );
LOGUSER("Writing CDM velocities");
the_output_plugin->write_dm_velocity(icoord, data_forIO);
@ -1105,7 +1106,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter velocity mode
add_constant_value( data_forIO, -counter_mode_amp[icoord] );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord] );
LOGUSER("Writing baryon velocities");
the_output_plugin->write_gas_velocity(icoord, data_forIO);
@ -1210,7 +1211,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
LOGUSER("Writing CDM displacements");
the_output_plugin->write_dm_position(icoord, data_forIO);
@ -1327,7 +1328,7 @@ int main(int argc, const char *argv[])
coarsen_density(rh_Poisson, data_forIO, false);
// add counter mode
add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
if( do_counter_mode ) add_constant_value( data_forIO, -counter_mode_amp[icoord]/cosmo.vfact );
LOGUSER("Writing baryon displacements");

10
src/poisson.cc
View file

@ -724,12 +724,16 @@ double fft_poisson_plugin::gradient(int dir, grid_hierarchy &u, grid_hierarchy &
}
#endif
/*double ktot = sqrt(ii*ii+jj*jj+k*k);
if( ktot >= nx/2 )//dir == 0 && i==nx/2 || dir == 1 && j==ny/2 || dir == 2 && k==nz/2 )
if( (dir == 0 && i==nx/2) || (dir == 1 && j==ny/2) || (dir == 2 && k==nz/2) )
{
#ifdef FFTW3
cdata[idx][0] = 0.0;
cdata[idx][1] = 0.0;
#else
cdata[idx].re = 0.0;
cdata[idx].im = 0.0;
}*/
#endif
}
}
RE(cdata[0]) = 0.0;