Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfix (ProbeVolume): Output IS (with natural ordering) to re-arrange sub-data (post-processing) #147

Merged
merged 2 commits into from
Feb 8, 2020
Merged

Bugfix (ProbeVolume): Output IS (with natural ordering) to re-arrange sub-data (post-processing) #147

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
28bde91
Bugfix: output IS (natural ordering) to re-arrange sub-data
mesnardo Feb 7, 2020
cb3a15c
Bump version to 0.5.1
mesnardo Feb 8, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 38 additions & 10 deletions include/petibm/probes.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,8 +144,11 @@ class ProbeVolume : public ProbeBase
/** \brief Limits of the volume. */
type::RealVec2D box;

/** \brief Index set for the grid points to monitor. */
IS is;
/** \brief Index set for the grid points to monitor (PETSc ordering). */
IS isPetsc;

/** \brief Index set for the grid points to monitor (Natural ordering). */
IS isNatural;

/** \brief Sub-vector of the region to monitor. */
Vec dvec;
Expand Down Expand Up @@ -182,8 +185,8 @@ class ProbeVolume : public ProbeBase
* \param box [in] Box area to monitor
* \return PetscErrorCode
*/
PetscErrorCode getSubMeshInfo(const type::Mesh &mesh,
const type::RealVec2D &box);
PetscErrorCode getInfo(const type::Mesh &mesh,
const type::RealVec2D &box);

/** \brief Create the index set for the points to monitor.
*
Expand All @@ -197,28 +200,53 @@ class ProbeVolume : public ProbeBase
* \param mesh [in] Cartesian mesh object
* \return PetscErrorCode
*/
PetscErrorCode createSubMesh(const type::Mesh &mesh);
PetscErrorCode createGrid(const type::Mesh &mesh);

/** \brief Write the sub mesh grid points into a file.
*
* Supported formats are HDF5 and ASCII.
*
* \param filePath [in] Path of the file to write in
* \return PetscErrorCode
*/
PetscErrorCode writeSubMesh(const std::string &filePath);
PetscErrorCode writeGrid(const std::string &filePath);

/** \brief Write the sub mesh into an ASCII file.
*
* \param filePath [in] Path of the file to write in
* \return PetscErrorCode
*/
PetscErrorCode writeSubMeshASCII(const std::string &filePath);
PetscErrorCode writeGrid_ASCII(const std::string &filePath);

/** \brief Write the sub mesh into a HDF5 file.
*
* \param filePath [in] Path of the file to write in
* \return PetscErrorCode
*/
PetscErrorCode writeSubMeshHDF5(const std::string &filePath);
PetscErrorCode writeGrid_HDF5(const std::string &filePath);

/** \brief Write index set (natural ordering) into a file.
*
* Supported formats are HDF5 and ASCII.
*
* \param filePath [in] Path of the file to write in
* \return PetscErrorCode
*/
PetscErrorCode writeIS(const std::string &filePath);

/** \brief Write index set (natural ordering) into a HDF5 file.
*
* \param filePath [in] Path of the file to write in
* \return PetscErrorCode
*/
PetscErrorCode writeIS_HDF5(const std::string &filePath);

/** \brief Write index set (natural ordering) into a HDF5 file.
*
* \param filePath [in] Path of the file to write in
* \return PetscErrorCode
*/
PetscErrorCode writeIS_ASCII(const std::string &filePath);

/** \copydoc ProbeBase::monitorVec() */
PetscErrorCode monitorVec(const DM &da,
Expand All @@ -242,15 +270,15 @@ class ProbeVolume : public ProbeBase
* \param t [in] Time
* \return PetscErrorCode
*/
PetscErrorCode writeVecASCII(const Vec &vec, const PetscReal &t);
PetscErrorCode writeVec_ASCII(const Vec &vec, const PetscReal &t);

/** \brief Output a PETSc Vec object to a HDF5 file.
*
* \param vec [in] PETSc Vec object to output
* \param t [in] Time
* \return PetscErrorCode
*/
PetscErrorCode writeVecHDF5(const Vec &vec, const PetscReal &t);
PetscErrorCode writeVec_HDF5(const Vec &vec, const PetscReal &t);

}; // ProbeVolume

