Begin DISORT integration
This PR will integrate the new scattering data system into DISORT.
Hi @riclarsson and @m-brath,
I put together a notebook to perform scattering simulations with the setup used by scatsolvercomp in v2.6
It currently fails because the propagation_matrix_scattering isn't setup properly in the agenda. Could you have a look at this and tell me what I am doing wrong?
From what I could gather it looks like we should be close to getting this to run. Great work on putting all the structure in place!
Made some progress. Notable code changes include:
- Fixed the order of the method arguments in
disort_settingsInit - ~~Removed test for polarization of propagation matrix because TRO phase matrices have polarization effects.~~ My bad, I got the format of propagation matrix wrong.
I fixed some more things in the calculation of the bulk scattering properties and checked that extinction, single-scattering albedo, and legendre coefficients agree with the ARTS 2.6. However, the Tbs for the cloudy sky simulations are way off.
@simonpf I rebased this on latest so I could use some of the changes. Please update on your side before you do further changes.
@olemke Can you help me here? I somehow pushed to my own branch isntead of here and when I tried to resolve it I think a lot of stuff got messed up.
My intent is that this PR should look like https://github.com/riclarsson/arts/tree/scattering_data, and it is not reporting any conflicts.
@riclarsson Since your branch is up to date with main and this one isn't, I assume that you did a rebase from main onto your branch and this causes merge issues with this one. I did a hard reset of this branch to the state of yours, essentially replacing this branch with your version. @simonpf Make sure you get a clean checkout of this branch before doing any further changes. If you have already done changes in your local branch, keep it around, check out a fresh copy of this branch under a different name, then cherry-pick the new commits you made in your original local branch into this one.
The disort.sht.ipynb test fails in the debug builds with the following assertion. I put it here since the assertion output is unfortunately not visible in the logs:
Assertion failed: (sht->get_n_zenith_angles() == n_za_scat_), function to_spectral, file phase_matrix.h, line 698.
The values are:
sht->get_n_zenith_angles(): 80
n_za_scat: 451
@olemke Thank you! I have three more patches to try and deal with the documentation failing. It builds for me now at least.
@simonpf One key thing here is that I also was forced to experience the slowness that was GCC on Linux. I fixed it. ARTS3 now compiles faster than ARTS2. At least on low GHz systems. (My laptop went from 2h to 20m compilation.) The main thing to avoid in the future is too many std::visit in header files. Or at all if the std::variant is larger than some number of types. Depending on how much we grow the types of scattering species we might have to create some layer inbetween to hide the variant. I moved a few parts of your code to cc files to avoid this also. The main slowdown before was not here though.
This updates the branch to current ARTS. New feature that is important: variants have their XML code automatically generated if they have unique names and all their components are XML-compatible. The XML code is now local to the current type, so I added a dummy XML IO for the scattering habit class. Note that it is also possible to save the functional classes in this new XML scheme, though it is experimental as intel Mac has problems remembering giving out names.
So I found the solution to why the cpp-disort code was misbehaving. A missing scale factor (1 - omega) in the source polynomial. With that added, the cpp-disort solver is now as close to cdisort when it is cloudy as when it is not.
So a few things had to be changed to make this pass all local tests (check and pyarts_docs targets).
@simonpf , I think it is mostly working now. There might be an edge on Linux doc/reldeb that requires some more patching, but I am sure those will be relatively easy to fix. I had to make a lot of changes to make this work, please have a look that we might merge this when all is OK:
- Fix the bug in my port of the DISORT code. I mean, it is required.
- I had to make the DISORT example an actual python file. This was because the data files are not found once the notebook is copied into the sphinx chain. @olemke , can you have a look and see if it is possible to add the test paths easily here? Otherwise we are effectively blocking notebooks that rely on data.
- I had to remove one of the assertions of the scattering species code. The code is
assert(sht->get_n_zenith_angles() == n_za_scat_);. The reason I had to remove it was because it fails in its related test, but the test still works on release mode. So I believe this check is either wrong, or has some logic in the code leading up to it not well-developed.
In addition to this I have some design questions for you @simonpf , some of which I changed already to get rid of warnings, which I think we must do because it is not nice to not be able to find real errors:
- I added default constructors and move operations to many classes. I think you are compiling without warnings on because my compiler was screaming at me that several default copy constructors were actually deleted. If you define
Foo() {}instead ofFoo() = default;, you cannot use any default stuff anymore. Is this OK or are some of these types not meant to be copied? - Do we really need all this
template <typename Scalar>stuff? Because it would be really nice if more of the scattering library compiled just once instead (and if small fixes here and there didn't trigger recompilation of the whole project) - You added a
get_bulk_scattering_properties_tro_spectralthat took a tolerance instead of a Legendre degree. This is not good in my opinion. The signature to call any of theget_bulk_scattering_properties_tro_spectralmust be the same (bar theconst-ness of theIndex). They must also all produce the same size output when they are called as well, based on the input. I do not see this being possible to maintain otherwise. If we need to compute the spectral matrices based on a tolerance level instead of a Legendre degree, I think we simply have to add another method for it. I changed your method to take anIndexfor now. I noticed that it doesn't actually care about the type or number being sent to it, so I think it will not respect sizes. This is problematic. Is there no way to just select the degrees it can, and leave things empty based on the value of the index passed to it?
This last commit is just a performance review. Having reviewed cdisort, I decided to just port their real-value eigenvector diagonalizer. So if you compile ARTS in GPL mode (i.e., the default) you now get a twice as fast diagonalize method as before.
@olemke , @m-brath , and @simonpf , I will go ahead and merge this now (after fixing whatever might fail the tests). There are too many changes now. And I think it fits the title "Begin DISORT integration".
More tests are welcome, and I will open an issue for Oliver to remember to have a look at turning the disort py-file back into a notebook when we figure out how to pass paths to it nicely.
Go ahead. I'll have a look when I'm back from vacation.