mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
commit
f65f60076f
22 changed files with 2738 additions and 1811 deletions
1
.gitignore
vendored
1
.gitignore
vendored
|
@ -53,3 +53,4 @@ src/CMakeCache.txt
|
|||
src/fastLPT
|
||||
src/input_powerspec.txt
|
||||
src/Makefile
|
||||
.DS_Store
|
||||
|
|
|
@ -1,18 +1,28 @@
|
|||
cmake_minimum_required(VERSION 3.9)
|
||||
set(PRGNAME fastLPT)
|
||||
project(fastLPT)
|
||||
set(PRGNAME monofonIC)
|
||||
project(monofonIC)
|
||||
|
||||
# include class submodule
|
||||
include(${CMAKE_CURRENT_SOURCE_DIR}/external/class.cmake)
|
||||
|
||||
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O0 -march=native -Wall -fno-omit-frame-pointer -g -fsanitize=undefined")
|
||||
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -march=native -Wall -pedantic")
|
||||
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O0 -march=native -Wall -fno-omit-frame-pointer -g -fsanitize=address")
|
||||
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -march=native -Wall -pedantic")
|
||||
find_package(PkgConfig REQUIRED)
|
||||
|
||||
set(CMAKE_MODULE_PATH
|
||||
"${CMAKE_MODULE_PATH};${PROJECT_SOURCE_DIR}")
|
||||
|
||||
|
||||
########################################################################################################################
|
||||
# Add a custom command that produces version.cc, plus
|
||||
# a dummy output that's not actually produced, in order
|
||||
# to force version.cmake to always be re-run before the build
|
||||
ADD_CUSTOM_COMMAND(
|
||||
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/version.cc
|
||||
${CMAKE_CURRENT_BINARY_DIR}/_version.cc
|
||||
COMMAND ${CMAKE_COMMAND} -P
|
||||
${CMAKE_CURRENT_SOURCE_DIR}/version.cmake)
|
||||
|
||||
########################################################################################################################
|
||||
# OpenMP
|
||||
find_package(OpenMP REQUIRED)
|
||||
|
@ -53,6 +63,8 @@ include_directories(${PROJECT_SOURCE_DIR}/include)
|
|||
file( GLOB SOURCES
|
||||
${PROJECT_SOURCE_DIR}/src/*.cc
|
||||
)
|
||||
# add the auto generated version file
|
||||
list (APPEND SOURCES "${CMAKE_CURRENT_BINARY_DIR}/version.cc")
|
||||
|
||||
# PLUGINS
|
||||
# get all the *.cc files in the plugin subfolder
|
||||
|
@ -82,6 +94,13 @@ if(FFTW_THREADS_FOUND)
|
|||
target_compile_options(${PRGNAME} PRIVATE "-DUSE_FFTW_THREADS")
|
||||
endif(FFTW_THREADS_FOUND)
|
||||
|
||||
if(HDF5_FOUND)
|
||||
# target_link_libraries(${PRGNAME} ${HDF5_C_LIBRARY_DIRS})
|
||||
target_link_libraries(${PRGNAME} ${HDF5_LIBRARIES})
|
||||
target_include_directories(${PRGNAME} PRIVATE ${HDF5_INCLUDE_DIRS})
|
||||
target_compile_options(${PRGNAME} PRIVATE "-DUSE_HDF5")
|
||||
endif(HDF5_FOUND)
|
||||
|
||||
target_link_libraries(${PRGNAME} ${FFTW_LIBRARIES})
|
||||
target_include_directories(${PRGNAME} PRIVATE ${FFTW_INCLUDE_DIR})
|
||||
|
||||
|
|
25
example.conf
25
example.conf
|
@ -1,29 +1,42 @@
|
|||
[setup]
|
||||
GridRes = 64
|
||||
BoxLength = 300
|
||||
BoxLength = 100
|
||||
zstart = 49.0
|
||||
LPTorder = 3
|
||||
LPTorder = 2
|
||||
SymplecticPT = no
|
||||
DoFixing = no
|
||||
BCClattice = no
|
||||
|
||||
[execution]
|
||||
NumThreads = 4
|
||||
|
||||
[output]
|
||||
fname_hdf5 = output.hdf5
|
||||
fname_hdf5 = output_sch.hdf5
|
||||
fbase_analysis = output
|
||||
format = gadget2
|
||||
filename = ics_gadget.dat
|
||||
#format = gadget2
|
||||
#filename = ics_gadget.dat
|
||||
#format = generic
|
||||
#filename = debug.hdf5
|
||||
#generic_out_eulerian = yes
|
||||
|
||||
format = grafic2
|
||||
filename = ics_ramses
|
||||
grafic_use_SPT = yes
|
||||
|
||||
[random]
|
||||
generator = NGENIC
|
||||
seed = 9001
|
||||
|
||||
[cosmology]
|
||||
transfer = CLASS #eisenstein
|
||||
#transfer = CLASS #eisenstein
|
||||
transfer = eisenstein
|
||||
Omega_m = 0.302
|
||||
Omega_b = 0.045
|
||||
Omega_L = 0.698
|
||||
H0 = 70.3
|
||||
sigma_8 = 0.811
|
||||
nspec = 0.961
|
||||
|
||||
[sch]
|
||||
hbar = 5.0 #0.5
|
||||
dt = 1.0
|
||||
|
|
|
@ -3,6 +3,7 @@
|
|||
#include <logger.hh>
|
||||
|
||||
#include <complex>
|
||||
#include <map>
|
||||
|
||||
#if defined(USE_MPI)
|
||||
#include <mpi.h>
|
||||
|
@ -21,6 +22,10 @@ using complex_t = fftw_complex;
|
|||
#define FFTW_PREFIX fftw
|
||||
#endif
|
||||
|
||||
enum class fluid_component { density, vx, vy, vz, dx, dy, dz };
|
||||
enum class cosmo_species { dm, baryon, neutrino };
|
||||
extern std::map<cosmo_species,std::string> cosmo_species_name;
|
||||
|
||||
using ccomplex_t = std::complex<real_t>;
|
||||
|
||||
#define FFTW_GEN_NAME_PRIM(a, b) a##_##b
|
||||
|
@ -102,6 +107,13 @@ MPI_Datatype GetMPIDatatype( void )
|
|||
#endif
|
||||
#endif
|
||||
|
||||
inline void multitask_sync_barrier( void )
|
||||
{
|
||||
#if defined(USE_MPI)
|
||||
MPI_Barrier( MPI_COMM_WORLD );
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
|
||||
namespace CONFIG
|
||||
|
@ -114,3 +126,13 @@ extern bool MPI_threads_ok;
|
|||
extern bool FFTW_threads_ok;
|
||||
extern int num_threads;
|
||||
} // namespace CONFIG
|
||||
|
||||
|
||||
// These variables are autogenerated and compiled
|
||||
// into the library by the version.cmake script
|
||||
extern "C"
|
||||
{
|
||||
extern const char* GIT_TAG;
|
||||
extern const char* GIT_REV;
|
||||
extern const char* GIT_BRANCH;
|
||||
}
|
|
@ -64,20 +64,26 @@ public:
|
|||
}
|
||||
|
||||
const Grid_FFT<data_t>* get_grid( size_t ilevel ) const { return this; }
|
||||
bool is_in_mask( size_t ilevel, size_t i, size_t j, size_t k ) const { return true; }
|
||||
bool is_refined( size_t ilevel, size_t i, size_t j, size_t k ) const { return false; }
|
||||
size_t levelmin() const {return 7;}
|
||||
size_t levelmax() const {return 7;}
|
||||
|
||||
void Setup();
|
||||
|
||||
//! return the (local) size of dimension i
|
||||
size_t size(size_t i) const { return sizes_[i]; }
|
||||
|
||||
//! return the (global) size of dimension i
|
||||
size_t global_size(size_t i) const { return n_[i]; }
|
||||
|
||||
|
||||
//! return locally stored number of elements of field
|
||||
size_t local_size( void ) const { return local_0_size_ * n_[1] * n_[2]; }
|
||||
|
||||
//! return a bounding box of the global extent of the field
|
||||
const bounding_box<size_t>& get_global_range( void ) const
|
||||
{
|
||||
return global_range_;
|
||||
}
|
||||
|
||||
//! set all field elements to zero
|
||||
void zero()
|
||||
{
|
||||
#pragma omp parallel for
|
||||
|
@ -85,6 +91,24 @@ public:
|
|||
data_[i] = 0.0;
|
||||
}
|
||||
|
||||
void copy_from( const Grid_FFT<data_t>& g ){
|
||||
// make sure the two fields are in the same space
|
||||
if( g.space_ != this->space_ ){
|
||||
if( this->space_ == kspace_id ) this->FourierTransformBackward(false);
|
||||
else this->FourierTransformForward(false);
|
||||
}
|
||||
|
||||
// make sure the two fields have the same dimensions
|
||||
assert( this->n_[0] == g.n_[0] );
|
||||
assert( this->n_[1] == g.n_[1] );
|
||||
assert( this->n_[2] == g.n_[2] );
|
||||
|
||||
// now we can copy all the data over
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < ntot_; ++i)
|
||||
data_[i] = g.data_[i];
|
||||
}
|
||||
|
||||
data_t& operator[]( size_t i ){
|
||||
return data_[i];
|
||||
}
|
||||
|
@ -135,6 +159,30 @@ public:
|
|||
return rr;
|
||||
}
|
||||
|
||||
template <typename ft>
|
||||
vec3<ft> get_unit_r(const size_t i, const size_t j, const size_t k) const
|
||||
{
|
||||
vec3<ft> rr;
|
||||
|
||||
rr[0] = real_t(i + local_0_start_) / real_t(n_[0]);
|
||||
rr[1] = real_t(j) / real_t(n_[1]);
|
||||
rr[2] = real_t(k) / real_t(n_[2]);
|
||||
|
||||
return rr;
|
||||
}
|
||||
|
||||
template <typename ft>
|
||||
vec3<ft> get_unit_r_staggered(const size_t i, const size_t j, const size_t k) const
|
||||
{
|
||||
vec3<ft> rr;
|
||||
|
||||
rr[0] = (real_t(i + local_0_start_)+0.5) / real_t(n_[0]);
|
||||
rr[1] = (real_t(j)+0.5) / real_t(n_[1]);
|
||||
rr[2] = (real_t(k)+0.5) / real_t(n_[2]);
|
||||
|
||||
return rr;
|
||||
}
|
||||
|
||||
void cell_pos( int ilevel, size_t i, size_t j, size_t k, double* x ) const {
|
||||
x[0] = double(i+local_0_start_)/size(0);
|
||||
x[1] = double(j)/size(1);
|
||||
|
@ -145,6 +193,14 @@ public:
|
|||
return n_[0]*n_[1]*n_[2];
|
||||
}
|
||||
|
||||
real_t get_dx( int idim ) const{
|
||||
return dx_[idim];
|
||||
}
|
||||
|
||||
const std::array<real_t, 3>& get_dx( void ) const{
|
||||
return dx_;
|
||||
}
|
||||
|
||||
template <typename ft>
|
||||
vec3<ft> get_k(const size_t i, const size_t j, const size_t k) const
|
||||
{
|
||||
|
@ -191,6 +247,16 @@ public:
|
|||
return *this;
|
||||
}
|
||||
|
||||
Grid_FFT<data_t>& apply_negative_Laplacian( void ){
|
||||
this->FourierTransformForward();
|
||||
this->apply_function_k_dep([&](auto x, auto k) {
|
||||
real_t kmod2 = k.norm_squared();
|
||||
return x*kmod2;
|
||||
});
|
||||
this->zero_DC_mode();
|
||||
return *this;
|
||||
}
|
||||
|
||||
Grid_FFT<data_t>& apply_InverseLaplacian( void ){
|
||||
this->FourierTransformForward();
|
||||
this->apply_function_k_dep([&](auto x, auto k) {
|
||||
|
@ -259,7 +325,9 @@ public:
|
|||
|
||||
double std(void)
|
||||
{
|
||||
real_t sum1{0.0}, sum2{0.0};
|
||||
double sum1{0.0}, sum2{0.0};
|
||||
size_t count{0};
|
||||
|
||||
#pragma omp parallel for reduction(+ : sum1, sum2)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
|
@ -273,15 +341,39 @@ public:
|
|||
}
|
||||
}
|
||||
}
|
||||
count = sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
sum2 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
#ifdef USE_MPI
|
||||
double globsum1{0.0}, globsum2{0.0};
|
||||
size_t globcount{0};
|
||||
|
||||
MPI_Allreduce( reinterpret_cast<const void*>(&sum1),
|
||||
reinterpret_cast<void*>(&globsum1),
|
||||
1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
|
||||
|
||||
MPI_Allreduce( reinterpret_cast<const void*>(&sum2),
|
||||
reinterpret_cast<void*>(&globsum2),
|
||||
1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
|
||||
|
||||
MPI_Allreduce( reinterpret_cast<const void*>(&count),
|
||||
reinterpret_cast<void*>(&globcount),
|
||||
1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD );
|
||||
|
||||
sum1 = globsum1;
|
||||
sum2 = globsum2;
|
||||
count = globcount;
|
||||
#endif
|
||||
sum1 /= count;
|
||||
sum2 /= count;
|
||||
|
||||
return std::sqrt(sum2 - sum1 * sum1);
|
||||
}
|
||||
|
||||
double mean(void)
|
||||
{
|
||||
real_t sum1{0.0};
|
||||
double sum1{0.0};
|
||||
size_t count{0};
|
||||
|
||||
#pragma omp parallel for reduction(+ : sum1)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
|
@ -294,8 +386,25 @@ public:
|
|||
}
|
||||
}
|
||||
}
|
||||
count = sizes_[0] * sizes_[1] * sizes_[2];
|
||||
|
||||
sum1 /= sizes_[0] * sizes_[1] * sizes_[2];
|
||||
#ifdef USE_MPI
|
||||
double globsum1{0.0};
|
||||
size_t globcount{0};
|
||||
|
||||
MPI_Allreduce( reinterpret_cast<const void*>(&sum1),
|
||||
reinterpret_cast<void*>(&globsum1),
|
||||
1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
|
||||
|
||||
MPI_Allreduce( reinterpret_cast<const void*>(&count),
|
||||
reinterpret_cast<void*>(&globcount),
|
||||
1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD );
|
||||
|
||||
sum1 = globsum1;
|
||||
count = globcount;
|
||||
#endif
|
||||
|
||||
sum1 /= count;
|
||||
|
||||
return sum1;
|
||||
}
|
||||
|
@ -321,6 +430,7 @@ public:
|
|||
}
|
||||
}
|
||||
|
||||
|
||||
template <typename functional, typename grid1_t, typename grid2_t>
|
||||
void assign_function_of_grids_r(const functional &f, const grid1_t &g1, const grid2_t &g2)
|
||||
{
|
||||
|
@ -373,6 +483,27 @@ public:
|
|||
}
|
||||
}
|
||||
|
||||
template <typename functional, typename grid_t>
|
||||
void assign_function_of_grids_k(const functional &f, const grid_t &g)
|
||||
{
|
||||
assert(g.size(0) == size(0) && g.size(1) == size(1)); // && g.size(2) == size(2) );
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
{
|
||||
for (size_t k = 0; k < sizes_[2]; ++k)
|
||||
{
|
||||
auto &elem = this->kelem(i, j, k);
|
||||
const auto &elemg = g.kelem(i, j, k);
|
||||
|
||||
elem = f(elemg);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename functional>
|
||||
void apply_function_k_dep(const functional &f)
|
||||
{
|
||||
|
@ -415,7 +546,7 @@ public:
|
|||
|
||||
void FillRandomReal(unsigned long int seed = 123456ul);
|
||||
|
||||
void Write_to_HDF5(std::string fname, std::string datasetname);
|
||||
void Write_to_HDF5(std::string fname, std::string datasetname) const;
|
||||
|
||||
void Write_PowerSpectrum( std::string ofname );
|
||||
|
||||
|
@ -423,6 +554,15 @@ public:
|
|||
|
||||
void Write_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
|
||||
|
||||
void stagger_field( void ){
|
||||
FourierTransformForward();
|
||||
apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
|
||||
real_t shift = k[0]*get_dx()[0] + k[1]*get_dx()[1] + k[2]*get_dx()[2];
|
||||
return x * std::exp(ccomplex_t(0.0,0.5*shift));
|
||||
});
|
||||
FourierTransformBackward();
|
||||
}
|
||||
|
||||
void zero_DC_mode(void)
|
||||
{
|
||||
if( space_ == kspace_id ){
|
||||
|
@ -432,7 +572,7 @@ public:
|
|||
cdata_[0] = (data_t)0.0;
|
||||
}else{
|
||||
data_t sum = 0.0;
|
||||
// #pragma omp parallel for reduction(+:sum)
|
||||
// #pragma omp parallel for reduction(+:sum)
|
||||
for (size_t i = 0; i < sizes_[0]; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < sizes_[1]; ++j)
|
||||
|
@ -464,4 +604,5 @@ public:
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
|
|
@ -1,129 +1,68 @@
|
|||
/*
|
||||
|
||||
output.hh - This file is part of MUSIC -
|
||||
output_plugin.hh - This file is part of MUSIC2 -
|
||||
a code to generate multi-scale initial conditions
|
||||
for cosmological simulations
|
||||
|
||||
Copyright (C) 2010 Oliver Hahn
|
||||
Copyright (C) 2019 Oliver Hahn
|
||||
|
||||
*/
|
||||
|
||||
#ifndef __OUTPUT_HH
|
||||
#define __OUTPUT_HH
|
||||
#pragma once
|
||||
|
||||
#include <string>
|
||||
#include <cstring>
|
||||
#include <map>
|
||||
|
||||
#include "general.hh"
|
||||
#include "grid_fft.hh"
|
||||
#include "config_file.hh"
|
||||
#include <particle_container.hh>
|
||||
#include <general.hh>
|
||||
#include <grid_fft.hh>
|
||||
#include <config_file.hh>
|
||||
|
||||
enum class output_type {particles,field_lagrangian,field_eulerian};
|
||||
|
||||
/*!
|
||||
* @class output_plugin
|
||||
* @brief abstract base class for output plug-ins
|
||||
*
|
||||
* This class provides the abstract base class for all output plug-ins.
|
||||
* All output plug-ins need to derive from it and implement the purely
|
||||
* virtual member functions.
|
||||
*/
|
||||
class output_plugin
|
||||
{
|
||||
|
||||
public:
|
||||
using grid_hierarchy = Grid_FFT<real_t>;
|
||||
|
||||
protected:
|
||||
|
||||
//! reference to the ConfigFile object that holds all configuration options
|
||||
ConfigFile& cf_;
|
||||
ConfigFile &cf_;
|
||||
|
||||
//! output file or directory name
|
||||
std::string fname_;
|
||||
|
||||
//! minimum refinement level
|
||||
// unsigned levelmin_;
|
||||
|
||||
//! maximum refinement level
|
||||
// unsigned levelmax_;
|
||||
|
||||
std::vector<unsigned>
|
||||
offx_, //!< vector describing the x-offset of each level
|
||||
offy_, //!< vector describing the y-offset of each level
|
||||
offz_, //!< vector describing the z-offset of each level
|
||||
sizex_, //!< vector describing the extent in x of each level
|
||||
sizey_, //!< vector describing the extent in y of each level
|
||||
sizez_; //!< vector describing the extent in z of each level
|
||||
|
||||
//! quick access function to query properties of the refinement grid from the configuration options
|
||||
/*! @param name name of the config property
|
||||
* @param icomp component index (0=x, 1=y, 2=z)
|
||||
* @param oit output iterator (e.g. std::back_inserter for vectors)
|
||||
*/
|
||||
template< typename output_iterator >
|
||||
void query_grid_prop( std::string name, int icomp, output_iterator oit )
|
||||
{
|
||||
char str[128];
|
||||
//for( unsigned i=levelmin_; i<=levelmax_; ++i )
|
||||
unsigned i=0;
|
||||
{
|
||||
sprintf( str, "%s(%u,%d)", name.c_str(), i, icomp );
|
||||
*oit = 0; //cf_.GetValue<unsigned>( "setup", str );
|
||||
++oit;
|
||||
}
|
||||
}
|
||||
|
||||
//! name of the output interface
|
||||
std::string interface_name_;
|
||||
public:
|
||||
|
||||
//! constructor
|
||||
explicit output_plugin( ConfigFile& cf )
|
||||
: cf_(cf)
|
||||
output_plugin(ConfigFile &cf, std::string interface_name )
|
||||
: cf_(cf), interface_name_(interface_name)
|
||||
{
|
||||
fname_ = cf.GetValue<std::string>("output","filename");
|
||||
// levelmin_ = cf.GetValue<unsigned>( "setup", "levelmin" );
|
||||
// levelmax_ = cf.GetValue<unsigned>( "setup", "levelmax" );
|
||||
|
||||
query_grid_prop( "offset", 0, std::back_inserter(offx_) );
|
||||
query_grid_prop( "offset", 1, std::back_inserter(offy_) );
|
||||
query_grid_prop( "offset", 2, std::back_inserter(offz_) );
|
||||
|
||||
query_grid_prop( "size", 0, std::back_inserter(sizex_) );
|
||||
query_grid_prop( "size", 1, std::back_inserter(sizey_) );
|
||||
query_grid_prop( "size", 2, std::back_inserter(sizez_) );
|
||||
fname_ = cf_.GetValue<std::string>("output", "filename");
|
||||
}
|
||||
|
||||
//! destructor
|
||||
virtual ~output_plugin()
|
||||
{ }
|
||||
//! virtual destructor
|
||||
virtual ~output_plugin(){}
|
||||
|
||||
//! purely virtual prototype to write the masses for each dark matter particle
|
||||
virtual void write_dm_mass( const grid_hierarchy& gh ) = 0;
|
||||
//! routine to write particle data for a species
|
||||
virtual void write_particle_data(const particle_container &pc, const cosmo_species &s ) {};
|
||||
|
||||
//! purely virtual prototype to write the dark matter density field
|
||||
virtual void write_dm_density( const grid_hierarchy& gh ) = 0;
|
||||
//! routine to write gridded fluid component data for a species
|
||||
virtual void write_grid_data(const Grid_FFT<real_t> &g, const cosmo_species &s, const fluid_component &c ) {};
|
||||
|
||||
//! purely virtual prototype to write the dark matter gravitational potential (from which displacements are computed in 1LPT)
|
||||
virtual void write_dm_potential( const grid_hierarchy& gh ) = 0;
|
||||
//! routine to query whether species is written as grid data
|
||||
virtual output_type write_species_as ( const cosmo_species &s ) const = 0;
|
||||
|
||||
//! purely virtual prototype to write dark matter particle velocities
|
||||
virtual void write_dm_velocity( int coord, const grid_hierarchy& gh ) = 0;
|
||||
//! routine to query whether species is written as grid data
|
||||
// virtual bool write_species_as_grid( const cosmo_species &s ) = 0;
|
||||
|
||||
//! purely virtual prototype to write dark matter particle positions
|
||||
virtual void write_dm_position( int coord, const grid_hierarchy& gh ) = 0;
|
||||
//! routine to query whether species is written as particle data
|
||||
// virtual bool write_species_as_particles( const cosmo_species &s ){ return !write_species_as_grid(s); }
|
||||
|
||||
//! purely virtual prototype to write the baryon velocities
|
||||
virtual void write_gas_velocity( int coord, const grid_hierarchy& gh ) = 0;
|
||||
//! routine to return a multiplicative factor that contains the desired position units for the output
|
||||
virtual real_t position_unit() const = 0;
|
||||
|
||||
//! purely virtual prototype to write the baryon coordinates
|
||||
virtual void write_gas_position( int coord, const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the baryon density field
|
||||
virtual void write_gas_density( const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype to write the baryon gravitational potential (from which displacements are computed in 1LPT)
|
||||
virtual void write_gas_potential( const grid_hierarchy& gh ) = 0;
|
||||
|
||||
//! purely virtual prototype for all things to be done at the very end
|
||||
virtual void finalize( void ) = 0;
|
||||
//! routine to return a multiplicative factor that contains the desired velocity units for the output
|
||||
virtual real_t velocity_unit() const = 0;
|
||||
};
|
||||
|
||||
/*!
|
||||
|
@ -132,14 +71,14 @@ public:
|
|||
struct output_plugin_creator
|
||||
{
|
||||
//! create an instance of a plug-in
|
||||
virtual std::unique_ptr<output_plugin> create( ConfigFile& cf ) const = 0;
|
||||
virtual std::unique_ptr<output_plugin> create(ConfigFile &cf) const = 0;
|
||||
|
||||
//! destroy an instance of a plug-in
|
||||
virtual ~output_plugin_creator() { }
|
||||
virtual ~output_plugin_creator() {}
|
||||
};
|
||||
|
||||
//! maps the name of a plug-in to a pointer of the factory pattern
|
||||
std::map< std::string, output_plugin_creator *>& get_output_plugin_map();
|
||||
std::map<std::string, output_plugin_creator *> &get_output_plugin_map();
|
||||
|
||||
//! print a list of all registered output plug-ins
|
||||
void print_output_plugins();
|
||||
|
@ -147,23 +86,22 @@ void print_output_plugins();
|
|||
/*!
|
||||
* @brief concrete factory pattern for output plug-ins
|
||||
*/
|
||||
template< class Derived >
|
||||
template <class Derived>
|
||||
struct output_plugin_creator_concrete : public output_plugin_creator
|
||||
{
|
||||
//! register the plug-in by its name
|
||||
output_plugin_creator_concrete( const std::string& plugin_name )
|
||||
output_plugin_creator_concrete(const std::string &plugin_name)
|
||||
{
|
||||
get_output_plugin_map()[ plugin_name ] = this;
|
||||
get_output_plugin_map()[plugin_name] = this;
|
||||
}
|
||||
|
||||
//! create an instance of the plug-in
|
||||
std::unique_ptr<output_plugin> create( ConfigFile& cf ) const
|
||||
std::unique_ptr<output_plugin> create(ConfigFile &cf) const
|
||||
{
|
||||
return std::make_unique<Derived>(cf);// Derived( cf );
|
||||
return std::make_unique<Derived>(cf); // Derived( cf );
|
||||
}
|
||||
};
|
||||
|
||||
//! failsafe version to select the output plug-in
|
||||
std::unique_ptr<output_plugin> select_output_plugin( ConfigFile& cf );
|
||||
std::unique_ptr<output_plugin> select_output_plugin(ConfigFile &cf);
|
||||
|
||||
#endif // __OUTPUT_HH
|
||||
|
|
108
include/particle_container.hh
Normal file
108
include/particle_container.hh
Normal file
|
@ -0,0 +1,108 @@
|
|||
#pragma once
|
||||
|
||||
#ifdef USE_MPI
|
||||
#include <mpi.h>
|
||||
#endif
|
||||
|
||||
#include <numeric>
|
||||
#include <vector>
|
||||
#include <general.hh>
|
||||
|
||||
class particle_container
|
||||
{
|
||||
public:
|
||||
std::vector<float> positions_, velocities_;
|
||||
std::vector<int> ids_;
|
||||
|
||||
particle_container()
|
||||
{
|
||||
}
|
||||
|
||||
particle_container(const particle_container &) = delete;
|
||||
|
||||
const void* get_pos_ptr() const{
|
||||
return reinterpret_cast<const void*>( &positions_[0] );
|
||||
}
|
||||
|
||||
const void* get_vel_ptr() const{
|
||||
return reinterpret_cast<const void*>( &velocities_[0] );
|
||||
}
|
||||
|
||||
const void* get_ids_ptr() const{
|
||||
return reinterpret_cast<const void*>( &ids_[0] );
|
||||
}
|
||||
|
||||
void allocate(size_t nump)
|
||||
{
|
||||
positions_.resize(3 * nump);
|
||||
velocities_.resize(3 * nump);
|
||||
ids_.resize(nump);
|
||||
}
|
||||
|
||||
void set_pos(size_t ipart, size_t idim, real_t p)
|
||||
{
|
||||
positions_[3 * ipart + idim] = p;
|
||||
}
|
||||
|
||||
void set_vel(size_t ipart, size_t idim, real_t p)
|
||||
{
|
||||
velocities_[3 * ipart + idim] = p;
|
||||
}
|
||||
|
||||
void set_id(size_t ipart, id_t id)
|
||||
{
|
||||
ids_[ipart] = id;
|
||||
}
|
||||
|
||||
size_t get_local_num_particles(void) const
|
||||
{
|
||||
return ids_.size();
|
||||
}
|
||||
|
||||
size_t get_global_num_particles(void) const
|
||||
{
|
||||
size_t local_nump = ids_.size(), global_nump;
|
||||
#ifdef USE_MPI
|
||||
MPI_Allreduce(reinterpret_cast<void *>(&local_nump), reinterpret_cast<void *>(&global_nump), 1,
|
||||
MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
|
||||
#else
|
||||
global_nump = local_nump;
|
||||
#endif
|
||||
return global_nump;
|
||||
}
|
||||
|
||||
size_t get_local_offset( void ) const
|
||||
{
|
||||
size_t this_offset = 0;
|
||||
|
||||
#ifdef USE_MPI
|
||||
int mpi_size, mpi_rank;
|
||||
MPI_Comm_size( MPI_COMM_WORLD, &mpi_size );
|
||||
MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank );
|
||||
|
||||
std::vector<size_t> nump_p_task(mpi_size,0), off_p_task;
|
||||
size_t num_p_this_task = this->get_local_num_particles();
|
||||
// MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
|
||||
// void *recvbuf, int recvcount, MPI_Datatype recvtype,
|
||||
// MPI_Comm comm)
|
||||
MPI_Allgather( reinterpret_cast<const void*>(&num_p_this_task), 1, MPI_UNSIGNED_LONG_LONG,
|
||||
reinterpret_cast<void*>(&nump_p_task[0]), 1, MPI_UNSIGNED_LONG_LONG, MPI_COMM_WORLD );
|
||||
|
||||
off_p_task.push_back( 0 );
|
||||
std::partial_sum(nump_p_task.begin(), nump_p_task.end(), std::back_inserter(off_p_task) );
|
||||
this_offset = nump_p_task.at(mpi_rank);
|
||||
#endif
|
||||
|
||||
return this_offset;
|
||||
}
|
||||
|
||||
void dump(void)
|
||||
{
|
||||
for (size_t i = 0; i < ids_.size(); ++i)
|
||||
{
|
||||
std::cout << positions_[3 * i + 0] << " " << positions_[3 * i + 1] << " " << positions_[3 * i + 2] << " "
|
||||
<< velocities_[3 * i + 0] << " " << velocities_[3 * i + 1] << " " << velocities_[3 * i + 2] << std::endl;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
|
@ -1,92 +1,139 @@
|
|||
#pragma once
|
||||
|
||||
#include <string>
|
||||
|
||||
#ifdef __APPLE__
|
||||
#include <sys/types.h>
|
||||
#include <sys/sysctl.h>
|
||||
#include <mach/mach.h>
|
||||
#include <mach/vm_statistics.h>
|
||||
#include <mach/mach_types.h>
|
||||
#include <mach/mach_init.h>
|
||||
#include <mach/mach_host.h>
|
||||
#include <unistd.h>
|
||||
#include <sys/types.h>
|
||||
#include <sys/sysctl.h>
|
||||
#include <mach/mach.h>
|
||||
#include <mach/vm_statistics.h>
|
||||
#include <mach/mach_types.h>
|
||||
#include <mach/mach_init.h>
|
||||
#include <mach/mach_host.h>
|
||||
#include <unistd.h>
|
||||
#elif __linux__
|
||||
#include <cstring>
|
||||
#include <cstdio>
|
||||
#include <strings.h>
|
||||
#include <cstring>
|
||||
#include <cstdio>
|
||||
#include <strings.h>
|
||||
#endif
|
||||
|
||||
namespace SystemStat
|
||||
{
|
||||
|
||||
namespace SystemStat {
|
||||
class Cpu
|
||||
{
|
||||
public:
|
||||
Cpu() {}
|
||||
|
||||
class Memory{
|
||||
private:
|
||||
size_t total;
|
||||
size_t avail;
|
||||
size_t used;
|
||||
public:
|
||||
Memory()
|
||||
std::string get_CPUstring() const
|
||||
{
|
||||
#ifdef __APPLE__
|
||||
char buffer[1024];
|
||||
size_t size=sizeof(buffer);
|
||||
if (sysctlbyname("machdep.cpu.brand_string", &buffer, &size, NULL, 0) < 0) {
|
||||
return "";
|
||||
}
|
||||
return std::string(buffer);
|
||||
#elif __linux__
|
||||
std::string str = "";
|
||||
FILE *cpuinfo = fopen("/proc/cpuinfo", "rb");
|
||||
char *arg = 0;
|
||||
size_t size = 0;
|
||||
while(getdelim(&arg, &size, '\n', cpuinfo) != -1)
|
||||
{
|
||||
if( strncmp(arg, "model name", 10) == 0 ){
|
||||
str = std::string(arg+13);
|
||||
break;
|
||||
}
|
||||
}
|
||||
free(arg);
|
||||
fclose(cpuinfo);
|
||||
return str;
|
||||
#endif
|
||||
}
|
||||
};
|
||||
|
||||
class Memory
|
||||
{
|
||||
private:
|
||||
size_t total;
|
||||
size_t avail;
|
||||
size_t used;
|
||||
|
||||
public:
|
||||
Memory()
|
||||
: total(0), avail(0), used(0)
|
||||
{
|
||||
this->get_statistics();
|
||||
}
|
||||
|
||||
size_t get_TotalMem() const { return this->total; }
|
||||
size_t get_AvailMem() const { return this->avail; }
|
||||
size_t get_UsedMem() const { return this->used; }
|
||||
void update() { this->get_statistics(); }
|
||||
|
||||
protected:
|
||||
int get_statistics(void)
|
||||
{
|
||||
#ifdef __APPLE__
|
||||
int64_t pagesize = int64_t(getpagesize());
|
||||
int mib[2] = {CTL_HW, HW_MEMSIZE};
|
||||
size_t length = sizeof(size_t);
|
||||
sysctl(mib, 2, &this->total, &length, nullptr, 0);
|
||||
|
||||
vm_statistics64 vmstat;
|
||||
natural_t mcount = HOST_VM_INFO64_COUNT;
|
||||
if (host_statistics64(mach_host_self(), HOST_VM_INFO64, reinterpret_cast<host_info64_t>(&vmstat), &mcount) == KERN_SUCCESS)
|
||||
{
|
||||
this->get_statistics();
|
||||
#if 1 // count inactive as available
|
||||
this->avail = (int64_t(vmstat.free_count) +
|
||||
int64_t(vmstat.inactive_count)) *
|
||||
pagesize;
|
||||
this->used = (int64_t(vmstat.active_count) +
|
||||
int64_t(vmstat.wire_count)) *
|
||||
pagesize;
|
||||
#else // count inactive as unavailable
|
||||
this->avail = int64_t(vmstat.free_count) * pagesize;
|
||||
this->used = (int64_t(vmstat.active_count) +
|
||||
int64_t(vmstat.inactive_count) +
|
||||
int64_t(vmstat.wire_count)) *
|
||||
pagesize;
|
||||
#endif
|
||||
}
|
||||
|
||||
size_t get_TotalMem() const { return this->total; }
|
||||
size_t get_AvailMem() const { return this->avail; }
|
||||
size_t get_UsedMem() const { return this->used; }
|
||||
void update(){ this->get_statistics(); }
|
||||
|
||||
protected:
|
||||
int get_statistics( void )
|
||||
#elif __linux__
|
||||
FILE *fd;
|
||||
char buf[1024];
|
||||
if ((fd = fopen("/proc/meminfo", "r")))
|
||||
{
|
||||
#ifdef __APPLE__
|
||||
int pagesize = getpagesize();
|
||||
int mib[2] = {CTL_HW, HW_MEMSIZE};
|
||||
size_t length = sizeof(size_t);
|
||||
sysctl(mib, 2, &this->total, &length, nullptr, 0);
|
||||
|
||||
mach_msg_type_number_t count = HOST_VM_INFO_COUNT;
|
||||
vm_statistics_data_t vmstat;
|
||||
if(KERN_SUCCESS != host_statistics(mach_host_self(), HOST_VM_INFO, (host_info_t)&vmstat, &count)) {
|
||||
return -2;
|
||||
}
|
||||
this->avail = (int64_t)vmstat.free_count * (int64_t)pagesize;
|
||||
this->used = ((int64_t)vmstat.active_count +
|
||||
(int64_t)vmstat.inactive_count +
|
||||
(int64_t)vmstat.wire_count) * (int64_t)pagesize;
|
||||
|
||||
#elif __linux__
|
||||
FILE *fd;
|
||||
char buf[1024];
|
||||
if((fd = fopen("/proc/meminfo", "r")))
|
||||
while (1)
|
||||
{
|
||||
while(1)
|
||||
if (fgets(buf, 500, fd) != buf)
|
||||
break;
|
||||
if (bcmp(buf, "MemTotal", 8) == 0)
|
||||
{
|
||||
if(fgets(buf, 500, fd) != buf) break;
|
||||
if(bcmp(buf, "MemTotal", 8) == 0)
|
||||
{
|
||||
this->total = atoll(buf + 10) * 1024; // in Mb
|
||||
}
|
||||
if(strncmp(buf, "Committed_AS", 12) == 0)
|
||||
{
|
||||
this->used = atoll(buf + 14) * 1024; // in Mb
|
||||
}
|
||||
// if(strncmp(buf, "SwapTotal", 9) == 0)
|
||||
// {
|
||||
// *SwapTotal = atoll(buf + 11);
|
||||
// }
|
||||
// if(strncmp(buf, "SwapFree", 8) == 0)
|
||||
// {
|
||||
// *SwapFree = atoll(buf + 10);
|
||||
// }
|
||||
this->total = atoll(buf + 10) * 1024; // in Mb
|
||||
}
|
||||
fclose(fd);
|
||||
if (strncmp(buf, "Committed_AS", 12) == 0)
|
||||
{
|
||||
this->used = atoll(buf + 14) * 1024; // in Mb
|
||||
}
|
||||
// if(strncmp(buf, "SwapTotal", 9) == 0)
|
||||
// {
|
||||
// *SwapTotal = atoll(buf + 11);
|
||||
// }
|
||||
// if(strncmp(buf, "SwapFree", 8) == 0)
|
||||
// {
|
||||
// *SwapFree = atoll(buf + 10);
|
||||
// }
|
||||
}
|
||||
this->avail = this->total - this->used;
|
||||
|
||||
#endif
|
||||
return 0;
|
||||
|
||||
fclose(fd);
|
||||
}
|
||||
};
|
||||
this->avail = this->total - this->used;
|
||||
|
||||
#endif
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
|
||||
} /* namespace SystemStat */
|
||||
|
|
17
include/testing.hh
Normal file
17
include/testing.hh
Normal file
|
@ -0,0 +1,17 @@
|
|||
#pragma once
|
||||
|
||||
#include <array>
|
||||
#include <general.hh>
|
||||
#include <config_file.hh>
|
||||
#include <grid_fft.hh>
|
||||
|
||||
namespace testing{
|
||||
void output_potentials_and_densities(
|
||||
ConfigFile& the_config,
|
||||
size_t ngrid, real_t boxlen,
|
||||
Grid_FFT<real_t>& phi,
|
||||
Grid_FFT<real_t>& phi2,
|
||||
Grid_FFT<real_t>& phi3a,
|
||||
Grid_FFT<real_t>& phi3b,
|
||||
std::array< Grid_FFT<real_t>*,3 >& A3 );
|
||||
}
|
0
src/cosmology_calculator.cc
Normal file
0
src/cosmology_calculator.cc
Normal file
|
@ -273,7 +273,7 @@ void create_hdf5(std::string Filename)
|
|||
}
|
||||
|
||||
template <typename data_t>
|
||||
void Grid_FFT<data_t>::Write_to_HDF5(std::string fname, std::string datasetname)
|
||||
void Grid_FFT<data_t>::Write_to_HDF5(std::string fname, std::string datasetname) const
|
||||
{
|
||||
hid_t file_id, dset_id; /* file and dataset identifiers */
|
||||
hid_t filespace, memspace; /* file and memory dataspace identifiers */
|
||||
|
|
|
@ -2,11 +2,19 @@
|
|||
#include <general.hh>
|
||||
#include <grid_fft.hh>
|
||||
#include <convolution.hh>
|
||||
#include <testing.hh>
|
||||
|
||||
#include <ic_generator.hh>
|
||||
|
||||
#include <unistd.h> // for unlink
|
||||
|
||||
std::map<cosmo_species,std::string> cosmo_species_name =
|
||||
{
|
||||
{cosmo_species::dm,"Dark matter"},
|
||||
{cosmo_species::baryon,"Baryons"},
|
||||
{cosmo_species::neutrino,"Neutrinos"}
|
||||
};
|
||||
|
||||
namespace ic_generator{
|
||||
|
||||
std::unique_ptr<RNG_plugin> the_random_number_generator;
|
||||
|
@ -32,6 +40,7 @@ int Run( ConfigFile& the_config )
|
|||
const real_t boxlen = the_config.GetValue<double>("setup", "BoxLength");
|
||||
const real_t zstart = the_config.GetValue<double>("setup", "zstart");
|
||||
int LPTorder = the_config.GetValueSafe<double>("setup","LPTorder",100);
|
||||
const bool initial_bcc_lattice = the_config.GetValueSafe<bool>("setup","BCClattice",false);
|
||||
const bool bSymplecticPT = the_config.GetValueSafe<bool>("setup","SymplecticPT",false);
|
||||
const real_t astart = 1.0/(1.0+zstart);
|
||||
const real_t volfac(std::pow(boxlen / ngrid / 2.0 / M_PI, 1.5));
|
||||
|
@ -53,11 +62,11 @@ int Run( ConfigFile& the_config )
|
|||
const real_t Dplus0 = the_cosmo_calc->CalcGrowthFactor(astart) / the_cosmo_calc->CalcGrowthFactor(1.0);
|
||||
const real_t vfac = the_cosmo_calc->CalcVFact(astart);
|
||||
|
||||
const double g1 = Dplus0;
|
||||
const double g2 = (LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0;
|
||||
const double g3a = (LPTorder>2)? -1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0;
|
||||
const double g3b = (LPTorder>2)? 10.0/21.*Dplus0*Dplus0*Dplus0 : 0.0;
|
||||
const double g3c = (LPTorder>2)? -1.0/7.0*Dplus0*Dplus0*Dplus0 : (bSymplecticPT)? g1*g2 : 0.0;
|
||||
const double g1 = -Dplus0;
|
||||
const double g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0);
|
||||
const double g3a = ((LPTorder>2)? -1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||
const double g3b = ((LPTorder>2)? 10.0/21.*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||
const double g3c = ((LPTorder>2)? -1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||
|
||||
const double vfac1 = vfac;
|
||||
const double vfac2 = 2*vfac1;
|
||||
|
@ -89,275 +98,420 @@ int Run( ConfigFile& the_config )
|
|||
auto add_twice_to = [](auto &g){return [&](auto i, auto v){ g[i] += 2*v; };};
|
||||
auto subtract_from = [](auto &g){return [&](auto i, auto v){ g[i] -= v; };};
|
||||
auto subtract_twice_from = [](auto &g){return [&](auto i, auto v){ g[i] -= 2*v; };};
|
||||
|
||||
//--------------------------------------------------------------------
|
||||
std::vector<cosmo_species> species_list;
|
||||
|
||||
//--------------------------------------------------------------------
|
||||
// Fill the grid with a Gaussian white noise field
|
||||
//--------------------------------------------------------------------
|
||||
the_random_number_generator->Fill_Grid( phi );
|
||||
species_list.push_back( cosmo_species::dm );
|
||||
species_list.push_back( cosmo_species::baryon );
|
||||
|
||||
//======================================================================
|
||||
//... compute 1LPT displacement potential ....
|
||||
//======================================================================
|
||||
// phi = - delta / k^2
|
||||
double wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(1) term" << std::flush;
|
||||
phi.FourierTransformForward();
|
||||
|
||||
phi.apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
|
||||
real_t kmod = k.norm();
|
||||
if( bDoFixing ) x = (std::abs(x)!=0.0)? x / std::abs(x) : x; //std::exp(ccomplex_t(0, iphase * PhaseRotation));
|
||||
ccomplex_t delta = x * the_cosmo_calc->GetAmplitude(kmod, total);
|
||||
return -delta / (kmod * kmod) / volfac;
|
||||
});
|
||||
for( auto& this_species : species_list )
|
||||
{
|
||||
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
|
||||
csoca::ilog << ">> Computing ICs for species \'" << cosmo_species_name[this_species] << "\'" << std::endl;
|
||||
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
|
||||
|
||||
phi.zero_DC_mode();
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
//======================================================================
|
||||
//... compute 1LPT displacement potential ....
|
||||
//======================================================================
|
||||
// phi = - delta / k^2
|
||||
double wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(1) term" << std::flush;
|
||||
|
||||
//======================================================================
|
||||
//... compute 2LPT displacement potential ....
|
||||
//======================================================================
|
||||
if( LPTorder > 1 || bSymplecticPT ){
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(2) term" << std::flush;
|
||||
phi2.FourierTransformForward(false);
|
||||
Conv.convolve_SumOfHessians( phi, {0,0}, phi, {1,1}, {2,2}, assign_to( phi2 ) );
|
||||
Conv.convolve_Hessians( phi, {1,1}, phi, {2,2}, add_to(phi2) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, subtract_from(phi2) );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, subtract_from(phi2) );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, subtract_from(phi2) );
|
||||
phi2.apply_InverseLaplacian();
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
}
|
||||
#if 1 // random ICs
|
||||
//--------------------------------------------------------------------
|
||||
// Fill the grid with a Gaussian white noise field
|
||||
//--------------------------------------------------------------------
|
||||
the_random_number_generator->Fill_Grid( phi );
|
||||
|
||||
//======================================================================
|
||||
//... compute 3LPT displacement potential
|
||||
//======================================================================
|
||||
if( LPTorder > 2 && !bSymplecticPT ){
|
||||
//... 3a term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(3a) term" << std::flush;
|
||||
phi3a.FourierTransformForward(false);
|
||||
Conv.convolve_Hessians( phi, {0,0}, phi, {1,1}, phi, {2,2}, assign_to(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,2}, phi, {1,2}, add_twice_to(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi, {0,0}, subtract_from(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, phi, {1,1}, subtract_from(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi, {2,2}, subtract_from(phi3a) );
|
||||
phi3a.apply_InverseLaplacian();
|
||||
phi.FourierTransformForward();
|
||||
|
||||
phi.apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
|
||||
real_t kmod = k.norm();
|
||||
if( bDoFixing ) x = (std::abs(x)!=0.0)? x / std::abs(x) : x;
|
||||
ccomplex_t delta = x * the_cosmo_calc->GetAmplitude(kmod, total);
|
||||
return -delta / (kmod * kmod) / volfac;
|
||||
});
|
||||
|
||||
phi.zero_DC_mode();
|
||||
#else // ICs with a given phi(1) potential function
|
||||
constexpr real_t twopi{2.0*M_PI};
|
||||
constexpr real_t epsilon_q1d{0.25};
|
||||
|
||||
constexpr real_t epsy{0.25};
|
||||
constexpr real_t epsz{0.0};//epsz{0.25};
|
||||
|
||||
phi.FourierTransformBackward(false);
|
||||
|
||||
phi.apply_function_r_dep([&](auto v, auto r) -> real_t {
|
||||
real_t q1 = r[0]-0.5*boxlen;//r[0]/boxlen * twopi - M_PI;
|
||||
real_t q2 = r[1]-0.5*boxlen;//r[1]/boxlen * twopi - M_PI;
|
||||
real_t q3 = r[2]-0.5*boxlen;//r[1]/boxlen * twopi - M_PI;
|
||||
|
||||
// std::cerr << q1 << " " << q2 << std::endl;
|
||||
|
||||
return -2.0*std::cos(q1+std::cos(q2));
|
||||
// return (-std::cos(q1) + epsilon_q1d * std::sin(q2));
|
||||
// return (-std::cos(q1) + epsy * std::sin(q2) + epsz * std::cos(q1) * std::sin(q3));
|
||||
});
|
||||
phi.FourierTransformForward();
|
||||
|
||||
|
||||
#endif
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
//... 3b term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(3b) term" << std::flush;
|
||||
phi3b.FourierTransformForward(false);
|
||||
Conv.convolve_SumOfHessians( phi, {0,0}, phi2, {1,1}, {2,2}, assign_to(phi3b) );
|
||||
Conv.convolve_SumOfHessians( phi, {1,1}, phi2, {2,2}, {0,0}, add_to(phi3b) );
|
||||
Conv.convolve_SumOfHessians( phi, {2,2}, phi2, {0,0}, {1,1}, add_to(phi3b) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi2, {0,1}, subtract_twice_from(phi3b) );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi2, {0,2}, subtract_twice_from(phi3b) );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi2, {1,2}, subtract_twice_from(phi3b) );
|
||||
phi3b.apply_InverseLaplacian();
|
||||
phi3b *= 0.5; // factor 1/2 from definition of phi(3b)!
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
//... transversal term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing A(3) term" << std::flush;
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
int idimp = (idim+1)%3, idimpp = (idim+2)%3;
|
||||
A3[idim]->FourierTransformForward(false);
|
||||
Conv.convolve_Hessians( phi2, {idim,idimp}, phi, {idim,idimpp}, assign_to(*A3[idim]) );
|
||||
Conv.convolve_Hessians( phi2, {idim,idimpp}, phi, {idim,idimp}, subtract_from(*A3[idim]) );
|
||||
Conv.convolve_DifferenceOfHessians( phi, {idimp,idimpp}, phi2,{idimp,idimp}, {idimpp,idimpp}, add_to(*A3[idim]) );
|
||||
Conv.convolve_DifferenceOfHessians( phi2,{idimp,idimpp}, phi, {idimp,idimp}, {idimpp,idimpp}, subtract_from(*A3[idim]) );
|
||||
A3[idim]->apply_InverseLaplacian();
|
||||
//======================================================================
|
||||
//... compute 2LPT displacement potential ....
|
||||
//======================================================================
|
||||
if( LPTorder > 1 || bSymplecticPT ){
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(2) term" << std::flush;
|
||||
phi2.FourierTransformForward(false);
|
||||
Conv.convolve_SumOfHessians( phi, {0,0}, phi, {1,1}, {2,2}, assign_to( phi2 ) );
|
||||
Conv.convolve_Hessians( phi, {1,1}, phi, {2,2}, add_to(phi2) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, subtract_from(phi2) );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, subtract_from(phi2) );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, subtract_from(phi2) );
|
||||
phi2.apply_InverseLaplacian();
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
}
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
}
|
||||
|
||||
if( bSymplecticPT ){
|
||||
//... transversal term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing vNLO(3) term" << std::flush;
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
A3[idim]->FourierTransformForward(false);
|
||||
Conv.convolve_Gradient_and_Hessian( phi, {0}, phi2, {idim,0}, assign_to(*A3[idim]) );
|
||||
Conv.convolve_Gradient_and_Hessian( phi, {1}, phi2, {idim,1}, add_to(*A3[idim]) );
|
||||
Conv.convolve_Gradient_and_Hessian( phi, {2}, phi2, {idim,2}, add_to(*A3[idim]) );
|
||||
//======================================================================
|
||||
//... compute 3LPT displacement potential
|
||||
//======================================================================
|
||||
if( LPTorder > 2 && !bSymplecticPT ){
|
||||
//... 3a term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(3a) term" << std::flush;
|
||||
phi3a.FourierTransformForward(false);
|
||||
Conv.convolve_Hessians( phi, {0,0}, phi, {1,1}, phi, {2,2}, assign_to(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,2}, phi, {1,2}, add_twice_to(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi, {0,0}, subtract_from(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, phi, {1,1}, subtract_from(phi3a) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi, {2,2}, subtract_from(phi3a) );
|
||||
phi3a.apply_InverseLaplacian();
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
//... 3b term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(3b) term" << std::flush;
|
||||
phi3b.FourierTransformForward(false);
|
||||
Conv.convolve_SumOfHessians( phi, {0,0}, phi2, {1,1}, {2,2}, assign_to(phi3b) );
|
||||
Conv.convolve_SumOfHessians( phi, {1,1}, phi2, {2,2}, {0,0}, add_to(phi3b) );
|
||||
Conv.convolve_SumOfHessians( phi, {2,2}, phi2, {0,0}, {1,1}, add_to(phi3b) );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi2, {0,1}, subtract_twice_from(phi3b) );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi2, {0,2}, subtract_twice_from(phi3b) );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi2, {1,2}, subtract_twice_from(phi3b) );
|
||||
phi3b.apply_InverseLaplacian();
|
||||
phi3b *= 0.5; // factor 1/2 from definition of phi(3b)!
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
//... transversal term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing A(3) term" << std::flush;
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
int idimp = (idim+1)%3, idimpp = (idim+2)%3;
|
||||
A3[idim]->FourierTransformForward(false);
|
||||
Conv.convolve_Hessians( phi2, {idim,idimp}, phi, {idim,idimpp}, assign_to(*A3[idim]) );
|
||||
Conv.convolve_Hessians( phi2, {idim,idimpp}, phi, {idim,idimp}, subtract_from(*A3[idim]) );
|
||||
Conv.convolve_DifferenceOfHessians( phi, {idimp,idimpp}, phi2,{idimp,idimp}, {idimpp,idimpp}, add_to(*A3[idim]) );
|
||||
Conv.convolve_DifferenceOfHessians( phi2,{idimp,idimpp}, phi, {idimp,idimp}, {idimpp,idimpp}, subtract_from(*A3[idim]) );
|
||||
A3[idim]->apply_InverseLaplacian();
|
||||
}
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
}
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
}
|
||||
if( bSymplecticPT ){
|
||||
//... transversal term ...
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing vNLO(3) term" << std::flush;
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
A3[idim]->FourierTransformForward(false);
|
||||
Conv.convolve_Gradient_and_Hessian( phi, {0}, phi2, {idim,0}, assign_to(*A3[idim]) );
|
||||
Conv.convolve_Gradient_and_Hessian( phi, {1}, phi2, {idim,1}, add_to(*A3[idim]) );
|
||||
Conv.convolve_Gradient_and_Hessian( phi, {2}, phi2, {idim,2}, add_to(*A3[idim]) );
|
||||
}
|
||||
csoca::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
///... scale all potentials with respective growth factors
|
||||
phi *= g1;
|
||||
phi2 *= g2;
|
||||
phi3a *= g3a;
|
||||
phi3b *= g3b;
|
||||
(*A3[0]) *= g3c;
|
||||
(*A3[1]) *= g3c;
|
||||
(*A3[2]) *= g3c;
|
||||
}
|
||||
|
||||
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
|
||||
///... scale all potentials with respective growth factors
|
||||
phi *= g1;
|
||||
phi2 *= g2;
|
||||
phi3a *= g3a;
|
||||
phi3b *= g3b;
|
||||
(*A3[0]) *= g3c;
|
||||
(*A3[1]) *= g3c;
|
||||
(*A3[2]) *= g3c;
|
||||
|
||||
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl;
|
||||
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
// we store the densities here if we compute them
|
||||
//======================================================================
|
||||
const bool compute_densities = false;
|
||||
if( compute_densities ){
|
||||
const std::string fname_hdf5 = the_config.GetValueSafe<std::string>("output", "fname_hdf5", "output.hdf5");
|
||||
const std::string fname_analysis = the_config.GetValueSafe<std::string>("output", "fbase_analysis", "output");
|
||||
|
||||
Grid_FFT<real_t> delta({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3a({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3b({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
delta.FourierTransformForward(false);
|
||||
delta2.FourierTransformForward(false);
|
||||
delta3a.FourierTransformForward(false);
|
||||
delta3b.FourierTransformForward(false);
|
||||
delta3.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i)
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
// we store the densities here if we compute them
|
||||
//======================================================================
|
||||
const bool testing_compute_densities = false;
|
||||
if( testing_compute_densities )
|
||||
{
|
||||
for (size_t j = 0; j < phi.size(1); ++j)
|
||||
testing::output_potentials_and_densities( the_config, ngrid, boxlen, phi, phi2, phi3a, phi3b, A3 );
|
||||
}
|
||||
else
|
||||
{
|
||||
// temporary storage of data
|
||||
Grid_FFT<real_t> tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
|
||||
//if( the_output_plugin->write_species_as( cosmo_species::dm ) == output_type::field_eulerian ){
|
||||
if( the_output_plugin->write_species_as(this_species) == output_type::field_eulerian )
|
||||
{
|
||||
for (size_t k = 0; k < phi.size(2); ++k)
|
||||
//======================================================================
|
||||
// use QPT to get density and velocity fields
|
||||
//======================================================================
|
||||
Grid_FFT<ccomplex_t> psi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> rho({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
//======================================================================
|
||||
// initialise psi = exp(i Phi(1)/hbar)
|
||||
//======================================================================
|
||||
phi.FourierTransformBackward();
|
||||
real_t std_phi1 = phi.std();
|
||||
|
||||
const real_t hbar = 2.0 * M_PI/ngrid * (2*std_phi1/Dplus0); //3sigma, but this might rather depend on gradients of phi...
|
||||
csoca::ilog << "Semiclassical PT : hbar = " << hbar << " from sigma(phi1) = " << std_phi1 << std::endl;
|
||||
|
||||
if( LPTorder == 1 ){
|
||||
psi.assign_function_of_grids_r([hbar,Dplus0]( real_t pphi ){
|
||||
return std::exp(ccomplex_t(0.0,1.0/hbar) * (pphi / Dplus0));
|
||||
}, phi );
|
||||
}else if( LPTorder >= 2 ){
|
||||
phi2.FourierTransformBackward();
|
||||
// we don't have a 1/2 in the Veff term because pre-factor is already 3/7
|
||||
psi.assign_function_of_grids_r([hbar,Dplus0]( real_t pphi, real_t pphi2 ){
|
||||
return std::exp(ccomplex_t(0.0,1.0/hbar) * (pphi + pphi2) / Dplus0);
|
||||
}, phi, phi2 );
|
||||
// phi2.FourierTransformBackward();
|
||||
}
|
||||
// phi.FourierTransformForward();
|
||||
|
||||
//======================================================================
|
||||
// evolve wave-function (one drift step) psi = psi *exp(-i hbar *k^2 dt / 2)
|
||||
//======================================================================
|
||||
psi.FourierTransformForward();
|
||||
psi.apply_function_k_dep([hbar,Dplus0]( auto epsi, auto k ){
|
||||
auto k2 = k.norm_squared();
|
||||
return epsi * std::exp( - ccomplex_t(0.0,0.5)*hbar* k2 * Dplus0);
|
||||
});
|
||||
psi.FourierTransformBackward();
|
||||
|
||||
if( LPTorder >= 2 ){
|
||||
psi.assign_function_of_grids_r([&](auto ppsi, auto pphi2) {
|
||||
return ppsi * std::exp(ccomplex_t(0.0,1.0/hbar) * (pphi2) / Dplus0);
|
||||
}, psi, phi2);
|
||||
}
|
||||
|
||||
//======================================================================
|
||||
// compute rho
|
||||
//======================================================================
|
||||
rho.assign_function_of_grids_r([&]( auto p ){
|
||||
auto pp = std::real(p)*std::real(p) + std::imag(p)*std::imag(p) - 1.0;
|
||||
return pp;
|
||||
}, psi);
|
||||
|
||||
the_output_plugin->write_grid_data( rho, this_species, fluid_component::density );
|
||||
rho.Write_PowerSpectrum("input_powerspec_sampled_evolved_semiclassical.txt");
|
||||
rho.FourierTransformBackward();
|
||||
|
||||
//======================================================================
|
||||
// compute v
|
||||
//======================================================================
|
||||
for( int idim=0; idim<3; ++idim )
|
||||
{
|
||||
auto kk = phi.get_k<real_t>(i,j,k);
|
||||
size_t idx = phi.get_idx(i,j,k);
|
||||
auto laplace = -kk.norm_squared();
|
||||
Grid_FFT<ccomplex_t> grad_psi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
grad_psi.copy_from(psi);
|
||||
grad_psi.FourierTransformForward();
|
||||
grad_psi.apply_function_k_dep([&](auto x, auto k) {
|
||||
return x * ccomplex_t(0.0,k[idim]);
|
||||
});
|
||||
grad_psi.FourierTransformBackward();
|
||||
psi.FourierTransformBackward();
|
||||
|
||||
// compute densities associated to respective potentials as well
|
||||
delta.kelem(idx) = laplace * phi.kelem(idx);
|
||||
delta2.kelem(idx) = laplace * phi2.kelem(idx);
|
||||
delta3a.kelem(idx) = laplace * phi3a.kelem(idx);
|
||||
delta3b.kelem(idx) = laplace * phi3b.kelem(idx);
|
||||
delta3.kelem(idx) = delta3a.kelem(idx) + delta3b.kelem(idx);
|
||||
tmp.assign_function_of_grids_r([&](auto ppsi, auto pgrad_psi, auto prho) {
|
||||
return std::real((std::conj(ppsi) * pgrad_psi - ppsi * std::conj(pgrad_psi)) / ccomplex_t(0.0, 2.0 / hbar)/(1.0+prho));
|
||||
}, psi, grad_psi, rho);
|
||||
|
||||
fluid_component fc = (idim==0)? fluid_component::vx : ((idim==1)? fluid_component::vy : fluid_component::vz );
|
||||
the_output_plugin->write_grid_data( tmp, this_species, fc );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
delta.Write_PowerSpectrum(fname_analysis+"_"+"power_delta1.txt");
|
||||
delta2.Write_PowerSpectrum(fname_analysis+"_"+"power_delta2.txt");
|
||||
delta3a.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3a.txt");
|
||||
delta3b.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3b.txt");
|
||||
delta3.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3.txt");
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles
|
||||
|| the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
|
||||
{
|
||||
//======================================================================
|
||||
// we store displacements and velocities here if we compute them
|
||||
//======================================================================
|
||||
size_t num_p_in_load = phi.local_size();
|
||||
particle_container particles;
|
||||
|
||||
phi.FourierTransformBackward();
|
||||
phi2.FourierTransformBackward();
|
||||
phi3a.FourierTransformBackward();
|
||||
phi3b.FourierTransformBackward();
|
||||
// if output plugin wants particles, then we need to store them, along with their IDs
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
// if particles occupy a bcc lattice, then there are 2 x N^3 of them...
|
||||
if( initial_bcc_lattice )
|
||||
particles.allocate( 2*num_p_in_load );
|
||||
else
|
||||
particles.allocate( num_p_in_load );
|
||||
|
||||
delta.FourierTransformBackward();
|
||||
delta2.FourierTransformBackward();
|
||||
delta3a.FourierTransformBackward();
|
||||
delta3b.FourierTransformBackward();
|
||||
delta3.FourierTransformBackward();
|
||||
|
||||
A3[0]->FourierTransformBackward();
|
||||
A3[1]->FourierTransformBackward();
|
||||
A3[2]->FourierTransformBackward();
|
||||
|
||||
#if defined(USE_MPI)
|
||||
if( CONFIG::MPI_task_rank == 0 )
|
||||
unlink(fname_hdf5.c_str());
|
||||
MPI_Barrier( MPI_COMM_WORLD );
|
||||
#else
|
||||
unlink(fname_hdf5.c_str());
|
||||
#endif
|
||||
|
||||
phi.Write_to_HDF5(fname_hdf5, "phi");
|
||||
phi2.Write_to_HDF5(fname_hdf5, "phi2");
|
||||
phi3a.Write_to_HDF5(fname_hdf5, "phi3a");
|
||||
phi3b.Write_to_HDF5(fname_hdf5, "phi3b");
|
||||
|
||||
delta.Write_to_HDF5(fname_hdf5, "delta");
|
||||
delta2.Write_to_HDF5(fname_hdf5, "delta2");
|
||||
delta3a.Write_to_HDF5(fname_hdf5, "delta3a");
|
||||
delta3b.Write_to_HDF5(fname_hdf5, "delta3b");
|
||||
delta3.Write_to_HDF5(fname_hdf5, "delta3");
|
||||
|
||||
A3[0]->Write_to_HDF5(fname_hdf5, "A3x");
|
||||
A3[1]->Write_to_HDF5(fname_hdf5, "A3y");
|
||||
A3[2]->Write_to_HDF5(fname_hdf5, "A3z");
|
||||
|
||||
}else{
|
||||
//======================================================================
|
||||
// we store displacements and velocities here if we compute them
|
||||
//======================================================================
|
||||
Grid_FFT<real_t> tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
// write out positions
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
int idimp = (idim+1)%3, idimpp = (idim+2)%3;
|
||||
|
||||
tmp.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i) {
|
||||
for (size_t j = 0; j < phi.size(1); ++j) {
|
||||
for (size_t k = 0; k < phi.size(2); ++k) {
|
||||
auto kk = phi.get_k<real_t>(i,j,k);
|
||||
size_t idx = phi.get_idx(i,j,k);
|
||||
|
||||
auto phitot = phi.kelem(idx) + phi2.kelem(idx) + phi3a.kelem(idx) + phi3b.kelem(idx);
|
||||
|
||||
// divide by Lbox, because displacement is in box units for output plugin
|
||||
tmp.kelem(idx) = ccomplex_t(0.0,1.0) * (kk[idim] * phitot + kk[idimp] * A3[idimpp]->kelem(idx) - kk[idimpp] * A3[idimp]->kelem(idx) ) / boxlen;
|
||||
}
|
||||
}
|
||||
}
|
||||
tmp.FourierTransformBackward();
|
||||
the_output_plugin->write_dm_position(idim, tmp );
|
||||
}
|
||||
|
||||
// write out velocities
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
int idimp = (idim+1)%3, idimpp = (idim+2)%3;
|
||||
|
||||
tmp.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i) {
|
||||
for (size_t j = 0; j < phi.size(1); ++j) {
|
||||
for (size_t k = 0; k < phi.size(2); ++k) {
|
||||
auto kk = phi.get_k<real_t>(i,j,k);
|
||||
size_t idx = phi.get_idx(i,j,k);
|
||||
|
||||
|
||||
// divide by Lbox, because displacement is in box units for output plugin
|
||||
if(!bSymplecticPT){
|
||||
auto phitot_v = vfac1 * phi.kelem(idx) + vfac2 * phi2.kelem(idx) + vfac3 * (phi3a.kelem(idx) + phi3b.kelem(idx));
|
||||
tmp.kelem(idx) = ccomplex_t(0.0,1.0) * (kk[idim] * phitot_v + vfac3 * (kk[idimp] * A3[idimpp]->kelem(idx) - kk[idimpp] * A3[idimp]->kelem(idx)) ) / boxlen;
|
||||
}else{
|
||||
auto phitot_v = vfac1 * phi.kelem(idx) + vfac2 * phi2.kelem(idx);
|
||||
tmp.kelem(idx) = ccomplex_t(0.0,1.0) * (kk[idim] * phitot_v) + vfac1 * A3[idim]->kelem(idx);
|
||||
// generate particle IDs
|
||||
auto ipcount0 = (initial_bcc_lattice? 2:1) * particles.get_local_offset();
|
||||
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
|
||||
for( size_t j=0; j<tmp.size(1); ++j){
|
||||
for( size_t k=0; k<tmp.size(2); ++k){
|
||||
particles.set_id( ipcount, ipcount+ipcount0 );
|
||||
++ipcount;
|
||||
if( initial_bcc_lattice ){
|
||||
particles.set_id( ipcount, ipcount+ipcount0 );
|
||||
++ipcount;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// write out positions
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
const int idimp = (idim+1)%3, idimpp = (idim+2)%3;
|
||||
const real_t lunit = the_output_plugin->position_unit();
|
||||
tmp.FourierTransformForward(false);
|
||||
|
||||
// combine the various LPT potentials into one and take gradient
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i) {
|
||||
for (size_t j = 0; j < phi.size(1); ++j) {
|
||||
for (size_t k = 0; k < phi.size(2); ++k) {
|
||||
auto kk = phi.get_k<real_t>(i,j,k);
|
||||
size_t idx = phi.get_idx(i,j,k);
|
||||
auto phitot = phi.kelem(idx) + phi2.kelem(idx) + phi3a.kelem(idx) + phi3b.kelem(idx);
|
||||
// divide by Lbox, because displacement is in box units for output plugin
|
||||
tmp.kelem(idx) = lunit * ccomplex_t(0.0,1.0) * (kk[idim] * phitot + kk[idimp] * A3[idimpp]->kelem(idx) - kk[idimpp] * A3[idimp]->kelem(idx) ) / boxlen;
|
||||
}
|
||||
}
|
||||
}
|
||||
tmp.FourierTransformBackward();
|
||||
|
||||
// if we write particle data, store particle data in particle structure
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
|
||||
for( size_t j=0; j<tmp.size(1); ++j){
|
||||
for( size_t k=0; k<tmp.size(2); ++k){
|
||||
auto pos = tmp.get_unit_r<float>(i,j,k);
|
||||
particles.set_pos( ipcount++, idim, pos[idim]*lunit + tmp.relem(i,j,k) );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if( initial_bcc_lattice ){
|
||||
tmp.stagger_field();
|
||||
auto ipcount0 = num_p_in_load;
|
||||
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
|
||||
for( size_t j=0; j<tmp.size(1); ++j){
|
||||
for( size_t k=0; k<tmp.size(2); ++k){
|
||||
auto pos = tmp.get_unit_r_staggered<float>(i,j,k);
|
||||
particles.set_pos( ipcount++, idim, pos[idim]*lunit + tmp.relem(i,j,k) );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// otherwise write out the grid data directly to the output plugin
|
||||
// else if( the_output_plugin->write_species_as( cosmo_species::dm ) == output_type::field_lagrangian )
|
||||
else if( the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
|
||||
{
|
||||
fluid_component fc = (idim==0)? fluid_component::dx : ((idim==1)? fluid_component::dy : fluid_component::dz );
|
||||
the_output_plugin->write_grid_data( tmp, this_species, fc );
|
||||
}
|
||||
}
|
||||
|
||||
// write out velocities
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
// cyclic rotations of indices
|
||||
int idimp = (idim+1)%3, idimpp = (idim+2)%3;
|
||||
const real_t vunit = the_output_plugin->velocity_unit();
|
||||
|
||||
tmp.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i) {
|
||||
for (size_t j = 0; j < phi.size(1); ++j) {
|
||||
for (size_t k = 0; k < phi.size(2); ++k) {
|
||||
auto kk = phi.get_k<real_t>(i,j,k);
|
||||
size_t idx = phi.get_idx(i,j,k);
|
||||
// divide by Lbox, because displacement is in box units for output plugin
|
||||
if(!bSymplecticPT){
|
||||
auto phitot_v = vfac1 * phi.kelem(idx) + vfac2 * phi2.kelem(idx) + vfac3 * (phi3a.kelem(idx) + phi3b.kelem(idx));
|
||||
tmp.kelem(idx) = vunit*ccomplex_t(0.0,1.0) * (kk[idim] * phitot_v + vfac3 * (kk[idimp] * A3[idimpp]->kelem(idx) - kk[idimpp] * A3[idimp]->kelem(idx)) ) / boxlen;
|
||||
}else{
|
||||
auto phitot_v = vfac1 * phi.kelem(idx) + vfac2 * phi2.kelem(idx);
|
||||
tmp.kelem(idx) = vunit*ccomplex_t(0.0,1.0) * (kk[idim] * phitot_v) + vfac1 * A3[idim]->kelem(idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
tmp.FourierTransformBackward();
|
||||
|
||||
// if we write particle data, store particle data in particle structure
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles ){
|
||||
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
|
||||
for( size_t j=0; j<tmp.size(1); ++j){
|
||||
for( size_t k=0; k<tmp.size(2); ++k){
|
||||
particles.set_vel( ipcount++, idim, tmp.relem(i,j,k) );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if( initial_bcc_lattice ){
|
||||
tmp.stagger_field();
|
||||
auto ipcount0 = num_p_in_load;
|
||||
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
|
||||
for( size_t j=0; j<tmp.size(1); ++j){
|
||||
for( size_t k=0; k<tmp.size(2); ++k){
|
||||
particles.set_vel( ipcount++, idim, tmp.relem(i,j,k) );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}// otherwise write out the grid data directly to the output plugin
|
||||
else if( the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
|
||||
{
|
||||
fluid_component fc = (idim==0)? fluid_component::vx : ((idim==1)? fluid_component::vy : fluid_component::vz );
|
||||
the_output_plugin->write_grid_data( tmp, this_species, fc );
|
||||
}
|
||||
}
|
||||
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
the_output_plugin->write_particle_data( particles, this_species );
|
||||
}
|
||||
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
|
||||
{
|
||||
// use density simply from 1st order SPT
|
||||
phi.FourierTransformForward();
|
||||
phi.apply_negative_Laplacian();
|
||||
phi.Write_PowerSpectrum("input_powerspec_sampled_SPT.txt");
|
||||
phi.FourierTransformBackward();
|
||||
the_output_plugin->write_grid_data( phi, this_species, fluid_component::density );
|
||||
}
|
||||
}
|
||||
tmp.FourierTransformBackward();
|
||||
the_output_plugin->write_dm_velocity(idim, tmp );
|
||||
|
||||
}
|
||||
|
||||
the_output_plugin->write_dm_mass(tmp);
|
||||
the_output_plugin->write_dm_density(tmp);
|
||||
|
||||
the_output_plugin->finalize();
|
||||
//delete the_output_plugin;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
} // end namespace ic_generator
|
||||
|
||||
|
|
23
src/main.cc
23
src/main.cc
|
@ -49,9 +49,25 @@ int main( int argc, char** argv )
|
|||
}
|
||||
#endif
|
||||
|
||||
csoca::ilog << "-----------------------------------------------------------------------------" << std::endl
|
||||
<< ">> FastLPT v0.1 >>" << std::endl
|
||||
<< "-----------------------------------------------------------------------------" << std::endl;
|
||||
csoca::ilog //<< " __ _____ _____ \n"
|
||||
//<< " / _| |_ _/ __ \\ \n"
|
||||
//<< " _ __ ___ ___ _ __ ___ | |_ ___ _ __ | | | / \\/ \n"
|
||||
//<< " | '_ ` _ \\ / _ \\| '_ \\ / _ \\| _/ _ \\| '_ \\ | | | | \n"
|
||||
//<< " | | | | | | (_) | | | | (_) | || (_) | | | || |_| \\__/\\ \n"
|
||||
//<< " |_| |_| |_|\\___/|_| |_|\\___/|_| \\___/|_| |_\\___/ \\____/ \n"
|
||||
//<< " \n"
|
||||
<< "\n\n"
|
||||
<< " .8888b dP a88888b. \n"
|
||||
<< " 88 \" 88 d8\' `88 \n"
|
||||
<< " 88d8b.d8b. .d8888b. 88d888b. .d8888b. 88aaa .d8888b. 88d888b. 88 88 \n"
|
||||
<< " 88\'`88\'`88 88\' `88 88\' `88 88\' `88 88 88\' `88 88\' `88 88 88 \n"
|
||||
<< " 88 88 88 88. .88 88 88 88. .88 88 88. .88 88 88 88 Y8. .88 \n"
|
||||
<< " dP dP dP `88888P\' dP dP `88888P\' dP `88888P\' dP dP dP Y88888P\' \n"
|
||||
<< " \n"
|
||||
<< "--------------------------------------------------------------------------------" << std::endl
|
||||
<< " version : v0.1a" << std::endl
|
||||
<< " git rev. : " << GIT_REV << ", tag: " << GIT_TAG << ", branch: " << GIT_BRANCH <<"\n"
|
||||
<< "--------------------------------------------------------------------------------" << std::endl;
|
||||
|
||||
|
||||
//------------------------------------------------------------------------------
|
||||
|
@ -107,6 +123,7 @@ int main( int argc, char** argv )
|
|||
//------------------------------------------------------------------------------
|
||||
// Write code configuration to screen
|
||||
//------------------------------------------------------------------------------
|
||||
csoca::ilog << std::setw(40) << std::left << "CPU vendor string" << " : " << SystemStat::Cpu().get_CPUstring() << std::endl;
|
||||
#if defined(USE_MPI)
|
||||
csoca::ilog << std::setw(40) << std::left << "MPI is enabled" << " : " << "yes (" << CONFIG::MPI_task_size << " tasks)" << std::endl;
|
||||
#else
|
||||
|
|
1085
src/plugins/HDF_IO.hh
Executable file
1085
src/plugins/HDF_IO.hh
Executable file
File diff suppressed because it is too large
Load diff
File diff suppressed because it is too large
Load diff
106
src/plugins/output_generic.cc
Normal file
106
src/plugins/output_generic.cc
Normal file
|
@ -0,0 +1,106 @@
|
|||
/*
|
||||
|
||||
output_generic.cc - This file is part of MUSIC2 - GPL
|
||||
Copyright (C) 2010-19 Oliver Hahn
|
||||
|
||||
*/
|
||||
|
||||
#ifdef USE_HDF5
|
||||
|
||||
#include <unistd.h> // for unlink
|
||||
|
||||
#include "HDF_IO.hh"
|
||||
#include <logger.hh>
|
||||
#include <output_plugin.hh>
|
||||
|
||||
class generic_output_plugin : public output_plugin
|
||||
{
|
||||
private:
|
||||
std::string get_field_name( const cosmo_species &s, const fluid_component &c );
|
||||
protected:
|
||||
bool out_eulerian_;
|
||||
public:
|
||||
//! constructor
|
||||
explicit generic_output_plugin(ConfigFile &cf )
|
||||
: output_plugin(cf, "Generic HDF5")
|
||||
{
|
||||
real_t astart = 1.0/(1.0+cf_.GetValue<double>("setup", "zstart"));
|
||||
real_t boxsize = cf_.GetValue<double>("setup", "BoxLength");
|
||||
|
||||
out_eulerian_ = cf_.GetValueSafe<bool>("output", "generic_out_eulerian",false);
|
||||
|
||||
if( CONFIG::MPI_task_rank == 0 )
|
||||
{
|
||||
unlink(fname_.c_str());
|
||||
HDFCreateFile( fname_ );
|
||||
HDFCreateGroup( fname_, "Header" );
|
||||
HDFWriteGroupAttribute<double>( fname_, "Header", "Boxsize", boxsize );
|
||||
HDFWriteGroupAttribute<double>( fname_, "Header", "astart", astart );
|
||||
}
|
||||
|
||||
#if defined(USE_MPI)
|
||||
MPI_Barrier( MPI_COMM_WORLD );
|
||||
#endif
|
||||
}
|
||||
|
||||
output_type write_species_as( const cosmo_species &s ) const
|
||||
{
|
||||
if( out_eulerian_ )
|
||||
return output_type::field_eulerian;
|
||||
return output_type::field_lagrangian;
|
||||
}
|
||||
|
||||
real_t position_unit() const { return 1.0; }
|
||||
|
||||
real_t velocity_unit() const { return 1.0; }
|
||||
|
||||
void write_grid_data(const Grid_FFT<real_t> &g, const cosmo_species &s, const fluid_component &c );
|
||||
};
|
||||
|
||||
|
||||
std::string generic_output_plugin::get_field_name( const cosmo_species &s, const fluid_component &c )
|
||||
{
|
||||
std::string field_name;
|
||||
switch( s ){
|
||||
case cosmo_species::dm:
|
||||
field_name += "DM"; break;
|
||||
case cosmo_species::baryon:
|
||||
field_name += "BA"; break;
|
||||
case cosmo_species::neutrino:
|
||||
field_name += "NU"; break;
|
||||
default: break;
|
||||
}
|
||||
field_name += "_";
|
||||
switch( c ){
|
||||
case fluid_component::density:
|
||||
field_name += "delta"; break;
|
||||
case fluid_component::vx:
|
||||
field_name += "vx"; break;
|
||||
case fluid_component::vy:
|
||||
field_name += "vy"; break;
|
||||
case fluid_component::vz:
|
||||
field_name += "vz"; break;
|
||||
case fluid_component::dx:
|
||||
field_name += "dx"; break;
|
||||
case fluid_component::dy:
|
||||
field_name += "dy"; break;
|
||||
case fluid_component::dz:
|
||||
field_name += "dz"; break;
|
||||
default: break;
|
||||
}
|
||||
return field_name;
|
||||
}
|
||||
|
||||
void generic_output_plugin::write_grid_data(const Grid_FFT<real_t> &g, const cosmo_species &s, const fluid_component &c )
|
||||
{
|
||||
std::string field_name = this->get_field_name( s, c );
|
||||
g.Write_to_HDF5(fname_, field_name);
|
||||
csoca::ilog << interface_name_ << " : Wrote field \'" << field_name << "\' to file \'" << fname_ << "\'" << std::endl;
|
||||
}
|
||||
|
||||
namespace
|
||||
{
|
||||
output_plugin_creator_concrete<generic_output_plugin> creator1("generic");
|
||||
} // namespace
|
||||
|
||||
#endif
|
284
src/plugins/output_grafic2.cc
Normal file
284
src/plugins/output_grafic2.cc
Normal file
|
@ -0,0 +1,284 @@
|
|||
#include <unistd.h> // for unlink
|
||||
#include <sys/types.h>
|
||||
#include <sys/stat.h>
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
|
||||
#include <general.hh>
|
||||
#include <output_plugin.hh>
|
||||
|
||||
//! Implementation of class grafic2_output_plugin
|
||||
/*!
|
||||
This class implements a grafic-2 (cf. Bertschinger 2001) compatible
|
||||
output format. With some RAMSES extras.
|
||||
*/
|
||||
class grafic2_output_plugin : public output_plugin
|
||||
{
|
||||
private:
|
||||
struct header
|
||||
{
|
||||
int n1, n2, n3;
|
||||
float dxini0;
|
||||
float xoff10, xoff20, xoff30;
|
||||
float astart0, omega_m0, omega_l0, h00;
|
||||
};
|
||||
|
||||
std::string get_file_name(const cosmo_species &s, const fluid_component &c) const;
|
||||
void write_ramses_namelist(void) const;
|
||||
|
||||
protected:
|
||||
header header_;
|
||||
real_t lunit_, vunit_;
|
||||
uint32_t levelmin_;
|
||||
bool bhavebaryons_;
|
||||
std::vector<float> data_buf_;
|
||||
std::string dirname_;
|
||||
bool bUseSPT_;
|
||||
|
||||
public:
|
||||
//! constructor
|
||||
explicit grafic2_output_plugin(ConfigFile &cf)
|
||||
: output_plugin(cf, "GRAFIC2/RAMSES")
|
||||
{
|
||||
lunit_ = 1.0;
|
||||
vunit_ = 1.0;
|
||||
|
||||
double
|
||||
boxlength = cf_.GetValue<double>("setup", "BoxLength"),
|
||||
H0 = cf_.GetValue<double>("cosmology", "H0"),
|
||||
zstart = cf_.GetValue<double>("setup", "zstart"),
|
||||
astart = 1.0 / (1.0 + zstart),
|
||||
omegam = cf_.GetValue<double>("cosmology", "Omega_m"),
|
||||
omegaL = cf_.GetValue<double>("cosmology", "Omega_L");
|
||||
uint32_t ngrid = cf_.GetValue<int>("setup", "GridRes");
|
||||
|
||||
bUseSPT_ = cf_.GetValueSafe<bool>("output", "grafic_use_SPT", false);
|
||||
levelmin_ = uint32_t(std::log2(double(ngrid)) + 1e-6);
|
||||
|
||||
if (std::abs(std::pow(2.0, levelmin_) - double(ngrid)) > 1e-4)
|
||||
{
|
||||
csoca::elog << interface_name_ << " plugin requires setup/GridRes to be power of 2!" << std::endl;
|
||||
abort();
|
||||
}
|
||||
|
||||
bhavebaryons_ = cf_.GetValueSafe<bool>("setup", "baryons", false);
|
||||
|
||||
header_.n1 = ngrid;
|
||||
header_.n2 = ngrid;
|
||||
header_.n3 = ngrid;
|
||||
header_.dxini0 = boxlength / (H0 * 0.01) / ngrid;
|
||||
header_.xoff10 = 0;
|
||||
header_.xoff20 = 0;
|
||||
header_.xoff30 = 0;
|
||||
header_.astart0 = astart;
|
||||
header_.omega_m0 = omegam;
|
||||
header_.omega_l0 = omegaL;
|
||||
header_.h00 = H0;
|
||||
|
||||
data_buf_.assign(ngrid * ngrid, 0.0f);
|
||||
|
||||
lunit_ = boxlength;
|
||||
vunit_ = boxlength;
|
||||
|
||||
// create directory structure
|
||||
dirname_ = this->fname_;
|
||||
remove(dirname_.c_str());
|
||||
mkdir(dirname_.c_str(), 0777);
|
||||
|
||||
// write RAMSES namelist file? if so only with one task
|
||||
if (cf_.GetValueSafe<bool>("output", "ramses_nml", true) && CONFIG::MPI_task_rank==0 )
|
||||
{
|
||||
write_ramses_namelist();
|
||||
}
|
||||
}
|
||||
|
||||
output_type write_species_as(const cosmo_species &s) const
|
||||
{
|
||||
if (s == cosmo_species::baryon && !bUseSPT_)
|
||||
return output_type::field_eulerian;
|
||||
return output_type::field_lagrangian;
|
||||
}
|
||||
|
||||
real_t position_unit() const { return lunit_; }
|
||||
|
||||
real_t velocity_unit() const { return vunit_; }
|
||||
|
||||
void write_grid_data(const Grid_FFT<real_t> &g, const cosmo_species &s, const fluid_component &c);
|
||||
};
|
||||
|
||||
std::string grafic2_output_plugin::get_file_name(const cosmo_species &s, const fluid_component &c) const
|
||||
{
|
||||
std::string file_name, species_str;
|
||||
|
||||
file_name = dirname_ + "/ic_";
|
||||
|
||||
switch (s)
|
||||
{
|
||||
case cosmo_species::dm:
|
||||
species_str = "c";
|
||||
break;
|
||||
case cosmo_species::baryon:
|
||||
species_str = "b";
|
||||
break;
|
||||
case cosmo_species::neutrino:
|
||||
species_str = "n";
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
switch (c)
|
||||
{
|
||||
case fluid_component::density:
|
||||
file_name += "delta" + species_str;
|
||||
break;
|
||||
case fluid_component::vx:
|
||||
file_name += "vel" + species_str + "z";
|
||||
break;
|
||||
case fluid_component::vy:
|
||||
file_name += "vel" + species_str + "y";
|
||||
break;
|
||||
case fluid_component::vz:
|
||||
file_name += "vel" + species_str + "x";
|
||||
break;
|
||||
case fluid_component::dx:
|
||||
file_name += "pos" + species_str + "z";
|
||||
break;
|
||||
case fluid_component::dy:
|
||||
file_name += "pos" + species_str + "y";
|
||||
break;
|
||||
case fluid_component::dz:
|
||||
file_name += "pos" + species_str + "x";
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
return file_name;
|
||||
}
|
||||
|
||||
void grafic2_output_plugin::write_grid_data(const Grid_FFT<real_t> &g, const cosmo_species &s, const fluid_component &c)
|
||||
{
|
||||
// ignore certain components
|
||||
if (s == cosmo_species::dm && c == fluid_component::density)
|
||||
return;
|
||||
if (s == cosmo_species::baryon && (c == fluid_component::dx || c == fluid_component::dy || c == fluid_component::dz))
|
||||
return;
|
||||
|
||||
// get file name based on species and fluid component type
|
||||
std::string file_name = this->get_file_name(s, c);
|
||||
|
||||
// serialize parallel write
|
||||
for (int write_rank = 0; write_rank < CONFIG::MPI_task_size; ++write_rank)
|
||||
{
|
||||
if (write_rank == CONFIG::MPI_task_rank)
|
||||
{
|
||||
if (write_rank == 0)
|
||||
{
|
||||
unlink(file_name.c_str());
|
||||
}
|
||||
std::ofstream ofs(file_name.c_str(), std::ios::binary|std::ios::app);
|
||||
|
||||
// write header or seek to end of file
|
||||
if (write_rank == 0)
|
||||
{
|
||||
uint32_t blocksz = sizeof(header);
|
||||
ofs.write(reinterpret_cast<const char *>(&blocksz), sizeof(int));
|
||||
ofs.write(reinterpret_cast<const char *>(&header_), blocksz);
|
||||
ofs.write(reinterpret_cast<const char *>(&blocksz), sizeof(int));
|
||||
}
|
||||
|
||||
// check field size against buffer size...
|
||||
uint32_t ngrid = cf_.GetValue<int>("setup", "GridRes");
|
||||
assert( g.global_size(0) == ngrid && g.global_size(1) == ngrid && g.global_size(2) == ngrid);
|
||||
assert( g.size(1) == ngrid && g.size(2) == ngrid);
|
||||
// write actual field slice by slice
|
||||
for (size_t i = 0; i < g.size(0); ++i)
|
||||
{
|
||||
for (unsigned j = 0; j < g.size(1); ++j)
|
||||
{
|
||||
for (unsigned k = 0; k < g.size(2); ++k)
|
||||
{
|
||||
data_buf_[j * ngrid + k] = g.relem(i, j, k);
|
||||
}
|
||||
}
|
||||
|
||||
uint32_t blocksz = ngrid * ngrid * sizeof(float);
|
||||
ofs.write(reinterpret_cast<const char *>(&blocksz), sizeof(uint32_t));
|
||||
ofs.write(reinterpret_cast<const char *>(&data_buf_[0]), blocksz);
|
||||
ofs.write(reinterpret_cast<const char *>(&blocksz), sizeof(uint32_t));
|
||||
}
|
||||
|
||||
ofs.close();
|
||||
}
|
||||
|
||||
multitask_sync_barrier();
|
||||
|
||||
} // end loop over write_rank
|
||||
|
||||
csoca::ilog << interface_name_ << " : Wrote field to file \'" << file_name << "\'" << std::endl;
|
||||
}
|
||||
|
||||
void grafic2_output_plugin::write_ramses_namelist(void) const
|
||||
{
|
||||
//... also write the refinement options to a dummy namelist file
|
||||
char ff[256];
|
||||
sprintf(ff, "%s/ramses.nml", dirname_.c_str());
|
||||
|
||||
std::ofstream ofst(ff, std::ios::trunc);
|
||||
|
||||
// -- RUN_PARAMS -- //
|
||||
ofst
|
||||
<< "&RUN_PARAMS\n"
|
||||
<< "cosmo=.true.\n"
|
||||
<< "pic=.true.\n"
|
||||
<< "poisson=.true.\n";
|
||||
|
||||
if (bhavebaryons_)
|
||||
ofst << "hydro=.true.\n";
|
||||
else
|
||||
ofst << "hydro=.false.\n";
|
||||
|
||||
ofst
|
||||
<< "nrestart=0\n"
|
||||
<< "nremap=1\n"
|
||||
<< "nsubcycle=1,2\n"
|
||||
<< "ncontrol=1\n"
|
||||
<< "verbose=.false.\n"
|
||||
<< "/\n\n";
|
||||
|
||||
// -- INIT_PARAMS -- //
|
||||
ofst
|
||||
<< "&INIT_PARAMS\n"
|
||||
<< "filetype=\'grafic\'\n"
|
||||
<< "initfile(1)=\'" << dirname_ << "\'\n"
|
||||
<< "/\n\n";
|
||||
|
||||
// initialize with settings for naddref additional levels of refinement
|
||||
unsigned naddref = 5;
|
||||
|
||||
// -- AMR_PARAMS -- //
|
||||
ofst << "&AMR_PARAMS\n"
|
||||
<< "levelmin=" << levelmin_ << "\n"
|
||||
<< "levelmax=" << levelmin_ + naddref << "\n"
|
||||
<< "nexpand=1";
|
||||
|
||||
ofst << "\n"
|
||||
<< "ngridmax=2000000\n"
|
||||
<< "npartmax=3000000\n"
|
||||
<< "/\n\n";
|
||||
|
||||
ofst << "&REFINE_PARAMS\n"
|
||||
<< "m_refine=" << 1 + naddref << "*8.,\n"
|
||||
<< "/\n";
|
||||
|
||||
csoca::ilog << interface_name_ << " wrote partial RAMSES namelist file \'" << fname_ << "\'" << std::endl;
|
||||
}
|
||||
|
||||
namespace
|
||||
{
|
||||
output_plugin_creator_concrete<grafic2_output_plugin> creator1("grafic2");
|
||||
} // namespace
|
|
@ -64,6 +64,7 @@ public:
|
|||
|
||||
void Fill_Grid(Grid_FFT<real_t> &g) const
|
||||
{
|
||||
g.zero();
|
||||
g.FourierTransformForward(false);
|
||||
|
||||
// transform is transposed!
|
||||
|
@ -88,7 +89,7 @@ public:
|
|||
|
||||
do {
|
||||
ampl = gsl_rng_uniform(pRandomGenerator_);
|
||||
} while (ampl == 0);
|
||||
} while (ampl == 0||ampl == 1);
|
||||
|
||||
if (i == nres_ / 2 || j == nres_ / 2 || k == nres_ / 2) continue;
|
||||
if (i == 0 && j == 0 && k == 0) continue;
|
||||
|
|
91
src/testing.cc
Normal file
91
src/testing.cc
Normal file
|
@ -0,0 +1,91 @@
|
|||
#include <testing.hh>
|
||||
#include <unistd.h> // for unlink
|
||||
|
||||
void output_potentials_and_densities(
|
||||
ConfigFile& the_config,
|
||||
size_t ngrid, real_t boxlen,
|
||||
Grid_FFT<real_t> &phi,
|
||||
Grid_FFT<real_t> &phi2,
|
||||
Grid_FFT<real_t> &phi3a,
|
||||
Grid_FFT<real_t> &phi3b,
|
||||
std::array<Grid_FFT<real_t> *, 3> &A3)
|
||||
{
|
||||
const std::string fname_hdf5 = the_config.GetValueSafe<std::string>("output", "fname_hdf5", "output.hdf5");
|
||||
const std::string fname_analysis = the_config.GetValueSafe<std::string>("output", "fbase_analysis", "output");
|
||||
|
||||
Grid_FFT<real_t> delta({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3a({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3b({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
Grid_FFT<real_t> delta3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
delta.FourierTransformForward(false);
|
||||
delta2.FourierTransformForward(false);
|
||||
delta3a.FourierTransformForward(false);
|
||||
delta3b.FourierTransformForward(false);
|
||||
delta3.FourierTransformForward(false);
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < phi.size(0); ++i)
|
||||
{
|
||||
for (size_t j = 0; j < phi.size(1); ++j)
|
||||
{
|
||||
for (size_t k = 0; k < phi.size(2); ++k)
|
||||
{
|
||||
auto kk = phi.get_k<real_t>(i, j, k);
|
||||
size_t idx = phi.get_idx(i, j, k);
|
||||
auto laplace = -kk.norm_squared();
|
||||
|
||||
// compute densities associated to respective potentials as well
|
||||
delta.kelem(idx) = laplace * phi.kelem(idx);
|
||||
delta2.kelem(idx) = laplace * phi2.kelem(idx);
|
||||
delta3a.kelem(idx) = laplace * phi3a.kelem(idx);
|
||||
delta3b.kelem(idx) = laplace * phi3b.kelem(idx);
|
||||
delta3.kelem(idx) = delta3a.kelem(idx) + delta3b.kelem(idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
delta.Write_PowerSpectrum(fname_analysis + "_" + "power_delta1.txt");
|
||||
delta2.Write_PowerSpectrum(fname_analysis + "_" + "power_delta2.txt");
|
||||
delta3a.Write_PowerSpectrum(fname_analysis + "_" + "power_delta3a.txt");
|
||||
delta3b.Write_PowerSpectrum(fname_analysis + "_" + "power_delta3b.txt");
|
||||
delta3.Write_PowerSpectrum(fname_analysis + "_" + "power_delta3.txt");
|
||||
|
||||
phi.FourierTransformBackward();
|
||||
phi2.FourierTransformBackward();
|
||||
phi3a.FourierTransformBackward();
|
||||
phi3b.FourierTransformBackward();
|
||||
|
||||
delta.FourierTransformBackward();
|
||||
delta2.FourierTransformBackward();
|
||||
delta3a.FourierTransformBackward();
|
||||
delta3b.FourierTransformBackward();
|
||||
delta3.FourierTransformBackward();
|
||||
|
||||
A3[0]->FourierTransformBackward();
|
||||
A3[1]->FourierTransformBackward();
|
||||
A3[2]->FourierTransformBackward();
|
||||
|
||||
#if defined(USE_MPI)
|
||||
if (CONFIG::MPI_task_rank == 0)
|
||||
unlink(fname_hdf5.c_str());
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
#else
|
||||
unlink(fname_hdf5.c_str());
|
||||
#endif
|
||||
|
||||
phi.Write_to_HDF5(fname_hdf5, "phi");
|
||||
phi2.Write_to_HDF5(fname_hdf5, "phi2");
|
||||
phi3a.Write_to_HDF5(fname_hdf5, "phi3a");
|
||||
phi3b.Write_to_HDF5(fname_hdf5, "phi3b");
|
||||
|
||||
delta.Write_to_HDF5(fname_hdf5, "delta");
|
||||
delta2.Write_to_HDF5(fname_hdf5, "delta2");
|
||||
delta3a.Write_to_HDF5(fname_hdf5, "delta3a");
|
||||
delta3b.Write_to_HDF5(fname_hdf5, "delta3b");
|
||||
delta3.Write_to_HDF5(fname_hdf5, "delta3");
|
||||
|
||||
A3[0]->Write_to_HDF5(fname_hdf5, "A3x");
|
||||
A3[1]->Write_to_HDF5(fname_hdf5, "A3y");
|
||||
A3[2]->Write_to_HDF5(fname_hdf5, "A3z");
|
||||
}
|
30
testing/RAMSES/ics_ramses.conf
Normal file
30
testing/RAMSES/ics_ramses.conf
Normal file
|
@ -0,0 +1,30 @@
|
|||
[setup]
|
||||
GridRes = 64
|
||||
BoxLength = 100
|
||||
zstart = 49.0
|
||||
LPTorder = 1
|
||||
SymplecticPT = no
|
||||
DoFixing = no
|
||||
BCClattice = no
|
||||
|
||||
[execution]
|
||||
NumThreads = 4
|
||||
|
||||
[output]
|
||||
format = grafic2
|
||||
filename = ics_ramses
|
||||
grafic_use_SPT = yes # set to no and LPTorder to 1 or 2 for semiclassical ICs
|
||||
|
||||
[random]
|
||||
generator = NGENIC
|
||||
seed = 9001
|
||||
|
||||
[cosmology]
|
||||
transfer = eisenstein
|
||||
Omega_m = 0.302
|
||||
Omega_b = 0.045
|
||||
Omega_L = 0.698
|
||||
H0 = 70.3
|
||||
sigma_8 = 0.811
|
||||
nspec = 0.961
|
||||
|
60
testing/RAMSES/ramses.nml
Normal file
60
testing/RAMSES/ramses.nml
Normal file
|
@ -0,0 +1,60 @@
|
|||
&RUN_PARAMS
|
||||
cosmo=.true.
|
||||
pic=.true.
|
||||
poisson=.true.
|
||||
hydro=.true.
|
||||
nrestart=0
|
||||
nremap=1
|
||||
nsubcycle=1,2
|
||||
ncontrol=1
|
||||
verbose=.false.
|
||||
/
|
||||
|
||||
&OUTPUT_PARAMS
|
||||
noutput=10
|
||||
aout=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0
|
||||
/
|
||||
|
||||
&INIT_PARAMS
|
||||
filetype='grafic'
|
||||
initfile(1)='ics_ramses'
|
||||
/
|
||||
|
||||
&AMR_PARAMS
|
||||
levelmin=6
|
||||
levelmax=11
|
||||
nexpand=1
|
||||
ngridmax=2000000
|
||||
npartmax=3000000
|
||||
/
|
||||
|
||||
&HYDRO_PARAMS
|
||||
gamma=1.6666667
|
||||
courant_factor=0.8
|
||||
slope_type=1
|
||||
pressure_fix=.true.
|
||||
scheme='muscl'
|
||||
/
|
||||
|
||||
&COOLING_PARAMS
|
||||
cooling=.false.
|
||||
metal=.false.
|
||||
/
|
||||
|
||||
&SF_PARAMS
|
||||
t_star=8.0
|
||||
n_star=0.1
|
||||
T2_star=1d4
|
||||
/
|
||||
|
||||
&FEEDBACK_PARAMS
|
||||
eta_sn=0.1
|
||||
yield=0.1
|
||||
/
|
||||
|
||||
&REFINE_PARAMS
|
||||
m_refine=7*8.
|
||||
interpol_var=1
|
||||
interpol_type=0
|
||||
/
|
||||
|
51
version.cmake
Normal file
51
version.cmake
Normal file
|
@ -0,0 +1,51 @@
|
|||
# CMake script to get git repository version at compile time
|
||||
# taken from Matt Keeter's blog and slightly adapted
|
||||
# https://www.mattkeeter.com/blog/2018-01-06-versioning/
|
||||
|
||||
execute_process(COMMAND git log --pretty=format:'%h' -n 1
|
||||
OUTPUT_VARIABLE GIT_REV
|
||||
ERROR_QUIET)
|
||||
|
||||
# Check whether we got any revision (which isn't
|
||||
# always the case, e.g. when someone downloaded a zip
|
||||
# file from Github instead of a checkout)
|
||||
if ("${GIT_REV}" STREQUAL "")
|
||||
set(GIT_REV "N/A")
|
||||
set(GIT_DIFF "")
|
||||
set(GIT_TAG "N/A")
|
||||
set(GIT_BRANCH "N/A")
|
||||
else()
|
||||
execute_process(
|
||||
COMMAND bash -c "git diff --quiet --exit-code || echo +"
|
||||
OUTPUT_VARIABLE GIT_DIFF)
|
||||
execute_process(
|
||||
COMMAND git describe --exact-match --tags
|
||||
OUTPUT_VARIABLE GIT_TAG ERROR_QUIET)
|
||||
execute_process(
|
||||
COMMAND git rev-parse --abbrev-ref HEAD
|
||||
OUTPUT_VARIABLE GIT_BRANCH)
|
||||
|
||||
string(STRIP "${GIT_REV}" GIT_REV)
|
||||
string(SUBSTRING "${GIT_REV}" 1 7 GIT_REV)
|
||||
string(STRIP "${GIT_DIFF}" GIT_DIFF)
|
||||
string(STRIP "${GIT_TAG}" GIT_TAG)
|
||||
string(STRIP "${GIT_BRANCH}" GIT_BRANCH)
|
||||
endif()
|
||||
|
||||
if("${GIT_TAG}" STREQUAL "")
|
||||
set(GIT_TAG "N/A")
|
||||
endif()
|
||||
|
||||
set(VERSION "const char* GIT_REV=\"${GIT_REV}${GIT_DIFF}\";
|
||||
const char* GIT_TAG=\"${GIT_TAG}\";
|
||||
const char* GIT_BRANCH=\"${GIT_BRANCH}\";")
|
||||
|
||||
if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/version.cc)
|
||||
file(READ ${CMAKE_CURRENT_SOURCE_DIR}/version.cc VERSION_)
|
||||
else()
|
||||
set(VERSION_ "")
|
||||
endif()
|
||||
|
||||
if (NOT "${VERSION}" STREQUAL "${VERSION_}")
|
||||
file(WRITE ${CMAKE_CURRENT_SOURCE_DIR}/version.cc "${VERSION}")
|
||||
endif()
|
Loading…
Reference in a new issue