diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index a37377ef8..b9ddde223 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -1,4 +1,4 @@ -name: "Test Builds" +name: "Build and Test" on: push: @@ -7,12 +7,13 @@ on: pull_request: branches: - main + - devel jobs: - pyexp: + exp: strategy: matrix: os: [ubuntu-latest] - cc: [gcc, mpicc] + cc: [gcc] name: "Test pyEXP Build" runs-on: ${{ matrix.os }} @@ -25,20 +26,21 @@ jobs: run: | sudo apt-get update sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev + sudo pip install numpy - name: Setup submodule and build run: | git submodule update --init --recursive mkdir -p build/install - - name: Compile pyEXP + - name: Compile EXP if: runner.os == 'Linux' env: CC: ${{ matrix.cc }} working-directory: ./build run: >- cmake - -DENABLE_NBODY=NO + -DENABLE_NBODY=YES -DENABLE_PYEXP=YES -DCMAKE_BUILD_TYPE=Release -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake @@ -48,48 +50,12 @@ jobs: - name: Make working-directory: ./build - run: make -j 2 - - # ----------------------------------------------------------------------------------- - - exp: - strategy: - matrix: - os: [ubuntu-latest] - cc: [gcc, mpicc] - - name: "Test Full EXP Build" - runs-on: ${{ matrix.os }} - steps: - - name: Checkout - uses: actions/checkout@v3 - - - name: Install core dependencies - ubuntu - if: runner.os == 'Linux' - run: | - sudo apt-get update - sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev - - - name: Setup submodule and build - run: | - git submodule update --init --recursive - mkdir -p build/install - - - name: Compile Full EXP - Linux - if: runner.os == 'Linux' - env: - CC: ${{ matrix.cc }} - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=YES - -DENABLE_PYEXP=NO - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -Wno-dev - .. + run: make -j 4 - - name: Make + - name: CTest Quick working-directory: ./build - run: make -j 2 + run: ctest -L quick + + #- name: CTest Long + #working-directory: ./build + #run: ctest -L long diff --git a/.github/workflows/buildfull.yml b/.github/workflows/buildfull.yml deleted file mode 100644 index eb9c86552..000000000 --- a/.github/workflows/buildfull.yml +++ /dev/null @@ -1,146 +0,0 @@ -name: "Test Builds" - -on: - -jobs: - pyexp: - strategy: - matrix: - os: [macos-latest, ubuntu-latest] - cc: [gcc, mpicc] - - name: "Test pyEXP Build" - runs-on: ${{ matrix.os }} - steps: - - name: Checkout - uses: actions/checkout@v3 - - - name: Install core dependencies - ubuntu - if: runner.os == 'Linux' - run: | - sudo apt-get update - sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev - - - name: Install core dependencies - mac - if: startsWith(matrix.os, 'mac') - run: | - brew update - brew reinstall gcc - brew install eigen fftw hdf5 open-mpi libomp - - - name: Setup submodule and build - run: | - git submodule update --init --recursive - mkdir -p build/install - - - name: Compile pyEXP - Linux - if: runner.os == 'Linux' - env: - CC: ${{ matrix.cc }} - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=NO - -DENABLE_PYEXP=YES - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -Wno-dev - .. - - # Note for future: The homebrew paths are for intel only. Once ARM macs are - # supported in here, we'll need to update to /opt/homebrew/... instead - - name: Compile pyEXP - Mac - if: startsWith(matrix.os, 'mac') - env: - CC: ${{ matrix.cc }} - LDFLAGS: -L/usr/local/opt/libomp/lib - CPPFLAGS: -I/usr/local/opt/libomp/include - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=NO - -DENABLE_PYEXP=YES - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/local/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -DOpenMP_CXX_INCLUDE_DIR=/usr/local/opt/libomp/include - -DOpenMP_C_INCLUDE_DIR=/usr/local/opt/libomp/include - -Wno-dev - .. - - - name: Make - working-directory: ./build - run: make -j 2 - - # ----------------------------------------------------------------------------------- - - exp: - strategy: - matrix: - os: [macos-latest, ubuntu-latest] - cc: [gcc, mpicc] - - name: "Test Full EXP Build" - runs-on: ${{ matrix.os }} - steps: - - name: Checkout - uses: actions/checkout@v3 - - - name: Install core dependencies - ubuntu - if: runner.os == 'Linux' - run: | - sudo apt-get update - sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev - - - name: Install core dependencies - mac - if: startsWith(matrix.os, 'mac') - run: | - brew update - brew reinstall gcc - brew install eigen fftw hdf5 open-mpi libomp - - - name: Setup submodule and build - run: | - git submodule update --init --recursive - mkdir -p build/install - - - name: Compile Full EXP - Linux - if: runner.os == 'Linux' - env: - CC: ${{ matrix.cc }} - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=YES - -DENABLE_PYEXP=NO - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -Wno-dev - .. - - # Note for future: The homebrew paths are for intel only. Once ARM macs are - # supported in here, we'll need to update to /opt/homebrew/... instead - - name: Compile Full EXP - Mac - if: startsWith(matrix.os, 'mac') - env: - CC: ${{ matrix.cc }} - LDFLAGS: -L/usr/local/opt/libomp/lib - CPPFLAGS: -I/usr/local/opt/libomp/include - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=YES - -DENABLE_PYEXP=NO - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/local/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -DOpenMP_CXX_INCLUDE_DIR=/usr/local/opt/libomp/include - -DOpenMP_C_INCLUDE_DIR=/usr/local/opt/libomp/include - -Wno-dev - .. - - - name: Make - working-directory: ./build - run: make -j 2 diff --git a/CMakeLists.txt b/CMakeLists.txt index 598bf4ad9..dbede3fe3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.21) # Needed for CUDA, MPI, and CTest features project( EXP - VERSION "7.7.30" + VERSION "7.7.99" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) @@ -228,7 +228,8 @@ add_subdirectory(extern/pybind11) # Set options for the HighFive git submodule in extern set(HIGHFIVE_EXAMPLES OFF CACHE BOOL "Do not build the examples") set(HIGHFIVE_BUILD_DOCS OFF CACHE BOOL "Do not build the documentation") -set(HIGHFIVE_USE_BOOST OFF CACHE BOOL "Do not use Boost in HighFIve") +set(HIGHFIVE_USE_BOOST OFF CACHE BOOL "Do not use Boost in HighFive") +set(HIGHFIVE_UNIT_TESTS OFF CACHE BOOL "Turn off internal testing for HighFIve") set(H5_USE_EIGEN TRUE CACHE BOOL "Eigen3 support in HighFive") add_subdirectory(extern/HighFive EXCLUDE_FROM_ALL) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 000000000..3e5c6ac52 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,240 @@ +# Contributing to EXP + +There are many ways to contribute to EXP. Here are some of them: + +- Blog about EXP. Cite the EXP published papers. Tell the world how + you're using EXP. This will help newcomers with more examples and + will help the EXP project to increase its visibility. +- Report bugs and request features in the [issue + tracker](https://github.com/EXP-code/EXP/issues), trying to follow + the guidelines detailed in **Reporting bugs** below. +- Submit new examples to the [EXP examples + repo](https://github.com/EXP-code/EXP-examples) or the [pyEXP + examples repo](https://github.com/EXP-code/pyEXP-examples). + +# Stable and development branches + +Our current branch policy changed as of May 1, 2023. Rather than a +`main` branch for a stable EXP/pyEXP release with development taking +place on the `devel` branch, the only official EXP branch will be +`main`. All development will take place through a [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +to `main`. All other branches are considered to be temporary. The +current stable release will be tagged in the GitHub release menu. The +HEAD of `main` will contain the latest new features and fixes. + +# Reporting bugs + +Well-written bug reports are very helpful, so keep in mind the following +guidelines when you're going to report a new bug. + +- check the [FAQ](https://exp-docs.readthedocs.io) first to see if + your issue is addressed in a well-known question +- check the [open issues](https://github.com/EXP-code/EXP/issues) + to see if the issue has already been reported. If it has, don't + dismiss the report, but check the ticket history and comments. If + you have additional useful information, please leave a comment, or + consider sending a pull request with a fix. +- write **complete, reproducible, specific bug reports**. The smaller + the test case, the better. Remember that other developers won't + have your project to reproduce the bug, so please include all + relevant files required to reproduce it. +- the most awesome way to provide a complete reproducible example is + to send a pull request which adds a failing test case to the EXP + testing suite. This is helpful even if you don't have an + intention to fix the issue yourselves. +- include the output of `exp -v` so developers working on your bug + know exactly which version and platform it occurred on, which is + often very helpful for reproducing it, or knowing if it was already + fixed. + +# Contributing to code development + +We are now using the [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +method for all EXP development, in including the code authors and +maintainers. In essence, a pull request is code patch that allows us +to easily review, test, and discuss the proposed change. The better a +patch is written, the higher the chances that it'll get accepted and +the sooner it will be merged. + +Well-written patches should: + +- contain the minimum amount of code required for the specific change. + Small patches are easier to review and merge. So, if you're doing + more than one change (or bug fix), please consider submitting one + patch per change. Do not collapse multiple changes into a single + patch. For big changes consider using a patch queue. +- pass all EXP basic example tests. See **Running tests** below. +- if you're contributing a feature, especially one that changes a + public (documented) API, please include the documentation changes + in the same patch. See **Documentation strategy** below. + +# Submitting patches + +The best way to submit a patch is to issue a [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +on GitHub. The work flow for this is: + +1. Fork the [EXP-code/EXP](https://github.com/EXP-code/EXP) repo in GitHub +2. Create a branch with the proposed change +3. Compile and test your changes +4. Once you are satisfied, create the [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +on the GitHub EXP code repo + +Remember to explain what was fixed or the new functionality (what it +is, why it's needed, etc). The more info you include, the easier will +be for core developers to understand and accept your patch. + +You can also discuss the new functionality (or bug fix) before creating +the patch, but it's always good to have a patch ready to illustrate +your arguments and show that you have put some additional thought into +the subject. A good starting point is to send a pull request on GitHub. +It can be simple enough to illustrate your idea, and leave +documentation/tests for later, after the idea has been validated and +proven useful. All functionality (including new features and bug fixes) +must include a test case to check that it works as expected, so please +include tests for your patches if you want them to get accepted sooner. + +There might be an existing pull request for the problem you'd like to +solve, so please do take a look at the existing pull request list. For +example, a pull request might be a good start, but changes are +requested by EXP maintainers, and the original pull request author +hasn't had time to address them. In this case consider picking up this +pull request: open a new pull request with all commits from the +original pull request, as well as additional changes to address the +raised issues. Doing so helps us a lot; it is not considered rude as +long as the original author is acknowledged by keeping his/her +commits. + +You can pull an existing pull request to a local branch by running +`git fetch https://github.com/EXP-code/code pull/$PR_NUMBER/head:$BRANCH_NAME_TO_CREATE` +(replace `$PR_NUMBER` with an ID of the pull request, and +`$BRANCH_NAME_TO_CREATE` with a name of the branch you want to create +locally). See also: +. + +When writing GitHub pull requests, try to keep titles short but +descriptive. Complete titles make it easy to skim through the issue +tracker. + +Finally, we all appreciate improving the code's readability and though +formatting and improved comments. But please, try to keep aesthetic +changes in separate commits from functional changes. This will make pull +requests easier to review and more likely to get merged. + +## Adding a Git Commit + +All code updates are done by pull request (PR) as described above. +The GitHub EXP-code includes Continuous Integration (CI) support which +can get very chatty and take up unnecessary compute resources. When +your changes only affect documentation (e.g. docstring additions or +software comments) or provide code snippets that you know still need +work, you may add a [no ci] in your commit message. For example: bash + +``` +> git commit -m "Updated some docstrings [no ci]" +``` + +When this commit is pushed out to your fork and branch associated with a +pull request, all CI will be skipped. + +# Documentation strategy + +EXP and pyEXP has three types of documentation: + +1. Inline Doxygen comments in C/C++ class headers. If you have not + used [Doxgyen](http://doxygen.org) in the past, we recommend that + you simply copy the style in existing headers in the `include` + directory in the EXP source tree. In short, Doxygen uses stylized + comments to extract documentation. Lines that start with `//!` or + blocks that start with `/**` and end with `*/` will end up + describing the immediately following class or member function. +2. Python wrapper code has embedded Python docstrings. For some + examples, check the C++ code in the `pyEXP`. +3. The ReadTheDocs manual that you are currently reading is designed to + provide an overview and tutorial for using EXP/pyEXP. You can fork + and issue pull requests against the `EXP-code/docs` repo just for + the EXP source code as described in + **Submitting-patches** above. + +## Contributing Documentation + +Our goal is a set of consistent documentation that provides users and +developers a shallow learning curve for using and adding to EXP. For end +users, we strive to write simple step-by-step instructions for common +tasks and give a clear description of the features and capabilities of +the EXP software. + +However, it *is* hard to know what everyone needs. As you work with this +package, we would love help with this and encourage all of your +contributions. + +This section is an attempt to provide some stylistic guidelines for +documentation writers and developers. For EXP and the documentation +overall, we hope that the existing documentation is a good starting +point. For internal Python documentation in pyEXP, we are trying to +follow the now familiar style of code documentation of both the Astropy +and Numpy projects. In particular, we have adopted the Numpy style +guidelines +[here](https://numpy.org/doc/devdocs/docs/howto_document.html). + +## Adding pyEXP docstrings + +pyEXP is a set of bindings implemented by +[pybind11](https://pybind11.readthedocs.io/) and a small amount of +native Python code. It has a full set of docstrings to provide users +with easy access to interactive tips and call signatures. If you would +like to contribute new code, please try to follow the following +guidelines: + +- Please write docstrings for all public classes, methods, and + functions +- Please consult the + [numpydoc](https://numpy.org/doc/devdocs/docs/howto_document.html) + format from the link above whenever possible +- We would like to encourage references to internal EXP project links + in docstrings when useful. This is still a work in progress. +- Examples and/or tutorials are strongly encouraged for typical + use-cases of a particular module or class. These can be included in + the [EXP examples + repository](https://github.com/EXP-code/EXP-examples) or the [pyEXP + examples repository](https://github.com/EXP-code/pyEXP-examples)) as + appropriate. + +# Writing examples + +We strongly encourage to contribute interesting examples and workflows +to one of the two example repositories: + +1. Use the [EXP examples + repo](https://github.com/EXP-code/EXP-examples) to illustrate + simulations. Check the existing ones for guidance. It's best if + your examples contain sample configuration files and phase-space + files, or possibly instructions for computing the phase-space files + if they are large. +2. Use the [pyEXP examples + repo](https://github.com/EXP-code/pyEXP-examples). Either documented + Python or Jupyter notebooks are ideal. + +# Running tests + +Please check any bug fixes and proposed new functionality works with +existing examples. Here are some suggested guidelines for EXP and pyEXP +changes, respectively: + +1. For EXP, clone the [EXP examples + repo](https://github.com/EXP-code/EXP-examples) to test changes to + the EXP simulation code. At the very least, please try the `Nbody` + example, but please try as many as you can to make sure that your + change will not break an existing use case for someone else. If your + change introduces a feature, please try to devise and contribute + example that demonstrates the new feature. That way, your new + feature will be tested by all future proposed changes and will help + others understand how to use your new feature. +2. The drill for pyEXP is similar. Clone [pyEXP examples + repo](https://github.com/EXP-code/pyEXP-examples) to test changes to + the pyEXP N-body analysis code. There are many work flows here, we + don't expect anyone to try all of them. But please use your best + judgment to try those affected by your proposed change. diff --git a/doc/exp.cfg b/doc/exp.cfg index b7ccb4716..32f7444ff 100644 --- a/doc/exp.cfg +++ b/doc/exp.cfg @@ -48,7 +48,7 @@ PROJECT_NAME = EXP # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 7.7.30 +PROJECT_NUMBER = 7.7.99 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/expui/CoefContainer.H b/expui/CoefContainer.H index 66a3c9951..0c8bf39d0 100644 --- a/expui/CoefContainer.H +++ b/expui/CoefContainer.H @@ -91,6 +91,7 @@ namespace MSSA void pack_sphere(); void unpack_sphere(); void restore_background_sphere(); + //@} //@{ //! Cylindrical coefficients @@ -99,6 +100,18 @@ namespace MSSA void restore_background_cylinder(); //@} + //@{ + //! Spherical field coefficients + void pack_sphfld(); + void unpack_sphfld(); + //@} + + //@{ + //! Cylindrical coefficients + void pack_cylfld(); + void unpack_cylfld(); + //@} + //@{ //! Slab coefficients void pack_slab(); diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index 03af71d92..d5deb13da 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include @@ -64,8 +65,14 @@ namespace MSSA else if (dynamic_cast(coefs.get())) { unpack_table(); } + else if (dynamic_cast(coefs.get())) { + unpack_sphfld(); + } + else if (dynamic_cast(coefs.get())) { + unpack_cylfld(); + } else { - throw std::runtime_error("CoefDB::unpack_channels(): can not reflect coefficient type"); + throw std::runtime_error(std::string("CoefDB::unpack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); } return coefs; @@ -79,8 +86,12 @@ namespace MSSA restore_background_cylinder(); else if (dynamic_cast(coefs.get())) { } // Do nothing + else if (dynamic_cast(coefs.get())) + { } // Do nothing + else if (dynamic_cast(coefs.get())) + { } // Do nothing else { - throw std::runtime_error("CoefDB::background(): can not reflect coefficient type"); + throw std::runtime_error(std::string("CoefDB::background(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); } } @@ -96,8 +107,12 @@ namespace MSSA pack_cube(); else if (dynamic_cast(coefs.get())) pack_table(); + else if (dynamic_cast(coefs.get())) + pack_sphfld(); + else if (dynamic_cast(coefs.get())) + pack_cylfld(); else { - throw std::runtime_error("CoefDB::pack_channels(): can not reflect coefficient type"); + throw std::runtime_error(std::string("CoefDB::pack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); } } @@ -108,7 +123,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int mmax = cf->mmax; int nmax = cf->nmax; @@ -165,7 +181,9 @@ namespace MSSA } for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { if (k[2]==0) data[k][t] = (*cf->coefs)(k[0], k[1]).real(); @@ -185,7 +203,8 @@ namespace MSSA void CoefDB::unpack_cylinder() { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -201,6 +220,89 @@ namespace MSSA // END time loop } + void CoefDB::pack_cylfld() + { + auto cur = dynamic_cast(coefs.get()); + + times = cur->Times(); + complexKey = true; + + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); + + int nfld = cf->nfld; + int mmax = cf->mmax; + int nmax = cf->nmax; + int ntimes = times.size(); + + // Promote desired keys into c/s pairs + // + keys.clear(); + for (auto v : keys0) { + // Sanity check rank + // + if (v.size() != 3) { + std::ostringstream sout; + sout << "CoefDB::pack_cylfld: key vector should have rank 3; " + << "found rank " << v.size() << " instead"; + throw std::runtime_error(sout.str()); + } + // Sanity check values + // + else if (v[0] >= 0 and v[0] < nfld and + v[1] >= 0 and v[1] <= mmax and + v[2] >= 0 and v[2] <= nmax ) { + auto c = v, s = v; + c.push_back(0); + s.push_back(1); + keys.push_back(c); + if (v[0]) keys.push_back(s); + } else { + throw std::runtime_error("CoefDB::pack_cylfld: key is out of bounds"); + } + } + + bkeys.clear(); // No background fields + + // Only pack the keys in the list + // + for (auto k : keys) { + data[k].resize(ntimes); + } + + for (int t=0; t + ( cur->getCoefStruct(times[t]).get() ); + + for (auto k : keys) { + if (k[3]==0) + data[k][t] = (*cf->coefs)(k[0], k[1], k[2]).real(); + else + data[k][t] = (*cf->coefs)(k[0], k[1], k[3]).imag(); + } + } + } + + void CoefDB::unpack_cylfld() + { + for (int i=0; i + ( coefs->getCoefStruct(times[i]).get() ); + + for (auto k : keys0) { + auto c = k, s = k; + c.push_back(0); s.push_back(1); + + int f = k[1], m = k[1], n = k[2]; + + if (m==0) (*cf->coefs)(f, m, n) = {data[c][i], 0.0}; + else (*cf->coefs)(f, m, n) = {data[c][i], data[s][i]}; + } + // END key loop + } + // END time loop + } + void CoefDB::pack_sphere() { auto cur = dynamic_cast(coefs.get()); @@ -208,7 +310,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int lmax = cf->lmax; int nmax = cf->nmax; @@ -272,7 +375,9 @@ namespace MSSA auto I = [](const Key& k) { return k[0]*(k[0]+1)/2 + k[1]; }; for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { auto c = (*cf->coefs)(I(k), k[2]); data[k][t] = c.real(); @@ -293,7 +398,8 @@ namespace MSSA for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -309,6 +415,99 @@ namespace MSSA } // END time loop } + + void CoefDB::pack_sphfld() + { + auto cur = dynamic_cast(coefs.get()); + + times = cur->Times(); + complexKey = true; + + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); + + int nfld = cf->nfld; + int lmax = cf->lmax; + int nmax = cf->nmax; + int ntimes = times.size(); + + // Make extended key list + // + keys.clear(); + for (auto k : keys0) { + // Sanity check rank + // + if (k.size() != 4) { + std::ostringstream sout; + sout << "CoefDB::pack_sphfld: key vector should have rank 4; " + << "found rank " << k.size() << " instead"; + throw std::runtime_error(sout.str()); + } + // Sanity check values + // + else if (k[0] < nfld and k[0] >= 0 and + k[1] <= lmax and k[1] >= 0 and + k[2] <= k[1] and k[2] >= 0 and + k[3] < nmax and k[3] >= 0 ) { + + auto v = k; + v.push_back(0); + keys.push_back(v); + data[v].resize(ntimes); + + if (k[2]>0) { + v[4] = 1; + keys.push_back(v); + data[v].resize(ntimes); + } + } + else { + throw std::runtime_error("CoefDB::pack_sphfld: key is out of bounds"); + } + } + + bkeys.clear(); // No background for field quantities + + auto I = [](const Key& k) { return k[1]*(k[1]+1)/2 + k[2]; }; + + for (int t=0; t + ( cur->getCoefStruct(times[t]).get() ); + + for (auto k : keys) { + auto c = (*cf->coefs)(k[0], I(k), k[3]); + data[k][t] = c.real(); + if (k[4]) data[k][t] = c.imag(); + } + } + } + + void CoefDB::unpack_sphfld() + { + auto I = [](const Key& k) { return k[1]*(k[1]+1)/2 + k[2]; }; + + for (int i=0; i + ( coefs->getCoefStruct(times[i]).get() ); + + for (auto k : keys0) { + + auto c = k, s = k; + + c.push_back(0); + s.push_back(1); + + int f = k[0], l = k[1], m = k[2], n = k[3]; + + if (m==0) (*cf->coefs)(f, I(k), n) = {data[c][i], 0.0 }; + else (*cf->coefs)(f, I(k), n) = {data[c][i], data[s][i]}; + } + // END key loop + } + // END time loop + } void CoefDB::pack_slab() { @@ -317,7 +516,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int nmaxx = cf->nmaxx; int nmaxy = cf->nmaxy; @@ -376,7 +576,9 @@ namespace MSSA } for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { auto c = (*cf->coefs)(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); @@ -398,7 +600,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int nmaxx = cf->nmaxx; int nmaxy = cf->nmaxy; @@ -457,7 +660,9 @@ namespace MSSA } for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { auto c = (*cf->coefs)(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); @@ -476,7 +681,8 @@ namespace MSSA { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -494,7 +700,8 @@ namespace MSSA { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -515,7 +722,8 @@ namespace MSSA times = cur->Times(); complexKey = false; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int cols = cf->cols; int ntimes = times.size(); @@ -552,7 +760,9 @@ namespace MSSA for (unsigned c=0; c( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + data[key][t] = (*cf->coefs)(c).real(); } } @@ -562,7 +772,8 @@ namespace MSSA { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); int cols = cf->cols; @@ -680,6 +891,7 @@ namespace MSSA auto I = [](const Key& k) { return k[0]*(k[0]+1)/2 + k[1]; }; for (int t=0; t (cur->getCoefStruct(times[t]).get()); @@ -700,6 +912,7 @@ namespace MSSA auto cur = dynamic_cast(coefs.get()); for (int t=0; t (cur->getCoefStruct(times[t]).get()); diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H index d405c3179..25ef5dbe3 100644 --- a/expui/FieldBasis.H +++ b/expui/FieldBasis.H @@ -53,9 +53,9 @@ namespace BasisClasses //@{ //! Parameters - std::string model, filename; + std::string model, modelname; int lmax, nmax; - double rmin, rmax, ascl, scale, delta; + double rmin, rmax, ascl, rmapping, delta; //@} //@{ @@ -116,8 +116,8 @@ namespace BasisClasses //! Register phase-space functionoid void addPSFunction(PSFunction func, std::vector& labels); - //! Prescaling factor - void set_scale(const double scl) { scale = scl; } + //! Coordinate mapping factor + void set_scale(const double scl) { rmapping = scl; } //! Zero out coefficients to prepare for a new expansion void reset_coefs(void); diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index eecdcc520..7c85e0aee 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -11,16 +11,23 @@ #include #endif -extern double Ylm01(int ll, int mm); extern double plgndr(int l, int m, double x); +static double Ylm(int l, int m) +{ + m = abs(m); + return sqrt( (2.0*l+1)/(4.0*M_PI) ) * + exp(0.5*(lgamma(1.0+l-m) - lgamma(1.0+l+m))); +} + + namespace BasisClasses { const std::set FieldBasis::valid_keys = { - "filename", + "modelname", "dof", - "scale", + "rmapping", "rmin", "rmax", "ascl", @@ -60,18 +67,18 @@ namespace BasisClasses // void FieldBasis::configure() { - nfld = 2; // Weight and density fields by default - lmax = 4; - nmax = 10; - rmin = 1.0e-4; - rmax = 2.0; - ascl = 0.01; - delta = 0.005; - scale = 0.05; - dof = 3; - model = "file"; - name = "field"; - filename = "SLGridSph.model"; + nfld = 2; // Weight and density fields by default + lmax = 4; + nmax = 10; + rmin = 1.0e-4; + rmax = 2.0; + ascl = 0.01; + delta = 0.005; + rmapping = 0.05; + dof = 3; + model = "file"; + name = "field"; + modelname = "SLGridSph.model"; initialize(); @@ -111,8 +118,8 @@ namespace BasisClasses // if (model == "file") { std::vector r, d; - std::ifstream in(filename); - if (not in) throw std::runtime_error("Error opening file: " + filename); + std::ifstream in(modelname); + if (not in) throw std::runtime_error("Error opening file: " + modelname); std::string line; while (std::getline(in, line)) { @@ -156,13 +163,30 @@ namespace BasisClasses // Generate the orthogonal function instance // - ortho = std::make_shared(nmax, densfunc, rmin, rmax, scale, dof); + ortho = std::make_shared(nmax, densfunc, rmin, rmax, rmapping, dof); // Initialize fieldlabels // fieldLabels.clear(); fieldLabels.push_back("weight"); fieldLabels.push_back("density"); + + // Debug + // + if (true) { + auto tst = ortho->testOrtho(); + double worst = 0.0; + for (int i=0; i(worst, fabs(1.0 - tst(i, j))); + else worst = std::max(worst, fabs(tst(i, j))); + } + } + if (myid==0) { + std::cout << "FieldBasis::orthoTest: worst=" << worst << std::endl; + ortho->dumpOrtho("fieldbasis_ortho.dump"); + } + } } void FieldBasis::allocateStore() @@ -187,17 +211,17 @@ namespace BasisClasses // Assign values from YAML // try { - if (conf["filename"]) filename = conf["filename"].as(); - if (conf["nfld" ]) nfld = conf["nfld" ].as(); - if (conf["lmax" ]) lmax = conf["lmax" ].as(); - if (conf["nmax" ]) nmax = conf["nmax" ].as(); - if (conf["dof" ]) dof = conf["dof" ].as(); - if (conf["rmin" ]) rmin = conf["rmin" ].as(); - if (conf["rmax" ]) rmax = conf["rmax" ].as(); - if (conf["ascl" ]) ascl = conf["ascl" ].as(); - if (conf["delta" ]) delta = conf["delta" ].as(); - if (conf["scale" ]) scale = conf["scale" ].as(); - if (conf["model" ]) model = conf["model" ].as(); + if (conf["modelname"]) modelname = conf["modelname"].as(); + if (conf["model" ]) model = conf["model" ].as(); + if (conf["nfld" ]) nfld = conf["nfld" ].as(); + if (conf["lmax" ]) lmax = conf["lmax" ].as(); + if (conf["nmax" ]) nmax = conf["nmax" ].as(); + if (conf["dof" ]) dof = conf["dof" ].as(); + if (conf["rmin" ]) rmin = conf["rmin" ].as(); + if (conf["rmax" ]) rmax = conf["rmax" ].as(); + if (conf["ascl" ]) ascl = conf["ascl" ].as(); + if (conf["delta" ]) delta = conf["delta" ].as(); + if (conf["rmapping" ]) rmapping = conf["rmapping" ].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in FieldBasis: " @@ -254,9 +278,33 @@ namespace BasisClasses if (dof==2) { auto p = dynamic_cast(c.get()); coefs[0] = p->coefs; + store[0] = p->store; + + // Sanity test dimensions + if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) { + std::ostringstream serr; + serr << "FieldBasis::set_coefs: dimension error! " + << " nfld [" << nfld << "!= " << p->nfld << "]" + << " mmax [" << lmax << "!= " << p->mmax << "]" + << " nmax [" << nmax << "!= " << p->nmax << "]"; + throw std::runtime_error(serr.str()); + } + } else { auto p = dynamic_cast(c.get()); coefs[0] = p->coefs; + store[0] = p->store; + + // Sanity test dimensions + if (nfld!=p->nfld || lmax!=p->lmax || nmax!=p->nmax) { + std::ostringstream serr; + serr << "FieldBasis::set_coefs: dimension error! " + << " nfld [" << nfld << "!= " << p->nfld << "]" + << " lmax [" << lmax << "!= " << p->lmax << "]" + << " nmax [" << nmax << "!= " << p->nmax << "]"; + throw std::runtime_error(serr.str()); + } + } } @@ -265,6 +313,8 @@ namespace BasisClasses double u, double v, double w) { constexpr std::complex I(0, 1); + constexpr double fac0 = 0.25*M_2_SQRTPI; + int tid = omp_get_thread_num(); PS3 pos{x, y, z}, vel{u, v, w}; @@ -287,12 +337,16 @@ namespace BasisClasses auto p = (*ortho)(R); - (*coefs[tid])(0, 0, 0) += mass*p(0); + (*coefs[tid])(0, 0, 0) += mass*p(0)*fac0; for (int m=0; m<=lmax; m++) { - std::complex P = std::exp(I*(phi*m)); + + std::complex P = std::exp(-I*(phi*m)); + for (int n=0; n P = - std::exp(I*(phi*m))*Ylm01(l, m)*plgndr(l, m, cth); - + std::exp(-I*(phi*m))*Ylm(l, m)*plgndr(l, m, cth) * s; + + s *= -1.0; // Flip sign for next evaluation + for (int n=0; n0 ? 2.0 : 1.0; for (int n=0; ndimensions(); + fout << "Dim size: " << d.size(); + for (int i=0; i ret(nfld, 0); auto p = (*ortho)(r); for (int i=0; i0 ? 2.0 : 1.0; for (int n=0; n::type>(type); - file.createAttribute("trendType", HighFive::DataSpace::From(trend)).write(trend); + file.createAttribute("trendType", HighFive::DataSpace::From(trend)).write(trend); // Save the key list // @@ -1536,7 +1536,18 @@ namespace MSSA { reconstructed = true; } - } catch (HighFive::Exception& err) { + } catch (HighFive::Exception& err) { // Number of channels + // + nkeys = mean.size(); + + // Make sure parameters are sane + // + if (numW<=0) numW = numT/2; + if (numW > numT/2) numW = numT/2; + + numK = numT - numW + 1; + + std::cerr << "**** Error opening or reading H5 file ****" << std::endl; throw; } @@ -1564,8 +1575,8 @@ namespace MSSA { // if (params["Traj"]) trajectory = params["Traj"].as(); - std::cout << "Trajectory is " << std::boolalpha << trajectory - << std::endl; + // std::cout << "Trajectory is " << std::boolalpha << trajectory + // << std::endl; // Eigen OpenMP reporting // diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 23d579c0f..e69be02aa 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -216,6 +216,128 @@ EmpCylSL::EmpCylSL(int nmax, int lmax, int mmax, int nord, maxSNR = 0.0; } +template +U getH5(std::string name, HighFive::File& file) +{ + if (file.hasAttribute(name)) { + U v; + HighFive::Attribute vv = file.getAttribute(name); + vv.read(v); + return v; + } else { + std::ostringstream sout; + sout << "EmpCylSL: could not find <" << name << ">"; + throw std::runtime_error(sout.str()); + } +} + +EmpCylSL::EmpCylSL(int mlim, std::string cachename) +{ + // Use default name? + // + if (cachename.size()==0) + throw std::runtime_error("EmpCylSL: you must specify a cachename"); + + cachefile = cachename; + + // Open and read the cache file to get the needed input parameters + // + + try { + // Silence the HDF5 error stack + // + HighFive::SilenceHDF5 quiet; + + // Open the hdf5 file + // + HighFive::File file(cachename, HighFive::File::ReadOnly); + + // For basis ID + std::string forceID("Cylinder"), geometry("cylinder"); + std::string model = EmpModelLabs[mtype]; + + // Version check + // + if (file.hasAttribute("Version")) { + auto ver = getH5(std::string("Version"), file); + if (ver.compare(Version)) + throw std::runtime_error("EmpCylSL: version mismatch"); + } else { + throw std::runtime_error("EmpCylSL: outdated cache"); + } + + NMAX = getH5("nmaxfid", file); + MMAX = getH5("mmax", file); + LMAX = getH5("lmaxfid", file); + NORDER = getH5("nmax", file); + CMAPR = getH5("cmapr", file); + CMAPZ = getH5("cmapz", file); + MMIN = 0; + MLIM = std::min(mlim, MMAX); + NMIN = 0; + NLIM = std::numeric_limits::max(); + Neven = getH5("neven", file); + Nodd = getH5("nodd", file); + ASCALE = getH5("ascl", file); + HSCALE = getH5("hscl", file); + } + catch (YAML::Exception& error) { + std::ostringstream sout; + sout << "EmpCylSL::getHeader: invalid cache file <" << cachefile << ">. "; + sout << "YAML error in getHeader: " << error.what() << std::endl; + throw GenericError(sout.str(), __FILE__, __LINE__, 1038, false); + } + + // Set EvenOdd if values seem sane + // + EvenOdd = false; + if (Nodd>=0 and Nodd<=NORDER and Nodd+Neven==NORDER) { + EvenOdd = true; + } + + pfac = 1.0/sqrt(ASCALE); + ffac = pfac/ASCALE; + dfac = ffac/ASCALE; + + EVEN_M = false; + + // Check whether MPI is initialized + // + int flag; + MPI_Initialized(&flag); + if (flag) use_mpi = true; + else use_mpi = false; + + // Enable MPI code in SLGridSph + // + if (use_mpi and numprocs>1) SLGridSph::mpi = 1; + + // Choose table dimension + // + MPItable = 4; + + // Initialize storage and values + // + coefs_made.resize(multistep+1); + std::fill(coefs_made.begin(), coefs_made.end(), false); + + eof_made = false; + + sampT = 1; + defSampT = 1; + tk_type = None; + + cylmass = 0.0; + cylmass1 = std::vector(nthrds); + cylmass_made = false; + + hallfile = ""; + minSNR = std::numeric_limits::max(); + maxSNR = 0.0; + + read_cache(); +} + void EmpCylSL::reset(int numr, int lmax, int mmax, int nord, double ascale, double hscale, int nodd, @@ -815,7 +937,7 @@ std::string compare_out(std::string str, U one, U two) return sout.str(); } -int EmpCylSL::cache_grid(int readwrite, string cachename) +int EmpCylSL::cache_grid(int readwrite, std::string cachename) { // Option to reset cache file name diff --git a/exputil/OrthoFunction.cc b/exputil/OrthoFunction.cc index 7a1badc20..e0bf4f1f9 100644 --- a/exputil/OrthoFunction.cc +++ b/exputil/OrthoFunction.cc @@ -1,3 +1,4 @@ +#include #include OrthoFunction::OrthoFunction @@ -93,6 +94,32 @@ Eigen::MatrixXd OrthoFunction::testOrtho() return ret; } +void OrthoFunction::dumpOrtho(const std::string& filename) +{ + std::ofstream fout(filename); + + if (fout) { + fout << "# OrthoFunction dump" << std::endl; + + const int number = 1000; + + for (int i=0; i mod, } +SLGridSph::SLGridSph(std::string cachename) +{ + if (cachename.size()) sph_cache_name = cachename; + else throw std::runtime_error("SLGridSph: you must specify a cachename"); + + tbdbg = false; + + int LMAX, NMAX, NUMR, CMAP, DIVERGE=0; + double RMIN, RMAX, RMAP, DFAC=1.0; + + try { + + auto node = getHeader(cachename); + + LMAX = node["lmax"].as(); + NMAX = node["nmax"].as(); + NUMR = node["numr"].as(); + CMAP = node["cmap"].as(); + RMIN = node["rmin"].as(); + RMAX = node["rmax"].as(); + RMAP = node["rmapping"].as(); + + model_file_name = node["model"].as(); + model = SphModTblPtr(new SphericalModelTable(model_file_name, diverge, dfac)); + } + catch (YAML::Exception& error) { + std::ostringstream sout; + sout << "SLGridMP2: error parsing parameters from getHeader: " + << error.what(); + throw GenericError(sout.str(), __FILE__, __LINE__, 1039, false); + } + + initialize(LMAX, NMAX, NUMR, RMIN, RMAX, false, CMAP, RMAP); +} + + std::map SLGridSph::cacheInfo(const std::string& cachefile, bool verbose) { diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 30a4a6587..233039ef7 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -501,6 +501,9 @@ public: double ascale, double hscale, int Nodd, std::string cachename); + //! Construct from cache file + EmpCylSL(int mlim, const std::string cache); + //! Destructor ~EmpCylSL(void); diff --git a/include/OrthoFunction.H b/include/OrthoFunction.H index b6f725ad6..6cc90ad5d 100644 --- a/include/OrthoFunction.H +++ b/include/OrthoFunction.H @@ -100,6 +100,9 @@ public: //! Test orthogonality Eigen::MatrixXd testOrtho(); + + //! Dump the orthogonal function table + void dumpOrtho(const std::string& filename); }; diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 95a9bf6b7..f83d98fe9 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -116,6 +116,10 @@ public: std::string cachename=".slgrid_sph_cache", bool Verbose=false); + //! Constructor (from cache file) + SLGridSph(std::string cachename); + + //! Destructor virtual ~SLGridSph(); @@ -216,6 +220,11 @@ public: double getRmax() { return rmax; } //@} + //@{ + //! Get expansion limits + int getLmax() { return lmax; } + int getNmax() { return nmax; } + //@} }; diff --git a/pyEXP/MSSAWrappers.cc b/pyEXP/MSSAWrappers.cc index 9fdc5658f..2b158ff2c 100644 --- a/pyEXP/MSSAWrappers.cc +++ b/pyEXP/MSSAWrappers.cc @@ -326,7 +326,6 @@ void MSSAtoolkitClasses(py::module &m) { dict({id: Coefs},...) reconstructed time series in the original coefficient form - Notes ----- The reconstructed data will overwrite the memory of the original coefficient @@ -542,7 +541,6 @@ void MSSAtoolkitClasses(py::module &m) { power value )"); - f.def("getRC", &expMSSA::getRC, R"( Access the detrended reconstructed channel series by internal key diff --git a/src/OutVel.H b/src/OutVel.H index 4fc907543..f2704120d 100644 --- a/src/OutVel.H +++ b/src/OutVel.H @@ -21,7 +21,7 @@ /** Dump velocity flow coefficients at each interval - @param filename of the coefficient file + @param modelname is the file specifying the density model @param name of the component to dump @@ -33,7 +33,7 @@ @param nmax is the maximum radial order - @param scale is the scale parameter for the density field + @param rmapping is the coordinate mapping parameter for the expansion @param rmin is the minimum radius for the density field @@ -50,7 +50,7 @@ class OutVel : public Output private: - std::string filename, model, outfile; + std::string modelname, model, outfile; double prev = -std::numeric_limits::max(); Component *tcomp; CoefClasses::CoefsPtr coefs; diff --git a/src/OutVel.cc b/src/OutVel.cc index 78225f9d2..cd877e4d3 100644 --- a/src/OutVel.cc +++ b/src/OutVel.cc @@ -9,12 +9,12 @@ const std::set OutVel::valid_keys = { - "filename", + "modelname", "nint", "nintsub", "name", "dof", - "scale", + "rmapping", "rmin", "rmax", "ascl", @@ -62,7 +62,10 @@ OutVel::OutVel(const YAML::Node& conf) : Output(conf) // Create the basis // - basis = std::make_shared(conf); + YAML::Node node; + node["parameters"] = conf; + + basis = std::make_shared(node); // Create the coefficient container based on the dimensionality. // Currently, these are spherical and polar grids. @@ -94,7 +97,7 @@ void OutVel::initialize() model = conf["model"].as(); else { std::string message = "OutVel: no model specified. Please specify " - "either 'file' with the model 'filename' or 'expon' for the\n" + "either 'file' with the model 'modelname' or 'expon' for the\n" "exponential disk model (i.e. Laguerre polynomials)"; throw std::runtime_error(message); } @@ -116,7 +119,7 @@ void OutVel::initialize() throw std::runtime_error(message); } - if (conf["filename"]) filename = conf["filename"].as(); + if (conf["modelname"]) modelname = conf["modelname"].as(); } catch (YAML::Exception & error) { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 803800330..9aee83328 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -30,14 +30,14 @@ if(ENABLE_PYEXP) WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) add_test(NAME removeCylCache - COMMAND ${CMAKE_COMMAND} -E remove .eof.cache.run0t test_adddisk_sl.0 + COMMAND ${CMAKE_COMMAND} -E remove .eof.cache.run0t WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Disk) set_tests_properties(removeSphCache PROPERTIES DEPENDS pyexpSphBasisTest REQUIRED_FILES ".slgrid_sph_cache") set_tests_properties(removeCylCache PROPERTIES DEPENDS pyexpCylBasisTest - REQUIRED_FILES ".eof.cache.run0t;test_adddisk_sl.0") + REQUIRED_FILES ".eof.cache.run0t") # Other tests for pyEXP go here ... @@ -56,7 +56,7 @@ if(ENABLE_NBODY) # Makes some spherical ICs using utils/ICs/gensph add_test(NAME makeICTest - COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 1000 -i SLGridSph.model + COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 10000 -i SLGridSph.model WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Runs those ICs using exp @@ -92,14 +92,14 @@ if(ENABLE_NBODY) # perhaps there is a better way? add_test(NAME removeTempFiles COMMAND ${CMAKE_COMMAND} -E remove - config.run0.yml current.processor.rates.run0 new.bods new.bods.0 + config.run0.yml current.processor.rates.run0 new.bods OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid - 'outcoef.dark halo.run0' + outcoef.halo.run0 SLGridSph.cache.run0 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Remove the temporary files set_tests_properties(removeTempFiles PROPERTIES DEPENDS expNbodyCheck2TW - REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;new.bods.0;OUTLOG.run0;run0.levels;SLGridSph.cache.run0;test.grid;" + REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;SLGridSph.cache.run0;test.grid;" ) # Makes some cube ICs using utils/ICs/cubeics @@ -137,7 +137,7 @@ if(ENABLE_NBODY) # Set labels for pyEXP tests set_tests_properties(expExecuteTest PROPERTIES LABELS "quick") set_tests_properties(makeICTest expNbodyTest expNbodyCheck2TW - removeTempFiles expCubeTest removeCubeFiles + removeTempFiles makeCubeICTest expCubeTest removeCubeFiles PROPERTIES LABELS "long") endif() diff --git a/tests/Cube/check.py b/tests/Cube/check.py index 1956e0ee6..b6b0a3eef 100644 --- a/tests/Cube/check.py +++ b/tests/Cube/check.py @@ -1,15 +1,26 @@ -import numpy as np +# open the output log file +file = open("OUTLOG.runS") -# Read the output log file -data = np.loadtxt("OUTLOG.runS", skiprows=6, delimiter="|") +n = 0 # Count lines +mean = [0.0, 0.0, 0.0] # Mean positions -# Columns 4, 5, 6 is mean position -x = np.mean(data[:,3]) -y = np.mean(data[:,4]) -z = np.mean(data[:,5]) +# Open the output log file +# +while (line := file.readline()) != "": + if n >= 6: # Skip the header stuff + v = [float(x) for x in line.split('|')] + mean[0] += v[3] # x pos + mean[1] += v[4] # y pos + mean[2] += v[5] # z pos + n = n + 1 # Count lines -# If the values are close to 0.5, assume it worked -if np.abs(x - 0.5) > 0.15 or np.abs(y - 0.5) > 0.15 or np.abs(z - 0.5) > 0.15: +if n>6: # Sanity check + for i in range(3): + mean[i] = mean[i]/(n-6) - 0.5 + +# If the squared values are close to 0.0, assume it worked +# +if mean[0]*mean[0] > 0.03 or mean[1]*mean[1] > 0.03 or mean[2]*mean[2] > 0.03: exit(1) else: exit(0) diff --git a/tests/Disk/cyl_basis.py b/tests/Disk/cyl_basis.py index 227b22f4a..2d2b6e0d3 100644 --- a/tests/Disk/cyl_basis.py +++ b/tests/Disk/cyl_basis.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # coding: utf-8 -import numpy as np +# import numpy as np import pyEXP # Make the disk basis config diff --git a/tests/Halo/check.py b/tests/Halo/check.py index 4299f403b..d05381c3d 100644 --- a/tests/Halo/check.py +++ b/tests/Halo/check.py @@ -1,14 +1,22 @@ -import numpy as np +# Open the output log file +file = open("OUTLOG.run0") -# Read the output log file -data = np.loadtxt("OUTLOG.run0", skiprows=6, delimiter="|") +n = 0 # Count lines +mean = 0.0 # Accumulate 2T/VC values -# Column 16 is -2T/VC. The mean should be 1 with a std dev < 0.03 -mean = np.mean(data[:,16]) -stdv = np.std (data[:,16]) +# Open the output log file +# +while (line := file.readline()) != "": + if n >= 6: # Skip the header stuff + v = [float(x) for x in line.split('|')] + mean += v[16] # This is the 2T/VC column + n = n + 1 # Count lines -# If the values are within 3 sigma, assume that the simulation worked -if np.abs(mean - 1.0) > 3.0*stdv: +if n>6: mean /= n-6 # Sanity check + +# Check closeness to 1.0 +# +if (mean-1.0)*(mean-1.0) > 0.003: exit(1) else: exit(0) diff --git a/tests/Halo/config.yml b/tests/Halo/config.yml index 3743a65e6..596d51f36 100644 --- a/tests/Halo/config.yml +++ b/tests/Halo/config.yml @@ -8,9 +8,9 @@ # ------------------------------------------------------------------------ Global: nthrds : 1 - dtime : 0.005 + dtime : 0.002 runtag : run0 - nsteps : 200 + nsteps : 500 multistep : 4 dynfracV : 0.01 dynfracA : 0.03 @@ -26,7 +26,7 @@ Global: # Each indented stanza beginning with '-' is a component # ------------------------------------------------------------------------ Components: - - name : dark halo + - name : halo parameters : {nlevel: 1, indexing: true} bodyfile : new.bods force : @@ -40,6 +40,7 @@ Components: rmapping : 0.0667 self_consistent: true modelname: SLGridSph.model + cachename: SLGridSph.cache.run0 # ------------------------------------------------------------------------ # This is a sequence of outputs @@ -48,7 +49,7 @@ Output: - id : outlog parameters : {nint: 10} - id : outcoef - parameters : {nint: 1, name: dark halo} + parameters : {nint: 1, name: halo} # ------------------------------------------------------------------------ # This is a sequence of external forces diff --git a/tests/Halo/readCoefs.py b/tests/Halo/readCoefs.py index ad7ab60e2..0004001f7 100644 --- a/tests/Halo/readCoefs.py +++ b/tests/Halo/readCoefs.py @@ -3,7 +3,7 @@ import pyEXP -coefs = pyEXP.coefs.Coefs.factory('outcoef.dark halo.run0') +coefs = pyEXP.coefs.Coefs.factory('outcoef.halo.run0') data = coefs.getAllCoefs() print(data.shape) print(coefs.getName()) diff --git a/tests/Halo/sph_basis.py b/tests/Halo/sph_basis.py index 1cee5d8ac..e72ca54e5 100644 --- a/tests/Halo/sph_basis.py +++ b/tests/Halo/sph_basis.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # coding: utf-8 -import numpy as np +# import numpy as np import pyEXP # Make the halo basis config diff --git a/utils/ICs/CMakeLists.txt b/utils/ICs/CMakeLists.txt index b0ec76de6..06cbb8b18 100644 --- a/utils/ICs/CMakeLists.txt +++ b/utils/ICs/CMakeLists.txt @@ -5,7 +5,8 @@ set(bin_PROGRAMS shrinkics gensph gendisk gendisk2d gsphere pstmod cubeics zangics slabics) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp exputil - ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) + ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES} + ${HDF5_CXX_LIBRARIES}) if(ENABLE_CUDA) list(APPEND common_LINKLIB CUDA::cudart CUDA::nvToolsExt) diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc index f685133c1..d921c883d 100644 --- a/utils/ICs/gensph.cc +++ b/utils/ICs/gensph.cc @@ -809,6 +809,7 @@ main(int argc, char **argv) if (myid==0) { std::ostringstream sout; sout << "cat " << OUTPS << ".* > " << OUTPS; + sout << "; rm " << OUTPS << ".*"; int ret = system(sout.str().c_str()); } diff --git a/utils/PhaseSpace/CMakeLists.txt b/utils/PhaseSpace/CMakeLists.txt index 3076d0b12..71d1c2141 100644 --- a/utils/PhaseSpace/CMakeLists.txt +++ b/utils/PhaseSpace/CMakeLists.txt @@ -14,7 +14,7 @@ if(HAVE_XDR) endif() set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp expui - exputil ${VTK_LIBRARIES} ${HDF5_LIBRARIES}) + exputil ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_CXX_LIBRARIES}) if(PNG_FOUND) list(APPEND common_LINKLIB PNG::PNG)