diff --git a/CMakeLists.txt b/CMakeLists.txt index c639586b7..0c5b629f3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.15) # This is the oldest version I've # preferred esp. for Cuda support project( EXP - VERSION "7.7.25" + VERSION "7.7.26" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) @@ -47,6 +47,7 @@ option(ENABLE_VTK "Configure VTK if available" FALSE) option(ENABLE_CUDA_SINGLE "Use real*4 instead of real*8 for CUDA" FALSE) option(ENABLE_DSMC "Enable DSCM module" FALSE) option(ENABLE_USER "Enable basic user modules" ON) +option(ENABLE_SLCHECK "Enable *careful* Sturm-Liouville solutions" TRUE) option(BUILD_SHARED_LIBS "Build using shared libraries" ON) option(BUILD_DOCS "Build documentation" OFF) @@ -146,6 +147,9 @@ endif() if(FFTW_FOUND) set(HAVE_FFTW TRUE) endif() +if(ENABLE_SLCHECK) + set(SLEDGE_THROW TRUE) +endif() if(PNG_FOUND AND ENABLE_PNG) set(HAVE_LIBPNGPP TRUE) endif() diff --git a/coefs/BasisFactory.H b/coefs/BasisFactory.H index 2b6f8f9c9..17a88c1be 100644 --- a/coefs/BasisFactory.H +++ b/coefs/BasisFactory.H @@ -101,6 +101,12 @@ namespace BasisClasses //! Using MPI bool use_mpi; + //! Return readable class name + virtual const std::string classname() = 0; + + //! Subspace index + virtual const std::string harmonic() = 0; + public: //! Constructor from YAML node @@ -198,7 +204,7 @@ namespace BasisClasses std::string model_file; int lmax, nmax, cmap, numr; - double rmin, rmax; + double rmin, rmax, rmap; bool NO_L0, NO_L1, EVEN_L, EVEN_M, M0_only; @@ -240,6 +246,12 @@ namespace BasisClasses //! Valid keys for YAML configurations static const std::set valid_keys; + //! Return readable class name + virtual const std::string classname() { return "SphericalSL";} + + //! Subspace index + virtual const std::string harmonic() { return "l";} + public: //! Constructor from YAML node @@ -291,7 +303,10 @@ namespace BasisClasses //! Compute the orthogonality of the basis by returning inner //! produce matrices - std::vector orthoCheck(int knots=40); + std::vector orthoCheck(int knots=40) + { + return sl->orthoCheck(knots); + } }; @@ -363,6 +378,12 @@ namespace BasisClasses //! Valid keys for YAML configurations static const std::set valid_keys; + //! Return readable class name + virtual const std::string classname() { return "FlatDisk";} + + //! Subspace index + virtual const std::string harmonic() { return "m";} + public: //! Constructor from YAML node @@ -493,6 +514,12 @@ namespace BasisClasses //! Valid keys for YAML configurations static const std::set valid_keys; + //! Return readable class name + virtual const std::string classname() { return "Cylindrical";} + + //! Subspace index + virtual const std::string harmonic() { return "m";} + public: //! Constructor from YAML node diff --git a/coefs/BasisFactory.cc b/coefs/BasisFactory.cc index 6b935bf23..48bb5817f 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -2,6 +2,7 @@ #include #include #include +#include #include #ifdef HAVE_FE_ENABLE @@ -92,7 +93,7 @@ namespace BasisClasses const std::set SphericalSL::valid_keys = { - "rs", + "rmapping", "cmap", "Lmax", "dof", @@ -172,7 +173,7 @@ namespace BasisClasses // Assign values from YAML // - double rs = 1.0; + double rmap = 1.0; try { if (conf["cmap"]) cmap = conf["cmap"].as(); @@ -180,10 +181,10 @@ namespace BasisClasses if (conf["nmax"]) nmax = conf["nmax"].as(); if (conf["modelname"]) model_file = conf["modelname"].as(); - if (conf["rs"]) - rs = conf["rs"].as(); + if (conf["rmapping"]) + rmap = conf["rmapping"].as(); else - rs = 1.0; + rmap = 1.0; if (conf["scale"]) scale = conf["scale"].as(); @@ -257,9 +258,12 @@ namespace BasisClasses // Finally, make the Sturm-Lioville basis... sl = std::make_shared - (model_file, lmax, nmax, numr, rmin, rmax, true, cmap, rs, + (model_file, lmax, nmax, numr, rmin, rmax, true, cmap, rmap, 0, 1, cachename); + // Test basis for consistency + orthoTest(orthoCheck(std::max(nmax*50, 200)), classname(), harmonic()); + // Number of possible threads int nthrds = omp_get_max_threads(); @@ -610,58 +614,6 @@ namespace BasisClasses } - std::vector SphericalSL::orthoCheck(int num) - { - // Gauss-Legendre knots and weights - LegeQuad lw(num); - - // Get the scaled coordinate limits - double ximin = sl->r_to_xi(rmin); - double ximax = sl->r_to_xi(rmax); - - // Initialize the return matrices - std::vector ret(lmax+1); - for (auto & v : ret) v.resize(nmax, nmax); - - // Do each harmonic order - for (int L=0; L<=lmax; L++) { - - // Unroll the loop for OpenMP parallelization -#pragma omp parallel for - for (int nn=0; nnxi_to_r(x); - - ans += r*r*sl->get_pot(x, L, n1, 0)* - sl->get_dens(x, L, n2, 0) / - sl->d_xi_to_r(x) * (ximax - ximin)*lw.weight(i); - - } - // END: inner product - - // Assign the matrix element - // - ret[L](n1, n2) = - ans; - // ^ - // | - // +--- Switch to normed scalar product rather - // that normed gravitational energy - } - // END: unrolled loop - } - // END: harmonic order loop - - return ret; - } - - SphericalSL::BasisArray SphericalSL::getBasis (double logxmin, double logxmax, int numgrid) { @@ -1009,14 +961,14 @@ namespace BasisClasses ncylnx = 256; ncylny = 128; ncylodd = 9; - ncylr = 200; + ncylr = 2000; eof_file = ".eof_cache_file"; Ignore = false; deproject = false; - rnum = 100; + rnum = 200; pnum = 1; - tnum = 40; + tnum = 80; ashift = 0.0; logarithmic = false; density = true; @@ -1089,7 +1041,6 @@ namespace BasisClasses if (conf["ashift" ]) ashift = conf["ashift" ].as(); if (conf["logr" ]) logarithmic = conf["logr" ].as(); - if (conf["density" ]) density = conf["density" ].as(); if (conf["EVEN_M" ]) EVEN_M = conf["EVEN_M" ].as(); if (conf["cmapr" ]) cmapR = conf["cmapr" ].as(); if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); @@ -1109,6 +1060,14 @@ namespace BasisClasses if (conf["dtype" ]) dtype = conf["dtype" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); + // Deprecation warning + if (conf["density" ]) { + if (myid==0) + std::cout << "Cylindrical: parameter 'density' is deprecated. " + << "The density field will be computed regardless." + << std::endl; + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing 'force' for Component <" @@ -1139,11 +1098,6 @@ namespace BasisClasses EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; - // For visualization - // - if (density) EmpCylSL::DENS = true; - - // Check for non-null cache file name. This must be specified // to prevent recomputation and unexpected behavior. // @@ -1275,7 +1229,12 @@ namespace BasisClasses sl->generate_eof(rnum, pnum, tnum, f); } + + // Orthogonality sanity check + // + orthoTest(orthoCheck(), classname(), harmonic()); } + void Cylindrical::getFields (double x, double y, double z, @@ -1553,6 +1512,10 @@ namespace BasisClasses // ortho = std::make_shared(conf); + // Orthogonality sanity check + // + orthoTest(orthoCheck(), classname(), harmonic()); + // Get max threads // int nthrds = omp_get_max_threads(); @@ -1881,6 +1844,7 @@ namespace BasisClasses return ortho->orthoCheck(); } + FlatDisk::BasisArray FlatDisk::getBasis (double logxmin, double logxmax, int numgrid) { diff --git a/config_cmake.h_in b/config_cmake.h_in index ab7e61f39..d1cf00db3 100644 --- a/config_cmake.h_in +++ b/config_cmake.h_in @@ -1,6 +1,9 @@ /* True if DSMC is enabled */ #cmakedefine DSMC_ENABLED @DSMC_ENABLED@ +/* Defined if ENABLE_SLCHECK is enabled */ +#cmakedefine SLEDGE_THROW @ENABLE_SLCHECK@ + /* "Have C++0x (set standard to cxx14 so legacy define)*/ #define HAVE_CXX0X diff --git a/doc/exp.cfg b/doc/exp.cfg index 0ecf59445..962a56ad3 100644 --- a/doc/exp.cfg +++ b/doc/exp.cfg @@ -38,7 +38,7 @@ PROJECT_NAME = "EXP" # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 7.7.24 +PROJECT_NUMBER = 7.7.26 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index b4bbeb7f3..006c8dce9 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -8,7 +8,7 @@ set(UTIL_SRC nrutil.cc elemfunc.cc euler.cc euler_slater.cc # Hankel.cc rotmatrix.cc wordSplit.cc FileUtils.cc BarrierWrapper.cc stack.cc localmpi.cc TableGrid.cc writePVD.cc libvars.cc TransformFFT.cc QDHT.cc YamlCheck.cc parseVersionString.cc EXPmath.cc laguerre_polynomial.cpp - YamlConfig.cc) + YamlConfig.cc orthoTest.cc) if(HAVE_VTK) list(APPEND UTIL_SRC VtkGrid.cc VtkPCA.cc) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 2e0c9514b..55fa4a0c0 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -48,7 +48,6 @@ using namespace __EXP__; // For reference to n-body globals #define TINY 1.0e-16 -bool EmpCylSL::DENS = false; bool EmpCylSL::PCAVAR = false; bool EmpCylSL::PCAVTK = false; bool EmpCylSL::PCAEOF = false; @@ -121,10 +120,7 @@ EmpCylSL::EmpCylSL() // Choose table dimension // - if (DENS) - MPItable = 4; - else - MPItable = 3; + MPItable = 4; // Initialize values // @@ -202,10 +198,7 @@ EmpCylSL::EmpCylSL(int nmax, int lmax, int mmax, int nord, // Choose table dimension // - if (DENS) - MPItable = 4; - else - MPItable = 3; + MPItable = 4; // Initialize storage and values // @@ -289,10 +282,7 @@ void EmpCylSL::reset(int numr, int lmax, int mmax, int nord, // Choose table dimension // - if (DENS) - MPItable = 4; - else - MPItable = 3; + MPItable = 4; // Initialize data and storage sampT = 1; @@ -607,12 +597,10 @@ void EmpCylSL::send_eof_grid() MPI_DOUBLE, 0, MPI_COMM_WORLD); blab("after", "zforceC", myid, m, v); - if (DENS) { - blab("before", "densC", myid, m, v); - MPI_Bcast(densC[m][v].data(), densC[m][v].size(), - MPI_DOUBLE, 0, MPI_COMM_WORLD); - blab("after", "densC", myid, m, v); - } + blab("before", "densC", myid, m, v); + MPI_Bcast(densC[m][v].data(), densC[m][v].size(), + MPI_DOUBLE, 0, MPI_COMM_WORLD); + blab("after", "densC", myid, m, v); } } @@ -637,13 +625,10 @@ void EmpCylSL::send_eof_grid() MPI_DOUBLE, 0, MPI_COMM_WORLD); blab("after", "zforceS", myid, m, v); - if (DENS) { - blab("before", "densS", myid, m, v); - MPI_Bcast(densS[m][v].data(), densS[m][v].size(), - MPI_DOUBLE, 0, MPI_COMM_WORLD); - blab("after", "densS", myid, m, v); - } - + blab("before", "densS", myid, m, v); + MPI_Bcast(densS[m][v].data(), densS[m][v].size(), + MPI_DOUBLE, 0, MPI_COMM_WORLD); + blab("after", "densS", myid, m, v); } } @@ -697,7 +682,6 @@ int EmpCylSL::read_eof_header(const std::string& eof_file) NUMY = node["numy" ].as(); NMAX = node["nmaxfid"].as(); NORDER = node["nmax" ].as(); - DENS = node["dens" ].as(); RMIN = node["rmin" ].as(); RMAX = node["rmax" ].as(); ASCALE = node["ascl" ].as(); @@ -727,7 +711,6 @@ int EmpCylSL::read_eof_header(const std::string& eof_file) in.read((char *)&NMAX, sizeof(int)); in.read((char *)&NORDER, sizeof(int)); in.read((char *)&tmp, sizeof(int)); - if (tmp) DENS = true; else DENS = false; in.read((char *)&CMAPR, sizeof(int)); in.read((char *)&RMIN, sizeof(double)); in.read((char *)&RMAX, sizeof(double)); @@ -744,7 +727,6 @@ int EmpCylSL::read_eof_header(const std::string& eof_file) cout << "NUMY=" << NUMY << endl; cout << "NMAX=" << NMAX << endl; cout << "NORDER=" << NORDER << endl; - cout << "DENS=" << DENS << endl; cout << "CMAPR=" << CMAPR << endl; cout << "CMAPZ=" << CMAPZ << endl; cout << "RMIN=" << RMIN << endl; @@ -884,14 +866,14 @@ int EmpCylSL::cache_grid(int readwrite, string cachename) if (ReadH5Cache()) return 1; else { if (myid==0) - std::cout << "---- EmpCylSL::cache_grid: HDF5 parameter mismatch" + std::cerr << "---- EmpCylSL::cache_grid: HDF5 parameter mismatch" << std::endl; return 0; } } catch (EXPException& error) { if (myid==0) { - std::cout << "---- EmpCylSL::cache_grid: file <" + std::cerr << "---- EmpCylSL::cache_grid: file <" << cachefile + ">, is not HDF5. Checking binary . . ." << std::endl; } @@ -997,7 +979,6 @@ int EmpCylSL::cache_grid(int readwrite, string cachename) (NUMY != numy ) | (NMAX != nmax ) | (NORDER != norder ) | - (DENS != dens ) | (CMAPR != cmapr ) | (fabs(rmin-RMIN)>1.0e-12 ) | (fabs(rmax-RMAX)>1.0e-12 ) | @@ -1019,7 +1000,6 @@ int EmpCylSL::cache_grid(int readwrite, string cachename) cout << compare_out("numy", NUMY, numy); cout << compare_out("nmax", NMAX, nmax); cout << compare_out("norder", NORDER, norder); - cout << compare_out("dens", DENS, dens); cout << compare_out("cmapr", CMAPR, cmapr); cout << compare_out("cmapz", CMAPZ, cmapz); cout << compare_out("rmin", RMIN, rmin); @@ -1050,12 +1030,9 @@ int EmpCylSL::cache_grid(int readwrite, string cachename) for (int iy=0; iy<=NUMY; iy++) in.read((char *)&zforceC[m][v](ix, iy), sizeof(double)); - if (DENS) { - for (int ix=0; ix<=NUMX; ix++) - for (int iy=0; iy<=NUMY; iy++) - in.read((char *)&densC[m][v](ix, iy), sizeof(double)); - - } + for (int ix=0; ix<=NUMX; ix++) + for (int iy=0; iy<=NUMY; iy++) + in.read((char *)&densC[m][v](ix, iy), sizeof(double)); } } @@ -1076,11 +1053,9 @@ int EmpCylSL::cache_grid(int readwrite, string cachename) for (int iy=0; iy<=NUMY; iy++) in.read((char *)&zforceS[m][v](ix, iy), sizeof(double)); - if (DENS) { - for (int ix=0; ix<=NUMX; ix++) - for (int iy=0; iy<=NUMY; iy++) - in.read((char *)&densS[m][v](ix, iy), sizeof(double)); - } + for (int ix=0; ix<=NUMX; ix++) + for (int iy=0; iy<=NUMY; iy++) + in.read((char *)&densS[m][v](ix, iy), sizeof(double)); } } @@ -1155,9 +1130,6 @@ YAML::Node EmpCylSL::getHeader_hdf5(const std::string& cachefile) return v; }; - bool dens = false; - if (getInt("idens")) dens = true; - node["mmax"] = getInt("mmax"); node["numx"] = getInt("numx"); node["numy"] = getInt("numy"); @@ -1165,7 +1137,6 @@ YAML::Node EmpCylSL::getHeader_hdf5(const std::string& cachefile) node["nmax"] = getInt("nmax"); node["neven"] = getInt("neven"); node["nodd"] = getInt("nodd"); - node["dens"] = dens; node["cmapr"] = getInt("cmapr"); node["cmapz"] = getInt("cmapz"); node["rmin"] = getDbl("rmin"); @@ -1269,10 +1240,10 @@ void EmpCylSL::receive_eof(int request_id, int MM) MPI_Recv(&mpi_double_buf2[MPIbufsz*(MPItable*n+2)], MPIbufsz, MPI_DOUBLE, current_source, 13 + MPItable*n+3, MPI_COMM_WORLD, &status); - if (DENS) - MPI_Recv(&mpi_double_buf2[MPIbufsz*(MPItable*n+3)], - MPIbufsz, MPI_DOUBLE, current_source, 13 + MPItable*n+4, - MPI_COMM_WORLD, &status); + + MPI_Recv(&mpi_double_buf2[MPIbufsz*(MPItable*n+3)], + MPIbufsz, MPI_DOUBLE, current_source, 13 + MPItable*n+4, + MPI_COMM_WORLD, &status); } @@ -1320,16 +1291,14 @@ void EmpCylSL::receive_eof(int request_id, int MM) zforceS[mm][n](ix, iy) = mpi_double_buf2[off+icnt++]; - if (DENS) { - off = MPIbufsz*(MPItable*n+3); - icnt = 0; - for (int ix=0; ix<=NUMX; ix++) - for (int iy=0; iy<=NUMY; iy++) - if (type) - densC[mm][n](ix, iy) = mpi_double_buf2[off+icnt++]; - else - densS[mm][n](ix, iy) = mpi_double_buf2[off+icnt++]; - } + off = MPIbufsz*(MPItable*n+3); + icnt = 0; + for (int ix=0; ix<=NUMX; ix++) + for (int iy=0; iy<=NUMY; iy++) + if (type) + densC[mm][n](ix, iy) = mpi_double_buf2[off+icnt++]; + else + densS[mm][n](ix, iy) = mpi_double_buf2[off+icnt++]; } if (VFLAG & 16) @@ -1359,7 +1328,7 @@ void EmpCylSL::compute_eof_grid(int request_id, int m) tpot[v].setZero(); trforce[v].setZero(); tzforce[v].setZero(); - if (DENS) tdens[v].setZero(); + tdens[v].setZero(); } for (int ix=0; ix<=NUMX; ix++) { @@ -1376,7 +1345,7 @@ void EmpCylSL::compute_eof_grid(int request_id, int m) ortho->get_pot(potd, rr/ASCALE); ortho->get_force(dpot, rr/ASCALE); - if (DENS) ortho->get_dens(dend, rr/ASCALE); + ortho->get_dens(dend, rr/ASCALE); double costh = z/rr; dlegendre_R(LMAX, costh, legs[0], dlegs[0]); @@ -1417,8 +1386,7 @@ void EmpCylSL::compute_eof_grid(int request_id, int m) tzforce[v](ix, iy) += -ef(nn, v) * (potr*z/rr + pott*r*r/(rr*rr*rr)); - if (DENS) - tdens[v](ix, iy) += ef(nn, v) * dens * 0.25/M_PI; + tdens[v](ix, iy) += ef(nn, v) * dens * 0.25/M_PI; } } } @@ -1483,22 +1451,20 @@ void EmpCylSL::compute_eof_grid(int request_id, int m) // Density // - if (DENS) { - off = MPIbufsz*(MPItable*n+3); - icnt = 0; - for (int ix=0; ix<=NUMX; ix++) - for (int iy=0; iy<=NUMY; iy++) - mpi_double_buf2[off + icnt++] = tdens[n](ix, iy); - - if (VFLAG & 16) - std::cerr << "Worker " << setw(4) << myid - << ": with request_id=" << request_id - << ", M=" << m << " sending Density" << std::endl; + off = MPIbufsz*(MPItable*n+3); + icnt = 0; + for (int ix=0; ix<=NUMX; ix++) + for (int iy=0; iy<=NUMY; iy++) + mpi_double_buf2[off + icnt++] = tdens[n](ix, iy); + + if (VFLAG & 16) + std::cerr << "Worker " << setw(4) << myid + << ": with request_id=" << request_id + << ", M=" << m << " sending Density" << std::endl; - MPI_Send(&mpi_double_buf2[off], MPIbufsz, MPI_DOUBLE, 0, - 13 + MPItable*n+4, MPI_COMM_WORLD); + MPI_Send(&mpi_double_buf2[off], MPIbufsz, MPI_DOUBLE, 0, + 13 + MPItable*n+4, MPI_COMM_WORLD); - } } } // END: MPI packing @@ -1511,12 +1477,12 @@ void EmpCylSL::compute_eof_grid(int request_id, int m) potC [m][n] = tpot[n]; rforceC[m][n] = trforce[n]; zforceC[m][n] = tzforce[n]; - if (DENS) densC[m][n] = tdens[n]; + densC [m][n] = tdens[n]; } else { potS [m][n] = tpot[n]; rforceS[m][n] = trforce[n]; zforceS[m][n] = trforce[n]; - if (DENS) densS[m][n] = tdens[n]; + densS [m][n] = tdens[n]; } } } @@ -1539,7 +1505,7 @@ void EmpCylSL::compute_even_odd(int request_id, int m) tpot[v].setZero(); trforce[v].setZero(); tzforce[v].setZero(); - if (DENS) tdens[v].setZero(); + tdens[v].setZero(); } for (int ix=0; ix<=NUMX; ix++) { @@ -1556,7 +1522,7 @@ void EmpCylSL::compute_even_odd(int request_id, int m) ortho->get_pot(potd, rr/ASCALE); ortho->get_force(dpot, rr/ASCALE); - if (DENS) ortho->get_dens(dend, rr/ASCALE); + ortho->get_dens(dend, rr/ASCALE); double costh = z/rr; dlegendre_R(LMAX, costh, legs[0], dlegs[0]); @@ -1601,8 +1567,7 @@ void EmpCylSL::compute_even_odd(int request_id, int m) tzforce[v](ix, iy) += -efE(nn, v) * (potr*z/rr + pott*r*r/(rr*rr*rr)); - if (DENS) - tdens[v](ix, iy) += efE(nn, v) * dens; + tdens[v](ix, iy) += efE(nn, v) * dens * 0.25/M_PI; } } } @@ -1647,8 +1612,7 @@ void EmpCylSL::compute_even_odd(int request_id, int m) tzforce[v](ix, iy) += -efO(nn, w) * (potr*z/rr + pott*r*r/(rr*rr*rr)); - if (DENS) - tdens[v](ix, iy) += efO(nn, w) * dens * 0.25/M_PI; + tdens[v](ix, iy) += efO(nn, w) * dens * 0.25/M_PI; } } } @@ -1714,21 +1678,19 @@ void EmpCylSL::compute_even_odd(int request_id, int m) // Density // - if (DENS) { - off = MPIbufsz*(MPItable*n+3); - icnt = 0; - for (int ix=0; ix<=NUMX; ix++) - for (int iy=0; iy<=NUMY; iy++) - mpi_double_buf2[off + icnt++] = tdens[n](ix, iy); - - if (VFLAG & 16) - std::cerr << "Worker " << setw(4) << myid - << ": with request_id=" << request_id - << ", M=" << m << " sending Density" << std::endl; - - MPI_Send(&mpi_double_buf2[off], MPIbufsz, MPI_DOUBLE, 0, - 13 + MPItable*n+4, MPI_COMM_WORLD); - } + off = MPIbufsz*(MPItable*n+3); + icnt = 0; + for (int ix=0; ix<=NUMX; ix++) + for (int iy=0; iy<=NUMY; iy++) + mpi_double_buf2[off + icnt++] = tdens[n](ix, iy); + + if (VFLAG & 16) + std::cerr << "Worker " << setw(4) << myid + << ": with request_id=" << request_id + << ", M=" << m << " sending Density" << std::endl; + + MPI_Send(&mpi_double_buf2[off], MPIbufsz, MPI_DOUBLE, 0, + 13 + MPItable*n+4, MPI_COMM_WORLD); } } // END: MPI packing @@ -1741,12 +1703,12 @@ void EmpCylSL::compute_even_odd(int request_id, int m) potC [m][n] = tpot[n]; rforceC[m][n] = trforce[n]; zforceC[m][n] = tzforce[n]; - if (DENS) densC[m][n] = tdens[n]; + densC [m][n] = tdens[n]; } else { potS [m][n] = tpot[n]; rforceS[m][n] = trforce[n]; zforceS[m][n] = trforce[n]; - if (DENS) densS[m][n] = tdens[n]; + densS [m][n] = tdens[n]; } } } @@ -1974,23 +1936,21 @@ void EmpCylSL::setup_table() rforceS .resize(MMAX+1); zforceS .resize(MMAX+1); - if (DENS) { - densC .resize(MMAX+1); - densS .resize(MMAX+1); - } + densC .resize(MMAX+1); + densS .resize(MMAX+1); for (int m=0; m<=MMAX; m++) { potC[m] .resize(rank3); rforceC[m].resize(rank3); zforceC[m].resize(rank3); - if (DENS) densC[m].resize(rank3); + densC[m] .resize(rank3); for (int v=0; v3) std::cerr << "Process " << myid << ": in EmpCylSL::accumlated_dens_eval, " @@ -5365,7 +5323,6 @@ void EmpCylSL::get_all(int mm, int nn, potC[mm][nn](ix+1, iy+1) * c11 ); - if (DENS) d += ccos * ( densC[mm][nn](ix , iy ) * c00 + @@ -5409,7 +5366,6 @@ void EmpCylSL::get_all(int mm, int nn, potS[mm][nn](ix+1, iy+1) * c11 ); - if (DENS) d += ssin * ( densS[mm][nn](ix , iy ) * c00 + @@ -5667,13 +5623,12 @@ void EmpCylSL::dump_basis(const string& name, int step, double Rmax) zforceC[mm][n](ix , iy+1) * c01 + zforceC[mm][n](ix+1, iy+1) * c11 ); - if (DENS) - outC << setw(15) - << fac*( - densC[mm][n](ix , iy ) * c00 + - densC[mm][n](ix+1, iy ) * c10 + - densC[mm][n](ix , iy+1) * c01 + - densC[mm][n](ix+1, iy+1) * c11 ); + outC << setw(15) + << fac*( + densC[mm][n](ix , iy ) * c00 + + densC[mm][n](ix+1, iy ) * c10 + + densC[mm][n](ix , iy+1) * c01 + + densC[mm][n](ix+1, iy+1) * c11 ); outC << endl; @@ -5699,13 +5654,12 @@ void EmpCylSL::dump_basis(const string& name, int step, double Rmax) zforceS[mm][n](ix , iy+1) * c01 + zforceS[mm][n](ix+1, iy+1) * c11 ); - if (DENS) - outS << setw(15) - << fac*( - densS[mm][n](ix , iy ) * c00 + - densS[mm][n](ix+1, iy ) * c10 + - densS[mm][n](ix , iy+1) * c01 + - densS[mm][n](ix+1, iy+1) * c11 ); + outS << setw(15) + << fac*( + densS[mm][n](ix , iy ) * c00 + + densS[mm][n](ix+1, iy ) * c10 + + densS[mm][n](ix , iy+1) * c01 + + densS[mm][n](ix+1, iy+1) * c11 ); outS << endl; } @@ -6322,7 +6276,6 @@ void EmpCylSL::dump_eof_file(const string& eof_file, const string& output) out << setw(20) << left << "NUMY" << " : " << numy << endl; out << setw(20) << left << "NMAXFID" << " : " << nmax << endl; out << setw(20) << left << "NMAX" << " : " << norder << endl; - out << setw(20) << left << "DENS" << " : " << std::boolalpha << dens << endl; out << setw(20) << left << "CMAPR" << " : " << cmap << endl; out << setw(20) << left << "RMIN" << " : " << rmin << endl; out << setw(20) << left << "RMAX" << " : " << rmax << endl; @@ -6339,8 +6292,7 @@ void EmpCylSL::dump_eof_file(const string& eof_file, const string& output) // Read table // - int nfield = 3; - if (DENS) nfield += 1; + int nfield = 4; Dynamic3dArray mat(nfield, numx+1, numy+1); @@ -6361,12 +6313,9 @@ void EmpCylSL::dump_eof_file(const string& eof_file, const string& output) for (int iy=0; iy<=numy; iy++) in.read((char *)&mat[2][ix][iy], sizeof(double)); - if (DENS) { - for (int ix=0; ix<=numx; ix++) - for (int iy=0; iy<=numy; iy++) - in.read((char *)&mat[3][ix][iy], sizeof(double)); - - } + for (int ix=0; ix<=numx; ix++) + for (int iy=0; iy<=numy; iy++) + in.read((char *)&mat[3][ix][iy], sizeof(double)); for (int ix=0; ixcur ? dif : cur; } - if (DENS) { - for (int ix=0; ix<=NUMX; ix++) - for (int iy=0; iy<=NUMY; iy++) { + for (int ix=0; ix<=NUMX; ix++) + for (int iy=0; iy<=NUMY; iy++) { - double one = densC[m][v](ix, iy); - double two = p->densC[m][v](ix, iy); + double one = densC[m][v](ix, iy); + double two = p->densC[m][v](ix, iy); - double cur = DBdif["densC"][m]; - double dif = fabs(one-two); - DBdif["densC"][m] = dif>cur ? dif : cur; - - cur = DBmax["densC"][m]; - dif = fabs(one); - DBmax["densC"][m] = dif>cur ? dif : cur; - } - } - + double cur = DBdif["densC"][m]; + double dif = fabs(one-two); + DBdif["densC"][m] = dif>cur ? dif : cur; + + cur = DBmax["densC"][m]; + dif = fabs(one); + DBmax["densC"][m] = dif>cur ? dif : cur; + } } } @@ -6718,22 +6662,19 @@ void EmpCylSL::compare_basis(const EmpCylSL *p) DBmax["zforceS"][m] = dif>cur ? dif : cur; } - if (DENS) { - for (int ix=0; ix<=NUMX; ix++) - for (int iy=0; iy<=NUMY; iy++) { - double one = densS[m][v](ix, iy); - double two = p->densS[m][v](ix, iy); + for (int ix=0; ix<=NUMX; ix++) + for (int iy=0; iy<=NUMY; iy++) { + double one = densS[m][v](ix, iy); + double two = p->densS[m][v](ix, iy); - double cur = DBdif["densS"][m]; - double dif = fabs(one-two); - DBdif["densS"][m] = dif>cur ? dif : cur; + double cur = DBdif["densS"][m]; + double dif = fabs(one-two); + DBdif["densS"][m] = dif>cur ? dif : cur; - cur = DBmax["densS"][m]; - dif = fabs(one); - DBmax["densS"][m] = dif>cur ? dif : cur; - } - } - + cur = DBmax["densS"][m]; + dif = fabs(one); + DBmax["densS"][m] = dif>cur ? dif : cur; + } } } @@ -6807,57 +6748,50 @@ double EmpCylSL::d_y_to_z(double y) // void EmpCylSL::ortho_check(std::ostream& out) { - if (DENS) { - - for (int mm=0; mm<=MMAX; mm++) { - // Header - // - out << std::string(60, '-') << std::endl - << " M=" << mm << std::endl - << std::string(60, '-') << std::endl; - - // Normalization: - // +--- Gravitational energy kernel - // | +--- Aximuthal - // | | - // v v - double fac = -4.0*M_PI * (2.0*M_PI) * dX * dY; - if (mm) fac *= 0.5; - - // Compute orthogonality matrix - // - for (int n1=0; n1 EmpCylSL::orthoCheck() { std::vector ret; - if (DENS) { - - ret.resize(MMAX+1); - for (auto & v : ret) v.resize(NORDER, NORDER); - - for (int mm=0; mm<=MMAX; mm++) { + ret.resize(MMAX+1); + for (auto & v : ret) v.resize(NORDER, NORDER); + + for (int mm=0; mm<=MMAX; mm++) { - // Normalization: - // +--- Gravitational energy kernel - // | +--- Aximuthal - // | | - // v v - double fac = -4.0*M_PI * (2.0*M_PI) * dX * dY; - if (mm) fac *= 0.5; + // Normalization: + // +--- Gravitational energy kernel + // | +--- Aximuthal + // | | + // v v + double fac = -4.0*M_PI * (2.0*M_PI) * dX * dY; + if (mm) fac *= 0.5; - // Compute orthogonality matrix - // + // Compute orthogonality matrix + // - // Unroll loop for OpenMP + // Unroll loop for OpenMP #pragma omp parallel for - for (int nn=0; nn0. + // + if (mm==0) + ret[mm](n1, n2) = sumC; + else ret[mm](n1, n2) = sqrt(0.5*(sumC*sumC + sumS*sumS)); - } } - } else { - std::cout << "EmpCylSL::orthoCheck: " - << "can not check orthogonality without density grid.\n" - << "Rerun EOF grid computation with DENS=true..." - << std::endl; } return ret; @@ -6933,8 +6868,6 @@ void EmpCylSL::getDensSC(int mm, int n, double R, double z, dC = 0.0; dS = 0.0; - if (not DENS) return; - if (R/ASCALE>Rtable or mm>MMAX or n>=rank3) return; double X = (r_to_xi(R) - XMIN)/dX; @@ -7059,11 +6992,6 @@ void EmpCylSL::WriteH5Cache() file.createAttribute("geometry", HighFive::DataSpace::From(geometry)).write(geometry); file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); - // HighFive boolean workarounds - // - int idens = 0; - if (DENS) idens = 1; - // Write the specific parameters // file.createAttribute("model", HighFive::DataSpace::From(model)). write(model); @@ -7075,7 +7003,6 @@ void EmpCylSL::WriteH5Cache() file.createAttribute ("nmaxfid", HighFive::DataSpace::From(NMAX)). write(NMAX); file.createAttribute ("neven", HighFive::DataSpace::From(Neven)). write(Neven); file.createAttribute ("nodd", HighFive::DataSpace::From(Nodd)). write(Nodd); - file.createAttribute ("idens", HighFive::DataSpace::From(idens)). write(idens); file.createAttribute ("cmapr", HighFive::DataSpace::From(CMAPR)). write(CMAPR); file.createAttribute ("cmapz", HighFive::DataSpace::From(CMAPZ)). write(CMAPZ); file.createAttribute ("rmin", HighFive::DataSpace::From(RMIN)). write(RMIN); @@ -7104,7 +7031,7 @@ void EmpCylSL::WriteH5Cache() order.createDataSet("potC", potC [m][n]); order.createDataSet("rforceC", rforceC[m][n]); order.createDataSet("zforceC", zforceC[m][n]); - if (DENS) order.createDataSet("densC", densC[m][n]); + order.createDataSet("densC", densC [m][n]); } } @@ -7128,7 +7055,7 @@ void EmpCylSL::WriteH5Cache() order.createDataSet("potS", potS [m][n]); order.createDataSet("rforceS", rforceS[m][n]); order.createDataSet("zforceS", zforceS[m][n]); - if (DENS) order.createDataSet("densS", densS[m][n]); + order.createDataSet("densS", densS [m][n]); } } @@ -7155,19 +7082,14 @@ bool EmpCylSL::ReadH5Cache() std::string forceID("Cylinder"), geometry("cylinder"); std::string model = EmpModelLabs[mtype]; - // HighFive boolean workarounds - // - int idens = 0; - if (DENS) idens = 1; - // Try checking the rest of the parameters before reading arrays // auto checkInt = [&file](int value, std::string name) { int v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); if (value == v) return true; - std::cout << "Parameter " << name << ": wanted " << value - << " found " << v << std::endl; + std::cout << "---- EmpCylSL cache parameter " << name << ": wanted " + << value << " found " << v << std::endl; return false; }; @@ -7175,8 +7097,8 @@ bool EmpCylSL::ReadH5Cache() { double v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); if (fabs(value - v) < 1.0e-16) return true; - std::cout << "Parameter " << name << ": wanted " << value - << " found " << v << std::endl; + std::cout << "---- EmpCylSL cache parameter " << name << ": wanted " + << value << " found " << v << std::endl; return false; }; @@ -7184,8 +7106,8 @@ bool EmpCylSL::ReadH5Cache() { std::string v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); if (value.compare(v)==0) return true; - std::cout << "Parameter " << name << ": wanted " << value - << " found " << v << std::endl; + std::cout << "--- EmpCylSL cache parameter " << name << ": wanted " + << value << " found " << v << std::endl; return false; }; @@ -7199,7 +7121,6 @@ bool EmpCylSL::ReadH5Cache() if (not checkInt(NMAX, "nmaxfid")) return false; if (not checkInt(Neven, "neven")) return false; if (not checkInt(Nodd, "nodd")) return false; - if (not checkInt(idens, "idens")) return false; if (not checkInt(CMAPR, "cmapr")) return false; if (not checkInt(CMAPZ, "cmapz")) return false; if (not checkDbl(RMIN, "rmin")) return false; @@ -7224,23 +7145,21 @@ bool EmpCylSL::ReadH5Cache() rforceS .resize(MMAX+1); zforceS .resize(MMAX+1); - if (DENS) { - densC .resize(MMAX+1); - densS .resize(MMAX+1); - } + densC .resize(MMAX+1); + densS .resize(MMAX+1); for (int m=0; m<=MMAX; m++) { potC[m] .resize(NORDER); rforceC[m].resize(NORDER); zforceC[m].resize(NORDER); - if (DENS) densC[m].resize(NORDER); + densC[m] .resize(NORDER); } for (int m=1; m<=MMAX; m++) { potS[m] .resize(NORDER); rforceS[m].resize(NORDER); zforceS[m].resize(NORDER); - if (DENS) densS[m].resize(NORDER); + densS[m] .resize(NORDER); } // Read arrays and data from H5 file @@ -7265,7 +7184,7 @@ bool EmpCylSL::ReadH5Cache() order.getDataSet("potC") .read(potC [m][n]); order.getDataSet("rforceC").read(rforceC[m][n]); order.getDataSet("zforceC").read(zforceC[m][n]); - if (DENS) order.getDataSet("densC").read(densC[m][n]); + order.getDataSet("densC") .read(densC[m][n]); } } @@ -7288,7 +7207,7 @@ bool EmpCylSL::ReadH5Cache() order.getDataSet("potS") .read(potS [m][n]); order.getDataSet("rforceS").read(rforceS[m][n]); order.getDataSet("zforceS").read(zforceS[m][n]); - if (DENS) order.getDataSet("densS").read(densS[m][n]); + order.getDataSet("densS") .read(densS[m][n]); } } diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 353c60a1d..37ed197f9 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -62,6 +62,30 @@ extern "C" { } +static +std::string sledge_error(int flag) +{ + if (flag==0) + return "reliable"; + else if (flag==-1) + return "too many levels for eigenvalues"; + else if (flag==-2) + return "too many levels for eigenvectors"; + else if (flag==1) + return "eigenvalue cluster?"; + else if (flag==2) + return "classification uncertainty"; + else if (flag>0) { + std::ostringstream sout; + sout << "unexpected warning: " << flag; + return sout.str(); + } else { + std::ostringstream sout; + sout << "unexpected fatal error: " << flag; + return sout.str(); + } +} + // Unit density exponential disk with scale length A class KuzminCyl : public CylModel @@ -177,7 +201,7 @@ void SLGridCyl::bomb(std::string oops) SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, double RMIN, double RMAX, double L, - bool CACHE, int CMAP, double SCALE, + bool CACHE, int CMAP, double RMAP, const std::string type, bool VERBOSE) { int m, k; @@ -193,7 +217,7 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, cache = CACHE; cmap = CMAP; - scale = SCALE; + rmap = RMAP; tbdbg = VERBOSE; @@ -229,6 +253,8 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, mpi_setup(); + int totbad = 0; + if (mpi_myid) { compute_table_worker(); @@ -286,12 +312,23 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, // // // - - MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, - MPI_ANY_TAG, MPI_COMM_WORLD, &status); - +#ifdef SLEDGE_THROW + int bad; // Get count of sledge failures + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_SOURCE, 10, + MPI_COMM_WORLD, &status); + totbad += bad; + // Get source node identity int retid = status.MPI_SOURCE; + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, retid, 11, + MPI_COMM_WORLD, &status); +#else + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, 11, + MPI_COMM_WORLD, &status); + + int retid = status.MPI_SOURCE; +#endif + mpi_unpack_table(); // @@ -323,9 +360,19 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, while (worker) { - +#ifdef SLEDGE_THROW + int bad; // Get count of sledge failures + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_SOURCE, 10, + MPI_COMM_WORLD, &status); + totbad += bad; + // Get source node identity + int retid = status.MPI_SOURCE; + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, retid, 11, + MPI_COMM_WORLD, &status); +#else MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status); +#endif mpi_unpack_table(); @@ -357,9 +404,31 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, } } - } // END Root +#ifdef SLEDGE_THROW + // Send total number of sledge failures to all nodes + MPI_Bcast(&totbad, 1, MPI_INT, 0, MPI_COMM_WORLD); + + // All nodes throw exception if any errors + if (totbad) { + std::ostringstream sout; + if (myid==0) { + sout << std::endl + << "SLGridCyl found " << totbad + << " tolerance errors in computing SL solutions" << std::endl + << "We suggest checking your model for sufficient smoothness and" + << std::endl + << "sufficiently many grid points that the relative difference" + << std::endl + << "between field quantities is <= 0.3"; + } else { + sout << std::endl << "sledge tolerance failure"; + } + + throw GenericError(sout.str(), __FILE__, __LINE__); + } +#endif } else { @@ -389,7 +458,7 @@ int SLGridCyl::read_cached_table(void) if (!in) return 0; int MMAX, NMAX, NUMR, NUMK, i, j, CMAP; - double RMIN, RMAX, L, AA, SCL; + double RMIN, RMAX, L, AA, RMAP; std::string MODEL; if (myid==0) @@ -423,23 +492,23 @@ int SLGridCyl::read_cached_table(void) std::ostringstream sout; sout << "YAML: error parsing <" << buf.get() << "> " << "in " << __FILE__ << ":" << __LINE__ << std::endl - << "YAML error: " << error.what() << std::endl; + << "YAML error: " << error.what(); throw GenericError(sout.str(), __FILE__, __LINE__, 1042, false); } // Get parameters // - MMAX = node["mmax" ].as(); - NMAX = node["nmax" ].as(); - NUMK = node["numk" ].as(); - NUMR = node["numr" ].as(); - CMAP = node["cmap" ].as(); - RMIN = node["rmin" ].as(); - RMAX = node["rmax" ].as(); - SCL = node["scale" ].as(); - L = node["L" ].as(); - AA = node["A" ].as(); - MODEL = node["model" ].as(); + MMAX = node["mmax" ].as(); + NMAX = node["nmax" ].as(); + NUMK = node["numk" ].as(); + NUMR = node["numr" ].as(); + CMAP = node["cmap" ].as(); + RMIN = node["rmin" ].as(); + RMAX = node["rmax" ].as(); + RMAP = node["rmapping"].as(); + L = node["L" ].as(); + AA = node["A" ].as(); + MODEL = node["model" ].as(); } else { std::cout << "---- SLGridCyl: bad magic number in cache file" << std::endl; return 0; @@ -495,10 +564,10 @@ int SLGridCyl::read_cached_table(void) return 0; } - if (SCL!=scale) { + if (RMAP!=rmap) { if (myid==0) - std::cout << "---- SLGridCyl::read_cached_table: found scale=" << SCL - << " wanted " << scale << std::endl; + std::cout << "---- SLGridCyl::read_cached_table: found rmapping=" << RMAP + << " wanted " << rmap << std::endl; return 0; } @@ -571,17 +640,17 @@ void SLGridCyl::write_cached_table(void) // content can be added as needed. YAML::Node node; - node["mmax" ] = mmax; - node["nmax" ] = nmax; - node["numk" ] = numk; - node["numr" ] = numr; - node["cmap" ] = cmap; - node["rmin" ] = rmin; - node["rmax" ] = rmax; - node["scale" ] = scale; - node["L" ] = l; - node["A" ] = A; - node["model" ] = cyl->ID(); + node["mmax" ] = mmax; + node["nmax" ] = nmax; + node["numk" ] = numk; + node["numr" ] = numr; + node["cmap" ] = cmap; + node["rmin" ] = rmin; + node["rmax" ] = rmax; + node["rmapping"] = rmap; + node["L" ] = l; + node["A" ] = A; + node["model" ] = cyl->ID(); // Serialize the node // @@ -645,7 +714,7 @@ double SLGridCyl::r_to_xi(double r) } if (cmap) { - return (r/scale-1.0)/(r/scale+1.0); + return (r/rmap-1.0)/(r/rmap+1.0); } else { return r; } @@ -666,7 +735,7 @@ double SLGridCyl::xi_to_r(double xi) bomb(ostr.str()); } - return (1.0+xi)/(1.0 - xi) * scale; + return (1.0+xi)/(1.0 - xi) * rmap; } else { return xi; } @@ -688,7 +757,7 @@ double SLGridCyl::d_xi_to_r(double xi) bomb(ostr.str()); } - return 0.5*(1.0-xi)*(1.0-xi)/scale; + return 0.5*(1.0-xi)*(1.0-xi)/rmap; } else { return 1.0; } @@ -1198,9 +1267,9 @@ void SLGridCyl::compute_table(struct TableCyl* table, int m, int k) { double cons[8] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - double tol[6] = {1.0e-4*scale,1.0e-5, - 1.0e-4*scale,1.0e-5, - 1.0e-4*scale,1.0e-5}; + double tol[6] = {1.0e-4*rmap,1.0e-5, + 1.0e-4*rmap,1.0e-5, + 1.0e-4*rmap,1.0e-5}; int VERBOSE=0; integer NUM, N, M; logical type[8]; @@ -1280,6 +1349,63 @@ void SLGridCyl::compute_table(struct TableCyl* table, int m, int k) sledge_(job, cons, endfin, invec, tol, type, ev, &NUM, xef, ef, pdef, t, rho, iflag, store); // + // Check for errors + // +#ifdef SLEDGE_THROW + unsigned bad = 0; // Count non-zero sledge flags + for (int i=0; i0) { + + if (myid==0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for m=" + << std::endl << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; iev.resize(N); for (int i=0; iev[i] = ev[i]; @@ -1347,8 +1472,8 @@ void SLGridCyl::init_table(void) d0.resize(numr); if (cmap) { - xmin = (rmin/scale - 1.0)/(rmin/scale + 1.0); - xmax = (rmax/scale - 1.0)/(rmax/scale + 1.0); + xmin = (rmin/rmap - 1.0)/(rmin/rmap + 1.0); + xmax = (rmax/rmap - 1.0)/(rmax/rmap + 1.0); dxi = (xmax-xmin)/(numr-1); } else { xmin = rmin; @@ -1370,9 +1495,9 @@ void SLGridCyl::compute_table_worker(void) { double cons[8] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - double tol[6] = {1.0e-4*scale,1.0e-5, - 1.0e-4*scale,1.0e-5, - 1.0e-4*scale,1.0e-5}; + double tol[6] = {1.0e-4*rmap,1.0e-5, + 1.0e-4*rmap,1.0e-5, + 1.0e-4*rmap,1.0e-5}; int i, j, VERBOSE=0; integer NUM; logical type[8]; @@ -1470,10 +1595,51 @@ void SLGridCyl::compute_table_worker(void) sledge_(job, cons, endfin, invec, tol, type, ev, &NUM, xef, ef, pdef, t, rho, iflag, store); + // - // Print results: + // Check for errors // +#ifdef SLEDGE_THROW + unsigned bad = 0; // Number of non-zero iflag values + for (int i=0; i0) { + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for m=" + << M << " K=" << K << std::endl << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; i mod, int LMAX, int NMAX, int NUMR, double RMIN, double RMAX, - bool CACHE, int CMAP, double SCALE, + bool CACHE, int CMAP, double RMAP, std::string cachename, bool VERBOSE) { mpi_buf = 0; @@ -1738,7 +1909,7 @@ SLGridSph::SLGridSph(std::shared_ptr mod, if (cachename.size()) sph_cache_name = cachename; else sph_cache_name = default_cache; - initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, SCALE); + initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, RMAP); } @@ -1768,7 +1939,7 @@ SLGridSph::cacheInfo(const std::string& cachefile, bool verbose) void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, double RMIN, double RMAX, - bool CACHE, int CMAP, double SCALE) + bool CACHE, int CMAP, double RMAP) { int l; @@ -1781,7 +1952,7 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, cache = CACHE; cmap = CMAP; - scale = SCALE; + rmap = RMAP; init_table(); @@ -1805,6 +1976,8 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, mpi_setup(); + int totbad = 0; // Count total number of sledge errors + if (mpi_myid) { // Begin workers compute_table_worker(); @@ -1854,10 +2027,21 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, // // +#ifdef SLEDGE_THROW + int bad; // Get the sledge error count from a + // worker + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_SOURCE, 10, + MPI_COMM_WORLD, &status); + totbad += bad; + int retid = status.MPI_SOURCE; + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, retid, 11, + MPI_COMM_WORLD, &status); +#else MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status); int retid = status.MPI_SOURCE; +#endif mpi_unpack_table(); @@ -1882,11 +2066,21 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, // while (worker) { - - - MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, - MPI_ANY_TAG, MPI_COMM_WORLD, &status); +#ifdef SLEDGE_THROW + int bad; // Get the sledge error count from a + // worker + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_SOURCE, 10, + MPI_COMM_WORLD, &status); + totbad += bad; + + int retid = status.MPI_SOURCE; + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, retid, 11, + MPI_COMM_WORLD, &status); +#else + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, 11, + MPI_COMM_WORLD, &status); +#endif mpi_unpack_table(); worker--; @@ -1913,6 +2107,31 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, } // END Root + +#ifdef SLEDGE_THROW + // Share the total sledge error count with all nodes + MPI_Bcast(&totbad, 1, MPI_INT, 0, MPI_COMM_WORLD); + + // Emit runtime exception on sledge errors + if (totbad) { + std::ostringstream sout; + + if (myid==0) { + sout << std::endl + << "SLGridSph found " << totbad + << " tolerance errors in computing SL solutions." << std::endl + << "We suggest checking your model file for smoothness and for" + << std::endl + << "a sufficient number grid points that the relative difference" + << std::endl + << "between field quantities is <= 0.3"; + } else { + sout << std::endl << "sledge tolerance failure"; + } + + throw GenericError(sout.str(), __FILE__, __LINE__); + } +#endif } // END MPI stanza, BEGIN single-process stanza else { @@ -2039,7 +2258,7 @@ bool SLGridSph::ReadH5Cache(void) if (not checkInt(cmap, "cmap")) return false; if (not checkDbl(rmin, "rmin")) return false; if (not checkDbl(rmax, "rmax")) return false; - if (not checkDbl(scale, "scale")) return false; + if (not checkDbl(rmap, "rmapping")) return false; if (not checkInt(diverge, "diverge")) return false; if (not checkDbl(dfac, "dfac")) return false; @@ -2063,7 +2282,7 @@ bool SLGridSph::ReadH5Cache(void) } if (myid==0) - std::cerr << "---- SLGridSph::ReadH5Cache: " + std::cout << "---- SLGridSph::ReadH5Cache: " << "successfully read basis cache <" << sph_cache_name << ">" << std::endl; @@ -2100,7 +2319,7 @@ void SLGridSph::WriteH5Cache(void) sout << "---- SLGridSph::WriteH5Cache write error: " << "what(): " << ex.what() << std::endl << "path1(): " << ex.path1() << std::endl - << "path2(): " << ex.path2() << std::endl; + << "path2(): " << ex.path2(); throw GenericError(sout.str(), __FILE__, __LINE__, 12, true); } @@ -2130,7 +2349,7 @@ void SLGridSph::WriteH5Cache(void) file.createAttribute ("cmap", HighFive::DataSpace::From(cmap)).write(cmap); file.createAttribute ("rmin", HighFive::DataSpace::From(rmin)).write(rmin); file.createAttribute ("rmax", HighFive::DataSpace::From(rmax)).write(rmax); - file.createAttribute ("scale", HighFive::DataSpace::From(scale)).write(scale); + file.createAttribute ("rmapping", HighFive::DataSpace::From(rmap)).write(rmap); file.createAttribute ("diverge", HighFive::DataSpace::From(diverge)).write(diverge); file.createAttribute ("dfac", HighFive::DataSpace::From(dfac)).write(dfac); @@ -2172,7 +2391,7 @@ double SLGridSph::r_to_xi(double r) if (cmap==1) { if (r<0.0) bomb("radius < 0!"); - ret = (r/scale-1.0)/(r/scale+1.0); + ret = (r/rmap-1.0)/(r/rmap+1.0); } else if (cmap==2) { if (r<=0.0) bomb("radius <= 0!"); ret = log(r); @@ -2191,7 +2410,7 @@ double SLGridSph::xi_to_r(double xi) if (xi<-1.0) bomb("xi < -1!"); if (xi>=1.0) bomb("xi >= 1!"); - ret =(1.0+xi)/(1.0 - xi) * scale; + ret =(1.0+xi)/(1.0 - xi) * rmap; } else if (cmap==2) { ret = exp(xi); } else { @@ -2211,7 +2430,7 @@ double SLGridSph::d_xi_to_r(double xi) if (xi<-1.0) bomb("xi < -1!"); if (xi>=1.0) bomb("xi >= 1!"); - ret = 0.5*(1.0-xi)*(1.0-xi)/scale; + ret = 0.5*(1.0-xi)*(1.0-xi)/rmap; } else if (cmap==2) { ret = exp(-xi); } else { @@ -2562,9 +2781,9 @@ void SLGridSph::compute_table(struct TableSph* table, int l) { double cons[8] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - double tol[6] = {1.0e-4*scale,1.0e-6, - 1.0e-4*scale,1.0e-6, - 1.0e-4*scale,1.0e-6}; + double tol [6] = {1.0e-4*rmap, 1.0e-6, + 1.0e-4*rmap, 1.0e-6, + 1.0e-4*rmap, 1.0e-6}; int VERBOSE=0; integer NUM, N; logical type[8] = {0, 0, 1, 0, 0, 0, 1, 0}; @@ -2638,6 +2857,66 @@ void SLGridSph::compute_table(struct TableSph* table, int l) sledge_(job, cons, endfin, invec, tol, type, ev, &NUM, xef, ef, pdef, t, rho, iflag, store); + + // + // Check for errors + // +#ifdef SLEDGE_THROW + // Accumulate error count + unsigned bad = 0; + for (int i=0; i0) { + + if (myid==0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for l=" << l + << std::endl << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; i0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for l=" << L + << std::endl << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; i. "; - sout << "YAML error in getHeader: " << error.what() << std::endl; + sout << "YAML error in getHeader: " << error.what(); throw GenericError(sout.str(), __FILE__, __LINE__, 1038, false); } @@ -3054,6 +3382,57 @@ YAML::Node SLGridSph::getHeader(const std::string& cachefile) } +std::vector SLGridSph::orthoCheck(int num) +{ + // Gauss-Legendre knots and weights + LegeQuad lw(num); + + // Get the scaled coordinate limits + double ximin = r_to_xi(rmin); + double ximax = r_to_xi(rmax); + + // Initialize the return matrices + std::vector ret(lmax+1); + for (auto & v : ret) v.resize(nmax, nmax); + + // Do each harmonic order + for (int L=0; L<=lmax; L++) { + + // Unroll the loop for OpenMP parallelization +#pragma omp parallel for + for (int nn=0; nn // - +#ifdef SLEDGE_THROW + int bad; // Get sledge error count + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_TAG, 10, + MPI_COMM_WORLD, &status); + totbad += bad; + // Get sledge comptuation result + int retid = status.MPI_SOURCE; + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, retid, 11, + MPI_COMM_WORLD, &status); +#else MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status); int retid = status.MPI_SOURCE; +#endif mpi_unpack_table(); @@ -3316,10 +3708,20 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, while (worker) { - - MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, - MPI_ANY_TAG, MPI_COMM_WORLD, &status); - +#ifdef SLEDGE_THROW + // Get sledge error count + int bad; + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_SOURCE, 10, MPI_COMM_WORLD, + &status); + totbad += bad; + // Get sledge computation result + int retid = status.MPI_SOURCE; + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, retid, 11, + MPI_COMM_WORLD, &status); +#else + MPI_Recv(mpi_buf, mpi_bufsz, MPI_PACKED, MPI_ANY_SOURCE, 11, + MPI_COMM_WORLD, &status); +#endif mpi_unpack_table(); worker--; @@ -3329,7 +3731,6 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, } - // // // @@ -3354,6 +3755,28 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, } // END Root +#ifdef SLEDGE_THROW + // Share total sledge error count with all nodes + MPI_Bcast(&totbad, 1, MPI_INT, 0, MPI_COMM_WORLD); + + // Throw runtime error if sledge comtputation has errors + if (totbad) { + std::ostringstream sout; + + if (myid==0) { + sout << std::endl << "SLGridSlab found " << totbad + << " tolerance errors in computing SL solutions." << std::endl + << "We suggest checking your model parameters to ensure a" + << std::endl + << "sufficient number of grid points that the relative difference" + << std::endl << "between field quantities is <= 0.3"; + } else { + sout << std::endl << "sledge tolerance failure"; + } + + throw GenericError(sout.str(), __FILE__, __LINE__); + } +#endif } else { @@ -3422,7 +3845,7 @@ int SLGridSlab::read_cached_table(void) std::ostringstream sout; sout << "YAML: error parsing <" << buf.get() << "> " << "in " << __FILE__ << ":" << __LINE__ << std::endl - << "YAML error: " << error.what() << std::endl; + << "YAML error: " << error.what(); throw GenericError(sout.str(), __FILE__, __LINE__, 1042, false); } @@ -4111,10 +4534,67 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) sledge_(job, cons, endfin, invec, tol, type, ev, &NUM, xef, ef, pdef, t, rho, iflag, store); + + // + // Check for errors + // +#ifdef SLEDGE_THROW + unsigned bad = 0; // Number of non-zero iflag values + for (int i=0; i0) { + + if (myid==0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", even" << std::endl << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; i0) { + + if (myid==0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", odd" << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; i0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", even" << std::endl << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; i0) { + + // Print sledge error info and throw runtime exception + if (myid==0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + + std::cout << std::endl + << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", odd" << std::endl << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" + << std::endl; + + for (int i=0; i +#include +#include +#include + +#include // For MatrixXd + +#include // For OpenMP API + +#include // MPI support +#include // For orthoTol + +void orthoTest(const std::vector& tests, + const std::string& classname, const std::string& indexname) +{ + // Number of possible threads + int nthrds = omp_get_max_threads(); + + // Worst so far + std::vector worst(nthrds), lworst(tests.size());; + + // Rank + int nmax = tests[0].rows(); + + // Test loop + for (int l=0; l(worst[tid], + fabs(1.0 - tests[l](n1, n2))); + else + worst[tid] = std::max(worst[tid], + fabs(tests[l](n1, n2))); + } + // END: unrolled loop + + lworst[l] = *std::max_element(worst.begin(), worst.end()); + } + // END: harmonic order loop + + double worst_ever = *std::max_element(lworst.begin(), lworst.end()); + + if (worst_ever > __EXP__::orthoTol) { + std::cout << classname << ": orthogonality failure" << std::endl + << std::right + << std::setw(4) << indexname + << std::setw(16) << "Worst" << std::endl; + for (int l=0; l AxiDiskPtr; - //! TRUE if density is computed (default: false) - static bool DENS; - //! TRUE if signal-to-noise methods are on (default: false) static bool PCAVAR; diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index dd2e3b2f1..2a04a4c44 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -58,7 +58,7 @@ private: double rmin, rmax, l; int cmap; - double scale; + double rmap; double dk; @@ -108,7 +108,7 @@ public: //! Constructor SLGridCyl(int mmax, int nmax, int numr, int numk, double rmin, double rmax, - double l, bool cache, int Cmap, double Scale, + double l, bool cache, int Cmap, double RMAP, const std::string type, bool Verbose=false); //! Destructor @@ -180,7 +180,7 @@ private: double rmin, rmax; int cmap, diverge; - double scale, dfac; + double rmap, dfac; double xmin, xmax, dxi; @@ -193,7 +193,7 @@ private: void initialize(int LMAX, int NMAX, int NUMR, double RMIN, double RMAX, - bool CACHE, int CMAP, double SCALE); + bool CACHE, int CMAP, double RMAP); void init_table(void); void compute_table(TableSph* table, int L); @@ -248,14 +248,14 @@ public: //! Constructor with model table SLGridSph(std::shared_ptr mod, int lmax, int nmax, int numr, double rmin, double rmax, - bool cache, int Cmap, double Scale, + bool cache, int Cmap, double RMAP, std::string cachename="", bool Verbose=false); //! Constructor (uses file *model_file_name* for file) SLGridSph(std::string modelname, int lmax, int nmax, int numr, double rmin, double rmax, - bool cache, int Cmap, double Scale, + bool cache, int Cmap, double RMAP, int DIVERGE, double DFAC, std::string cachename="", bool Verbose=false); @@ -303,6 +303,10 @@ public: static std::map cacheInfo(const std::string& cachefile, bool verbose=true); + //! Compute the orthogonality of the basis by returning inner + //! produce matrices + std::vector orthoCheck(int knots=40); + #if HAVE_LIBCUDA==1 void initialize_cuda(std::vector& cuArray, thrust::host_vector& tex); @@ -311,7 +315,7 @@ public: { cudaMappingConstants ret; - ret.rscale = scale; + ret.rscale = rmap; ret.hscale = 0.0; ret.xmin = xmin; ret.xmax = xmax; diff --git a/include/exputils.H b/include/exputils.H new file mode 100644 index 000000000..d638535eb --- /dev/null +++ b/include/exputils.H @@ -0,0 +1,12 @@ +// exputils signatures + +#ifndef _EXPUTILS_H +#define _EXPUTILS_H + +#include +#include + +void orthoTest(const std::vector& tests, + const std::string& classname, const std::string& indexname); + +#endif diff --git a/include/libvars.H b/include/libvars.H index ee9e4f2f3..12dfd45d3 100644 --- a/include/libvars.H +++ b/include/libvars.H @@ -39,6 +39,9 @@ namespace __EXP__ //! Sign convention element rank for eigenfunctions and eigenvectors extern int nevsign; + //! Sanity tolerance for orthogonality + extern double orthoTol; + }; #endif // END _LIBVARS_H diff --git a/include/testembed.cc b/include/testembed.cc deleted file mode 100644 index e69de29bb..000000000 diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 4b6ce281a..2b71ff366 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -134,6 +134,14 @@ void BasisFactoryClasses(py::module &m) { PYBIND11_OVERRIDE_PURE(void, Basis, set_coefs, coefs); } + const std::string classname() override { + PYBIND11_OVERRIDE_PURE(std::string, Basis, classname); + } + + const std::string harmonic() override { + PYBIND11_OVERRIDE_PURE(std::string, Basis, harmonic); + } + public: // Inherit the constructors using BasisClasses::Basis::Basis; @@ -183,6 +191,14 @@ void BasisFactoryClasses(py::module &m) { PYBIND11_OVERRIDE(void, SphericalSL, set_coefs, coefs); } + const std::string classname() override { + PYBIND11_OVERRIDE(std::string, SphericalSL, classname); + } + + const std::string harmonic() override { + PYBIND11_OVERRIDE(std::string, SphericalSL, harmonic); + } + public: // Inherit the constructors @@ -232,6 +248,14 @@ void BasisFactoryClasses(py::module &m) { PYBIND11_OVERRIDE(void, Cylindrical, set_coefs, coefs); } + const std::string classname() override { + PYBIND11_OVERRIDE(std::string, Cylindrical, classname); + } + + const std::string harmonic() override { + PYBIND11_OVERRIDE(std::string, Cylindrical, harmonic); + } + public: // Inherit the constructors @@ -283,6 +307,14 @@ void BasisFactoryClasses(py::module &m) { PYBIND11_OVERRIDE(void, FlatDisk, set_coefs, coefs); } + const std::string classname() override { + PYBIND11_OVERRIDE(std::string, FlatDisk, classname); + } + + const std::string harmonic() override { + PYBIND11_OVERRIDE(std::string, FlatDisk, harmonic); + } + public: // Inherit the constructors diff --git a/src/Cylinder.H b/src/Cylinder.H index 064ce30a6..b2c5ed665 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -94,8 +94,6 @@ class MixtureBasis; @param try_cache false suppresses cache reading on restart (default: true) - @param density builds grid with density fields (e.g. for post processing) - @param EVEN_M true uses even harmonic orders only @param cmap is the coordinate mapping type (deprecated but kept for backward consistency) @@ -136,7 +134,7 @@ private: double hcyl, hexp, snr, rem; int nmax, ncylodd, ncylrecomp, npca, npca0, nvtk, cmapR, cmapZ; string eof_file; - bool self_consistent, logarithmic, density, pcavar, pcainit, pcavtk, pcadiag, pcaeof, eof_over; + bool self_consistent, logarithmic, pcavar, pcainit, pcavtk, pcadiag, pcaeof, eof_over; bool try_cache, firstime, dump_basis, compute, firstime_coef; // These should be ok for all derived classes, hence declared private @@ -357,7 +355,6 @@ public: //! \param pcavtk set to true dumps PCA functions in VTK format for diagnostics (default: false) //! \param try_cache set to true means try to read basis from cache (default: true) //! \param dump_basis set to true outputs basis into file - //! \param density set to true means compute density basis in addition to potential (default: false) //! \param cmapR selects the radial coordinate mapping (default: 2, power mapping) //! \param cmapZ selects the vertical coordinate mapping (default: 1, sinh mapping) //! diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 6c3ab77f0..5ee27894f 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -11,11 +11,13 @@ #include #include #include +#include Timer timer_debug; double EXPSCALE=1.0, HSCALE=1.0, ASHIFT=0.25; + //@{ //! These are for testing exclusively (should be set false for production) static bool cudaAccumOverride = false; @@ -166,7 +168,6 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : subsamp = false; nvtk = 1; pcainit = true; - density = false; coef_dump = true; try_cache = true; dump_basis = false; @@ -208,11 +209,6 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : // if (eof_file.size()) cachename = eof_file; - // For debugging; no use by force algorithm - // - if (density) EmpCylSL::DENS = true; - - // Make the empirical orthogonal basis instance // ortho = std::make_shared @@ -379,6 +375,14 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : } #endif + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) + // + std::cout << "---- "; + orthoTest(ortho->orthoCheck(), "Cylinder", "m"); + + // Initialize internal variables + // ncompcyl = 0; pos.resize(nthrds); @@ -443,13 +447,20 @@ void Cylinder::initialize() if (conf["pcadiag" ]) pcadiag = conf["pcadiag" ].as(); if (conf["subsamp" ]) subsamp = conf["subsamp" ].as(); if (conf["try_cache" ]) try_cache = conf["try_cache" ].as(); - if (conf["density" ]) density = conf["density" ].as(); if (conf["EVEN_M" ]) EVEN_M = conf["EVEN_M" ].as(); if (conf["cmap" ]) cmapR = conf["cmap" ].as(); if (conf["cmapr" ]) cmapR = conf["cmapr" ].as(); if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); + // Deprecation warning + if (conf["density" ]) { + if (myid==0) + std::cout << "Cylinder: parameter 'density' is deprecated. " + << "The density field will be computed regardless." + << std::endl; + } + if (not conf["ncylodd"]) { ncylodd = nmax/4; } @@ -1477,10 +1488,7 @@ void Cylinder::determine_fields_at_point_cyl(double r, double z, double phi, *tpotr *= -1.0; *tpotz *= -1.0; *tpotp *= -1.0; - if (density) - *tdens = ortho->accumulated_dens_eval(r, z, phi, *tdens0); - else - *tdens = 0.0; + *tdens = ortho->accumulated_dens_eval(r, z, phi, *tdens0); } // Dump coefficients to a file diff --git a/src/Sphere.H b/src/Sphere.H index e6ade7f32..1178b996b 100644 --- a/src/Sphere.H +++ b/src/Sphere.H @@ -19,7 +19,7 @@ typedef std::shared_ptr SLGridSphPtr; \param rmin is the minimum value in the table for the radial basis functions (default is 0.001) - \param rs is the halo scale length (default is 0.067*rmax). This + \param rmapping is the halo scale length (default is 0.067*rmax). This is used to for the coordinate mapping (see below). \param cmap should be 1 for algebraic coordinate scaling, 2 for log @@ -88,7 +88,7 @@ private: // Parameters double rsphSL; - double rs; + double rmap; double tnext, dtime; int numr; int nums; @@ -120,7 +120,7 @@ public: @param m allows the spherical basis to be used for multiple center expansions Input line parameters include: - @param rs is the radius for coordinate scaling + @param rmapping is the radius for coordinate scaling @param numr is the number of radial grid points @param cmap set to true for scaling coordinates from the semi-infinite to finite segment @param diverge set to true means assume a cuspy profile diff --git a/src/Sphere.cc b/src/Sphere.cc index 3029959c8..51ffe5fa1 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -7,10 +7,11 @@ #include #include #include +#include const std::set Sphere::valid_keys = { - "rs", + "rmapping", "numr", "nums", "cmap", @@ -28,7 +29,7 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : { id = "Sphere SL"; // Defaults - rs = 0.067*rmax; + rmap = 0.067*rmax; numr = 2000; nums = 2000; cmap = 1; @@ -60,7 +61,7 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : // Generate Sturm-Liouville grid ortho = std::make_shared(modelname, Lmax, nmax, numr, rmin, rmax, true, - cmap, rs, diverge, dfac, cachename); + cmap, rmap, diverge, dfac, cachename); // Get the min and max expansion radii rmin = ortho->getRmin(); @@ -71,6 +72,7 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : std::cout << "---- Sphere parameters: " << std::endl << sep << "lmax=" << Lmax << std::endl << sep << "nmax=" << nmax + << std::endl << sep << "rmapping=" << rmap << std::endl << sep << "cmap=" << cmap << std::endl << sep << "rmin=" << rmin << std::endl << sep << "rmax=" << rmax @@ -85,6 +87,12 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : } + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) + // + std::cout << "---- "; + orthoTest(ortho->orthoCheck(std::max(nmax*50, 200)), "Sphere", "l"); + setup(); } @@ -98,7 +106,7 @@ void Sphere::initialize() // Assign values from YAML // try { - if (conf["rs"]) rs = conf["rs"].as(); + if (conf["rmapping"]) rmap = conf["rmapping"].as(); if (conf["numr"]) numr = conf["numr"].as(); if (conf["nums"]) nums = conf["nums"].as(); if (conf["cmap"]) cmap = conf["cmap"].as(); @@ -319,6 +327,12 @@ void Sphere::make_model_bin() std::string cachename = outdir + cache_file + "." + runtag; ortho = std::make_shared(mod, Lmax, nmax, numR, Rmin, Rmax, false, 1, 1.0, cachename); + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) + // + std::cout << "---- "; + orthoTest(ortho->orthoCheck(std::max(nmax*50, 200)), "Sphere", "l"); + // Update time trigger // tnext = tnow + dtime; @@ -459,6 +473,12 @@ void Sphere::make_model_plummer() std::string cachename = outdir + cache_file + "." + runtag; ortho = std::make_shared(mod, Lmax, nmax, numr, Rmin, Rmax, false, 1, 1.0, cachename); + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) + // + std::cout << "---- "; + orthoTest(ortho->orthoCheck(std::max(nmax*50, 200)), "Sphere", "l"); + // Update time trigger // tnext = tnow + dtime; diff --git a/utils/Analysis/KDcyltest.cc b/utils/Analysis/KDcyltest.cc index ac91693e0..b272ccdbb 100644 --- a/utils/Analysis/KDcyltest.cc +++ b/utils/Analysis/KDcyltest.cc @@ -153,7 +153,6 @@ main(int argc, char **argv) int mmax, numx, numy, norder, cmapr, cmapz, nodd=-1; double rcylmin, rcylmax, vscale; - bool DENS; if (not ignore) { @@ -210,7 +209,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); NMAX = node["nmax" ].as(); norder = node["norder"].as(); - DENS = node["dens" ].as(); if (node["nodd"]) nodd = node["nodd" ].as(); if (node["cmap"]) @@ -242,10 +240,6 @@ main(int argc, char **argv) in.read((char *)&rcylmax, sizeof(double)); in.read((char *)&rscale, sizeof(double)); in.read((char *)&vscale, sizeof(double)); - - std::cout << "idens=" << idens << std::endl; - - if (idens) DENS = true; } } } @@ -257,7 +251,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = true; - EmpCylSL::DENS = DENS; // Create expansion // diff --git a/utils/Analysis/KL_cyl.cc b/utils/Analysis/KL_cyl.cc index 97e2586b4..c77a9416b 100644 --- a/utils/Analysis/KL_cyl.cc +++ b/utils/Analysis/KL_cyl.cc @@ -228,7 +228,6 @@ main(int argc, char **argv) int mmax, numx, numy, norder, cmapr, cmapz, nodd=-1; double rcylmin, rcylmax, vscale; - bool DENS; if (not ignore) { @@ -285,7 +284,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); NMAX = node["nmax" ].as(); norder = node["norder"].as(); - DENS = node["dens" ].as(); if (node["nodd"]) nodd = node["nodd" ].as(); if (node["cmap"]) @@ -317,8 +315,6 @@ main(int argc, char **argv) in.read((char *)&rcylmax, sizeof(double)); in.read((char *)&rscale, sizeof(double)); in.read((char *)&vscale, sizeof(double)); - - if (idens) DENS = true; } } } @@ -330,7 +326,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = LOG; - EmpCylSL::DENS = DENS; EmpCylSL::PCAVAR = true; EmpCylSL::PCADRY = true; diff --git a/utils/Analysis/cross_validation_cyl.cc b/utils/Analysis/cross_validation_cyl.cc index b2c6ca446..b5264ef58 100644 --- a/utils/Analysis/cross_validation_cyl.cc +++ b/utils/Analysis/cross_validation_cyl.cc @@ -228,7 +228,6 @@ main(int argc, char **argv) int mmax, numx, numy, norder, cmapr, cmapz, nodd=-1; double rcylmin, rcylmax, vscale; - bool DENS; if (not ignore) { @@ -285,7 +284,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); NMAX = node["nmax" ].as(); norder = node["norder"].as(); - DENS = node["dens" ].as(); if (node["odd"]) nodd = node["odd" ].as(); if (node["cmap"]) @@ -317,8 +315,6 @@ main(int argc, char **argv) in.read((char *)&rcylmax, sizeof(double)); in.read((char *)&rscale, sizeof(double)); in.read((char *)&vscale, sizeof(double)); - - if (idens) DENS = true; } } } @@ -330,7 +326,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = LOG; - EmpCylSL::DENS = DENS; EmpCylSL::PCAVAR = true; // Create expansion diff --git a/utils/Analysis/cross_validation_cyl2.cc b/utils/Analysis/cross_validation_cyl2.cc index 60a75268e..05ff3d5cd 100644 --- a/utils/Analysis/cross_validation_cyl2.cc +++ b/utils/Analysis/cross_validation_cyl2.cc @@ -286,7 +286,6 @@ main(int argc, char **argv) int mmax, numx, numy, norder, cmapr, cmapz, nodd=-1; double rcylmin, rcylmax, vscale; - bool DENS; if (not ignore) { @@ -343,7 +342,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); NMAX = node["nmax" ].as(); norder = node["norder"].as(); - DENS = node["dens" ].as(); if (node["nodd"]) nodd = node["nodd" ].as(); if (node["cmap"]) @@ -375,8 +373,6 @@ main(int argc, char **argv) in.read((char *)&rcylmax, sizeof(double)); in.read((char *)&rscale, sizeof(double)); in.read((char *)&vscale, sizeof(double)); - - if (idens) DENS = true; } } } @@ -388,7 +384,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = LOG; - EmpCylSL::DENS = DENS; EmpCylSL::PCAVAR = true; EmpCylSL::PCADRY = true; diff --git a/utils/Analysis/diskeof.cc b/utils/Analysis/diskeof.cc index 6696eed6e..df36500a8 100644 --- a/utils/Analysis/diskeof.cc +++ b/utils/Analysis/diskeof.cc @@ -251,8 +251,6 @@ main(int argc, char **argv) sleep(20); #endif - bool DENS = false; - // ================================================== // Read basis cache // ================================================== @@ -308,7 +306,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); nmax = node["nmax" ].as(); norder = node["norder"].as(); - DENS = node["dens" ].as(); if (node["nodd"]) nodd = node["nodd" ].as(); if (node["cmap"]) @@ -340,8 +337,6 @@ main(int argc, char **argv) in.read((char *)&rcylmax, sizeof(double)); in.read((char *)&rscale, sizeof(double)); in.read((char *)&vscale, sizeof(double)); - - if (tmp) DENS = true; } } @@ -352,7 +347,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = true; - EmpCylSL::DENS = DENS; // Create expansion instance // diff --git a/utils/Analysis/diskfreqs.cc b/utils/Analysis/diskfreqs.cc index bace01f41..5f7fb72bb 100644 --- a/utils/Analysis/diskfreqs.cc +++ b/utils/Analysis/diskfreqs.cc @@ -62,7 +62,7 @@ main(int argc, char **argv) { int numx, numy, lmax=36, mmax, nmax, norder, cmapr, cmapz, numr, nodd=-1; double rcylmin, rcylmax, rscale, vscale, RMAX, Tmin, Tmax, dT, H, eps; - bool DENS, verbose = false, mask = false, ignore, logl; + bool verbose = false, mask = false, ignore, logl; std::string CACHEFILE, COEFFILE, COEFFILE2, MODEL, OUTFILE, fileType, filePrefix; // ================================================== @@ -200,7 +200,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); nmax = node["nmax" ].as(); norder = node["norder"].as(); - DENS = node["dens" ].as(); if (node["nodd"]) nodd = node["nodd" ].as(); if (node["cmap"]) @@ -226,7 +225,7 @@ main(int argc, char **argv) in.read((char *)&numy, sizeof(int)); in.read((char *)&nmax, sizeof(int)); in.read((char *)&norder, sizeof(int)); - in.read((char *)&DENS, sizeof(int)); + in.read((char *)&tmp, sizeof(int)); in.read((char *)&cmapr, sizeof(int)); in.read((char *)&rcylmin, sizeof(double)); in.read((char *)&rcylmax, sizeof(double)); @@ -242,7 +241,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = logl; - EmpCylSL::DENS = DENS; // Create expansion // diff --git a/utils/Analysis/diskprof.cc b/utils/Analysis/diskprof.cc index c01d02e21..254e48cd0 100644 --- a/utils/Analysis/diskprof.cc +++ b/utils/Analysis/diskprof.cc @@ -842,7 +842,7 @@ main(int argc, char **argv) int nice, numx, numy, lmax, mmax, nmax, norder, m1, m2, n1, n2, nodd=-1; int initc, partc, init, cmapr, cmapz; double rcylmin, rcylmax, rscale, vscale, snr, Hexp=4.0; - bool DENS, PCA, PVD, verbose = false, mask = false, ignore, logl; + bool PCA, PVD, verbose = false, mask = false, ignore, logl; std::string CACHEFILE, COEFFILE, cname, dir("."), fileType, filePrefix, config; std::string psfiles, delim; @@ -916,8 +916,6 @@ main(int argc, char **argv) ("C,center", "Accumulation center", cxxopts::value >(c0)) ("diff", "render the difference between the trimmed and untrimmed basis") - ("density", "compute density", - cxxopts::value(DENS)->default_value("true")) ("compname", "train on Component (default=stars)", cxxopts::value(cname)->default_value("stars")) ("init", "fiducial index", @@ -1125,7 +1123,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); nmax = node["nmax" ].as(); norder = node["norder"].as(); - DENS = node["dens" ].as(); if (node["nodd"]) nodd = node["nodd" ].as(); if (node["cmap"]) @@ -1151,11 +1148,7 @@ main(int argc, char **argv) in.read((char *)&numy, sizeof(int)); in.read((char *)&nmax, sizeof(int)); in.read((char *)&norder, sizeof(int)); - in.read((char *)&tmp, sizeof(int)); - if (tmp) DENS = true; - else DENS = false; - in.read((char *)&cmapr, sizeof(int)); in.read((char *)&rcylmin, sizeof(double)); in.read((char *)&rcylmax, sizeof(double)); @@ -1172,7 +1165,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = logl; - EmpCylSL::DENS = DENS; // Create expansion // diff --git a/utils/Analysis/diskprof_coef.cc b/utils/Analysis/diskprof_coef.cc index 096c9402f..08b945e6f 100644 --- a/utils/Analysis/diskprof_coef.cc +++ b/utils/Analysis/diskprof_coef.cc @@ -525,7 +525,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = true; - EmpCylSL::DENS = dens; // Create expansion // diff --git a/utils/Analysis/mssaprof_disk.cc b/utils/Analysis/mssaprof_disk.cc index 11d869795..a73f56bcf 100644 --- a/utils/Analysis/mssaprof_disk.cc +++ b/utils/Analysis/mssaprof_disk.cc @@ -401,7 +401,7 @@ main(int argc, char **argv) int lmax=36, stride=1; double rcylmin, rcylmax, rscale, vscale; - bool DENS, verbose = false, mask = false; + bool verbose = false, mask = false; std::string CACHEFILE, coeffile, config; cxxopts::Options options(argv[0], overview); @@ -509,7 +509,6 @@ main(int argc, char **argv) // ================================================== int mmax, numx, numy, nmax, norder, cmapr=1, cmapz=1, tmp, nodd=-1; - bool dens=false; double rmin, rmax, ascl, hscl; // Open EOF cachefile @@ -563,7 +562,6 @@ main(int argc, char **argv) numy = node["numy" ].as(); nmax = node["nmax" ].as(); norder = node["norder"].as(); - dens = node["dens" ].as(); if (node["nodd"]) nodd = node["nodd" ].as(); if (node["cmap"]) @@ -587,7 +585,7 @@ main(int argc, char **argv) in.read((char *)&numy, sizeof(int)); in.read((char *)&nmax, sizeof(int)); in.read((char *)&norder, sizeof(int)); - in.read((char *)&dens, sizeof(int)); if (tmp) dens = true; + in.read((char *)&tmp, sizeof(int)); in.read((char *)&cmapr, sizeof(int)); in.read((char *)&rmin, sizeof(double)); in.read((char *)&rmax, sizeof(double)); @@ -602,7 +600,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = cmapr; EmpCylSL::CMAPZ = cmapz; EmpCylSL::logarithmic = true; - EmpCylSL::DENS = dens; // Create expansion // diff --git a/utils/Analysis/pcatest.cc b/utils/Analysis/pcatest.cc index 9fdc97805..0bb9a6933 100644 --- a/utils/Analysis/pcatest.cc +++ b/utils/Analysis/pcatest.cc @@ -69,7 +69,7 @@ main(int argc, char **argv) int OUTR, OUTZ, lmax, mmax, nmax, norder, numx, numy, nodd; std::string CACHEFILE, OUTTAG; double rcylmin, rcylmax, rscale, vscale; - bool DENS, LOGSC; + bool LOGSC; // // Parse Command line @@ -114,8 +114,6 @@ main(int argc, char **argv) cxxopts::value(CACHEFILE)->default_value(".eof.cache.file")) ("outtag", "outtag for basis files", cxxopts::value(OUTTAG)->default_value("basis")) - ("density", "compute density", - cxxopts::value(DENS)->default_value("true")) ("logscale", "logscale for output basis", cxxopts::value(LOGSC)->default_value("false")) ; @@ -154,7 +152,6 @@ main(int argc, char **argv) EmpCylSL::CMAPR = 1; EmpCylSL::CMAPZ = 1; EmpCylSL::logarithmic = true; - EmpCylSL::DENS = DENS; EmpCylSL::NOUT = norder; // Create expansion diff --git a/utils/ICs/check_coefs.cc b/utils/ICs/check_coefs.cc index 8dba7e2ca..a9faf3af5 100644 --- a/utils/ICs/check_coefs.cc +++ b/utils/ICs/check_coefs.cc @@ -331,7 +331,6 @@ main(int ac, char **av) double ppower; bool SVD; int NINT; - bool DENS; bool ignore; bool orthotst; string cachefile; @@ -397,8 +396,6 @@ main(int ac, char **av) cxxopts::value(VFLAG)->default_value("0")) ("expcond", "EOF conditioning from analytic density profile", cxxopts::value(expcond)->default_value("true")) - ("DENS", "Compute density basis functions", - cxxopts::value(DENS)->default_value("true")) ("ignore", "Ignore parameters in cachefile and recompute basis", cxxopts::value(ignore)->default_value("false")) ("orthotst", "Compute orthgonality of basis functions", @@ -604,7 +601,6 @@ main(int ac, char **av) NUMY = node["numy" ].as(); NMAX = node["nmax" ].as(); NORDER = node["norder"].as(); - DENS = node["dens" ].as(); // RMIN = node["rmin" ].as(); // RMAX = node["rmax" ].as(); ASCALE = node["ascl" ].as(); @@ -632,11 +628,7 @@ main(int ac, char **av) in.read((char *)&NUMY, sizeof(int)); in.read((char *)&NMAX, sizeof(int)); in.read((char *)&NORDER, sizeof(int)); - in.read((char *)&tmp, sizeof(int)); - if (tmp) DENS = true; - else DENS = false; - in.read((char *)&CMTYPE, sizeof(int)); in.read((char *)&RCYLMIN, sizeof(double)); in.read((char *)&RCYLMAX, sizeof(double)); @@ -663,7 +655,6 @@ main(int ac, char **av) EmpCylSL::CMAPZ = CMAPZ; EmpCylSL::VFLAG = VFLAG; EmpCylSL::logarithmic = LOGR; - EmpCylSL::DENS = DENS; // For DiskDens AA = ASCALE; diff --git a/utils/ICs/check_coefs2.cc b/utils/ICs/check_coefs2.cc index f13f0cccf..fdee055b9 100644 --- a/utils/ICs/check_coefs2.cc +++ b/utils/ICs/check_coefs2.cc @@ -324,7 +324,6 @@ main(int ac, char **av) double rtrunc; double rwidth; int NINT; - bool DENS; bool ignore; bool orthotst; bool use_progress; @@ -394,8 +393,6 @@ main(int ac, char **av) cxxopts::value(VFLAG)->default_value("0")) ("expcond", "EOF conditioning from analytic density profile", cxxopts::value(expcond)->default_value("true")) - ("DENS", "Compute density basis functions", - cxxopts::value(DENS)->default_value("true")) ("ignore", "Ignore parameters in cachefile and recompute basis", cxxopts::value(ignore)->default_value("false")) ("orthotst", "Compute orthgonality of basis functions", @@ -612,7 +609,6 @@ main(int ac, char **av) NUMY = node["numy" ].as(); NMAX = node["nmax" ].as(); NORDER = node["norder"].as(); - DENS = node["dens" ].as(); // RMIN = node["rmin" ].as(); // RMAX = node["rmax" ].as(); ASCALE = node["ascl" ].as(); @@ -643,11 +639,7 @@ main(int ac, char **av) in.read((char *)&NUMY, sizeof(int)); in.read((char *)&NMAX, sizeof(int)); in.read((char *)&NORDER, sizeof(int)); - in.read((char *)&tmp, sizeof(int)); - if (tmp) DENS = true; - else DENS = false; - in.read((char *)&CMTYPE, sizeof(int)); in.read((char *)&RCYLMIN, sizeof(double)); in.read((char *)&RCYLMAX, sizeof(double)); @@ -674,7 +666,6 @@ main(int ac, char **av) EmpCylSL::CMAPZ = CMAPZ; EmpCylSL::VFLAG = VFLAG; EmpCylSL::logarithmic = LOGR; - EmpCylSL::DENS = DENS; // For DiskDens AA = ASCALE; diff --git a/utils/ICs/cylcache.cc b/utils/ICs/cylcache.cc index 34e89b902..9a2422c6f 100644 --- a/utils/ICs/cylcache.cc +++ b/utils/ICs/cylcache.cc @@ -304,7 +304,6 @@ main(int ac, char **av) int NMAX; int NCYLODD; double PPower; - bool DENS; std::string cachefile; std::string config; std::string dtype; @@ -349,8 +348,6 @@ main(int ac, char **av) cxxopts::value(expcond)->default_value("true")) ("LOGR", "Logarithmic scaling for model table in EmpCylSL", cxxopts::value(LOGR)->default_value("true")) - ("DENS", "Compute and cache basis density field", - cxxopts::value(DENS)->default_value("true")) ("RCYLMIN", "Minimum disk radius for EmpCylSL", cxxopts::value(RCYLMIN)->default_value("0.001")) ("RCYLMAX", "Maximum disk radius for EmpCylSL", @@ -546,7 +543,6 @@ main(int ac, char **av) EmpCylSL::CMAPZ = CMAPZ; EmpCylSL::VFLAG = VFLAG; EmpCylSL::logarithmic = LOGR; - EmpCylSL::DENS = DENS; // Create expansion only if needed . . . std::shared_ptr expandd; diff --git a/utils/ICs/eof_compare.cc b/utils/ICs/eof_compare.cc index 57ae8ffbc..936f19d38 100644 --- a/utils/ICs/eof_compare.cc +++ b/utils/ICs/eof_compare.cc @@ -94,7 +94,6 @@ main(int argc, char **argv) EmpCylSL::NUMY = 64; EmpCylSL::CMAPR = 1; EmpCylSL::CMAPZ = 1; - EmpCylSL::DENS = true; EmpCylSL::VFLAG = 26; EmpCylSL test1(nmax, lmax, mmax, norder, acyl, hcyl, nodd, eof1); diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index ef6cf2b2f..3287754a3 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -368,7 +368,7 @@ main(int ac, char **av) double PPower, R_DF, DR_DF; double Hratio, scale_height, scale_length, scale_lenfkN; double disk_mass, gas_mass, gscal_length, ToomreQ, Temp, Tmin; - bool const_height, images, multi, SVD, DENS, basis, zeropos, zerovel; + bool const_height, images, multi, SVD, basis, zeropos, zerovel; bool report, ignore, evolved, diskmodel; int nhalo, ndisk, ngas, ngparam; std::string hbods, dbods, gbods, suffix, centerfile, halofile1, halofile2; @@ -467,8 +467,6 @@ main(int ac, char **av) cxxopts::value(LOGR)->default_value("true")) ("CHEBY", "(boolean) Use Chebyshev smoothing for disc velocities", cxxopts::value(CHEBY)->default_value("false")) - ("DENS", "(boolean) Output density cylindrical basis tables", - cxxopts::value(DENS)->default_value("false")) ("zeropos", "(boolean) Zero center of mass", cxxopts::value(zeropos)->default_value("true")) ("zerovel", "(boolean) Zero center of velocity", @@ -872,7 +870,6 @@ main(int ac, char **av) EmpCylSL::CMAPZ = CMAPZ; EmpCylSL::VFLAG = VFLAG; EmpCylSL::logarithmic = LOGR; - EmpCylSL::DENS = DENS; EmpCylSL::PCAVAR = SELECT; // Create expansion only if needed . . . diff --git a/utils/SL/slcheck.cc b/utils/SL/slcheck.cc index bc85b6935..85a1cea9c 100644 --- a/utils/SL/slcheck.cc +++ b/utils/SL/slcheck.cc @@ -6,11 +6,12 @@ #include -#include -#include +#include #include -#include +#include +#include #include +#include int main(int argc, char** argv) { @@ -100,23 +101,39 @@ int main(int argc, char** argv) if (rmin<0.0) rmin = model.get_min_radius(); if (rmax<0.0) rmax = model.get_max_radius(); - // Generate Sturm-Liouville grid - auto ortho = std::make_shared(filename, Lmax, nmax, numr, rmin, rmax, - true, cmap, rs, 0, 1.0, cachefile, true); - // ^ ^ ^ - // | | | - // Use cache file------------------------+ | | - // | | - // Cusp extrapolation------------------------------------+ | - // | - // Turn on diagnostic output in SL creation---------------------------------+ + bool bad = false; + std::shared_ptr ortho; + + // Generate Sturm-Liouville grid + // + try { + ortho = std::make_shared(filename, Lmax, nmax, numr, rmin, rmax, + true, cmap, rs, 0, 1.0, cachefile, false); + // ^ ^ ^ + // | | | + // Use cache file----------------------+ | | + // | | + // Cusp extrapolation----------------------------------+ | + // | + // Turn on diagnostic output in SL creation-------------------------------+ + } + catch (const EXPException& err) { + if (myid==0) { + std::cerr << err.what() << std::endl; + } + bad = true; + } // Slaves exit if (use_mpi && myid>0) { MPI_Finalize(); exit(0); } + if (bad) { + MPI_Finalize(); + exit(0); + } // Do what? while (1) { bool done=false;