Expand Down
131 changes: 100 additions & 31 deletions src/misc/probes.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,8 @@ PetscErrorCode ProbeVolume::init(const MPI_Comm &comm,
// data are added together and we write the time averaged data
n_sum = node["n_sum"].as<PetscInt>(0);

is = PETSC_NULL;
isPetsc = PETSC_NULL;
isNatural = PETSC_NULL;
dvec = PETSC_NULL;

// store information about the sub-volume to monitor
Expand All @@ -198,13 +199,11 @@ PetscErrorCode ProbeVolume::init(const MPI_Comm &comm,
startIdxDir.resize(3, 0);

// get information about the location of the sub-mesh
ierr = getSubMeshInfo(mesh, box); CHKERRQ(ierr);
ierr = getInfo(mesh, box); CHKERRQ(ierr);
// create gridline coordinates for the sub-mesh
ierr = createSubMesh(mesh); CHKERRQ(ierr);
ierr = createGrid(mesh); CHKERRQ(ierr);
// write the sub-mesh to file
ierr = writeSubMesh(path); CHKERRQ(ierr);
// create a PETSc Index Set object to easily grab a sub-vector
ierr = createIS(mesh); CHKERRQ(ierr);
ierr = writeGrid(path); CHKERRQ(ierr);

// create a PETSc Viewer to output the data
ierr = PetscViewerCreate(comm, &viewer); CHKERRQ(ierr);
Expand All @@ -214,6 +213,11 @@ PetscErrorCode ProbeVolume::init(const MPI_Comm &comm,
ierr = PetscViewerFileSetMode(viewer, FILE_MODE_APPEND); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, path.c_str()); CHKERRQ(ierr);

// create a PETSc Index Set objects to the sub-vector of interest
ierr = createIS(mesh); CHKERRQ(ierr);
// write Index Set (natural ordering) to file
ierr = writeIS(path); CHKERRQ(ierr);

PetscFunctionReturn(0);
} // ProbeVolume::init

Expand All @@ -225,15 +229,16 @@ PetscErrorCode ProbeVolume::destroy()
PetscFunctionBeginUser;

ierr = ProbeBase::destroy(); CHKERRQ(ierr);
if (is != PETSC_NULL) {ierr = ISDestroy(&is); CHKERRQ(ierr);}
if (isPetsc != PETSC_NULL) {ierr = ISDestroy(&isPetsc); CHKERRQ(ierr);}
if (isNatural != PETSC_NULL) {ierr = ISDestroy(&isNatural); CHKERRQ(ierr);}
if (dvec != PETSC_NULL) {ierr = VecDestroy(&dvec); CHKERRQ(ierr);}

PetscFunctionReturn(0);
} // ProbeVolume::destroy

// Get information about the sub-mesh area to monitor.
PetscErrorCode ProbeVolume::getSubMeshInfo(const type::Mesh &mesh,
const type::RealVec2D &box)
PetscErrorCode ProbeVolume::getInfo(const type::Mesh &mesh,
const type::RealVec2D &box)
{
std::vector<PetscReal>::iterator low, up;

Expand All @@ -256,9 +261,9 @@ PetscErrorCode ProbeVolume::getSubMeshInfo(const type::Mesh &mesh,
1, std::multiplies<PetscInt>());

PetscFunctionReturn(0);
} // ProbeVolume::getSubMeshInfo
} // ProbeVolume::getInfo

// Create the index set for the points to monitor.
// Create the index sets (PETSc and natural ordering) of the points to monitor.
PetscErrorCode ProbeVolume::createIS(const type::Mesh &mesh)
{
PetscErrorCode ierr;
Expand Down Expand Up @@ -287,17 +292,25 @@ PetscErrorCode ProbeVolume::createIS(const type::Mesh &mesh)
}
indices.resize(count);

// AO will map `indices` from natural ordering to PETSc ordering.
// However, we need to keep the vector with natural-ordering indices
// so we can post-process the solution in the sub-volume.
std::vector<PetscInt> indices2(indices);

AO ao;
ierr = DMDAGetAO(mesh->da[field], &ao); CHKERRQ(ierr);
ierr = AOApplicationToPetsc(ao, count, &indices[0]); CHKERRQ(ierr);
ierr = ISCreateGeneral(
comm, count, &indices[0], PETSC_COPY_VALUES, &is); CHKERRQ(ierr);
ierr = ISCreateGeneral(comm, count, &indices[0],
PETSC_COPY_VALUES, &isPetsc); CHKERRQ(ierr);
// Create IS containing the (natural )index of the points in the sub-volume.
ierr = ISCreateGeneral(comm, count, &indices2[0],
PETSC_COPY_VALUES, &isNatural); CHKERRQ(ierr);

