From 2f11d9269f1f9660278813f6a149ce258ea50dfc Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 28 Sep 2023 10:02:22 -0400 Subject: [PATCH 01/30] Allow accurate SL solutions only --- CMakeLists.txt | 4 + config_cmake.h_in | 3 + exputil/SLGridMP2.cc | 458 ++++++++++++++++++++++++++++++++++++++++- include/EXPException.H | 4 +- utils/SL/slcheck.cc | 40 ++-- 5 files changed, 484 insertions(+), 25 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c639586b7..46f059537 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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_SL_EXCEPTIONS "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_SL_EXCEPTIONS) + set(SLEDGE_THROW TRUE) +endif() if(PNG_FOUND AND ENABLE_PNG) set(HAVE_LIBPNGPP TRUE) endif() diff --git a/config_cmake.h_in b/config_cmake.h_in index ab7e61f39..8563fbe09 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_SL_EXCEPTIONS is enabled */ +#cmakedefine SLEDGE_THROW @ENABLE_SL_EXCEPTIONS@ + /* "Have C++0x (set standard to cxx14 so legacy define)*/ #define HAVE_CXX0X diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 353c60a1d..5d49cea6e 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 @@ -423,7 +447,7 @@ 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); } @@ -1280,6 +1304,60 @@ 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 iflag + for (int i=0; i0) { + + if (myid==0) { + + std::cout.precision(6); + std::cout.setf(ios::scientific); + std::cout << std::left; + + std::cout << "Tolerance errors in Sturm-Liouville solver for m=" + << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << 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]; @@ -1470,10 +1547,58 @@ 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 << "Tolerance errors in Sturm-Liouville solver for m=" + << M << " K=" << K << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << 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 << "Tolerance errors in Sturm-Liouville solver for l=" << l + << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << 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 << "Tolerance errors in Sturm-Liouville solver for l=" << L + << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << 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); } @@ -3422,7 +3654,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 +4343,63 @@ 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 << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", even" << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << 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 << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", odd" << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << 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 << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", even" << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << 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 << "Tolerance errors in Sturm-Liouville solver for Kx=" << KX + << " Ky=" << KY << ", odd" << std::endl; + + std::cout << std::setw(15) << "order" + << std::setw(15) << "ev" + << std::setw(40) << "error" + << std::endl + << std::setw(15) << "-----" + << std::setw(15) << "-----" + << std::setw(40) << "---" + << std::endl; + + for (int i=0; i -#include -#include +#include #include -#include +#include +#include #include +#include int main(int argc, char** argv) { @@ -100,23 +101,36 @@ 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, true); + // ^ ^ ^ + // | | | + // 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) exit(-1); // Do what? while (1) { bool done=false; From c57f3d37b40b5a157b8ea1c058729521abc1f181 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 28 Sep 2023 11:42:59 -0400 Subject: [PATCH 02/30] Rename the ENABLE flag for mnemonic consistency --- CMakeLists.txt | 4 ++-- config_cmake.h_in | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 46f059537..ba0cbb012 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,7 +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_SL_EXCEPTIONS "Enable *careful* Sturm-Liouville solutions" TRUE) +option(ENABLE_SLCHECK "Enable *careful* Sturm-Liouville solutions" TRUE) option(BUILD_SHARED_LIBS "Build using shared libraries" ON) option(BUILD_DOCS "Build documentation" OFF) @@ -147,7 +147,7 @@ endif() if(FFTW_FOUND) set(HAVE_FFTW TRUE) endif() -if(ENABLE_SL_EXCEPTIONS) +if(ENABLE_SLCHECK) set(SLEDGE_THROW TRUE) endif() if(PNG_FOUND AND ENABLE_PNG) diff --git a/config_cmake.h_in b/config_cmake.h_in index 8563fbe09..d1cf00db3 100644 --- a/config_cmake.h_in +++ b/config_cmake.h_in @@ -1,8 +1,8 @@ /* True if DSMC is enabled */ #cmakedefine DSMC_ENABLED @DSMC_ENABLED@ -/* Defined if ENABLE_SL_EXCEPTIONS is enabled */ -#cmakedefine SLEDGE_THROW @ENABLE_SL_EXCEPTIONS@ +/* Defined if ENABLE_SLCHECK is enabled */ +#cmakedefine SLEDGE_THROW @ENABLE_SLCHECK@ /* "Have C++0x (set standard to cxx14 so legacy define)*/ #define HAVE_CXX0X From 84be4aa7d44f0e9139274879340cd078b39bc5f2 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 28 Sep 2023 12:23:36 -0400 Subject: [PATCH 03/30] Improve output readability for user --- exputil/SLGridMP2.cc | 207 +++++++++++++++++++++++-------------------- 1 file changed, 110 insertions(+), 97 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 5d49cea6e..7133983f6 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -1321,17 +1321,18 @@ void SLGridCyl::compute_table(struct TableCyl* table, int m, int k) std::cout.precision(6); std::cout.setf(ios::scientific); std::cout << std::left; - - std::cout << "Tolerance errors in Sturm-Liouville solver for m=" - << std::endl; + + 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) << "ev" - << std::setw(40) << "error" + << std::setw(15) << "eigenvalue" + << std::setw(40) << "condition" << std::endl << std::setw(15) << "-----" - << std::setw(15) << "-----" - << std::setw(40) << "-----" + << std::setw(15) << "----------" + << std::setw(40) << "---------" << std::endl; for (int i=0; i Date: Thu, 28 Sep 2023 16:15:25 -0400 Subject: [PATCH 04/30] Updated to cleanly terminate when computing SLGrid* with MPI --- exputil/SLGridMP2.cc | 198 +++++++++++++++++++++++++++++++++---------- utils/SL/slcheck.cc | 7 +- 2 files changed, 156 insertions(+), 49 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 7133983f6..a8d6aa142 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -253,6 +253,8 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, mpi_setup(); + int totbad = 0; + if (mpi_myid) { compute_table_worker(); @@ -310,12 +312,21 @@ 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; + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_SOURCE, 10, + MPI_COMM_WORLD, &status); int retid = status.MPI_SOURCE; + totbad += bad; + 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(); // @@ -347,9 +358,18 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, while (worker) { - +#ifdef SLEDGE_THROW + int bad; + 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); +#endif mpi_unpack_table(); @@ -381,9 +401,29 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, } } - } // END Root +#ifdef SLEDGE_THROW + MPI_Bcast(&totbad, 1, MPI_INT, 0, MPI_COMM_WORLD); + + 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 { @@ -1586,17 +1626,6 @@ void SLGridCyl::compute_table_worker(void) << std::endl; } std::cout << std::endl; - - sout << std::endl - << "SLGridCyl found " << bad - << " 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"; - - throw GenericError(sout.str(), __FILE__, __LINE__); } #endif @@ -1681,6 +1710,10 @@ void SLGridCyl::compute_table_worker(void) table.m = M; table.k = K; +#ifdef SLEDGE_THROW + MPI_Send(&bad, 1, MPI_INT, 0, 10, MPI_COMM_WORLD); +#endif + int position = mpi_pack_table(&table, M, K); MPI_Send(mpi_buf, position, MPI_PACKED, 0, 11, MPI_COMM_WORLD); @@ -1934,6 +1967,8 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, mpi_setup(); + int totbad = 0; + if (mpi_myid) { // Begin workers compute_table_worker(); @@ -1983,10 +2018,20 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, // // +#ifdef SLEDGE_THROW + int bad; + 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(); @@ -2011,11 +2056,20 @@ 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; + 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--; @@ -2042,6 +2096,29 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, } // END Root + +#ifdef SLEDGE_THROW + MPI_Bcast(&totbad, 1, MPI_INT, 0, MPI_COMM_WORLD); + + 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 { @@ -3064,16 +3141,6 @@ void SLGridSph::compute_table_worker(void) } std::cout << std::endl; - sout << std::endl - << "SLGridSph found " << bad - << " 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"; - - throw GenericError(sout.str(), __FILE__, __LINE__); } #endif @@ -3132,6 +3199,10 @@ void SLGridSph::compute_table_worker(void) table.l = L; +#ifdef SLEDGE_THROW + MPI_Send(&bad, 1, MPI_INT, 0, 10, MPI_COMM_WORLD); +#endif + int position = mpi_pack_table(&table, L); MPI_Send(mpi_buf, position, MPI_PACKED, 0, 11, MPI_COMM_WORLD); @@ -3462,7 +3533,10 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, mpi_setup(); + int totbad = 0; + if (mpi_myid) { + compute_table_worker(); // @@ -3520,11 +3594,20 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, // // // - +#ifdef SLEDGE_THROW + int bad; + MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_TAG, 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(); @@ -3556,10 +3639,18 @@ 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 + int bad; + 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--; @@ -3569,7 +3660,6 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, } - // // // @@ -3594,6 +3684,26 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, } // END Root +#ifdef SLEDGE_THROW + MPI_Bcast(&totbad, 1, MPI_INT, 0, MPI_COMM_WORLD); + + 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 { @@ -4759,16 +4869,7 @@ void SLGridSlab::compute_table_worker(void) << std::endl; } std::cout << std::endl; - - sout << std::endl << "SLGridSlab found " << bad - << " 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"; - - throw GenericError(sout.str(), __FILE__, __LINE__); - } + } #endif if (tbdbg) { @@ -4932,6 +5033,9 @@ void SLGridSlab::compute_table_worker(void) table.kx = KX; table.ky = KY; +#ifdef SLEDGE_THROW + MPI_Send(&bad, 1, MPI_INT, 0, 10, MPI_COMM_WORLD); +#endif int position = mpi_pack_table(&table, KX, KY); MPI_Send(mpi_buf, position, MPI_PACKED, 0, 11, MPI_COMM_WORLD); diff --git a/utils/SL/slcheck.cc b/utils/SL/slcheck.cc index 6607d8d6b..85a1cea9c 100644 --- a/utils/SL/slcheck.cc +++ b/utils/SL/slcheck.cc @@ -109,7 +109,7 @@ int main(int argc, char** argv) // try { ortho = std::make_shared(filename, Lmax, nmax, numr, rmin, rmax, - true, cmap, rs, 0, 1.0, cachefile, true); + true, cmap, rs, 0, 1.0, cachefile, false); // ^ ^ ^ // | | | // Use cache file----------------------+ | | @@ -130,7 +130,10 @@ int main(int argc, char** argv) exit(0); } - if (bad) exit(-1); + if (bad) { + MPI_Finalize(); + exit(0); + } // Do what? while (1) { bool done=false; From 8995d7031e0cca46a1ca98f57db1ba6a8c0533df Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 28 Sep 2023 16:37:04 -0400 Subject: [PATCH 05/30] Added more comments to new code --- exputil/SLGridMP2.cc | 57 +++++++++++++++++++++++++++++++++----------- 1 file changed, 43 insertions(+), 14 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index a8d6aa142..e63003f8d 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -313,11 +313,13 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, // // #ifdef SLEDGE_THROW - int bad; + int bad; // Get count of sledge failures MPI_Recv(&bad, 1, MPI_INT, MPI_ANY_SOURCE, 10, MPI_COMM_WORLD, &status); - int retid = status.MPI_SOURCE; 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 @@ -359,10 +361,11 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, while (worker) { #ifdef SLEDGE_THROW - int bad; + 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); @@ -404,8 +407,10 @@ 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) { @@ -1347,13 +1352,14 @@ void SLGridCyl::compute_table(struct TableCyl* table, int m, int k) // Check for errors // #ifdef SLEDGE_THROW - unsigned bad = 0; // Count non-zero iflag + unsigned bad = 0; // Count non-zero sledge flags for (int i=0; i0) { if (myid==0) { @@ -1601,6 +1607,8 @@ void SLGridCyl::compute_table_worker(void) std::ostringstream sout; // Runtime exception message + // Print info if we have errors. Number of errors will be sent to + // and accumulated by the root node. if (bad>0) { std::cout.precision(6); std::cout.setf(ios::scientific); @@ -1711,6 +1719,7 @@ void SLGridCyl::compute_table_worker(void) table.k = K; #ifdef SLEDGE_THROW + // Send the failure code to the root MPI_Send(&bad, 1, MPI_INT, 0, 10, MPI_COMM_WORLD); #endif @@ -1967,7 +1976,7 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, mpi_setup(); - int totbad = 0; + int totbad = 0; // Count total number of sledge errors if (mpi_myid) { // Begin workers compute_table_worker(); @@ -2019,7 +2028,8 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, // #ifdef SLEDGE_THROW - int bad; + 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; @@ -2057,7 +2067,8 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, while (worker) { #ifdef SLEDGE_THROW - int bad; + 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; @@ -2098,8 +2109,10 @@ 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; @@ -2849,14 +2862,15 @@ void SLGridSph::compute_table(struct TableSph* table, int l) // Check for errors // #ifdef SLEDGE_THROW + // Accumulate error count unsigned bad = 0; - for (int i=0; i0) { if (myid==0) { @@ -3106,14 +3120,16 @@ void SLGridSph::compute_table_worker(void) // Check for errors // #ifdef SLEDGE_THROW + // Get sledge error count unsigned bad = 0; - for (int i=0; i0) { std::cout.precision(6); @@ -3200,9 +3216,10 @@ void SLGridSph::compute_table_worker(void) table.l = L; #ifdef SLEDGE_THROW + // Send sledge error count to root MPI_Send(&bad, 1, MPI_INT, 0, 10, MPI_COMM_WORLD); #endif - + // Send sledge comptuation to root int position = mpi_pack_table(&table, L); MPI_Send(mpi_buf, position, MPI_PACKED, 0, 11, MPI_COMM_WORLD); @@ -3533,7 +3550,7 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, mpi_setup(); - int totbad = 0; + int totbad = 0; // Accumulate total sledge error count if (mpi_myid) { @@ -3595,10 +3612,11 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, // // #ifdef SLEDGE_THROW - int bad; + 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); @@ -3640,10 +3658,12 @@ SLGridSlab::SLGridSlab(int NUMK, int NMAX, int NUMZ, double ZMAX, while (worker) { #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); @@ -3685,8 +3705,10 @@ 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; @@ -4473,6 +4495,8 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) std::ostringstream sout; // Runtime error message + // Print info if we have sledge errors and throw a runtime + // exception if (bad>0) { if (myid==0) { @@ -4578,13 +4602,14 @@ void SLGridSlab::compute_table(struct TableSlab* table, int KX, int KY) // Check for errors // #ifdef SLEDGE_THROW - bad = 0; + bad = 0; // Accumulate sledge error count for (int i=0; i0) { if (myid==0) { @@ -4843,6 +4868,8 @@ void SLGridSlab::compute_table_worker(void) std::ostringstream sout; // Runtime error message + // Print info if we have errors. Number of errors will be sent to + // and accumulated by the root node. if (bad>0) { std::cout.precision(6); @@ -4931,8 +4958,8 @@ void SLGridSlab::compute_table_worker(void) // Check for errors // #ifdef SLEDGE_THROW + // Accumulate sledge error count bad = 0; - for (int i=0; i0) { + // Print sledge error info and throw runtime exception if (myid==0) { std::cout.precision(6); @@ -5034,6 +5062,7 @@ void SLGridSlab::compute_table_worker(void) table.ky = KY; #ifdef SLEDGE_THROW + // Send sledge error count to root MPI_Send(&bad, 1, MPI_INT, 0, 10, MPI_COMM_WORLD); #endif int position = mpi_pack_table(&table, KX, KY); From 44a8c28b478b5886131e24467aaeabe52583cdc2 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 30 Sep 2023 18:35:57 -0400 Subject: [PATCH 06/30] Preliminary implementation of enforced orthogonality checks --- coefs/BasisFactory.H | 30 ++++++++++++++++ coefs/BasisFactory.cc | 77 ++++++++++++++++++++++++++++++++++++++++-- pyEXP/BasisWrappers.cc | 32 ++++++++++++++++++ 3 files changed, 137 insertions(+), 2 deletions(-) diff --git a/coefs/BasisFactory.H b/coefs/BasisFactory.H index 2b6f8f9c9..eb902f9e3 100644 --- a/coefs/BasisFactory.H +++ b/coefs/BasisFactory.H @@ -101,6 +101,18 @@ namespace BasisClasses //! Using MPI bool use_mpi; + //! Test the orthogonality of the basis (for all derived classes) + void orthoTest(const std::vector&); + + //! Sanity tolerance for orthogonality + static constexpr double orthoTol = 1.0e-2; + + //! Return readable class name + virtual const std::string classname() = 0; + + //! Subspace index + virtual const std::string harmonic() = 0; + public: //! Constructor from YAML node @@ -240,6 +252,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 @@ -363,6 +381,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 +517,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..83f5f2423 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -90,6 +90,60 @@ namespace BasisClasses } + void Basis::orthoTest(const std::vector& tests) + { + // 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 > orthoTol) { + std::cout << classname() << ": orthogonality failure" << std::endl + << std::right + << std::setw(4) << harmonic() << std::setw(16) << "Worst" << std::endl; + for (int l=0; l SphericalSL::valid_keys = { "rs", @@ -260,6 +314,9 @@ namespace BasisClasses (model_file, lmax, nmax, numr, rmin, rmax, true, cmap, rs, 0, 1, cachename); + // Test basis for consistency + orthoTest(orthoCheck(std::max(nmax*5, 100))); + // Number of possible threads int nthrds = omp_get_max_threads(); @@ -661,7 +718,6 @@ namespace BasisClasses return ret; } - SphericalSL::BasisArray SphericalSL::getBasis (double logxmin, double logxmax, int numgrid) { @@ -1089,7 +1145,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 +1164,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 <" @@ -1275,7 +1338,12 @@ namespace BasisClasses sl->generate_eof(rnum, pnum, tnum, f); } + + // Orthogonality sanity check + // + orthoTest(orthoCheck()); } + void Cylindrical::getFields (double x, double y, double z, @@ -1553,6 +1621,10 @@ namespace BasisClasses // ortho = std::make_shared(conf); + // Orthogonality sanity check + // + orthoTest(orthoCheck()); + // Get max threads // int nthrds = omp_get_max_threads(); @@ -1881,6 +1953,7 @@ namespace BasisClasses return ortho->orthoCheck(); } + FlatDisk::BasisArray FlatDisk::getBasis (double logxmin, double logxmax, int numgrid) { 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 From 83cd83d5e422e0093c01ab3496e232d0670cf0bf Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 30 Sep 2023 19:43:16 -0400 Subject: [PATCH 07/30] Fix missing factor of 1/(4*pi) in density field; correct return matrix so the diagonal is 1 for both m=0 and m>0 --- exputil/EmpCylSL.cc | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 2e0c9514b..6992cf6c6 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -1602,7 +1602,7 @@ void EmpCylSL::compute_even_odd(int request_id, int m) -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; } } } @@ -6912,7 +6912,10 @@ std::vector EmpCylSL::orthoCheck() // Combine sines and cosines // - ret[mm](n1, n2) = sqrt(0.5*(sumC*sumC + sumS*sumS)); + if (mm==0) + ret[mm](n1, n2) = sumC; + else + ret[mm](n1, n2) = sqrt(0.5*(sumC*sumC + sumS*sumS)); } } } else { From 72bc2dbdd29009ee26631bcb128f3fb02abe8195 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 30 Sep 2023 20:03:55 -0400 Subject: [PATCH 08/30] Just an additional comment before I forget --- exputil/EmpCylSL.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 6992cf6c6..3ce3cd36b 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -6910,7 +6910,9 @@ std::vector EmpCylSL::orthoCheck() } } - // Combine sines and cosines + // Combine sines and cosines. sumC and sumS should each by + // unity or zero so the returns below combine sumC and sumS + // for m>0. // if (mm==0) ret[mm](n1, n2) = sumC; From 6defe9e3a5668bf636977ad048f816aaf42e9616 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 09:46:22 -0400 Subject: [PATCH 09/30] Remove DENS flag altogether; density functions will *always* be computed --- exputil/EmpCylSL.cc | 488 ++++++++++++++++++-------------------------- include/EmpCylSL.H | 3 - 2 files changed, 199 insertions(+), 292 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 3ce3cd36b..b8975fe3d 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; @@ -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 * 0.25/M_PI; + 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 cur = DBdif["densC"][m]; - double dif = fabs(one-two); - DBdif["densC"][m] = dif>cur ? dif : cur; + double one = densC[m][v](ix, iy); + double two = p->densC[m][v](ix, iy); - 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)); } + + // Combine sines and cosines. sumC and sumS should each by + // unity or zero so the returns below combine sumC and sumS + // for m>0. + // + 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; @@ -6938,8 +6864,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; @@ -7064,11 +6988,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); @@ -7080,7 +6999,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); @@ -7109,7 +7027,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]); } } @@ -7133,7 +7051,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]); } } @@ -7160,11 +7078,6 @@ 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) @@ -7204,7 +7117,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; @@ -7229,23 +7141,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 @@ -7270,7 +7180,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]); } } @@ -7293,7 +7203,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/include/EmpCylSL.H b/include/EmpCylSL.H index d4692f5a0..9c33136ff 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -387,9 +387,6 @@ public: typedef std::shared_ptr AxiDiskPtr; - //! TRUE if density is computed (default: false) - static bool DENS; - //! TRUE if signal-to-noise methods are on (default: false) static bool PCAVAR; From d80ea383772e78f62e02201944409010f2ef7cc6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 09:46:41 -0400 Subject: [PATCH 10/30] Deprecate 'density' flag; density will always be computed --- src/Cylinder.cc | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 6c3ab77f0..6d30082c9 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -166,7 +166,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 +207,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 @@ -443,13 +437,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 +1478,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 From 1cb0a15e1ffcf281bccd89caef89c47756bd0231 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 09:46:54 -0400 Subject: [PATCH 11/30] Deprecate 'density' flag; density will always be computed --- src/Cylinder.H | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) 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) //! From 8fa362c5c38d56818c87ede06f10883dbaf47620 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 09:47:35 -0400 Subject: [PATCH 12/30] Ignore 'density' since DENS is gone from EmpCylSL --- coefs/BasisFactory.cc | 5 ----- 1 file changed, 5 deletions(-) diff --git a/coefs/BasisFactory.cc b/coefs/BasisFactory.cc index 83f5f2423..9f8d84b14 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -1202,11 +1202,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. // From 927d73d553cd5eafc747d386ac7192dadede17a7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 09:48:10 -0400 Subject: [PATCH 13/30] Remove DENS parameter --- utils/Analysis/KDcyltest.cc | 7 ------- utils/Analysis/KL_cyl.cc | 5 ----- utils/Analysis/cross_validation_cyl.cc | 5 ----- utils/Analysis/cross_validation_cyl2.cc | 5 ----- utils/Analysis/diskeof.cc | 6 ------ utils/Analysis/diskfreqs.cc | 6 ++---- utils/Analysis/diskprof.cc | 10 +--------- utils/Analysis/diskprof_coef.cc | 1 - utils/Analysis/mssaprof_disk.cc | 7 ++----- utils/Analysis/pcatest.cc | 5 +---- utils/ICs/check_coefs.cc | 9 --------- utils/ICs/check_coefs2.cc | 9 --------- utils/ICs/cylcache.cc | 4 ---- utils/ICs/eof_compare.cc | 1 - utils/ICs/initial.cc | 5 +---- 15 files changed, 7 insertions(+), 78 deletions(-) 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 . . . From 25238d3785ba3fae5ca6f66c5e903433935807b0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 10:58:47 -0400 Subject: [PATCH 14/30] Move orthoTest() function to exputils --- coefs/BasisFactory.H | 11 ++-- coefs/BasisFactory.cc | 115 +++--------------------------------------- exputil/orthoTest.cc | 70 +++++++++++++++++++++++++ 3 files changed, 81 insertions(+), 115 deletions(-) create mode 100644 exputil/orthoTest.cc diff --git a/coefs/BasisFactory.H b/coefs/BasisFactory.H index eb902f9e3..8ff10bd86 100644 --- a/coefs/BasisFactory.H +++ b/coefs/BasisFactory.H @@ -101,12 +101,6 @@ namespace BasisClasses //! Using MPI bool use_mpi; - //! Test the orthogonality of the basis (for all derived classes) - void orthoTest(const std::vector&); - - //! Sanity tolerance for orthogonality - static constexpr double orthoTol = 1.0e-2; - //! Return readable class name virtual const std::string classname() = 0; @@ -309,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); + } }; diff --git a/coefs/BasisFactory.cc b/coefs/BasisFactory.cc index 9f8d84b14..3835de2ff 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -8,6 +8,10 @@ #include #endif +extern void orthoTest(const std::vector& tests, + const std::string& classname, + const std::string& indexname); + namespace BasisClasses { @@ -90,60 +94,6 @@ namespace BasisClasses } - void Basis::orthoTest(const std::vector& tests) - { - // 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 > orthoTol) { - std::cout << classname() << ": orthogonality failure" << std::endl - << std::right - << std::setw(4) << harmonic() << std::setw(16) << "Worst" << std::endl; - for (int l=0; l SphericalSL::valid_keys = { "rs", @@ -315,7 +265,7 @@ namespace BasisClasses 0, 1, cachename); // Test basis for consistency - orthoTest(orthoCheck(std::max(nmax*5, 100))); + orthoTest(orthoCheck(std::max(nmax*5, 100)), classname(), harmonic()); // Number of possible threads int nthrds = omp_get_max_threads(); @@ -667,57 +617,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) { @@ -1336,7 +1235,7 @@ namespace BasisClasses // Orthogonality sanity check // - orthoTest(orthoCheck()); + orthoTest(orthoCheck(), classname(), harmonic()); } @@ -1618,7 +1517,7 @@ namespace BasisClasses // Orthogonality sanity check // - orthoTest(orthoCheck()); + orthoTest(orthoCheck(), classname(), harmonic()); // Get max threads // diff --git a/exputil/orthoTest.cc b/exputil/orthoTest.cc new file mode 100644 index 000000000..a0eef6ff7 --- /dev/null +++ b/exputil/orthoTest.cc @@ -0,0 +1,70 @@ +// Helper routine to check an orthogonality matrix and throw an +// exception if tolerance is not met + +#include +#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 Date: Sun, 1 Oct 2023 10:59:23 -0400 Subject: [PATCH 15/30] Move orthoTol to EXP libvars --- exputil/libvars.cc | 2 ++ include/libvars.H | 3 +++ 2 files changed, 5 insertions(+) diff --git a/exputil/libvars.cc b/exputil/libvars.cc index a2c160300..8ad8da230 100644 --- a/exputil/libvars.cc +++ b/exputil/libvars.cc @@ -34,6 +34,8 @@ namespace __EXP__ //! Sign convention element rank for eigenfunctions and eigenvectors int nevsign = 4; + //! Sanity tolerance for orthogonality + double orthoTol = 1.0e-2; }; 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 From 31094621b5d43c1422e9ecaf832e3c5315671198 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 10:59:59 -0400 Subject: [PATCH 16/30] Add ported orthtoTest to CMake --- exputil/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) From fe8abec214757b5dd75131bef86d1c8da8c86cb7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 11:00:21 -0400 Subject: [PATCH 17/30] Add orthoCheck from BasisFactory to SLGridMP --- exputil/SLGridMP2.cc | 53 +++++++++++++++++++++++++++++++++++++++++++- include/SLGridMP2.H | 4 ++++ 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index e63003f8d..6d27408b8 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -3382,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 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); From 0cbf0ef30122de2de66156877af624ed3cb40347 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 11:00:45 -0400 Subject: [PATCH 18/30] Add orthoTest to EXP side --- src/Cylinder.cc | 8 ++++++++ src/Sphere.cc | 8 ++++++++ 2 files changed, 16 insertions(+) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 6d30082c9..286c7cfd8 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -16,6 +16,10 @@ Timer timer_debug; double EXPSCALE=1.0, HSCALE=1.0, ASHIFT=0.25; +extern void orthoTest(const std::vector& tests, + const std::string& classname, + const std::string& indexname); + //@{ //! These are for testing exclusively (should be set false for production) static bool cudaAccumOverride = false; @@ -212,6 +216,10 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : ortho = std::make_shared (nmaxfid, lmaxfid, mmax, nmax, acyl, hcyl, ncylodd, cachename); + // Test for basis consistency + // + orthoTest(ortho->orthoCheck(), "Cylinder", "m"); + // Set azimuthal harmonic order restriction? // if (mlim>=0) ortho->set_mlim(mlim); diff --git a/src/Sphere.cc b/src/Sphere.cc index 3029959c8..20d6674eb 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -8,6 +8,10 @@ #include #include +extern void orthoTest(const std::vector& tests, + const std::string& classname, + const std::string& indexname); + const std::set Sphere::valid_keys = { "rs", @@ -319,6 +323,10 @@ 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 + // + orthoTest(ortho->orthoCheck(std::max(nmax*5, 100)), "Sphere", "l"); + // Update time trigger // tnext = tnow + dtime; From 4e580dfc734c85a8441c04ba2213a818edb4f3cb Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 11:08:01 -0400 Subject: [PATCH 19/30] Remove this crufty empty file --- include/testembed.cc | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 include/testembed.cc diff --git a/include/testembed.cc b/include/testembed.cc deleted file mode 100644 index e69de29bb..000000000 From 13bebf11433c86d27bbd99d5a7ca93287fbb6fec Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 11:13:40 -0400 Subject: [PATCH 20/30] Move checkTest() signature to a header file --- coefs/BasisFactory.cc | 5 +---- src/Cylinder.cc | 4 +--- src/Sphere.cc | 5 +---- 3 files changed, 3 insertions(+), 11 deletions(-) diff --git a/coefs/BasisFactory.cc b/coefs/BasisFactory.cc index 3835de2ff..120e113e1 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -2,16 +2,13 @@ #include #include #include +#include #include #ifdef HAVE_FE_ENABLE #include #endif -extern void orthoTest(const std::vector& tests, - const std::string& classname, - const std::string& indexname); - namespace BasisClasses { diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 286c7cfd8..d28352b9b 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -11,14 +11,12 @@ #include #include #include +#include Timer timer_debug; double EXPSCALE=1.0, HSCALE=1.0, ASHIFT=0.25; -extern void orthoTest(const std::vector& tests, - const std::string& classname, - const std::string& indexname); //@{ //! These are for testing exclusively (should be set false for production) diff --git a/src/Sphere.cc b/src/Sphere.cc index 20d6674eb..d79c48c6e 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -7,10 +7,7 @@ #include #include #include - -extern void orthoTest(const std::vector& tests, - const std::string& classname, - const std::string& indexname); +#include const std::set Sphere::valid_keys = { From 4979ca4b64e1307ff6178ed008a7e433b3989042 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 11:46:00 -0400 Subject: [PATCH 21/30] Add biorthogonal check to all reset locations --- src/Cylinder.cc | 11 +++++++---- src/Sphere.cc | 11 +++++++++++ 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index d28352b9b..4c3f82af1 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -214,10 +214,6 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : ortho = std::make_shared (nmaxfid, lmaxfid, mmax, nmax, acyl, hcyl, ncylodd, cachename); - // Test for basis consistency - // - orthoTest(ortho->orthoCheck(), "Cylinder", "m"); - // Set azimuthal harmonic order restriction? // if (mlim>=0) ortho->set_mlim(mlim); @@ -379,6 +375,13 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : } #endif + // Test for basis consistency + // + std::cout << "---- "; + orthoTest(ortho->orthoCheck(), "Cylinder", "m"); + + // Initialize internal variables + // ncompcyl = 0; pos.resize(nthrds); diff --git a/src/Sphere.cc b/src/Sphere.cc index d79c48c6e..a5e016bdb 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -86,6 +86,11 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : } + // Test for basis consistency + // + std::cout << "---- "; + orthoTest(ortho->orthoCheck(std::max(nmax*5, 100)), "Sphere", "l"); + setup(); } @@ -322,6 +327,7 @@ void Sphere::make_model_bin() // Test for basis consistency // + std::cout << "---- "; orthoTest(ortho->orthoCheck(std::max(nmax*5, 100)), "Sphere", "l"); // Update time trigger @@ -464,6 +470,11 @@ 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 + // + std::cout << "---- "; + orthoTest(ortho->orthoCheck(std::max(nmax*5, 100)), "Sphere", "l"); + // Update time trigger // tnext = tnow + dtime; From 8e27bd62c17764c21eaea61526f8e56a0512404c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 11:57:47 -0400 Subject: [PATCH 22/30] Missing new header file --- include/exputils.H | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 include/exputils.H 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 From 0f49664a5d754490dac0e1978033e08ed94ab1cc Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 13:24:47 -0400 Subject: [PATCH 23/30] Switched to trapezoidal rule for meridional grid quadrature instead of centered rectangles - This did not help much - Since we are linearly interpolating, doing any higher order thing will not help --- exputil/EmpCylSL.cc | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index b8975fe3d..ddeb4ac3b 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -6826,19 +6826,23 @@ std::vector EmpCylSL::orthoCheck() double sumC = 0.0, sumS = 0.0; for (int ix=0; ix<=NUMX; ix++) { - double x = XMIN + dX*ix; - double r = xi_to_r(x); + double x = XMIN + dX*ix; + double r = xi_to_r(x); + double fx = 1.0; + // Trapezoidal rule + if (ix ==0 or ix==NUMX) fx = 0.5; for (int iy=0; iy<=NUMY; iy++) { double y = YMIN + dY*iy; + double fy = 1.0; + // Trapezoidal rule + if (iy ==0 or iy==NUMX) fy = 0.5; - sumC += fac * r/d_xi_to_r(x) * d_y_to_z(y) * - potC[mm][n1](ix, iy) * densC[mm][n2](ix, iy); - - if (mm) - sumS += fac * r/d_xi_to_r(x) * d_y_to_z(y) * - potS[mm][n1](ix, iy) * densS[mm][n2](ix, iy); + double jac = fac * r/d_xi_to_r(x) * d_y_to_z(y) * fx * fy; + + sumC += jac * potC[mm][n1](ix, iy) * densC[mm][n2](ix, iy); + if (mm) sumS += jac * potS[mm][n1](ix, iy) * densS[mm][n2](ix, iy); } } From 1239586470a0352a6c26786eebd7a9b82435f756 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 15:33:18 -0400 Subject: [PATCH 24/30] Revised the default parameters --- coefs/BasisFactory.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/coefs/BasisFactory.cc b/coefs/BasisFactory.cc index 120e113e1..af31eb374 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -961,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; From 835f27443e3329588bab5ba1679b388addd85326 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 16:43:34 -0400 Subject: [PATCH 25/30] Added some comments and increased default knots for orthoTest() in Sphere --- coefs/BasisFactory.cc | 2 +- src/Cylinder.cc | 3 ++- src/Sphere.cc | 15 +++++++++------ 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/coefs/BasisFactory.cc b/coefs/BasisFactory.cc index af31eb374..120223711 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -262,7 +262,7 @@ namespace BasisClasses 0, 1, cachename); // Test basis for consistency - orthoTest(orthoCheck(std::max(nmax*5, 100)), classname(), harmonic()); + orthoTest(orthoCheck(std::max(nmax*50, 200)), classname(), harmonic()); // Number of possible threads int nthrds = omp_get_max_threads(); diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 4c3f82af1..5ee27894f 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -375,7 +375,8 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : } #endif - // Test for basis consistency + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) // std::cout << "---- "; orthoTest(ortho->orthoCheck(), "Cylinder", "m"); diff --git a/src/Sphere.cc b/src/Sphere.cc index a5e016bdb..e672f9b69 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -86,10 +86,11 @@ Sphere::Sphere(Component* c0, const YAML::Node& conf, MixtureBasis* m) : } - // Test for basis consistency + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) // std::cout << "---- "; - orthoTest(ortho->orthoCheck(std::max(nmax*5, 100)), "Sphere", "l"); + orthoTest(ortho->orthoCheck(std::max(nmax*50, 200)), "Sphere", "l"); setup(); } @@ -325,10 +326,11 @@ 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 + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) // std::cout << "---- "; - orthoTest(ortho->orthoCheck(std::max(nmax*5, 100)), "Sphere", "l"); + orthoTest(ortho->orthoCheck(std::max(nmax*50, 200)), "Sphere", "l"); // Update time trigger // @@ -470,10 +472,11 @@ 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 + // Test for basis consistency (will generate an exception if maximum + // error is out of tolerance) // std::cout << "---- "; - orthoTest(ortho->orthoCheck(std::max(nmax*5, 100)), "Sphere", "l"); + orthoTest(ortho->orthoCheck(std::max(nmax*50, 200)), "Sphere", "l"); // Update time trigger // From 6dd42efb134eaeede5078554e038a88b3e88e26f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 17:00:50 -0400 Subject: [PATCH 26/30] Make use of 'rs' and 'scale' consistent between SLGridSph and SphericalSL in pyEXP or Sphere in EXP --- exputil/SLGridMP2.cc | 78 ++++++++++++++++++++++---------------------- include/SLGridMP2.H | 14 ++++---- 2 files changed, 46 insertions(+), 46 deletions(-) diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 6d27408b8..b9c0735e7 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -201,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 RS, const std::string type, bool VERBOSE) { int m, k; @@ -217,7 +217,7 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, cache = CACHE; cmap = CMAP; - scale = SCALE; + rs = RS; tbdbg = VERBOSE; @@ -458,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, RS; std::string MODEL; if (myid==0) @@ -505,7 +505,7 @@ int SLGridCyl::read_cached_table(void) CMAP = node["cmap" ].as(); RMIN = node["rmin" ].as(); RMAX = node["rmax" ].as(); - SCL = node["scale" ].as(); + RS = node["rs" ].as(); L = node["L" ].as(); AA = node["A" ].as(); MODEL = node["model" ].as(); @@ -564,10 +564,10 @@ int SLGridCyl::read_cached_table(void) return 0; } - if (SCL!=scale) { + if (RS!=rs) { if (myid==0) - std::cout << "---- SLGridCyl::read_cached_table: found scale=" << SCL - << " wanted " << scale << std::endl; + std::cout << "---- SLGridCyl::read_cached_table: found rs=" << RS + << " wanted " << rs << std::endl; return 0; } @@ -647,7 +647,7 @@ void SLGridCyl::write_cached_table(void) node["cmap" ] = cmap; node["rmin" ] = rmin; node["rmax" ] = rmax; - node["scale" ] = scale; + node["rs" ] = rs; node["L" ] = l; node["A" ] = A; node["model" ] = cyl->ID(); @@ -714,7 +714,7 @@ double SLGridCyl::r_to_xi(double r) } if (cmap) { - return (r/scale-1.0)/(r/scale+1.0); + return (r/rs-1.0)/(r/rs+1.0); } else { return r; } @@ -735,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) * rs; } else { return xi; } @@ -757,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)/rs; } else { return 1.0; } @@ -1267,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*rs,1.0e-5, + 1.0e-4*rs,1.0e-5, + 1.0e-4*rs,1.0e-5}; int VERBOSE=0; integer NUM, N, M; logical type[8]; @@ -1472,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/rs - 1.0)/(rmin/rs + 1.0); + xmax = (rmax/rs - 1.0)/(rmax/rs + 1.0); dxi = (xmax-xmin)/(numr-1); } else { xmin = rmin; @@ -1495,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*rs,1.0e-5, + 1.0e-4*rs,1.0e-5, + 1.0e-4*rs,1.0e-5}; int i, j, VERBOSE=0; integer NUM; logical type[8]; @@ -1876,7 +1876,7 @@ void SLGridSph::bomb(string oops) SLGridSph::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 RS, int DIVERGE, double DFAC, std::string cachename, bool VERBOSE) { @@ -1892,12 +1892,12 @@ SLGridSph::SLGridSph(std::string modelname, diverge = DIVERGE; dfac = DFAC; - initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, SCALE); + initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, RS); } SLGridSph::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 RS, std::string cachename, bool VERBOSE) { mpi_buf = 0; @@ -1909,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, RS); } @@ -1939,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 RS) { int l; @@ -1952,7 +1952,7 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, cache = CACHE; cmap = CMAP; - scale = SCALE; + rs = RS; init_table(); @@ -2258,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(rs, "rs")) return false; if (not checkInt(diverge, "diverge")) return false; if (not checkDbl(dfac, "dfac")) return false; @@ -2349,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 ("rs", HighFive::DataSpace::From(rs)).write(rs); file.createAttribute ("diverge", HighFive::DataSpace::From(diverge)).write(diverge); file.createAttribute ("dfac", HighFive::DataSpace::From(dfac)).write(dfac); @@ -2391,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/rs-1.0)/(r/rs+1.0); } else if (cmap==2) { if (r<=0.0) bomb("radius <= 0!"); ret = log(r); @@ -2410,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) * rs; } else if (cmap==2) { ret = exp(xi); } else { @@ -2430,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)/rs; } else if (cmap==2) { ret = exp(-xi); } else { @@ -2781,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*rs, 1.0e-6, + 1.0e-4*rs, 1.0e-6, + 1.0e-4*rs, 1.0e-6}; int VERBOSE=0; integer NUM, N; logical type[8] = {0, 0, 1, 0, 0, 0, 1, 0}; @@ -2992,8 +2992,8 @@ void SLGridSph::init_table(void) d0.resize(numr); if (cmap==1) { - xmin = (rmin/scale - 1.0)/(rmin/scale + 1.0); - xmax = (rmax/scale - 1.0)/(rmax/scale + 1.0); + xmin = (rmin/rs - 1.0)/(rmin/rs + 1.0); + xmax = (rmax/rs - 1.0)/(rmax/rs + 1.0); } else if (cmap==2) { xmin = log(rmin); @@ -3018,9 +3018,9 @@ void SLGridSph::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-1*scale,1.0e-6, - 1.0e-1*scale,1.0e-6, - 1.0e-1*scale,1.0e-6}; + double tol [6] = {1.0e-1*rs, 1.0e-6, + 1.0e-1*rs, 1.0e-6, + 1.0e-1*rs, 1.0e-6}; int VERBOSE=0; integer NUM; @@ -3367,7 +3367,7 @@ YAML::Node SLGridSph::getHeader(const std::string& cachefile) node["cmap"] = getInt("cmap"); node["rmin"] = getDbl("rmin"); node["rmax"] = getDbl("rmax"); - node["scale"] = getDbl("scale"); + node["rs" ] = getDbl("rs"); node["diverge"] = getInt("diverge"); node["dfac"] = getDbl("dfac"); } diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 2fd11c127..1967fc415 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -58,7 +58,7 @@ private: double rmin, rmax, l; int cmap; - double scale; + double rs; 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 RS, const std::string type, bool Verbose=false); //! Destructor @@ -180,7 +180,7 @@ private: double rmin, rmax; int cmap, diverge; - double scale, dfac; + double rs, 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 RS); 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 RS, 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 RS, int DIVERGE, double DFAC, std::string cachename="", bool Verbose=false); @@ -315,7 +315,7 @@ public: { cudaMappingConstants ret; - ret.rscale = scale; + ret.rscale = rs; ret.hscale = 0.0; ret.xmin = xmin; ret.xmax = xmax; From d1782ca43a1114c0c8f6b09f638b978748e3e563 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 17:16:38 -0400 Subject: [PATCH 27/30] Version bump --- CMakeLists.txt | 2 +- doc/exp.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ba0cbb012..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) 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 From b6f58cb18023766b7a57235d4d059c50c5274c4b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 18:53:22 -0400 Subject: [PATCH 28/30] Identify source of parameter mistmatches as EmpCylSL for readablility --- exputil/EmpCylSL.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index ddeb4ac3b..6adce4ed1 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -7088,8 +7088,8 @@ bool EmpCylSL::ReadH5Cache() { 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: parameter " << name << ": wanted " + << value << " found " << v << std::endl; return false; }; @@ -7097,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: parameter " << name << ": wanted " + << value << " found " << v << std::endl; return false; }; @@ -7106,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: parameter " << name << ": wanted " + << value << " found " << v << std::endl; return false; }; From ff115598142beec050791f9566d3d13178f04e16 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 1 Oct 2023 19:06:03 -0400 Subject: [PATCH 29/30] A few minor stdout/err tweaks for readability --- exputil/EmpCylSL.cc | 10 +++++----- exputil/SLGridMP2.cc | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 6adce4ed1..55fa4a0c0 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -866,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; } @@ -7088,7 +7088,7 @@ bool EmpCylSL::ReadH5Cache() { int v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); if (value == v) return true; - std::cout << "---- EmpCylSL: parameter " << name << ": wanted " + std::cout << "---- EmpCylSL cache parameter " << name << ": wanted " << value << " found " << v << std::endl; return false; }; @@ -7097,7 +7097,7 @@ 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 << "---- EmpCylSL: parameter " << name << ": wanted " + std::cout << "---- EmpCylSL cache parameter " << name << ": wanted " << value << " found " << v << std::endl; return false; }; @@ -7106,7 +7106,7 @@ bool EmpCylSL::ReadH5Cache() { std::string v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v); if (value.compare(v)==0) return true; - std::cout << "--- EmpCylSL: parameter " << name << ": wanted " + std::cout << "--- EmpCylSL cache parameter " << name << ": wanted " << value << " found " << v << std::endl; return false; }; diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index b9c0735e7..41adae9d3 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -2282,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; From 862c60c406ab8edcd65d5d0e54cc01d3327323ff Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 2 Oct 2023 09:36:25 -0400 Subject: [PATCH 30/30] Rename 'rs' to 'rmapping' in pyEXP and EXP --- coefs/BasisFactory.H | 2 +- coefs/BasisFactory.cc | 12 ++-- exputil/SLGridMP2.cc | 136 +++++++++++++++++++++--------------------- include/SLGridMP2.H | 14 ++--- src/Sphere.H | 6 +- src/Sphere.cc | 9 +-- 6 files changed, 90 insertions(+), 89 deletions(-) diff --git a/coefs/BasisFactory.H b/coefs/BasisFactory.H index 8ff10bd86..17a88c1be 100644 --- a/coefs/BasisFactory.H +++ b/coefs/BasisFactory.H @@ -204,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; diff --git a/coefs/BasisFactory.cc b/coefs/BasisFactory.cc index 120223711..48bb5817f 100644 --- a/coefs/BasisFactory.cc +++ b/coefs/BasisFactory.cc @@ -93,7 +93,7 @@ namespace BasisClasses const std::set SphericalSL::valid_keys = { - "rs", + "rmapping", "cmap", "Lmax", "dof", @@ -173,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(); @@ -181,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(); @@ -258,7 +258,7 @@ 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 diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 41adae9d3..37ed197f9 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -201,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 RS, + bool CACHE, int CMAP, double RMAP, const std::string type, bool VERBOSE) { int m, k; @@ -217,7 +217,7 @@ SLGridCyl::SLGridCyl(int MMAX, int NMAX, int NUMR, int NUMK, cache = CACHE; cmap = CMAP; - rs = RS; + rmap = RMAP; tbdbg = VERBOSE; @@ -458,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, RS; + double RMIN, RMAX, L, AA, RMAP; std::string MODEL; if (myid==0) @@ -498,17 +498,17 @@ int SLGridCyl::read_cached_table(void) // 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(); - RS = node["rs" ].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; @@ -564,10 +564,10 @@ int SLGridCyl::read_cached_table(void) return 0; } - if (RS!=rs) { + if (RMAP!=rmap) { if (myid==0) - std::cout << "---- SLGridCyl::read_cached_table: found rs=" << RS - << " wanted " << rs << std::endl; + std::cout << "---- SLGridCyl::read_cached_table: found rmapping=" << RMAP + << " wanted " << rmap << std::endl; return 0; } @@ -640,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["rs" ] = rs; - 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 // @@ -714,7 +714,7 @@ double SLGridCyl::r_to_xi(double r) } if (cmap) { - return (r/rs-1.0)/(r/rs+1.0); + return (r/rmap-1.0)/(r/rmap+1.0); } else { return r; } @@ -735,7 +735,7 @@ double SLGridCyl::xi_to_r(double xi) bomb(ostr.str()); } - return (1.0+xi)/(1.0 - xi) * rs; + return (1.0+xi)/(1.0 - xi) * rmap; } else { return xi; } @@ -757,7 +757,7 @@ double SLGridCyl::d_xi_to_r(double xi) bomb(ostr.str()); } - return 0.5*(1.0-xi)*(1.0-xi)/rs; + return 0.5*(1.0-xi)*(1.0-xi)/rmap; } else { return 1.0; } @@ -1267,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*rs,1.0e-5, - 1.0e-4*rs,1.0e-5, - 1.0e-4*rs,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]; @@ -1472,8 +1472,8 @@ void SLGridCyl::init_table(void) d0.resize(numr); if (cmap) { - xmin = (rmin/rs - 1.0)/(rmin/rs + 1.0); - xmax = (rmax/rs - 1.0)/(rmax/rs + 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; @@ -1495,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*rs,1.0e-5, - 1.0e-4*rs,1.0e-5, - 1.0e-4*rs,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]; @@ -1876,7 +1876,7 @@ void SLGridSph::bomb(string oops) SLGridSph::SLGridSph(std::string modelname, int LMAX, int NMAX, int NUMR, double RMIN, double RMAX, - bool CACHE, int CMAP, double RS, + bool CACHE, int CMAP, double RMAP, int DIVERGE, double DFAC, std::string cachename, bool VERBOSE) { @@ -1892,12 +1892,12 @@ SLGridSph::SLGridSph(std::string modelname, diverge = DIVERGE; dfac = DFAC; - initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, RS); + initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, RMAP); } SLGridSph::SLGridSph(std::shared_ptr mod, int LMAX, int NMAX, int NUMR, double RMIN, double RMAX, - bool CACHE, int CMAP, double RS, + bool CACHE, int CMAP, double RMAP, std::string cachename, bool VERBOSE) { mpi_buf = 0; @@ -1909,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, RS); + initialize(LMAX, NMAX, NUMR, RMIN, RMAX, CACHE, CMAP, RMAP); } @@ -1939,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 RS) + bool CACHE, int CMAP, double RMAP) { int l; @@ -1952,7 +1952,7 @@ void SLGridSph::initialize(int LMAX, int NMAX, int NUMR, cache = CACHE; cmap = CMAP; - rs = RS; + rmap = RMAP; init_table(); @@ -2258,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(rs, "rs")) return false; + if (not checkDbl(rmap, "rmapping")) return false; if (not checkInt(diverge, "diverge")) return false; if (not checkDbl(dfac, "dfac")) return false; @@ -2349,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 ("rs", HighFive::DataSpace::From(rs)).write(rs); + 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); @@ -2391,7 +2391,7 @@ double SLGridSph::r_to_xi(double r) if (cmap==1) { if (r<0.0) bomb("radius < 0!"); - ret = (r/rs-1.0)/(r/rs+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); @@ -2410,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) * rs; + ret =(1.0+xi)/(1.0 - xi) * rmap; } else if (cmap==2) { ret = exp(xi); } else { @@ -2430,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)/rs; + ret = 0.5*(1.0-xi)*(1.0-xi)/rmap; } else if (cmap==2) { ret = exp(-xi); } else { @@ -2781,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*rs, 1.0e-6, - 1.0e-4*rs, 1.0e-6, - 1.0e-4*rs, 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}; @@ -2992,8 +2992,8 @@ void SLGridSph::init_table(void) d0.resize(numr); if (cmap==1) { - xmin = (rmin/rs - 1.0)/(rmin/rs + 1.0); - xmax = (rmax/rs - 1.0)/(rmax/rs + 1.0); + xmin = (rmin/rmap - 1.0)/(rmin/rmap + 1.0); + xmax = (rmax/rmap - 1.0)/(rmax/rmap + 1.0); } else if (cmap==2) { xmin = log(rmin); @@ -3018,9 +3018,9 @@ void SLGridSph::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-1*rs, 1.0e-6, - 1.0e-1*rs, 1.0e-6, - 1.0e-1*rs, 1.0e-6}; + double tol [6] = {1.0e-1*rmap, 1.0e-6, + 1.0e-1*rmap, 1.0e-6, + 1.0e-1*rmap, 1.0e-6}; int VERBOSE=0; integer NUM; @@ -3360,16 +3360,16 @@ YAML::Node SLGridSph::getHeader(const std::string& cachefile) return v; }; - node["model"] = getStr("model"); - node["lmax"] = getInt("lmax"); - node["nmax"] = getInt("nmax"); - node["numr"] = getInt("numr"); - node["cmap"] = getInt("cmap"); - node["rmin"] = getDbl("rmin"); - node["rmax"] = getDbl("rmax"); - node["rs" ] = getDbl("rs"); - node["diverge"] = getInt("diverge"); - node["dfac"] = getDbl("dfac"); + node["model"] = getStr("model"); + node["lmax"] = getInt("lmax"); + node["nmax"] = getInt("nmax"); + node["numr"] = getInt("numr"); + node["cmap"] = getInt("cmap"); + node["rmin"] = getDbl("rmin"); + node["rmax"] = getDbl("rmax"); + node["rmapping"] = getDbl("rmapping"); + node["diverge"] = getInt("diverge"); + node["dfac"] = getDbl("dfac"); } catch (YAML::Exception& error) { std::ostringstream sout; diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 1967fc415..2a04a4c44 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -58,7 +58,7 @@ private: double rmin, rmax, l; int cmap; - double rs; + 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 RS, + 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 rs, 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 RS); + 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 RS, + 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 RS, + bool cache, int Cmap, double RMAP, int DIVERGE, double DFAC, std::string cachename="", bool Verbose=false); @@ -315,7 +315,7 @@ public: { cudaMappingConstants ret; - ret.rscale = rs; + ret.rscale = rmap; ret.hscale = 0.0; ret.xmin = xmin; ret.xmax = xmax; 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 e672f9b69..51ffe5fa1 100644 --- a/src/Sphere.cc +++ b/src/Sphere.cc @@ -11,7 +11,7 @@ const std::set Sphere::valid_keys = { - "rs", + "rmapping", "numr", "nums", "cmap", @@ -29,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; @@ -61,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(); @@ -72,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 @@ -105,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();