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

introduced some stl arrays

Browse source
This commit is contained in:
Oliver Hahn 2023-02-04 16:24:33 -08:00
parent 5bd8c45341
commit 33d784a137
9 changed files with 1485 additions and 644 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/general.hh
View file

@ -54,6 +54,8 @@
typedef double real_t;
#endif
#include <array>
using vec3_t = std::array<real_t,3>;
#ifdef FFTW3
#define RE(x) ((x)[0])

35
src/plugins/convex_hull.hh
View file

@ -11,6 +11,9 @@
#include "log.hh"
#include <array>
/***** Slow but short convex hull Implementation ******/
/*
* Finds the convex hull of a set of data points.
@ -25,14 +28,16 @@
template< typename real_t >
struct convex_hull{
typedef const real_t *cpr_;
using vec3_t = std::array<real_t,3>;
size_t npoints_;
std::vector<int> faceidx_L_, faceidx_U_;
std::vector<real_t> normals_L_, normals_U_;
std::vector<real_t> x0_L_, x0_U_;
real_t centroid_[3], volume_;
real_t left_[3], right_[3];
real_t anchor_pt_[3];
vec3_t centroid_;
real_t volume_;
vec3_t left_, right_;
vec3_t anchor_pt_;
inline double turn( cpr_ p, cpr_ q, cpr_ r ) const
{ return (q[0]-p[0])*(r[1]-p[1]) - (r[0]-p[0])*(q[1]-p[1]); }
@ -224,32 +229,33 @@ struct convex_hull{
dist *= -1.0;
// take care of possible periodic boundaries
T x[3];
std::array<T,3> x;
for( size_t p=0; p<3; ++p )
{
T d = xp[p] - anchor_pt_[p];
real_t d = xp[p] - anchor_pt_[p];
if( d>0.5 ) x[p] = xp[p]-1.0; else if ( d<-0.5 ) x[p] = xp[p]+1.0; else x[p] = xp[p];
}
// check for inside vs. outside
for( size_t i=0; i<normals_L_.size()/3; ++i )
{
double xp[3] = {x[0]-x0_L_[3*i+0],x[1]-x0_L_[3*i+1],x[2]-x0_L_[3*i+2]};
if( dot( xp, &normals_L_[3*i]) < dist ) return false;
std::array<T,3> xp{{x[0]-x0_L_[3*i+0],x[1]-x0_L_[3*i+1],x[2]-x0_L_[3*i+2]}};
if( dot( &xp[0], &normals_L_[3*i]) < dist ) return false;
}
for( size_t i=0; i<normals_U_.size()/3; ++i )
{
double xp[3] = {x[0]-x0_U_[3*i+0],x[1]-x0_U_[3*i+1],x[2]-x0_U_[3*i+2]};
if( dot( xp, &normals_U_[3*i]) < dist ) return false;
std::array<T,3> xp{{x[0]-x0_U_[3*i+0],x[1]-x0_U_[3*i+1],x[2]-x0_U_[3*i+2]}};
if( dot( &xp[0], &normals_U_[3*i]) < dist ) return false;
}
return true;
}
void expand_vector_from_centroid( real_t *v, double dr )
void expand_vector_from_centroid( real_t* v, double dr )
{
double dx[3], d = 0.0;
vec3_t dx;
real_t d = 0.0;
for( int i=0; i<3; ++i )
{
dx[i] = v[i]-centroid_[i];
@ -267,9 +273,8 @@ struct convex_hull{
for( size_t i=0; i<normals_U_.size(); i+=3 )
expand_vector_from_centroid( &x0_U_[i], dr );
expand_vector_from_centroid( left_, dr );
expand_vector_from_centroid( right_, dr );
expand_vector_from_centroid( &left_[0], dr );
expand_vector_from_centroid( &right_[0], dr );
}
void get_defining_indices( std::set<int>& unique ) const
@ -282,7 +287,7 @@ struct convex_hull{
unique.insert( faceidx_U_[i] );
}
convex_hull( cpr_ points, size_t npoints, const double* anchor )
convex_hull( cpr_ points, size_t npoints, const vec3_t& anchor )
: npoints_( npoints )
{
anchor_pt_[0] = anchor[0];

2
src/plugins/output_arepo.cc
View file

@ -563,7 +563,7 @@ public:
if (levelmin_ != levelmax_)
{
double lxref[3], x0ref[3], x1ref[3];
vec3_t lxref, x0ref, x1ref;
double pmgrid_new;
the_region_generator->get_AABB(x0ref, x1ref, levelmax_); // generalized beyond box

809
src/plugins/output_swift.cc Normal file
View file

@ -0,0 +1,809 @@
/*
* output_arepo.cc - This file is part of MUSIC -
* a code to generate multi-scale initial conditions
* for cosmological simulations
*
* Copyright (C) 2010 Oliver Hahn
*
* Plugin: Dylan Nelson (dnelson@cfa.harvard.edu)
*/
#ifdef HAVE_HDF5
#define GAS_PARTTYPE 0
#define HIGHRES_DM_PARTTYPE 1
#define COARSE_DM_DEFAULT_PARTTYPE 2
#define STAR_PARTTYPE 4
#define NTYPES 6
#include <sstream>
#include <string>
#include <algorithm>
#include "output.hh"
#include "HDF_IO.hh"
class arepo_output_plugin : public output_plugin
{
protected:
// header/config
std::vector<std::vector<unsigned int>> nPart;
std::vector<long long> nPartTotal;
std::vector<double> massTable;
double time, redshift, boxSize;
unsigned int numFiles, coarsePartType;
double omega0, omega_L, hubbleParam;
// configuration
double UnitLength_in_cm, UnitMass_in_g, UnitVelocity_in_cm_per_s;
double omega_b, rhoCrit;
double posFac, velFac;
long long nPartTotAllTypes;
bool doBaryons, useLongIDs, doublePrec;
size_t npfine, npart, npcoarse;
std::vector<size_t> levelcounts;
// parameter file hints
int pmgrid, gridboost;
float softening, Tini;
using output_plugin::cf_;
// Nx1 vector (e.g. masses,particleids)
template <typename T>
void writeHDF5_a(std::string fieldName, int partTypeNum, const std::vector<T> &data)
{
hid_t HDF_FileID, HDF_GroupID, HDF_DatasetID, HDF_DataspaceID, HDF_Type;
hsize_t HDF_Dims, offset = 0;
std::stringstream GrpName;
GrpName << "PartType" << partTypeNum;
for (unsigned i = 0; i < numFiles; i++)
{
std::string filename = fname_;
HDF_Dims = data.size();
// modify local filename and write size
if (numFiles > 1)
{
std::stringstream s;
s << "." << i << ".hdf5";
filename.replace(filename.find(".hdf5"), 5, s.str());
HDF_Dims = ceil(data.size() / numFiles);
if (i == numFiles - 1)
HDF_Dims = data.size() - offset;
}
HDF_FileID = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
HDF_GroupID = H5Gopen(HDF_FileID, GrpName.str().c_str());
HDF_Type = GetDataType<T>();
HDF_DataspaceID = H5Screate_simple(1, &HDF_Dims, NULL);
HDF_DatasetID = H5Dcreate(HDF_GroupID, fieldName.c_str(), HDF_Type, HDF_DataspaceID, H5P_DEFAULT);
// write and close
H5Dwrite(HDF_DatasetID, HDF_Type, H5S_ALL, H5S_ALL, H5P_DEFAULT, &data[offset]);
H5Dclose(HDF_DatasetID);
H5Sclose(HDF_DataspaceID);
H5Gclose(HDF_GroupID);
H5Fclose(HDF_FileID);
offset += HDF_Dims;
}
}
// Nx3 vector (e.g. pos,vel), where coord = index of the second dimension (written one at a time)
template <typename T>
void writeHDF5_b(std::string fieldName, int coord, int partTypeNum, std::vector<T> &data, bool readFlag = false)
{
hid_t HDF_FileID, HDF_GroupID, HDF_DatasetID, HDF_DataspaceID, HDF_Type;
hsize_t HDF_Dims[2], HDF_DimsMem[2], w_offset = 0;
std::stringstream GrpName;
GrpName << "PartType" << partTypeNum;
for (unsigned i = 0; i < numFiles; i++)
{
std::string filename = fname_;
HDF_Dims[0] = data.size();
// modify local filename and write size
if (numFiles > 1)
{
std::stringstream s;
s << "." << i << ".hdf5";
filename.replace(filename.find(".hdf5"), 5, s.str());
HDF_Dims[0] = ceil(data.size() / numFiles);
if (i == numFiles - 1)
HDF_Dims[0] = data.size() - w_offset;
}
HDF_FileID = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
HDF_GroupID = H5Gopen(HDF_FileID, GrpName.str().c_str());
HDF_Type = GetDataType<T>();
HDF_Dims[1] = 3;
// if dataset does not yet exist, create it (on first coord call)
if (!(H5Lexists(HDF_GroupID, fieldName.c_str(), H5P_DEFAULT)))
{
HDF_DataspaceID = H5Screate_simple(2, HDF_Dims, NULL);
HDF_DatasetID = H5Dcreate(HDF_GroupID, fieldName.c_str(), HDF_Type, HDF_DataspaceID, H5P_DEFAULT);
H5Sclose(HDF_DataspaceID);
H5Dclose(HDF_DatasetID);
}
// make memory space (just indicates the size/shape of data)
HDF_DimsMem[0] = HDF_Dims[0];
HDF_DimsMem[1] = 1;
hid_t HDF_MemoryspaceID = H5Screate_simple(2, HDF_DimsMem, NULL);
// open hyperslab
hsize_t count[2] = {1, 1}, stride[2] = {1, 1}, offset[2] = {0, 0};
offset[1] = coord; // set where in the second dimension to write
count[0] = HDF_Dims[0]; // set size in the first dimension (num particles of this type)
HDF_DatasetID = H5Dopen(HDF_GroupID, fieldName.c_str());
HDF_DataspaceID = H5Dget_space(HDF_DatasetID);
H5Sselect_hyperslab(HDF_DataspaceID, H5S_SELECT_SET, offset, stride, count, NULL);
// write (or read) and close
if (readFlag)
H5Dread(HDF_DatasetID, HDF_Type, HDF_MemoryspaceID, HDF_DataspaceID, H5P_DEFAULT,
&data[w_offset]);
else
H5Dwrite(HDF_DatasetID, HDF_Type, HDF_MemoryspaceID, HDF_DataspaceID, H5P_DEFAULT,
&data[w_offset]);
H5Dclose(HDF_DatasetID);
H5Gclose(HDF_GroupID);
H5Fclose(HDF_FileID);
w_offset += HDF_Dims[0];
}
}
// called from finalize()
void generateAndWriteIDs(void)
{
long long offset = 1; // don't use ID==0
nPartTotAllTypes = 0;
for (size_t i = 0; i < nPartTotal.size(); i++)
{
if (!nPartTotal[i])
continue;
nPartTotAllTypes += nPartTotal[i];
if (!useLongIDs)
{
std::vector<int> ids = std::vector<int>(nPartTotal[i]);
for (int j = 0; j < nPartTotal[i]; j++)
ids[j] = offset + j;
writeHDF5_a("ParticleIDs", i, ids);
}
else
{
std::vector<long long> ids = std::vector<long long>(nPartTotal[i]);
for (long long j = 0; j < nPartTotal[i]; j++)
ids[j] = offset + j;
writeHDF5_a("ParticleIDs", i, ids);
}
// make IDs of all particle types sequential (unique) = unnecessary, but consistent with gadget output format
offset += nPartTotal[i];
}
}
void countLeafCells(const grid_hierarchy &gh)
{
npfine = 0;
npart = 0;
npcoarse = 0;
npfine = gh.count_leaf_cells(gh.levelmax(), gh.levelmax());
npart = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());
if (levelmax_ != levelmin_) // multimass
npcoarse = gh.count_leaf_cells(gh.levelmin(), gh.levelmax() - 1);
}
template <typename T>
void __write_dm_mass(const grid_hierarchy &gh)
{
countLeafCells(gh);
// fill levelcount for header
levelcounts = std::vector<size_t>(levelmax_ - levelmin_ + 1);
for (int ilevel = gh.levelmax(); ilevel >= (int)gh.levelmin(); --ilevel)
levelcounts[gh.levelmax() - ilevel] = gh.count_leaf_cells(ilevel, ilevel);
if (levelmax_ > levelmin_ + 1) // morethan2bnd
{
// DM particles will have variable masses
size_t count = 0;
std::vector<T> data(npcoarse);
for (int ilevel = gh.levelmax() - 1; ilevel >= (int)gh.levelmin(); --ilevel)
{
// baryon particles live only on finest grid, these particles here are total matter particles
T pmass = omega0 * rhoCrit * pow(boxSize * posFac, 3.0) / pow(2, 3 * ilevel);
for (unsigned i = 0; i < gh.get_grid(ilevel)->size(0); ++i)
for (unsigned j = 0; j < gh.get_grid(ilevel)->size(1); ++j)
for (unsigned k = 0; k < gh.get_grid(ilevel)->size(2); ++k)
if (gh.is_in_mask(ilevel, i, j, k) && !gh.is_refined(ilevel, i, j, k))
{
data[count++] = pmass;
}
}
if (count != npcoarse)
throw std::runtime_error("Internal consistency error while writing masses");
writeHDF5_a("Masses", coarsePartType, data); // write DM
}
else
{
// DM particles will all have the same mass, just write to massTable
if (levelmax_ != levelmin_) // multimass
massTable[coarsePartType] = omega0 * rhoCrit * pow(boxSize * posFac, 3.0) / pow(2, 3 * levelmin_);
}
}
template <typename T>
void __write_dm_position(int coord, const grid_hierarchy &gh)
{
countLeafCells(gh);
// update header
hsize_t offset_fine = 0, offset_coarse = 0;
for (unsigned i = 0; i < numFiles; i++)
{
hsize_t dims_fine = ceil(npfine / numFiles);
hsize_t dims_coarse = ceil(npcoarse / numFiles);
if (i == numFiles - 1)
{
dims_fine = npfine - offset_fine;
dims_coarse = npcoarse - offset_coarse;
}
nPart[i][HIGHRES_DM_PARTTYPE] = dims_fine;
nPart[i][coarsePartType] = dims_coarse;
offset_fine += dims_fine;
offset_coarse += dims_coarse;
}
nPartTotal[HIGHRES_DM_PARTTYPE] = npfine;
nPartTotal[coarsePartType] = npcoarse;
// FINE: collect displacements and convert to absolute coordinates with correct units
int ilevel = gh.levelmax();
std::vector<T> data(npfine);
size_t count = 0;
for (unsigned i = 0; i < gh.get_grid(ilevel)->size(0); ++i)
for (unsigned j = 0; j < gh.get_grid(ilevel)->size(1); ++j)
for (unsigned k = 0; k < gh.get_grid(ilevel)->size(2); ++k)
if (gh.is_in_mask(ilevel, i, j, k) && !gh.is_refined(ilevel, i, j, k))
{
double xx[3];
gh.cell_pos(ilevel, i, j, k, xx);
xx[coord] = (xx[coord] + (*gh.get_grid(ilevel))(i, j, k)) * boxSize;
xx[coord] = fmod(xx[coord] + boxSize, boxSize);
data[count++] = (T)(xx[coord] * posFac);
}
writeHDF5_b("Coordinates", coord, HIGHRES_DM_PARTTYPE, data); // write fine DM
if (count != npfine)
throw std::runtime_error("Internal consistency error while writing fine DM pos");
// COARSE: collect displacements and convert to absolute coordinates with correct units
if (levelmax_ != levelmin_) // multimass
{
data = std::vector<T>(npcoarse, 0.0);
count = 0;
for (int ilevel = gh.levelmax() - 1; ilevel >= (int)gh.levelmin(); --ilevel)
for (unsigned i = 0; i < gh.get_grid(ilevel)->size(0); ++i)
for (unsigned j = 0; j < gh.get_grid(ilevel)->size(1); ++j)
for (unsigned k = 0; k < gh.get_grid(ilevel)->size(2); ++k)
if (gh.is_in_mask(ilevel, i, j, k) && !gh.is_refined(ilevel, i, j, k))
{
double xx[3];
gh.cell_pos(ilevel, i, j, k, xx);
xx[coord] = (xx[coord] + (*gh.get_grid(ilevel))(i, j, k)) * boxSize;
if (!doBaryons) // if so, we will handle the mod in write_gas_position
xx[coord] = fmod(xx[coord] + boxSize, boxSize) * posFac;
data[count++] = (T)xx[coord];
}
if (count != npcoarse)
throw std::runtime_error("Internal consistency error while writing coarse DM pos");
writeHDF5_b("Coordinates", coord, coarsePartType, data); // write coarse DM
}
}
template <typename T>
void __write_dm_velocity(int coord, const grid_hierarchy &gh)
{
countLeafCells(gh);
// FINE: collect velocities and convert to correct units
int ilevel = gh.levelmax();
std::vector<T> data(npfine);
size_t count = 0;
for (unsigned i = 0; i < gh.get_grid(ilevel)->size(0); ++i)
for (unsigned j = 0; j < gh.get_grid(ilevel)->size(1); ++j)
for (unsigned k = 0; k < gh.get_grid(ilevel)->size(2); ++k)
if (gh.is_in_mask(ilevel, i, j, k) && !gh.is_refined(ilevel, i, j, k))
{
data[count++] = (T)(*gh.get_grid(ilevel))(i, j, k) * velFac;
}
writeHDF5_b("Velocities", coord, HIGHRES_DM_PARTTYPE, data); // write fine DM
if (count != npfine)
throw std::runtime_error("Internal consistency error while writing fine DM pos");
// COARSE: collect velocities and convert to correct units
if (levelmax_ != levelmin_) // multimass
{
data = std::vector<T>(npcoarse, 0.0);
count = 0;
for (int ilevel = gh.levelmax() - 1; ilevel >= (int)gh.levelmin(); --ilevel)
for (unsigned i = 0; i < gh.get_grid(ilevel)->size(0); ++i)
for (unsigned j = 0; j < gh.get_grid(ilevel)->size(1); ++j)
for (unsigned k = 0; k < gh.get_grid(ilevel)->size(2); ++k)
if (gh.is_in_mask(ilevel, i, j, k) && !gh.is_refined(ilevel, i, j, k))
{
data[count++] = (T)(*gh.get_grid(ilevel))(i, j, k) * velFac;
}
if (count != npcoarse)
throw std::runtime_error("Internal consistency error while writing coarse DM pos");
writeHDF5_b("Velocities", coord, coarsePartType, data); // write coarse DM
}
}
template <typename T>
void __write_gas_velocity(int coord, const grid_hierarchy &gh)
{
countLeafCells(gh);
std::vector<T> gas_data(npart); // read/write gas at all levels from the gh
size_t count = 0;
for (int ilevel = levelmax_; ilevel >= (int)levelmin_; --ilevel)
for (unsigned i = 0; i < gh.get_grid(ilevel)->size(0); ++i)
for (unsigned j = 0; j < gh.get_grid(ilevel)->size(1); ++j)
for (unsigned k = 0; k < gh.get_grid(ilevel)->size(2); ++k)
if (gh.is_in_mask(ilevel, i, j, k) && !gh.is_refined(ilevel, i, j, k))
{
gas_data[count++] = (T)(*gh.get_grid(ilevel))(i, j, k) * velFac;
}
if (count != npart)
throw std::runtime_error("Internal consistency error while writing GAS pos");
// calculate modified DM velocities if: multimass and baryons present
if (doBaryons && npcoarse)
{
double facb = omega_b / omega0;
double facc = (omega0 - omega_b) / omega0;
std::vector<T> dm_data(npcoarse);
writeHDF5_b("Velocities", coord, coarsePartType, dm_data, true); // read coarse DM vels
// overwrite
for (size_t i = 0; i < npcoarse; i++)
dm_data[i] = facc * dm_data[i] + facb * gas_data[npfine + i];
writeHDF5_b("Velocities", coord, coarsePartType, dm_data); // overwrite coarse DM vels
} // dm_data deallocated
// restrict gas_data to fine only and request write
std::vector<T> data(gas_data.begin() + 0, gas_data.begin() + npfine);
std::vector<T>().swap(gas_data); // deallocate
writeHDF5_b("Velocities", coord, GAS_PARTTYPE, data); // write highres gas
}
template <typename T>
void __write_gas_position(int coord, const grid_hierarchy &gh)
{
countLeafCells(gh);
// update header (will actually write only gas at levelmax)
hsize_t offset = 0;
for (unsigned i = 0; i < numFiles; i++)
{
hsize_t dims = ceil(npfine / numFiles);
if (i == numFiles - 1)
dims = npfine - offset;
nPart[i][GAS_PARTTYPE] = dims;
offset += dims;
}
nPartTotal[GAS_PARTTYPE] = npfine;
std::vector<double> gas_data(npart); // read/write gas at all levels from the gh
size_t count = 0;
double h = 1.0 / (1ul << gh.levelmax());
for (int ilevel = gh.levelmax(); ilevel >= (int)gh.levelmin(); --ilevel)
for (unsigned i = 0; i < gh.get_grid(ilevel)->size(0); ++i)
for (unsigned j = 0; j < gh.get_grid(ilevel)->size(1); ++j)
for (unsigned k = 0; k < gh.get_grid(ilevel)->size(2); ++k)
if (gh.is_in_mask(ilevel, i, j, k) && !gh.is_refined(ilevel, i, j, k))
{
double xx[3];
gh.cell_pos(ilevel, i, j, k, xx);
// shift particle positions (this has to be done as the same shift
// is used when computing the convolution kernel for SPH baryons)
xx[coord] += 0.5 * h;
xx[coord] = (xx[coord] + (*gh.get_grid(ilevel))(i, j, k)) * boxSize;
gas_data[count++] = xx[coord];
}
if (count != npart)
throw std::runtime_error("Internal consistency error while writing coarse DM pos");
// calculate modified DM coordinates if: multimass and baryons present
if (doBaryons && npcoarse)
{
double facb = omega_b / omega0;
double facc = (omega0 - omega_b) / omega0;
std::vector<T> dm_data(npcoarse);
writeHDF5_b("Coordinates", coord, coarsePartType, dm_data, true); // read coarse DM vels
// overwrite
for (size_t i = 0; i < npcoarse; i++)
{
dm_data[i] = facc * dm_data[i] + facb * gas_data[npfine + i];
dm_data[i] = fmod(dm_data[i] + boxSize, boxSize) * posFac;
}
writeHDF5_b("Coordinates", coord, coarsePartType, dm_data); // overwrite coarse DM vels
}
// restrict gas_data to fine only and request write
//std::vector<float> data( gas_data.begin() + 0, gas_data.begin() + npfine );
std::vector<T> data(npfine);
for (size_t i = 0; i < npfine; i++)
data[i] = (T)(fmod(gas_data[i] + boxSize, boxSize) * posFac);
std::vector<double>().swap(gas_data); // deallocate
writeHDF5_b("Coordinates", coord, GAS_PARTTYPE, data); // write highres gas
}
public:
arepo_output_plugin(config_file &cf) : output_plugin(cf)
{
// ensure that everyone knows we want to do SPH, implies: bsph=1, bbshift=1, decic_baryons=1
// -> instead of just writing gas densities (which are here ignored), the gas displacements are also written
cf.insertValue("setup", "do_SPH", "yes");
// init header and config parameters
nPartTotal = std::vector<long long>(NTYPES, 0);
massTable = std::vector<double>(NTYPES, 0.0);
coarsePartType = cf.getValueSafe<unsigned>("output", "arepo_coarsetype", COARSE_DM_DEFAULT_PARTTYPE);
UnitLength_in_cm = cf.getValueSafe<double>("output", "arepo_unitlength", 3.085678e21); // 1.0 kpc
UnitMass_in_g = cf.getValueSafe<double>("output", "arepo_unitmass", 1.989e43); // 1.0e10 solar masses
UnitVelocity_in_cm_per_s = cf.getValueSafe<double>("output", "arepo_unitvel", 1e5); // 1 km/sec
omega0 = cf.getValue<double>("cosmology", "Omega_m");
omega_b = cf.getValue<double>("cosmology", "Omega_b");
omega_L = cf.getValue<double>("cosmology", "Omega_L");
redshift = cf.getValue<double>("setup", "zstart");
boxSize = cf.getValue<double>("setup", "boxlength");
doBaryons = cf.getValueSafe<bool>("setup", "baryons", false);
useLongIDs = cf.getValueSafe<bool>("output", "arepo_longids", false);
numFiles = cf.getValueSafe<unsigned>("output", "arepo_num_files", 1);
doublePrec = cf.getValueSafe<bool>("output", "arepo_doubleprec", 0);
for (unsigned i = 0; i < numFiles; i++)
nPart.push_back(std::vector<unsigned int>(NTYPES, 0));
// factors which multiply positions and velocities
time = 1.0 / (1.0 + redshift);
posFac = 3.085678e24 / UnitLength_in_cm; // MUSIC uses Mpc internally, i.e. posFac=1e3 for kpc output
velFac = (1.0f / sqrt(time)) * boxSize;
// critical density
rhoCrit = 27.7519737e-9; // in h^2 1e10 M_sol / kpc^3
rhoCrit *= pow(UnitLength_in_cm / 3.085678e21, 3.0);
rhoCrit *= (1.989e43 / UnitMass_in_g);
// calculate PMGRID suggestion
pmgrid = pow(2, levelmin_) * 2; // unigrid
gridboost = 1;
if (levelmin_ != levelmax_)
{
vec3_t lxref, x0ref, x1ref;
double pmgrid_new;
the_region_generator->get_AABB(x0ref, x1ref, levelmax_); // generalized beyond box
for (int i = 0; i < 3; i++)
lxref[i] = x1ref[i] - x0ref[i];
// fraction box length of the zoom region
lxref[0] = pow((lxref[0] * lxref[1] * lxref[2]), 0.333);
pmgrid_new = pow(2, levelmax_) * 2; // to cover entire box at highest resolution
pmgrid_new *= lxref[0]; // only need to cover a fraction
if ((gridboost = round(pmgrid_new / pmgrid)) > 1)
gridboost = pow(2, ceil(log(gridboost) / log(2.0))); // round to nearest, higher power of 2
if (gridboost == 0)
gridboost = 1;
}
// calculate Tini for gas
hubbleParam = cf.getValue<double>("cosmology", "H0") / 100.0;
double astart = 1.0 / (1.0 + redshift);
double h2 = hubbleParam * hubbleParam;
double adec = 1.0 / (160.0 * pow(omega_b * h2 / 0.022, 2.0 / 5.0));
double Tcmb0 = 2.726;
Tini = astart < adec ? Tcmb0 / astart : Tcmb0 / astart / astart * adec;
// calculate softening suggestion
softening = (boxSize * posFac) / pow(2, levelmax_) / 40.0;
// header and sanity checks
if (!doBaryons)
massTable[HIGHRES_DM_PARTTYPE] = omega0 * rhoCrit * pow(boxSize * posFac, 3.0) / pow(2, 3 * levelmax_);
else
massTable[HIGHRES_DM_PARTTYPE] = (omega0 - omega_b) * rhoCrit * pow(boxSize * posFac, 3.0) / pow(2, 3 * levelmax_);
if (coarsePartType == GAS_PARTTYPE || coarsePartType == HIGHRES_DM_PARTTYPE)
throw std::runtime_error("Error: Specified illegal Arepo particle type for coarse particles.");
if (coarsePartType == STAR_PARTTYPE)
LOGWARN("WARNING: Specified coarse particle type will collide with stars if USE_SFR enabled.");
// create file(s)
for (unsigned i = 0; i < numFiles; i++)
{
std::string filename = fname_;
if (numFiles > 1)
{
size_t pos = filename.find(".hdf5");
if (pos != filename.length() - 5)
throw std::runtime_error("Error: Unexpected output filename (doesn't end in .hdf5).");
std::stringstream s;
s << "." << i << ".hdf5";
filename.replace(pos, 5, s.str());
}
HDFCreateFile(filename);
// create particle type groups
std::stringstream GrpName;
GrpName << "PartType" << HIGHRES_DM_PARTTYPE;
HDFCreateGroup(filename, GrpName.str().c_str()); // highres or unigrid DM
if (doBaryons)
{
GrpName.str("");
GrpName << "PartType" << GAS_PARTTYPE;
HDFCreateGroup(filename, GrpName.str().c_str()); // gas
}
if (levelmax_ != levelmin_) // multimass
{
GrpName.str("");
GrpName << "PartType" << coarsePartType;
HDFCreateGroup(filename, GrpName.str().c_str()); // coarse DM
}
}
}
~arepo_output_plugin()
{
}
/* ------------------------------------------------------------------------------- */
void write_dm_mass(const grid_hierarchy &gh)
{
if (!doublePrec)
__write_dm_mass<float>(gh);
else
__write_dm_mass<double>(gh);
}
void write_dm_position(int coord, const grid_hierarchy &gh)
{
if (!doublePrec)
__write_dm_position<float>(coord, gh);
else
__write_dm_position<double>(coord, gh);
}
void write_dm_velocity(int coord, const grid_hierarchy &gh)
{
if (!doublePrec)
__write_dm_velocity<float>(coord, gh);
else
__write_dm_velocity<double>(coord, gh);
}
void write_dm_density(const grid_hierarchy &gh)
{ /* skip */
}
void write_dm_potential(const grid_hierarchy &gh)
{ /* skip */
}
/* ------------------------------------------------------------------------------- */
void write_gas_velocity(int coord, const grid_hierarchy &gh)
{
if (!doublePrec)
__write_gas_velocity<float>(coord, gh);
else
__write_gas_velocity<double>(coord, gh);
}
void write_gas_position(int coord, const grid_hierarchy &gh)
{
if (!doublePrec)
__write_gas_position<float>(coord, gh);
else
__write_gas_position<double>(coord, gh);
}
void write_gas_density(const grid_hierarchy &gh)
{
// if only saving highres gas, then all gas cells have the same initial mass
// do not write out densities as we write out displacements
if (doBaryons)
massTable[GAS_PARTTYPE] = omega_b * rhoCrit * pow(boxSize * posFac, 3.0) / pow(2, 3 * levelmax_);
}
void write_gas_potential(const grid_hierarchy &gh)
{ /* skip */
}
void finalize(void)
{
// generate and add contiguous IDs for each particle type we have written
generateAndWriteIDs();
std::vector<unsigned int> nPartTotalLW(nPartTotal.size());
std::vector<unsigned int> nPartTotalHW(nPartTotal.size());
for (size_t i = 0; i < nPartTotalHW.size(); i++)
{
nPartTotalHW[i] = (unsigned)((size_t)nPartTotal[i] >> 32);
nPartTotalLW[i] = (unsigned int)((size_t)nPartTotal[i]);
}
// output particle counts
std::cout << " - Arepo : wrote " << nPartTotAllTypes << " particles..." << std::endl;
for (size_t i = 0; i < nPartTotal.size(); i++)
std::cout << " type [" << i << "] : " << std::setw(12) << nPartTotal[i] << std::endl;
std::cout << std::endl;
// write final header (some of these fields are required, others are extra info)
for (unsigned i = 0; i < numFiles; i++)
{
std::string filename = fname_;
if (numFiles > 1)
{
std::stringstream s;
s << "." << i << ".hdf5";
filename.replace(filename.find(".hdf5"), 5, s.str());
std::cout << " " << filename;
for (size_t j = 0; j < nPart[i].size(); j++)
std::cout << " " << std::setw(10) << nPart[i][j];
std::cout << std::endl;
}
HDFCreateGroup(filename, "Header");
HDFWriteGroupAttribute(filename, "Header", "NumPart_ThisFile", nPart[i]);
HDFWriteGroupAttribute(filename, "Header", "NumPart_Total", nPartTotalLW);
HDFWriteGroupAttribute(filename, "Header", "NumPart_Total_HighWord", nPartTotalHW);
HDFWriteGroupAttribute(filename, "Header", "MassTable", massTable);
HDFWriteGroupAttribute(filename, "Header", "BoxSize", boxSize * posFac);
HDFWriteGroupAttribute(filename, "Header", "NumFilesPerSnapshot", numFiles);
HDFWriteGroupAttribute(filename, "Header", "Time", time);
HDFWriteGroupAttribute(filename, "Header", "Redshift", redshift);
HDFWriteGroupAttribute(filename, "Header", "Omega0", omega0);
HDFWriteGroupAttribute(filename, "Header", "OmegaLambda", omega_L);
HDFWriteGroupAttribute(filename, "Header", "OmegaBaryon", omega_b);
HDFWriteGroupAttribute(filename, "Header", "HubbleParam", hubbleParam);
HDFWriteGroupAttribute(filename, "Header", "Flag_Sfr", 0);
HDFWriteGroupAttribute(filename, "Header", "Flag_Cooling", 0);
HDFWriteGroupAttribute(filename, "Header", "Flag_StellarAge", 0);
HDFWriteGroupAttribute(filename, "Header", "Flag_Metals", 0);
HDFWriteGroupAttribute(filename, "Header", "Flag_Feedback", 0);
HDFWriteGroupAttribute(filename, "Header", "Flag_DoublePrecision", (int)doublePrec);
HDFWriteGroupAttribute(filename, "Header", "Music_levelmin", levelmin_);
HDFWriteGroupAttribute(filename, "Header", "Music_levelmax", levelmax_);
HDFWriteGroupAttribute(filename, "Header", "Music_levelcounts", levelcounts);
HDFWriteGroupAttribute(filename, "Header", "haveBaryons", (int)doBaryons);
HDFWriteGroupAttribute(filename, "Header", "longIDs", (int)useLongIDs);
HDFWriteGroupAttribute(filename, "Header", "suggested_pmgrid", pmgrid);
HDFWriteGroupAttribute(filename, "Header", "suggested_gridboost", gridboost);
HDFWriteGroupAttribute(filename, "Header", "suggested_highressoft", softening);
HDFWriteGroupAttribute(filename, "Header", "suggested_gas_Tinit", Tini);
}
// give config/parameter file hints
if (useLongIDs)
std::cout << " - Arepo: Wrote 64bit IDs, enable LONGIDS." << std::endl;
if (doublePrec)
std::cout << " - Arepo: Double precision ICs, set INPUT_IN_DOUBLEPRECISION." << std::endl;
if (NTYPES != 6)
std::cout << " - Arepo: Using [" << NTYPES << "] particle types, set NTYPES to match." << std::endl;
if (doBaryons)
std::cout << " - Arepo: Wrote high-res gas (only), set REFINEMENT_HIGH_RES_GAS and GENERATE_GAS_IN_ICS with "
<< "SPLIT_PARTICLE_TYPE=" << pow(2, coarsePartType) << "." << std::endl;
if (levelmax_ != levelmin_)
std::cout << " - Arepo: Have zoom type ICs, set PLACEHIGHRESREGION=" << pow(2, HIGHRES_DM_PARTTYPE)
<< " (suggest PMGRID=" << pmgrid << " with GRIDBOOST=" << gridboost << ")." << std::endl;
else
std::cout << " - Arepo: Have unigrid type ICs (suggest PMGRID=" << pmgrid << ")." << std::endl;
if (levelmax_ > levelmin_ + 1)
std::cout << " - Arepo: More than one coarse DM mass using same type, set INDIVIDUAL_GRAVITY_SOFTENING="
<< pow(2, coarsePartType) << " (+" << pow(2, STAR_PARTTYPE) << " if including stars)." << std::endl;
if (doBaryons)
std::cout << " - Arepo: Set initial gas temperature to " << std::fixed << std::setprecision(3) << Tini << " K." << std::endl;
std::cout << " - Arepo: Suggest grav softening = " << std::setprecision(3) << softening << " for high res DM." << std::endl;
}
};
namespace
{
output_plugin_creator_concrete<arepo_output_plugin> creator("arepo");
}
#endif // HAVE_HDF5