PetscFunctionReturn(0);
} // ProbeVolume::createIS

// Get the sub-mesh area to monitor.
PetscErrorCode ProbeVolume::createSubMesh(const type::Mesh &mesh)
PetscErrorCode ProbeVolume::createGrid(const type::Mesh &mesh)
{
PetscFunctionBeginUser;

Expand All @@ -313,10 +326,10 @@ PetscErrorCode ProbeVolume::createSubMesh(const type::Mesh &mesh)
}

PetscFunctionReturn(0);
} // ProbeVolume::createSubMesh
} // ProbeVolume::createGrid

// Write the sub-mesh grid points into a file.
PetscErrorCode ProbeVolume::writeSubMesh(const std::string &filePath)
PetscErrorCode ProbeVolume::writeGrid(const std::string &filePath)
{
PetscErrorCode ierr;

Expand All @@ -325,22 +338,22 @@ PetscErrorCode ProbeVolume::writeSubMesh(const std::string &filePath)
// only ASCII and HDF5 formats are currently supported
if (std::strcmp(viewerType, "ascii") == 0)
{
ierr = writeSubMeshASCII(filePath); CHKERRQ(ierr);
ierr = writeGrid_ASCII(filePath); CHKERRQ(ierr);
}
else if (std::strcmp(viewerType, "hdf5") == 0)
{
ierr = writeSubMeshHDF5(filePath); CHKERRQ(ierr);
ierr = writeGrid_HDF5(filePath); CHKERRQ(ierr);
}
else
SETERRQ(comm, PETSC_ERR_ARG_UNKNOWN_TYPE,
"Unsupported PetscViewerType. Supported types are:\n"
"\t PETSCVIEWERASCII and PETSCVIEWERHDF5");

PetscFunctionReturn(0);
} // ProbeVolume::writeSubMesh
} // ProbeVolume::writeGrid

// Write the sub mesh into a HDF5 file.
PetscErrorCode ProbeVolume::writeSubMeshHDF5(const std::string &filePath)
PetscErrorCode ProbeVolume::writeGrid_HDF5(const std::string &filePath)
{
PetscErrorCode ierr;

Expand Down Expand Up @@ -373,10 +386,10 @@ PetscErrorCode ProbeVolume::writeSubMeshHDF5(const std::string &filePath)
}

PetscFunctionReturn(0);
} // ProbeVolume::writeSubMeshHDF5
} // ProbeVolume::writeGrid_HDF5

// Write the sub mesh into an ASCII file.
PetscErrorCode ProbeVolume::writeSubMeshASCII(const std::string &filePath)
PetscErrorCode ProbeVolume::writeGrid_ASCII(const std::string &filePath)
{
PetscErrorCode ierr;

Expand Down Expand Up @@ -412,7 +425,63 @@ PetscErrorCode ProbeVolume::writeSubMeshASCII(const std::string &filePath)
ierr = MPI_Barrier(comm); CHKERRQ(ierr);

PetscFunctionReturn(0);
} // ProbeVolume::writeSubMeshASCII
} // ProbeVolume::writeGrid_ASCII

// Write the index set (natural ordering) into a file.
PetscErrorCode ProbeVolume::writeIS(const std::string &filePath)
{
PetscErrorCode ierr;

PetscFunctionBeginUser;

// only ASCII and HDF5 formats are currently supported
if (std::strcmp(viewerType, "ascii") == 0)
{
ierr = writeIS_ASCII(filePath); CHKERRQ(ierr);
}
else if (std::strcmp(viewerType, "hdf5") == 0)
{
ierr = writeIS_HDF5(filePath); CHKERRQ(ierr);
}
else
SETERRQ(comm, PETSC_ERR_ARG_UNKNOWN_TYPE,
"Unsupported PetscViewerType. Supported types are:\n"
"\t PETSCVIEWERASCII and PETSCVIEWERHDF5");

PetscFunctionReturn(0);
} // ProbeVolume::writeIS

PetscErrorCode ProbeVolume::writeIS_HDF5(const std::string &filePath)
{
PetscErrorCode ierr;

PetscFunctionBeginUser;

ierr = PetscViewerFileSetMode(viewer, FILE_MODE_APPEND); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, filePath.c_str()); CHKERRQ(ierr);

