Skip to content

Commit

Permalink
Implementation of a point set evaluator for FieldGenerator
Browse files Browse the repository at this point in the history
  • Loading branch information
Martin D. Weinberg committed Aug 7, 2024
1 parent 7a0e463 commit 76cdbd5
Show file tree
Hide file tree
Showing 3 changed files with 273 additions and 14 deletions.
23 changes: 22 additions & 1 deletion expui/FieldGenerator.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ namespace Field

std::vector<double> times, pmin, pmax;
std::vector<int> grid;
Eigen::MatrixXd mesh;

//! Sanity check time vector with coefficient DB
void check_times(CoefClasses::CoefsPtr coefs);
Expand All @@ -35,12 +36,32 @@ namespace Field

public:

//! Constructor
//! Constructor for a rectangular grid
FieldGenerator(const std::vector<double> &time,
const std::vector<double> &pmin,
const std::vector<double> &pmax,
const std::vector<int> &grid);

//! Constructor for an arbitrary point mesh. The mesh should be
//! an Nx3 array
FieldGenerator(const std::vector<double> &time,
const Eigen::MatrixXd &mesh);

/** Get field quantities at user defined points
For example:
.
.
// Generate the fields for all coefficients in 'coefs'
auto db = points(basis, coefs);
// Get fields evaluated at Time=3.14 and for all points in
// your supplied mesh for density ("dens")
Eigen::MatrixXf points = db["3.14"]["dens"];
*/
std::map<double, std::map<std::string, Eigen::VectorXf>>
points(BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs);