16
src/plugins/region_convex_hull.cc
View file

@ -31,7 +31,7 @@ private:
double vfac_;
bool do_extra_padding_;
double anchor_pt_[3];
convex_hull<double>::vec3_t anchor_pt_;
std::vector<float> level_dist_;
@ -123,7 +123,7 @@ public:
delete phull_;
}
void get_AABB( double *left, double *right, unsigned level )
void get_AABB( vec3_t& left, vec3_t& right, unsigned level ) const
{
for( int i=0; i<3; ++i )
{
@ -145,26 +145,26 @@ public:
}
void update_AABB( double *left, double *right )
void update_AABB( vec3_t& left, vec3_t& right )
{
// we ignore this, the grid generator must have generated a grid that contains the ellipsoid
// it might have enlarged it, but who cares...
}
bool query_point( double *x, int ilevel )
{ return phull_->check_point( x, level_dist_[ilevel] ); }
bool query_point( const vec3_t& x, int ilevel ) const
{ return phull_->check_point( &x[0], level_dist_[ilevel] ); }
bool is_grid_dim_forced( size_t* ndims )
bool is_grid_dim_forced( index3_t& ndims ) const
{ return false; }
void get_center( double *xc )
void get_center( vec3_t& xc ) const
{
xc[0] = phull_->centroid_[0];
xc[1] = phull_->centroid_[1];
xc[2] = phull_->centroid_[2];
}
void get_center_unshifted( double *xc )
void get_center_unshifted( vec3_t& xc ) const
{
double dx = 1.0/(1<<shift_level);
double c[3] = { phull_->centroid_[0], phull_->centroid_[1], phull_->centroid_[2] };

718
src/plugins/region_ellipsoid.cc
View file

File diff suppressed because it is too large Load diff

18
src/plugins/region_multibox.cc
View file

@ -35,12 +35,12 @@ typedef std::vector<coord> region;
class region_multibox_plugin : public region_generator_plugin{
private:
double cen[3];
vec3_t cen;
unsigned res;
double vfac_;
grid refgrid;
region where(unsigned level)
region where(unsigned level) const
{
slice curslice;
col curcol;
@ -169,7 +169,7 @@ public:
}
void get_AABB( double *left, double *right, unsigned level )
void get_AABB( vec3_t& left, vec3_t &right, unsigned level ) const
{
left[0] = left[1] = left[2] = 1.0;
right[0] = right[1] = right[2] = 0.0;
@ -216,12 +216,12 @@ public:
fclose(gridfile);
}
void update_AABB( double *left, double *right )
void update_AABB( vec3_t& left, vec3_t& right )
{
//no need for this I think?
}
bool query_point( double *x, int level )
bool query_point( const vec3_t& x, int level ) const
{
for(int i=0; i<3; ++i)
{
@ -230,14 +230,14 @@ public:
return (level <= int(refgrid[(x[0])*res][(x[1])*res][(x[2])*res]));
}
bool is_grid_dim_forced( size_t* ndims )
bool is_grid_dim_forced( index3_t& ndims ) const
{
return false; //is this true?
}
void get_center( double *xc )
void get_center( vec3_t& xc ) const
{
double xc2[3];
vec3_t xc2;
get_AABB(xc, xc2, levelmax_);
for(int i=0; i<3; ++i)
{
@ -246,7 +246,7 @@ public:
}
}
void get_center_unshifted( double *xc )
void get_center_unshifted( vec3_t& xc ) const
{
get_center(xc);
}

247
src/region_generator.cc
View file

@ -1,51 +1,49 @@
#include <algorithm>
#include "region_generator.hh"
std::map< std::string, region_generator_plugin_creator *>&
std::map<std::string, region_generator_plugin_creator *> &
get_region_generator_plugin_map()
{
static std::map< std::string, region_generator_plugin_creator* > region_generator_plugin_map;
return region_generator_plugin_map;
static std::map<std::string, region_generator_plugin_creator *> region_generator_plugin_map;
return region_generator_plugin_map;
}
void print_region_generator_plugins()
{
std::map< std::string, region_generator_plugin_creator *>& m = get_region_generator_plugin_map();
std::map< std::string, region_generator_plugin_creator *>::iterator it;
it = m.begin();
std::cout << " - Available region generator plug-ins:\n";
while( it!=m.end() )
{
if( (*it).second )
std::cout << "\t\'" << (*it).first << "\'\n";
++it;
}
std::map<std::string, region_generator_plugin_creator *> &m = get_region_generator_plugin_map();
std::map<std::string, region_generator_plugin_creator *>::iterator it;
it = m.begin();
std::cout << " - Available region generator plug-ins:\n";
while (it != m.end())
{
if ((*it).second)
std::cout << "\t\'" << (*it).first << "\'\n";
++it;
}
}
region_generator_plugin *select_region_generator_plugin( config_file& cf )
region_generator_plugin *select_region_generator_plugin(config_file &cf)
{
std::string rgname = cf.getValueSafe<std::string>( "setup", "region", "box" );
std::string rgname = cf.getValueSafe<std::string>("setup", "region", "box");
region_generator_plugin_creator *the_region_generator_plugin_creator
= get_region_generator_plugin_map()[ rgname ];
region_generator_plugin_creator *the_region_generator_plugin_creator = get_region_generator_plugin_map()[rgname];
if( !the_region_generator_plugin_creator )
{
std::cerr << " - Error: region generator plug-in \'" << rgname << "\' not found." << std::endl;
LOGERR("Invalid/Unregistered region generator plug-in encountered : %s",rgname.c_str() );
print_region_generator_plugins();
throw std::runtime_error("Unknown region generator plug-in");
if (!the_region_generator_plugin_creator)
{
std::cerr << " - Error: region generator plug-in \'" << rgname << "\' not found." << std::endl;
LOGERR("Invalid/Unregistered region generator plug-in encountered : %s", rgname.c_str());
print_region_generator_plugins();
throw std::runtime_error("Unknown region generator plug-in");
}
else
{
std::cout << " - Selecting region generator plug-in \'" << rgname << "\'..." << std::endl;
LOGUSER("Selecting region generator plug-in : %s", rgname.c_str());
}
}else
{
std::cout << " - Selecting region generator plug-in \'" << rgname << "\'..." << std::endl;
LOGUSER("Selecting region generator plug-in : %s",rgname.c_str() );
}
region_generator_plugin *the_region_generator_plugin = the_region_generator_plugin_creator->create(cf);
region_generator_plugin *the_region_generator_plugin
= the_region_generator_plugin_creator->create( cf );
return the_region_generator_plugin;
return the_region_generator_plugin;
}
/*******************************************************************************/
@ -54,11 +52,12 @@ region_generator_plugin *select_region_generator_plugin( config_file& cf )
#include <cmath>
class region_box_plugin : public region_generator_plugin{
class region_box_plugin : public region_generator_plugin
{
private:
double
x0ref_[3], //!< coordinates of refinement region origin (in [0..1[)
lxref_[3], //!< extent of refinement region (int [0..1[)
x0ref_[3], //!< coordinates of refinement region origin (in [0..1[)
lxref_[3], //!< extent of refinement region (int [0..1[)
xcref_[3];
size_t lnref_[3];
bool bhave_nref_;
@ -68,84 +67,91 @@ private:
double padding_fine_;
public:
region_box_plugin( config_file& cf )
: region_generator_plugin( cf )
region_box_plugin(config_file &cf)
: region_generator_plugin(cf)
{
levelmin_ = pcf_->getValue<unsigned>("setup","levelmin");
levelmax_ = pcf_->getValue<unsigned>("setup","levelmax");
levelmin_ = pcf_->getValue<unsigned>("setup", "levelmin");
levelmax_ = pcf_->getValue<unsigned>("setup", "levelmax");
if( levelmin_ != levelmax_ )
if (levelmin_ != levelmax_)
{
padding_ = cf.getValue<int>("setup","padding");
padding_ = cf.getValue<int>("setup", "padding");
std::string temp;
if( !pcf_->containsKey("setup","ref_offset") && !pcf_->containsKey("setup","ref_center") )
if (!pcf_->containsKey("setup", "ref_offset") && !pcf_->containsKey("setup", "ref_center"))
{
LOGERR("Found levelmin!=levelmax but neither ref_offset nor ref_center was specified.");
throw std::runtime_error("Found levelmin!=levelmax but neither ref_offset nor ref_center was specified.");
}
if( !pcf_->containsKey("setup","ref_extent") && !pcf_->containsKey("setup","ref_dims") )
if (!pcf_->containsKey("setup", "ref_extent") && !pcf_->containsKey("setup", "ref_dims"))
{
LOGERR("Found levelmin!=levelmax but neither ref_extent nor ref_dims was specified.");
throw std::runtime_error("Found levelmin!=levelmax but neither ref_extent nor ref_dims was specified.");
}
if( pcf_->containsKey("setup","ref_extent") )
if (pcf_->containsKey("setup", "ref_extent"))
{
temp = pcf_->getValue<std::string>( "setup", "ref_extent" );
std::remove_if(temp.begin(),temp.end(),isspace);
if(sscanf( temp.c_str(), "%lf,%lf,%lf", &lxref_[0],&lxref_[1],&lxref_[2] )!=3){
LOGERR("Error parsing triple for ref_extent");
throw std::runtime_error("Error parsing triple for ref_extent");
}
temp = pcf_->getValue<std::string>("setup", "ref_extent");
std::remove_if(temp.begin(), temp.end(), isspace);
if (sscanf(temp.c_str(), "%lf,%lf,%lf", &lxref_[0], &lxref_[1], &lxref_[2]) != 3)
{
LOGERR("Error parsing triple for ref_extent");
throw std::runtime_error("Error parsing triple for ref_extent");
}
bhave_nref_ = false;
}else if( pcf_->containsKey("setup","ref_dims") ){
temp = pcf_->getValue<std::string>("setup","ref_dims");
std::remove_if(temp.begin(),temp.end(),isspace);
if(sscanf( temp.c_str(), "%lu,%lu,%lu", &lnref_[0],&lnref_[1],&lnref_[2] )!=3){
LOGERR("Error parsing triple for ref_dims");
throw std::runtime_error("Error parsing triple for ref_dims");
}
}
else if (pcf_->containsKey("setup", "ref_dims"))
{
temp = pcf_->getValue<std::string>("setup", "ref_dims");
std::remove_if(temp.begin(), temp.end(), isspace);
if (sscanf(temp.c_str(), "%lu,%lu,%lu", &lnref_[0], &lnref_[1], &lnref_[2]) != 3)
{
LOGERR("Error parsing triple for ref_dims");
throw std::runtime_error("Error parsing triple for ref_dims");
}
bhave_nref_ = true;
lxref_[0] = lnref_[0] * 1.0/(double)(1<<levelmax_);
lxref_[1] = lnref_[1] * 1.0/(double)(1<<levelmax_);
lxref_[2] = lnref_[2] * 1.0/(double)(1<<levelmax_);
lxref_[0] = lnref_[0] * 1.0 / (double)(1 << levelmax_);
lxref_[1] = lnref_[1] * 1.0 / (double)(1 << levelmax_);
lxref_[2] = lnref_[2] * 1.0 / (double)(1 << levelmax_);
}
if( pcf_->containsKey("setup","ref_center") )
if (pcf_->containsKey("setup", "ref_center"))
{
temp = pcf_->getValue<std::string>( "setup", "ref_center" );
std::remove_if(temp.begin(),temp.end(),isspace);
if(sscanf( temp.c_str(), "%lf,%lf,%lf", &xcref_[0], &xcref_[1], &xcref_[2] )!=3){
LOGERR("Error parsing triple for ref_center");
throw std::runtime_error("Error parsing triple for ref_center");
}
x0ref_[0] = fmod( xcref_[0]-0.5*lxref_[0]+1.0,1.0);
x0ref_[1] = fmod( xcref_[1]-0.5*lxref_[1]+1.0,1.0);
x0ref_[2] = fmod( xcref_[2]-0.5*lxref_[2]+1.0,1.0);
}else if( pcf_->containsKey("setup","ref_offset") ){
temp = pcf_->getValue<std::string>( "setup", "ref_offset" );
std::remove_if(temp.begin(),temp.end(),isspace);
if(sscanf( temp.c_str(), "%lf,%lf,%lf", &x0ref_[0], &x0ref_[1], &x0ref_[2] )!=3){
LOGERR("Error parsing triple for ref_offset");
throw std::runtime_error("Error parsing triple for ref_offset");
}
xcref_[0] = fmod( x0ref_[0]+0.5*lxref_[0], 1.0 );
xcref_[1] = fmod( x0ref_[1]+0.5*lxref_[1], 1.0 );
xcref_[2] = fmod( x0ref_[2]+0.5*lxref_[2], 1.0 );
temp = pcf_->getValue<std::string>("setup", "ref_center");
std::remove_if(temp.begin(), temp.end(), isspace);
if (sscanf(temp.c_str(), "%lf,%lf,%lf", &xcref_[0], &xcref_[1], &xcref_[2]) != 3)
{
LOGERR("Error parsing triple for ref_center");
throw std::runtime_error("Error parsing triple for ref_center");
}
x0ref_[0] = fmod(xcref_[0] - 0.5 * lxref_[0] + 1.0, 1.0);
x0ref_[1] = fmod(xcref_[1] - 0.5 * lxref_[1] + 1.0, 1.0);
x0ref_[2] = fmod(xcref_[2] - 0.5 * lxref_[2] + 1.0, 1.0);
}
else if (pcf_->containsKey("setup", "ref_offset"))
{
temp = pcf_->getValue<std::string>("setup", "ref_offset");
std::remove_if(temp.begin(), temp.end(), isspace);
if (sscanf(temp.c_str(), "%lf,%lf,%lf", &x0ref_[0], &x0ref_[1], &x0ref_[2]) != 3)
{
LOGERR("Error parsing triple for ref_offset");
throw std::runtime_error("Error parsing triple for ref_offset");
}
xcref_[0] = fmod(x0ref_[0] + 0.5 * lxref_[0], 1.0);
xcref_[1] = fmod(x0ref_[1] + 0.5 * lxref_[1], 1.0);
xcref_[2] = fmod(x0ref_[2] + 0.5 * lxref_[2], 1.0);
}
// conditions should be added here
{
do_extra_padding_ = false;
std::string output_plugin = cf.getValue<std::string>("output","format");
if( output_plugin == std::string("grafic2") )
std::string output_plugin = cf.getValue<std::string>("output", "format");
if (output_plugin == std::string("grafic2"))
do_extra_padding_ = true;
padding_fine_ = 0.0;
if( do_extra_padding_ )
padding_fine_ = (double)(padding_+1) * 1.0/(1ul<<levelmax_);
if (do_extra_padding_)
padding_fine_ = (double)(padding_ + 1) * 1.0 / (1ul << levelmax_);
}
}
else
@ -153,78 +159,81 @@ public:
x0ref_[0] = x0ref_[1] = x0ref_[2] = 0.0;
lxref_[0] = lxref_[1] = lxref_[2] = 1.0;
xcref_[0] = xcref_[1] = xcref_[2] = 0.5;
}
}
void get_AABB( double *left, double *right, unsigned level )
void get_AABB(vec3_t &left, vec3_t &right, unsigned level) const
{
double dx = 1.0/(1ul<<level);
double pad = (double)(padding_+1) * dx;
double dx = 1.0 / (1ul << level);
double pad = (double)(padding_ + 1) * dx;
if( ! do_extra_padding_ ) pad = 0.0;
if (!do_extra_padding_)
pad = 0.0;
for( int i=0; i<3; ++i )
for (int i = 0; i < 3; ++i)
{
left[i] = x0ref_[i] - pad;
right[i] = x0ref_[i] + lxref_[i] + pad;
}
}
void update_AABB( double *left, double *right )
void update_AABB(vec3_t &left, vec3_t &right)
{
for( int i=0; i<3; ++i )
{
double dx = right[i] - left[i];
if( dx < -0.5 ) dx += 1.0; else if (dx > 0.5 ) dx -= 1.0;
x0ref_[i] = left[i];
lxref_[i] = dx;
xcref_[i] = left[i] + 0.5 * dx;
}
//fprintf(stderr,"left = %f,%f,%f - right = %f,%f,%f\n",left[0],left[1],left[2],right[0],right[1],right[2]);
for (int i = 0; i < 3; ++i)
{
double dx = right[i] - left[i];
if (dx < -0.5)
dx += 1.0;
else if (dx > 0.5)
dx -= 1.0;
x0ref_[i] = left[i];
lxref_[i] = dx;
xcref_[i] = left[i] + 0.5 * dx;
}
// fprintf(stderr,"left = %f,%f,%f - right = %f,%f,%f\n",left[0],left[1],left[2],right[0],right[1],right[2]);
}
bool query_point( double *x, int ilevel )
bool query_point( const vec3_t &x, int ilevel) const
{
if( !do_extra_padding_ )
if (!do_extra_padding_)
return true;
bool check = true;
double dx;
for( int i=0; i<3; ++i )
for (int i = 0; i < 3; ++i)
{
dx = x[i] - x0ref_[i];
if( dx < -0.5 ) dx += 1.0;
else if (dx > 0.5 ) dx -= 1.0;
if (dx < -0.5)
dx += 1.0;
else if (dx > 0.5)
dx -= 1.0;
check &= ((dx >= padding_fine_) & (dx <= lxref_[i]-padding_fine_));
check &= ((dx >= padding_fine_) & (dx <= lxref_[i] - padding_fine_));
}
return check;
}
bool is_grid_dim_forced( size_t* ndims )
bool is_grid_dim_forced(index3_t& ndims) const
{
for( int i=0; i<3; ++i )
for (int i = 0; i < 3; ++i)
ndims[i] = lnref_[i];
return bhave_nref_;
}
void get_center( double *xc )
void get_center(vec3_t &xc) const
{
xc[0] = xcref_[0];
xc[1] = xcref_[1];
xc[2] = xcref_[2];
}
void get_center_unshifted( double *xc )
{
get_center( xc );
}
void get_center_unshifted(vec3_t &xc) const
{
get_center(xc);
}
};
namespace{
region_generator_plugin_creator_concrete< region_box_plugin > creator("box");
namespace
{
region_generator_plugin_creator_concrete<region_box_plugin> creator("box");
}

16
src/region_generator.hh
View file

@ -4,6 +4,10 @@
#include <vector>
#include "config_file.hh"
#include <array>
using vec3_t = std::array<double,3>;
using index3_t = std::array<ptrdiff_t,3>;
//! Abstract base class for region generators
/*!
This class implements a purely virtual interface that can be
@ -26,22 +30,22 @@ public:
virtual ~region_generator_plugin() { };
//! compute the bounding box of the region
virtual void get_AABB( double *left, double *right, unsigned level) = 0;
virtual void get_AABB( vec3_t& left, vec3_t& right, unsigned level) const = 0;
//! query whether a point intersects the region
virtual bool query_point( double *x, int level ) = 0;
virtual bool query_point( const vec3_t& x, int level ) const = 0;
//! query whether the region generator explicitly forces the grid dimensions
virtual bool is_grid_dim_forced( size_t *ndims ) = 0;
virtual bool is_grid_dim_forced( index3_t& ndims ) const = 0;
//! get the center of the region
virtual void get_center( double *xc ) = 0;
virtual void get_center( vec3_t& xc ) const = 0;
//! get the center of the region with a possible re-centering unapplied
virtual void get_center_unshifted( double *xc ) = 0;
virtual void get_center_unshifted( vec3_t& xc ) const = 0;
//! update the highres bounding box to what the grid generator actually uses
virtual void update_AABB( double *left, double *right ) = 0;
virtual void update_AABB( vec3_t& left, vec3_t& right ) = 0;
};
//! Implements abstract factory design pattern for region generator plug-ins