Skip to content

Commit

Permalink
Fixed AMR segfault
Browse files Browse the repository at this point in the history 
Had to move variables from velocity_mesh_amr.h to
velocity_mesh_parameters.h in order to get rid of segfaults.

Changed how velocity mesh max refinement level is written out,
vlsvextract now correctly converts multi-population AMR meshes.
  • Loading branch information
@sandroos
sandroos committed Jun 5, 2015
1 parent 47c3ac1 commit 3e4de6b
Show file tree
Hide file tree
Showing 12 changed files with 218 additions and 141 deletions.
10 changes: 5 additions & 5 deletions Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ DEPS_VLSVMOVER = ${DEPS_CELL} vlasovsolver/vlasovmover.cpp vlasovsolver/cpu_acc_

DEPS_VLSVMOVER_AMR = ${DEPS_CELL} vlasovsolver_amr/vlasovmover.cpp vlasovsolver_amr/cpu_acc_map.hpp vlasovsolver_amr/cpu_acc_intersections.hpp \
vlasovsolver_amr/cpu_acc_intersections.hpp vlasovsolver_amr/cpu_acc_semilag.hpp vlasovsolver_amr/cpu_acc_transform.hpp \
vlasovsolver_amr/cpu_moments.h vlasovsolver_amr/cpu_trans_map.hpp velocity_blocks.h
vlasovsolver/cpu_moments.h vlasovsolver_amr/cpu_trans_map.hpp velocity_blocks.h

#DEPS_PROJECTS = projects/project.h projects/project.cpp \
# projects/MultiPeak/MultiPeak.h projects/MultiPeak/MultiPeak.cpp ${DEPS_CELL}
Expand All @@ -204,7 +204,7 @@ OBJS = version.o memoryallocation.o backgroundfield.o quadr.o dipole.o linedipo

# Add Vlasov solver objects (depend on mesh: AMR or non-AMR)
ifeq ($(MESH),AMR)

OBJS += cpu_moments.o
else
OBJS += cpu_acc_intersections.o cpu_acc_map.o cpu_acc_sort_blocks.o cpu_acc_semilag.o cpu_acc_transform.o \
cpu_moments.o cpu_trans_map.o
Expand Down Expand Up @@ -419,16 +419,16 @@ cpu_acc_sort_blocks.o: ${DEPS_CPU_ACC_SORT_BLOCKS}
cpu_acc_transform.o: ${DEPS_CPU_ACC_TRANSFORM}
${CMP} ${CXXFLAGS} ${FLAG_OPENMP} ${MATHFLAGS} ${FLAGS} -c vlasovsolver/cpu_acc_transform.cpp ${INC_EIGEN} ${INC_DCCRG} ${INC_PROFILE} ${INC_ZOLTAN}

cpu_moments.o: ${DEPS_CPU_MOMENTS}
${CMP} ${CXXFLAGS} ${FLAG_OPENMP} ${MATHFLAGS} ${FLAGS} -c vlasovsolver/cpu_moments.cpp ${INC_DCCRG} ${INC_ZOLTAN} ${INC_PROFILE}

cpu_trans_map.o: ${DEPS_CPU_TRANS_MAP}
${CMP} ${CXXFLAGS} ${FLAG_OPENMP} ${MATHFLAGS} ${FLAGS} -c vlasovsolver/cpu_trans_map.cpp ${INC_EIGEN} ${INC_DCCRG} ${INC_PROFILE} ${INC_VECTORCLASS} ${INC_ZOLTAN} ${INC_VLSV}

vlasovmover.o: ${DEPS_VLSVMOVER}
${CMP} ${CXXFLAGS} ${FLAG_OPENMP} ${MATHFLAGS} ${FLAGS} -c vlasovsolver/vlasovmover.cpp -I$(CURDIR) ${INC_BOOST} ${INC_EIGEN} ${INC_DCCRG} ${INC_ZOLTAN} ${INC_PROFILE} ${INC_VECTORCLASS} ${INC_EIGEN} ${INC_VLSV}
endif

