diff --git a/CMakeLists.txt b/CMakeLists.txt index f61876725..db6280fbd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -62,7 +62,7 @@ set(CMAKE_FIND_PACKAGE_SORT_DIRECTION DEC) find_package(MPI REQUIRED COMPONENTS C CXX) find_package(OpenMP) -find_package(FFTW) +find_package(FFTW REQUIRED) find_package(HDF5 COMPONENTS C CXX HL REQUIRED) find_package(TIRPC) # Check for alternative Sun rpc support find_package(Eigen3 REQUIRED) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 9625e4f2d..bcca340bf 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -368,7 +368,7 @@ namespace BasisClasses { // Sanity check on derived class type // - if (typeid(*coef) != typeid(CoefClasses::SphStruct)) + if (not dynamic_cast(coef.get())) throw std::runtime_error("Spherical::set_coefs: you must pass a CoefClasses::SphStruct"); // Sanity check on dimensionality @@ -1431,7 +1431,7 @@ namespace BasisClasses void Cylindrical::set_coefs(CoefClasses::CoefStrPtr coef) { - if (typeid(*coef) != typeid(CoefClasses::CylStruct)) + if (not dynamic_cast(coef.get())) throw std::runtime_error("Cylindrical::set_coefs: you must pass a CoefClasses::CylStruct"); CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); @@ -1708,7 +1708,7 @@ namespace BasisClasses { // Sanity check on derived class type // - if (typeid(*coef) != typeid(CoefClasses::CylStruct)) + if (not dynamic_cast(coef.get())) throw std::runtime_error("FlatDisk::set_coefs: you must pass a CoefClasses::CylStruct"); // Sanity check on dimensionality @@ -2140,7 +2140,7 @@ namespace BasisClasses { // Sanity check on derived class type // - if (typeid(*coef) != typeid(CoefClasses::SlabStruct)) + if (not dynamic_cast(coef.get())) throw std::runtime_error("Slab::set_coefs: you must pass a CoefClasses::SlabStruct"); // Sanity check on dimensionality @@ -2572,7 +2572,7 @@ namespace BasisClasses { // Sanity check on derived class type // - if (typeid(*coef) != typeid(CoefClasses::CubeStruct)) + if (not dynamic_cast(coef.get())) throw std::runtime_error("Cube::set_coefs: you must pass a CoefClasses::CubeStruct"); // Sanity check on dimensionality diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index d5deb13da..0d48b17fe 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -72,7 +72,7 @@ namespace MSSA unpack_cylfld(); } else { - throw std::runtime_error(std::string("CoefDB::unpack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); + throw std::runtime_error(std::string("CoefDB::unpack_channels(): can not reflect coefficient type=") + typeid(coefs.get()).name()); } return coefs; @@ -91,7 +91,7 @@ namespace MSSA else if (dynamic_cast(coefs.get())) { } // Do nothing else { - throw std::runtime_error(std::string("CoefDB::background(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); + throw std::runtime_error(std::string("CoefDB::background(): can not reflect coefficient type=") + typeid(coefs.get()).name()); } } @@ -112,7 +112,7 @@ namespace MSSA else if (dynamic_cast(coefs.get())) pack_cylfld(); else { - throw std::runtime_error(std::string("CoefDB::pack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); + throw std::runtime_error(std::string("CoefDB::pack_channels(): can not reflect coefficient type=") + typeid(coefs.get()).name()); } } 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/exputil/BarrierWrapper.cc b/exputil/BarrierWrapper.cc index 03755f62d..514aa4e05 100644 --- a/exputil/BarrierWrapper.cc +++ b/exputil/BarrierWrapper.cc @@ -101,12 +101,13 @@ void BarrierWrapper::light_operator(const string& label, // Compare adjacent strings in the list if (localid==0) { - char tmp[cbufsz]; // Working buffer + // Working buffer + auto tmp = std::make_unique(cbufsz); bool firstime = true; - string one, two = strncpy(tmp, &bufferT[0], cbufsz); + string one, two = strncpy(tmp.get(), &bufferT[0], cbufsz); for (int n=1; n(Number); for (int j=0; j #include #include +#include #include #include #include @@ -530,11 +531,11 @@ namespace PR { int rank = dataspace.getSimpleExtentNdims(); - hsize_t dims[rank]; + auto dims = std::make_unique(rank); - int ndims = dataspace.getSimpleExtentDims(dims, NULL); + int ndims = dataspace.getSimpleExtentDims(dims.get(), NULL); - H5::DataSpace mspace(rank, dims); + H5::DataSpace mspace(rank, dims.get()); std::vector masses(dims[0]); dataset.read(&masses[0], H5::PredType::NATIVE_FLOAT, mspace, dataspace ); @@ -569,11 +570,11 @@ namespace PR { int rank = dataspace.getSimpleExtentNdims(); - hsize_t dims[rank]; + auto dims = std::make_unique(rank); - int ndims = dataspace.getSimpleExtentDims( dims, NULL); + int ndims = dataspace.getSimpleExtentDims( dims.get(), NULL); - H5::DataSpace mspace(rank, dims); + H5::DataSpace mspace(rank, dims.get()); std::vector seq(dims[0]); dataset.read(&seq[0], H5::PredType::NATIVE_UINT32, mspace, dataspace ); diff --git a/exputil/TransformFFT.cc b/exputil/TransformFFT.cc index 08bfdac56..8b0bb54a7 100644 --- a/exputil/TransformFFT.cc +++ b/exputil/TransformFFT.cc @@ -12,9 +12,11 @@ TransformFFT::TransformFFT(double DR, std::vector& Y) dk = 2.0*M_PI/dr/N; in = Y; - out = new fftw_complex [N/2+1]; + out.resize(N/2+1); - p = fftw_plan_dft_r2c_1d(N, in.data(), out, FFTW_ESTIMATE); + p = fftw_plan_dft_r2c_1d(N, in.data(), + reinterpret_cast(out.data()), + FFTW_ESTIMATE); fftw_execute(p); } @@ -30,9 +32,11 @@ TransformFFT::TransformFFT(double DR, Eigen::VectorXd& Y) in.resize(N); for (int j=0; j(out.data()), + FFTW_ESTIMATE); fftw_execute(p); @@ -41,7 +45,6 @@ TransformFFT::TransformFFT(double DR, Eigen::VectorXd& Y) TransformFFT::~TransformFFT() { - delete [] out; fftw_destroy_plan(p); } @@ -53,15 +56,15 @@ void TransformFFT::Power(Eigen::VectorXd& F, Eigen::VectorXd& P) double d2 = dr*dr; F(0) = 0.0; - P(0) = d2 * (out[0][0]*out[0][0] + out[0][1]*out[0][1]); + P(0) = d2 * std::norm(out[0]); for (int j=1; j& F, std::vector& P) @@ -72,15 +75,15 @@ void TransformFFT::Power(std::vector& F, std::vector& P) double d2 = dr*dr; F[0] = 0.0; - P[0] = d2 * ( out[0][0]*out[0][0] + out[0][1]*out[0][1] ); + P[0] = d2 * std::norm(out[0]); for (int j=1; j(out[j][0], out[j][1]) * dr; + W[j] = out[j] * dr; } } @@ -105,8 +108,8 @@ void TransformFFT::Inverse(std::vector& F, for (int j=0; j data(count); data[0] = new char [11]; strcpy(data[0], "LoadConfig"); // Emulate the caller name diff --git a/include/EXPini.H b/include/EXPini.H index 79ac22e59..0690eea98 100644 --- a/include/EXPini.H +++ b/include/EXPini.H @@ -88,7 +88,7 @@ cxxopts::ParseResult LoadConfig(cxxopts::Options& options, YAML::Node conf = YAML::LoadFile(config); const int count = conf.size()*2+1; - char* data[count]; + std::vector data(count); int cnt = 1; data[0] = new char [11]; @@ -131,7 +131,7 @@ cxxopts::ParseResult LoadConfig(cxxopts::Options& options, auto vm = options.parse(cnt, &data[0]); - for (int i=0; i +#include + #include #include -#include -#include class TransformFFT { @@ -15,7 +16,7 @@ class TransformFFT fftw_plan p; std::vector in; - fftw_complex* out; + std::vector> out; public: diff --git a/pyEXP/FieldWrappers.cc b/pyEXP/FieldWrappers.cc index 5e0c1dd4e..165e58daf 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,6 +248,7 @@ void FieldGeneratorClasses(py::module &m) { See also -------- + 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 @@ -282,6 +335,10 @@ void FieldGeneratorClasses(py::module &m) { See also -------- +<<<<<<< HEAD +======= + points : generate fields at an array of mesh points +>>>>>>> pointMesh slices : generate fields in a surface slice given by the initializtion grid volumes : generate fields in volume given by the initializtion grid @@ -316,6 +373,10 @@ void FieldGeneratorClasses(py::module &m) { See also -------- +<<<<<<< HEAD +======= + points : generate fields at an array of mesh points +>>>>>>> pointMesh slices : generate fields in a surface slice given by the initializtion grid volumes : generate fields in volume given by the initializtion grid @@ -360,6 +421,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 )"); diff --git a/utils/Analysis/gas2dcyl.cc b/utils/Analysis/gas2dcyl.cc index f37d5aae7..1fd2e6722 100644 --- a/utils/Analysis/gas2dcyl.cc +++ b/utils/Analysis/gas2dcyl.cc @@ -171,9 +171,9 @@ main(int argc, char **argv) char *c = const_cast(files[n].c_str()); MPI_Bcast(c, sz+1, MPI_CHAR, 0, MPI_COMM_WORLD); } else { - char l[sz+1]; - MPI_Bcast(&l[0], sz+1, MPI_CHAR, 0, MPI_COMM_WORLD); - files.push_back(l); + auto l = std::make_unique(sz+1); + MPI_Bcast(l.get(), sz+1, MPI_CHAR, 0, MPI_COMM_WORLD); + files.push_back(l.get()); } MPI_Barrier(MPI_COMM_WORLD); } diff --git a/utils/Analysis/psphisto.cc b/utils/Analysis/psphisto.cc index 01a57c70f..cf04001b7 100644 --- a/utils/Analysis/psphisto.cc +++ b/utils/Analysis/psphisto.cc @@ -34,8 +34,9 @@ #include #include #include -#include +#include #include +#include // System libs #include @@ -177,9 +178,9 @@ main(int argc, char **argv) char *c = const_cast(files[n].c_str()); MPI_Bcast(c, sz+1, MPI_CHAR, 0, MPI_COMM_WORLD); } else { - char l[sz+1]; - MPI_Bcast(&l[0], sz+1, MPI_CHAR, 0, MPI_COMM_WORLD); - files.push_back(l); + auto l = std::make_unique(sz+1); + MPI_Bcast(l.get(), sz+1, MPI_CHAR, 0, MPI_COMM_WORLD); + files.push_back(l.get()); } MPI_Barrier(MPI_COMM_WORLD); } diff --git a/utils/ICs/SphericalSL.cc b/utils/ICs/SphericalSL.cc index 04171b51c..a4390b5c2 100644 --- a/utils/ICs/SphericalSL.cc +++ b/utils/ICs/SphericalSL.cc @@ -245,7 +245,7 @@ void SphericalSL::compute_coefficients_single(vector &part) void SphericalSL::compute_coefficients_thread(vector& part) { - std::thread t[nthrds]; + std::vector t(nthrds); // Launch the threads for (int id=0; id #ifdef HAVE_FFTW +#include #include #endif @@ -355,13 +356,21 @@ main(int argc, char **argv) #ifdef HAVE_FFTW if (SMOOTH>0.0) { - double a[NUMR], b[NUMR], c[NUMR]; - fftw_complex A[NUMR], B[NUMR], C[NUMR]; + std::vector a(NUMR), b(NUMR), c(NUMR); + std::vector> A(NUMR), B(NUMR), C(NUMR); fftw_plan pa, pb, pinv; - pa = fftw_plan_dft_r2c_1d(NUMR, a, A, FFTW_ESTIMATE); - pb = fftw_plan_dft_r2c_1d(NUMR, b, B, FFTW_ESTIMATE); - pinv = fftw_plan_dft_c2r_1d(NUMR, C, c, FFTW_ESTIMATE); + pa = fftw_plan_dft_r2c_1d(NUMR, a.data(), + reinterpret_cast(A.data()), + FFTW_ESTIMATE); + + pb = fftw_plan_dft_r2c_1d(NUMR, b.data(), + reinterpret_cast(B.data()), + FFTW_ESTIMATE); + + pinv = fftw_plan_dft_c2r_1d(NUMR, + reinterpret_cast(C.data()), + c.data(), FFTW_ESTIMATE); double xmin=rmin, xmax=rmax; @@ -370,7 +379,7 @@ main(int argc, char **argv) double scale = dk / NUMR; int indx; - ofstream tin, tout; + std::ofstream tin, tout; if (myid==0) { tin.open("ebar_fft.input"); @@ -403,8 +412,7 @@ main(int argc, char **argv) fftw_execute(pb); for (int i=0; i rr(NUMR), mm(NUMR); + std::vector rr(NUMR), mm(NUMR); double fac; for (int i=0; i(buf_size); std::vector or_time, or_c[3]; while (!in.eof()) { - in.getline(buf, buf_size); + in.getline(buf.get(), buf_size); if (in.eof()) break; - istringstream ins(buf); + istringstream ins(buf.get()); ins >> val; or_time.push_back(val); for (int k=0; k<5; k++) ins >> val;