From 76cdbd551708481e4a4805e547b82637bbc1d5f4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 7 Aug 2024 17:10:57 -0400 Subject: [PATCH] Implementation of a point set evaluator for FieldGenerator --- expui/FieldGenerator.H | 23 +++++- expui/FieldGenerator.cc | 178 ++++++++++++++++++++++++++++++++++++++++ pyEXP/FieldWrappers.cc | 86 ++++++++++++++++--- 3 files changed, 273 insertions(+), 14 deletions(-) diff --git a/expui/FieldGenerator.H b/expui/FieldGenerator.H index 3b4873ade..4ab573aef 100644 --- a/expui/FieldGenerator.H +++ b/expui/FieldGenerator.H @@ -20,6 +20,7 @@ namespace Field std::vector times, pmin, pmax; std::vector grid; + Eigen::MatrixXd mesh; //! Sanity check time vector with coefficient DB void check_times(CoefClasses::CoefsPtr coefs); @@ -35,12 +36,32 @@ namespace Field public: - //! Constructor + //! Constructor for a rectangular grid FieldGenerator(const std::vector &time, const std::vector &pmin, const std::vector &pmax, const std::vector &grid); + //! Constructor for an arbitrary point mesh. The mesh should be + //! an Nx3 array + FieldGenerator(const std::vector &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> + points(BasisClasses::BasisPtr basis, CoefClasses::CoefsPtr coefs); + /** Get a field slices as a map in time and type For example: diff --git a/expui/FieldGenerator.cc b/expui/FieldGenerator.cc index 5e00e1c82..73a1595f7 100644 --- a/expui/FieldGenerator.cc +++ b/expui/FieldGenerator.cc @@ -41,6 +41,33 @@ namespace Field } } + FieldGenerator::FieldGenerator(const std::vector &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 ctimes = coefs->Times(); @@ -867,5 +894,156 @@ namespace Field } // END histogram1d + std::map> + 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> 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 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 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 bf(9); + + for (int n=1; n 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 diff --git a/pyEXP/FieldWrappers.cc b/pyEXP/FieldWrappers.cc index fb79bef5c..995f91358 100644 --- a/pyEXP/FieldWrappers.cc +++ b/pyEXP/FieldWrappers.cc @@ -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 ------------------ @@ -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 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 @@ -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")); @@ -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"), @@ -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 )", @@ -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 )", @@ -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 )");