cpu_moments.o: ${DEPS_CPU_MOMENTS}
${CMP} ${CXXFLAGS} ${FLAG_OPENMP} ${MATHFLAGS} ${FLAGS} -c vlasovsolver/cpu_moments.cpp ${INC_DCCRG} ${INC_ZOLTAN} ${INC_PROFILE}

derivatives.o: ${DEPS_FSOLVER} fieldsolver/fs_limiters.h fieldsolver/fs_limiters.cpp fieldsolver/derivatives.hpp fieldsolver/derivatives.cpp
${CMP} ${CXXFLAGS} ${FLAGS} -c fieldsolver/derivatives.cpp -I$(CURDIR) ${INC_BOOST} ${INC_EIGEN} ${INC_DCCRG} ${INC_PROFILE} ${INC_ZOLTAN}

Expand Down
2 changes: 1 addition & 1 deletion grid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -959,7 +959,7 @@ bool validateMesh(dccrg::Dccrg<SpatialCell,dccrg::Cartesian_Geometry>& mpiGrid,c
// min value, add the block to remove list
#pragma omp for
for (size_t b=0; b<newBlocks[c].size(); ++b) {
if (getObjectWrapper().project->setVelocityBlock(cell,newBlocks[c][b].second) <= Parameters::sparseMinValue) {
if (getObjectWrapper().project->setVelocityBlock(cell,newBlocks[c][b].second,popID) <= Parameters::sparseMinValue) {
threadRemBlocks[tid].push_back(newBlocks[c][b].first);
++counter[tid];
}
Expand Down
8 changes: 7 additions & 1 deletion iowrite.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,6 +157,12 @@ bool writeVelocityDistributionData(const int& popID,Writer& vlsvWriter,
attribs.clear();
attribs["mesh"] = getObjectWrapper().particleSpecies[popID].name;
attribs["type"] = vlsv::mesh::STRING_UCD_AMR;

// stringstream is necessary here to correctly convert refLevelMaxAllowed into a string
stringstream ss;
ss << static_cast<unsigned int>(getObjectWrapper().velocityMeshes[meshID].refLevelMaxAllowed);
attribs["max_velocity_ref_level"] = ss.str();

if (mpiGrid.get_rank() == MASTER_RANK) {
if (vlsvWriter.writeArray("MESH_BBOX",attribs,6,1,bbox) == false) success = false;

Expand Down Expand Up @@ -426,7 +432,7 @@ bool writeCommonGridData(
//if( vlsvWriter.writeParameter("vyblocks_ini", &P::vyblocks_ini) == false ) { return false; }
//if( vlsvWriter.writeParameter("vzblocks_ini", &P::vzblocks_ini) == false ) { return false; }
//if( vlsvWriter.writeParameter("max_velocity_ref_level", &P::amrMaxVelocityRefLevel) == false) {return false;}

//Mark the new version:
float version = 1.00;
if( vlsvWriter.writeParameter( "version", &version ) == false ) { return false; }
Expand Down
129 changes: 63 additions & 66 deletions projects/project.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ struct VelocityMeshParams {
vector<double> vx_length;
vector<double> vy_length;
vector<double> vz_length;
vector<unsigned int> maxRefLevels;

void resize(const size_t& size) {
name.resize(1);
Expand All @@ -64,6 +65,7 @@ struct VelocityMeshParams {
vx_length.resize(1);
vy_length.resize(1);
vz_length.resize(1);
maxRefLevels.resize(1);
}
};

Expand Down Expand Up @@ -123,6 +125,7 @@ namespace projects {
RP::addComposing("velocitymesh.vx_length","Initial number of velocity blocks in vx-direction.");
RP::addComposing("velocitymesh.vy_length","Initial number of velocity blocks in vy-direction.");
RP::addComposing("velocitymesh.vz_length","Initial number of velocity blocks in vz-direction.");
RP::addComposing("velocitymesh.max_refinement_level","Maximum allowed mesh refinement level.");

// These parameters are only read if the 'velocitymesh.' parameters are not defined
// in order to support older configuration files.
Expand Down Expand Up @@ -159,6 +162,7 @@ namespace projects {
RP::get("velocitymesh.vx_length",velMeshParams->vx_length);
RP::get("velocitymesh.vy_length",velMeshParams->vy_length);
RP::get("velocitymesh.vz_length",velMeshParams->vz_length);
RP::get("velocitymesh.max_refinement_level",velMeshParams->maxRefLevels);
}

bool Project::initialize() {
Expand Down Expand Up @@ -190,6 +194,7 @@ namespace projects {
if (velMeshParams->vx_length.size() != velMeshParams->name.size()) success = false;
if (velMeshParams->vy_length.size() != velMeshParams->name.size()) success = false;
if (velMeshParams->vz_length.size() != velMeshParams->name.size()) success = false;
if (velMeshParams->maxRefLevels.size() != velMeshParams->name.size()) success = false;
if (success == false) {
stringstream ss;
ss << "(PROJECT) ERROR in configuration file velocity mesh definitions at ";
Expand All @@ -215,6 +220,7 @@ namespace projects {
RP::get("gridbuilder.vx_length",velMeshParams->vx_length[0]);
RP::get("gridbuilder.vy_length",velMeshParams->vy_length[0]);
RP::get("gridbuilder.vz_length",velMeshParams->vz_length[0]);
velMeshParams->maxRefLevels[0] = 0;
}

// Store velocity mesh parameters
Expand Down Expand Up @@ -249,6 +255,7 @@ namespace projects {
meshParams.blockLength[0] = WID;
meshParams.blockLength[1] = WID;
meshParams.blockLength[2] = WID;
meshParams.refLevelMaxAllowed = velMeshParams->maxRefLevels[m];
owrapper.velocityMeshes.push_back(meshParams);
}

Expand Down Expand Up @@ -415,26 +422,23 @@ namespace projects {
logFile << write;
}

void Project::setVelocitySpace(const int& popID,SpatialCell* cell) const {
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& vmesh = cell->get_velocity_mesh(popID);

vector<vmesh::GlobalID> blocksToInitialize = this->findBlocksToInitialize(cell,popID);
Real* parameters = cell->get_block_parameters(popID);

creal x = cell->parameters[CellParams::XCRD];
creal y = cell->parameters[CellParams::YCRD];
creal z = cell->parameters[CellParams::ZCRD];
creal dx = cell->parameters[CellParams::DX];
creal dy = cell->parameters[CellParams::DY];
creal dz = cell->parameters[CellParams::DZ];
Realf* data = cell->get_data(popID);

/** Calculate the volume averages of distribution function for the
* given particle population in the given spatial cell. The velocity block
* is defined by its local ID. The function returns the maximum value of the
* distribution function within the velocity block. If it is below the sparse
* min value for the population, this block should be removed or marked
* as a no-content block.
* @param cell Spatial cell.
* @param blockLID Velocity block local ID within the spatial cell.
* @param popID Population ID.
* @return Maximum value of the calculated distribution function.*/
Real Project::setVelocityBlock(spatial_cell::SpatialCell* cell,const vmesh::LocalID& blockLID,const int& popID) const {
// If simulation doesn't use one or more velocity coordinates,
// only calculate the distribution function for one layer of cells.
uint WID_VX = WID;
uint WID_VY = WID;
uint WID_VZ = WID;
switch (Parameters::geometry) {
switch (Parameters::geometry) {
case geometry::XY4D:
WID_VZ=1;
break;
Expand All @@ -445,6 +449,48 @@ namespace projects {
break;
}

// Fetch spatial cell coordinates and size
creal x = cell->parameters[CellParams::XCRD];
creal y = cell->parameters[CellParams::YCRD];
creal z = cell->parameters[CellParams::ZCRD];
creal dx = cell->parameters[CellParams::DX];
creal dy = cell->parameters[CellParams::DY];
creal dz = cell->parameters[CellParams::DZ];

const Real* parameters = cell->get_block_parameters(popID);
Realf* data = cell->get_data(popID);

creal vxBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VXCRD];
creal vyBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VYCRD];
creal vzBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VZCRD];
creal dvxCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVX];
creal dvyCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVY];
creal dvzCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVZ];

// Calculate volume average of distribution function for each phase-space cell in the block.
Real maxValue = 0.0;
for (uint kc=0; kc<WID_VZ; ++kc) for (uint jc=0; jc<WID_VY; ++jc) for (uint ic=0; ic<WID_VX; ++ic) {
creal vxCell = vxBlock + ic*dvxCell;
creal vyCell = vyBlock + jc*dvyCell;
creal vzCell = vzBlock + kc*dvzCell;
creal average =
calcPhaseSpaceDensity(
x, y, z, dx, dy, dz,
vxCell,vyCell,vzCell,
dvxCell,dvyCell,dvzCell,popID);
if (average != 0.0) {
data[blockLID*SIZE_VELBLOCK+cellIndex(ic,jc,kc)] = average;
maxValue = max(maxValue,average);
}
}

return maxValue;
}

void Project::setVelocitySpace(const int& popID,SpatialCell* cell) const {
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& vmesh = cell->get_velocity_mesh(popID);

vector<vmesh::GlobalID> blocksToInitialize = this->findBlocksToInitialize(cell,popID);
vector<vmesh::GlobalID> removeList;
for (uint i=0; i<blocksToInitialize.size(); ++i) {
const vmesh::GlobalID blockGID = blocksToInitialize[i];
Expand All @@ -454,31 +500,7 @@ namespace projects {
exit(1);
}

creal vxBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VXCRD];
creal vyBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VYCRD];
creal vzBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VZCRD];
creal dvxCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVX];
creal dvyCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVY];
creal dvzCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVZ];

// Calculate volume average of distrib. function for each cell in the block.
Real maxValue = 0.0;
for (uint kc=0; kc<WID_VZ; ++kc) for (uint jc=0; jc<WID_VY; ++jc) for (uint ic=0; ic<WID_VX; ++ic) {
creal vxCell = vxBlock + ic*dvxCell;
creal vyCell = vyBlock + jc*dvyCell;
creal vzCell = vzBlock + kc*dvzCell;
creal average =
calcPhaseSpaceDensity(
x, y, z, dx, dy, dz,
vxCell,vyCell,vzCell,
dvxCell,dvyCell,dvzCell,popID);

if (average != 0.0) {
data[blockLID*SIZE_VELBLOCK+cellIndex(ic,jc,kc)] = average;
maxValue = max(maxValue,average);
}
}

const Real maxValue = setVelocityBlock(cell,blockLID,popID);
if (maxValue < getObjectWrapper().particleSpecies[popID].sparseMinValue) removeList.push_back(blockGID);
}

Expand Down Expand Up @@ -533,32 +555,7 @@ namespace projects {
for (map<vmesh::GlobalID,vmesh::LocalID>::const_iterator it=insertedBlocks.begin(); it!=insertedBlocks.end(); ++it) {
const vmesh::GlobalID blockGID = it->first;
const vmesh::LocalID blockLID = it->second;
parameters = cell->get_block_parameters(popID);
creal vxBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VXCRD];
creal vyBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VYCRD];
creal vzBlock = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::VZCRD];
creal dvxCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVX];
creal dvyCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVY];
creal dvzCell = parameters[blockLID*BlockParams::N_VELOCITY_BLOCK_PARAMS + BlockParams::DVZ];