ierr = PetscViewerHDF5PushGroup(viewer, "mesh"); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) isNatural, "IS"); CHKERRQ(ierr);
ierr = ISView(isNatural, viewer); CHKERRQ(ierr);

PetscFunctionReturn(0);
} // ProbeVolume::writeIS_HDF5

// Write index set (natural ordering) into HDF5 file.
PetscErrorCode ProbeVolume::writeIS_ASCII(const std::string &filePath)
{
PetscErrorCode ierr;

PetscFunctionBeginUser;

ierr = PetscViewerPushFormat(
viewer, PETSC_VIEWER_ASCII_SYMMODU); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) isNatural, "IS"); CHKERRQ(ierr);
ierr = ISView(isNatural, viewer); CHKERRQ(ierr);
ierr = PetscViewerPopFormat(viewer); CHKERRQ(ierr);

PetscFunctionReturn(0);
} // ProbeVolume::writeIS_ASCII

// Monitor a sub-region of the full-domain PETSc Vec object.
PetscErrorCode ProbeVolume::monitorVec(const DM &da,
Expand All @@ -426,7 +495,7 @@ PetscErrorCode ProbeVolume::monitorVec(const DM &da,

// grab the part of the vector that corresponds to the sub-volume
Vec svec;
ierr = VecGetSubVector(fvec, is, &svec); CHKERRQ(ierr);
ierr = VecGetSubVector(fvec, isPetsc, &svec); CHKERRQ(ierr);
if (n_sum != 0) // we accumulate the data over the time-steps
{
if (dvec == PETSC_NULL)
Expand All @@ -449,7 +518,7 @@ PetscErrorCode ProbeVolume::monitorVec(const DM &da,
{
ierr = writeVec(svec, t); CHKERRQ(ierr);
}
ierr = VecRestoreSubVector(fvec, is, &svec); CHKERRQ(ierr);
ierr = VecRestoreSubVector(fvec, isPetsc, &svec); CHKERRQ(ierr);

PetscFunctionReturn(0);
} // ProbeVolume::monitorVec
Expand All @@ -464,11 +533,11 @@ PetscErrorCode ProbeVolume::writeVec(const Vec &vec, const PetscReal &t)
// only ASCII and HDF5 formats are currently supported
if (std::strcmp(viewerType, "ascii") == 0)
{
ierr = writeVecASCII(vec, t); CHKERRQ(ierr);
ierr = writeVec_ASCII(vec, t); CHKERRQ(ierr);
}
else if (std::strcmp(viewerType, "hdf5") == 0)
{
ierr = writeVecHDF5(vec, t); CHKERRQ(ierr);
ierr = writeVec_HDF5(vec, t); CHKERRQ(ierr);
}
else
SETERRQ(comm, PETSC_ERR_ARG_UNKNOWN_TYPE,
Expand All @@ -479,7 +548,7 @@ PetscErrorCode ProbeVolume::writeVec(const Vec &vec, const PetscReal &t)
} // ProbeVolume::writeVec

// Write vector data data into a HDF5 file.
PetscErrorCode ProbeVolume::writeVecHDF5(const Vec &vec, const PetscReal &t)
PetscErrorCode ProbeVolume::writeVec_HDF5(const Vec &vec, const PetscReal &t)
{
PetscErrorCode ierr;

Expand All @@ -498,10 +567,10 @@ PetscErrorCode ProbeVolume::writeVecHDF5(const Vec &vec, const PetscReal &t)
}

PetscFunctionReturn(0);
} // ProbeVolume::writeVecHDF5
} // ProbeVolume::writeVec_HDF5

// Write vector data into an ASCII file.
PetscErrorCode ProbeVolume::writeVecASCII(const Vec &vec, const PetscReal &t)
PetscErrorCode ProbeVolume::writeVec_ASCII(const Vec &vec, const PetscReal &t)
{
PetscErrorCode ierr;

Expand All @@ -517,7 +586,7 @@ PetscErrorCode ProbeVolume::writeVecASCII(const Vec &vec, const PetscReal &t)
ierr = PetscViewerPopFormat(viewer); CHKERRQ(ierr);

PetscFunctionReturn(0);
} // ProbeVolume::writeVecASCII
} // ProbeVolume::writeVec_ASCII

//***************************************************************************//
//************************* ProbePoint *************************//
Expand Down