/** Get a field slices as a map in time and type
For example:
Expand Down
178 changes: 178 additions & 0 deletions expui/FieldGenerator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,33 @@ namespace Field
}
}

FieldGenerator::FieldGenerator(const std::vector<double> &time,
const Eigen::MatrixXd &mesh) :

times(time), mesh(mesh)
{
// Sanity check on mesh
//
if (mesh.cols() != 3)
throw std::runtime_error("FieldGenerator: bad mesh specification. The mesh must be an Nx3 array where the columns are Cartesian points");

// Check whether MPI is initialized
//
int flag;
MPI_Initialized(&flag);
if (flag) use_mpi = true;
else use_mpi = false;


// Fall back sanity (works for me but this needs to be fixed
// generally)
//
if (use_mpi) {
MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
}
}

void FieldGenerator::check_times(CoefClasses::CoefsPtr coefs)
{
std::vector<double> ctimes = coefs->Times();
Expand Down Expand Up @@ -867,5 +894,156 @@ namespace Field
}
// END histogram1d

std::map<double, std::map<std::string, Eigen::VectorXf>>
FieldGenerator::points(BasisClasses::BasisPtr basis,
CoefClasses::CoefsPtr coefs)
{
// Sanity check on mesh
//
if (mesh.size() == 0)
throw std::runtime_error("FieldGenerator::points: bad mesh specification. Did you call the mesh constructor?");

// Set midplane evaluation parameters
//
basis->setMidplane(midplane);
basis->setColumnHeight(colheight);

// Check
//
check_times(coefs);

std::map<double, std::map<std::string, Eigen::VectorXf>> ret;

// Now get the desired coordinate type
//
auto ctype = basis->coordinates;

// Field labels (force field labels added below)
//
auto labels = basis->getFieldLabels(ctype);

int ncnt = 0; // Process counter for MPI

std::map<std::string, Eigen::VectorXf> frame;
for (auto label : labels) {
frame[label].resize(mesh.rows());
}

for (auto T : times) {

if (ncnt++ % numprocs > 0) continue;

if (not coefs->getCoefStruct(T)) {
std::cout << "Could not find time=" << T << ", continuing" << std::endl;
continue;
}

basis->set_coefs(coefs->getCoefStruct(T));

#pragma omp parallel for
for (int k=0; k<mesh.rows(); k++) {

// Cartesian to spherical for all_eval
//
double x = mesh(k, 0);
double y = mesh(k, 1);
double z = mesh(k, 2);

// Coordinate values
double r, costh, phi, R;

// Return values
double p0, p1, d0, d1, f1, f2, f3;
std::vector<double> v;

if (ctype == BasisClasses::Basis::Coord::Spherical) {
r = sqrt(x*x + y*y + z*z) + 1.0e-18;
costh = z/r;
phi = atan2(y, x);
v = (*basis)(r, costh, phi, ctype);
} else if (ctype == BasisClasses::Basis::Coord::Cylindrical) {
R = sqrt(x*x + y*y) + 1.0e-18;
phi = atan2(y, x);
v = (*basis)(R, z, phi, ctype);
} else {
v = (*basis)(x, y, z, BasisClasses::Basis::Coord::Cartesian);
}

// Pack the frame structure
//
for (int n=0; n<labels.size(); n++)
frame[labels[n]](k) = v[n];
}

ret[T] = frame;
}

if (use_mpi) {

std::vector<char> bf(9);

for (int n=1; n<numprocs; n++) {

if (myid==n) {
int sz = ret.size();
MPI_Send(&sz, 1, MPI_INT, 0, 102, MPI_COMM_WORLD);
for (auto & v : ret) {
MPI_Send(&v.first, 1, MPI_DOUBLE, 0, 103, MPI_COMM_WORLD);
int fsz = v.second.size();
MPI_Send(&fsz, 1, MPI_INT, 0, 104, MPI_COMM_WORLD);

for (auto & f : v.second) {
MPI_Send(f.first.c_str(), f.first.length(), MPI_CHAR, 0, 105,
MPI_COMM_WORLD);

MPI_Send(f.second.data(), f.second.size(),
MPI_FLOAT, 0, 106, MPI_COMM_WORLD);
}
}
}

if (myid==0) {
MPI_Status status;
int sz, fsz, l;
double T;

MPI_Recv(&sz, 1, MPI_INT, n, 102, MPI_COMM_WORLD, &status);
for (int i=0; i<sz; i++) {
MPI_Recv(&T, 1, MPI_DOUBLE, n, 103, MPI_COMM_WORLD, &status);
MPI_Recv(&fsz, 1, MPI_INT, n, 104, MPI_COMM_WORLD, &status);

for (int j=0; j<fsz; j++) {
// Get the field name
//
MPI_Probe(n, 105, MPI_COMM_WORLD, &status);
MPI_Get_count(&status, MPI_CHAR, &l);
MPI_Recv(bf.data(), l, MPI_CHAR, n, 105, MPI_COMM_WORLD, &status);
std::string s(bf.data(), l);

// Sanity check
//
if (frame.find(s) == frame.end()) {
std::cerr << "Error finding <" << s << "> in field map"
<< std::endl;
}

// Get the data
//
MPI_Recv(frame[s].data(), frame[s].size(), MPI_FLOAT, n, 106,
MPI_COMM_WORLD, &status);
}

ret[T] = frame;
}
}
}
}

// Toggle off midplane evaluation
basis->setMidplane(false);

return ret;
}

}
// END namespace Field
86 changes: 73 additions & 13 deletions pyEXP/FieldWrappers.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,22 +24,24 @@ void FieldGeneratorClasses(py::module &m) {
a list of upper bounds, and a list of knots per dimension. These
lists all have rank 3 for (x, y, z). For a two-dimensional surface,
one of the knot array values must be zero. The member functions
lines, slices and volumes, called with the basis and coefficient
objects, return a numpy.ndarray containing the field evaluations.
Each of these functions returns a dictionary of times to a dictionary
of field names to numpy.ndarrays at each time. There are also members
which will write these generated fields to files. The linear probe
members, 'lines' and 'file_lines', evaluate 'num' field points along
a user-specified segment between the 3d points 'beg' and 'end'. See
help(pyEXP.basis) and help(pyEXP.coefs) for info on the basis and
coefficient objects.
lines, slices, an arbitrary point-set, and volumes, called with the
basis and coefficient objects, return a numpy.ndarray containing the
field evaluations. Each of these functions returns a dictionary of
times to a dictionary of field names to numpy.ndarrays at each time.
There are also members which will write these generated fields to files.
The linear probe members, 'lines' and 'file_lines', evaluate 'num'
field points along a user-specified segment between the 3d points 'beg'
and 'end'. See help(pyEXP.basis) and help(pyEXP.coefs) for info on
the basis and coefficient objects.
Data packing
------------
All slices and volumes are returned as numpy.ndarrays in row-major
order. That is, the first index is x, the second is y, and the
third is z. The ranks are specified by the 'gridsize' array with
(nx, ny, nz) as input to the FieldGenerator constructor.
(nx, ny, nz) as input to the FieldGenerator constructor. A point
mesh is rows of (x, y, z) Cartesian coordinates. The output field
values returned by points is in the same order as the input array.
Coordinate systems
------------------
Expand Down Expand Up @@ -116,6 +118,24 @@ void FieldGeneratorClasses(py::module &m) {
py::arg("times"), py::arg("lower"), py::arg("upper"),
py::arg("gridsize"));

f.def(py::init<const std::vector<double>, const Eigen::MatrixXd>(),
R"(
Create fields for given times and and provided point set
Parameters
----------
times : list(float,...)
list of evaluation times
mesh : ndarray of Nx3 floats
point set (x, y, z) for field evaluation
Returns
-------
FieldGenerator
new object
)",
py::arg("times"), py::arg("mesh"));

f.def("setMidplane", &Field::FieldGenerator::setMidplane,
R"(
Set the field generator to generate midplane fields
Expand Down Expand Up @@ -160,6 +180,38 @@ void FieldGeneratorClasses(py::module &m) {
See also
--------
points : generate fields at an array of mesh points
lines : generate fields along a line given by its end points
volumes : generate fields in volume given by the initializtion grid
)", py::arg("basis"), py::arg("coefs"));

f.def("points", &Field::FieldGenerator::points,
R"(
Return a dictionary of arrays (1d numpy arrays) indexed by time
and field type corresponding to the mesh points
Parameters
----------
basis : Basis
basis instance of any geometry; geometry will be deduced by the generator
coefs : Coefs
coefficient container instance
Returns
-------
dict({time: {field-name: numpy.ndarray})
dictionary of times, field names, and data arrays
Notes
-----
Each data array in the dictionary is a 1-dimensional array with the
same order as the mesh points in the constructor
routines.
See also
--------
slices : generate fields in a surface slice given by the
initializtion grid
lines : generate fields along a line given by its end points
volumes : generate fields in volume given by the initializtion grid
)", py::arg("basis"), py::arg("coefs"));
Expand Down Expand Up @@ -196,7 +248,9 @@ void FieldGeneratorClasses(py::module &m) {
See also
--------
slices : generate fields in a surface slice given by the initializtion grid
points : generate fields at an array of mesh points
slices : generate fields in a surface slice given by the
initializtion grid
volumes : generate fields in volume given by the initializtion grid
)",
py::arg("basis"), py::arg("coefs"),
Expand Down Expand Up @@ -281,7 +335,9 @@ void FieldGeneratorClasses(py::module &m) {
See also
--------
slices : generate fields in a surface slice given by the initializtion grid
points : generate fields at an array of mesh points
slices : generate fields in a surface slice given by the
initializtion grid
volumes : generate fields in volume given by the initializtion grid
lines : generate fields along a line given by its end points
)",
Expand Down Expand Up @@ -314,7 +370,9 @@ void FieldGeneratorClasses(py::module &m) {
See also
--------
slices : generate fields in a surface slice given by the initializtion grid
points : generate fields at an array of mesh points
slices : generate fields in a surface slice given by the
initializtion grid
volumes : generate fields in volume given by the initializtion grid
lines : generate fields along a line given by its end points
)",
Expand Down Expand Up @@ -357,6 +415,8 @@ void FieldGeneratorClasses(py::module &m) {
See also
--------
points : generate fields at an array of mesh points
lines : generate fields along a line given by its end points
slices : generate fields in a slice given by the initializtion grid
)");
Expand Down

0 comments on commit 76cdbd5

Please sign in to comment.