Real maxValue = 0.0;
for (uint kc=0; kc<WID; ++kc) {
for (uint jc=0; jc<WID; ++jc) {
for (uint ic=0; ic<WID; ++ic) {
creal vxCell = vxBlock + ic*dvxCell;
creal vyCell = vyBlock + jc*dvyCell;
creal vzCell = vzBlock + kc*dvzCell;
creal average =
calcPhaseSpaceDensity(
x, y, z, dx, dy, dz,
vxCell,vyCell,vzCell,
dvxCell,dvyCell,dvzCell,popID);
cell->get_data(popID)[blockLID*SIZE_VELBLOCK + kc*WID2+jc*WID+ic] = average;
maxValue = max(maxValue,average);
}
}
}

const Real maxValue = setVelocityBlock(cell,blockLID,popID);
if (maxValue <= getObjectWrapper().particleSpecies[popID].sparseMinValue)
removeList.push_back(it->first);
}
Expand Down
2 changes: 1 addition & 1 deletion projects/project.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ namespace projects {
*/
void setCell(spatial_cell::SpatialCell* cell);

Real setVelocityBlock(spatial_cell::SpatialCell* cell,const vmesh::LocalID& blockLID);
Real setVelocityBlock(spatial_cell::SpatialCell* cell,const vmesh::LocalID& blockLID,const int& popID) const;

protected:
/*! \brief Returns a list of blocks to loop through when initialising.
Expand Down
24 changes: 24 additions & 0 deletions readparameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,6 +315,30 @@ bool Readparameters::get(const std::string& name,std::vector<int>& value) {
return ret;
}

/** Get the value of the given parameter added with addComposing(). This may be called after having called Parse, and it may be called
* by any process, in any order.
* @param name The name of the parameter.
* @param value A variable where the value of the parameter is written.
* @return If true, the given parameter was found and its value was written to value.
*/
bool Readparameters::get(const std::string& name,std::vector<unsigned int>& value) {
vector<string> stringValue;
bool ret;
using boost::lexical_cast;
ret=Readparameters::get(name,stringValue);
if (ret) {
for (vector<string>::iterator i = stringValue.begin(); i!=stringValue.end(); ++i) {
try {
value.push_back(lexical_cast<unsigned int>(*i));
}
catch (...){
if(Readparameters::rank==0) cerr << "Problems casting value to unsigned int " << name <<" = " << *i <<endl;
return false;
}
}
}
return ret;
}


/** Get the value of the given parameter added with addComposing(). This may be called after having called Parse, and it may be called
Expand Down
1 change: 1 addition & 0 deletions readparameters.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ struct Readparameters {
static bool addComposing(const std::string& name,const std::string& desc);
static bool get(const std::string& name,std::vector<std::string>& value);
static bool get(const std::string& name,std::vector<int>& value);
static bool get(const std::string& name,std::vector<unsigned int>& value);
static bool get(const std::string& name,std::vector<float>& value);
static bool get(const std::string& name,std::vector<double>& value);

Expand Down
6 changes: 4 additions & 2 deletions tools/vlsvextract.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1190,8 +1190,10 @@ bool setVelocityMeshVariables(vlsv::Reader& vlsvReader,CellStructure& cellStruct
// was actually given as a parameter.
uint32_t dummyUInt;
cellStruct.maxVelRefLevel = 0;
if (vlsvReader.readParameter("max_velocity_ref_level",dummyUInt) == true) {
cellStruct.maxVelRefLevel = dummyUInt;
map<string,string> attribsOut;
vlsvReader.getArrayAttributes("MESH_BBOX",attribs,attribsOut);
if (attribsOut.find("max_velocity_ref_level") != attribsOut.end()) {
cellStruct.maxVelRefLevel = atoi(attribsOut["max_velocity_ref_level"].c_str());
}

return true;
Expand Down
Loading

0 comments on commit 3e4de6b

Please sign in to comment.