mctools / ncrystal

NCrystal : a library for thermal neutron transport in crystals and other materials
https://mctools.github.io/ncrystal/
Other
38 stars 17 forks source link

Manual alpha, beta grid points with loadKernel(vdoslux=5) command #85

Open chapmancw opened 2 years ago

chapmancw commented 2 years ago

Is your feature request related to a problem? Please describe. When generating an S(alpha,beta) object using the loadKernel(vdoslux=5) command, an automated alpha,beta grid is generated. For the application I am interested in, this grid has an insufficient density of alpha,beta points.

Describe the solution you'd like I would like to be able to manually define the alpha,beta grid points to calculate the S(alpha,beta) object on. It could either be added onto the pre-defined grid or replace it all together. If possible, it would be preferable if this could be made accessible in the Python API.

Describe alternatives you've considered I am unfortunately not familiar enough with NCrystal to know of any other potential work arounds for this problem. We have previously used the LEAPR module of NJOY to generate an S(alpha,beta) object on a given alpha,beta grid, but NCrystal has additional functionalities not currently implemented in NJOY that we would like to use.

tkittel commented 2 years ago

Thanks @chapmancw for the request, I think I can cook up something. I'll ask some more pointed follow-up questions later, but for now a few initial thoughts:

  1. I am pretty surprised that vdoslux=5 does not give you enough points, at all the examples I investigated vdoslux=5 was completely overkill.
  2. Looking at the code, I think beta-grids used for VDOS->S(alpha,beta) expansion will for technical reasons always need to have more points at beta<0 and for the positive and negative grids to be otherwise symmetric around beta=0. I.e.: [ bmin, ..., -bmax, ..., 0, .... bmax] where the part [-bmax, ..., 0, .... bmax] is completely symmetrical around beta=0. However one can of course have a beta grid request which does not strictly adhere to this and then just have NCrystal internally add any missing points. Let me know what you think.
  3. I assume you are not just satisfied with specifying a beta-grid or an alpha-grid, but might also like to control how many orders of VDOS expansion to do -- albeit indirectly by specifying Emax. Emax is typically 5eV for ENDF kernels, but if you are fitting data in a low-E region of phase-space and need to constantly redo the VDOS expansion, I could imagine you would be happy to be able to lower this number. Thus, I am imagining that it might be better to have the VDOS->SAB functionality available as a completely standalone C++ and Python function, which could then take whatever parameters (vdos, atomdata, betagrid, Emax, ...). Crucially the result of such a standalone function would for simplicity not participate in NCrystal's internal caching scheme. Of course, such a function should be able to work with VDOS loaded from .ncmat files as well as specified manually.
  4. Such a function would clearly be experts-only, with all the caveats ("here be dragons") that would imply :-)
  5. If the intent is also to use the created kernels for simulations with NCrystal, it is probably a good idea to somehow merge your desired grid points with those suggested by NCrystal (for any reasonable vdoslux), since some of the values are added with the intent of minimising artifacts introduced by the simulation code (in the ideal world we should not have such artifacts, but I haven't had time to make the world that ideal yet ;-) ).
tkittel commented 2 years ago

I am adding @marquezj and @dddijulio as well, since in the past they also had a similar request. Guys, can you also review the request from @chapmancw discussed here, and provide any feedback. Then I can hopefully figure out what is needed to make everyone happy all in one go :-)

chapmancw commented 2 years ago
  1. The energy range I'm specifically interested in is between 1-1000 meV, and when I set vdoslux=5, I only get 445/3200 of the total beta grid in that energy range, with 386 of those points being between 100-1000 meV, leaving only 59 points from 1-100 meV, where I previously used ~440 points.
  2. Without knowing much about the exact inner workings of NCrystal, would it be possible to just add the requested points to the NCrystal generated alpha,beta grids? I could just run this with a smaller vdoslux (e.g., vdoslux=4) to ensure the alpha,beta grids to blow up beyond any memory/size/etc. limits within NCrystal.
  3. I actually don't know if that would be necessary; NCrystal currently generates the S(alpha,beta) fast enough for our application & unless modifying the VDOS expansion order would significantly impact the energy region below 1000 meV, I would think that would only potentially hinder us (though I might be mistaken there).
  4. Agreed. This is a fairly niche edge case, so I would think this would be better as an optional flag/variable/etc. instead of a required one.
  5. I think that's more or less what I mentioned in 2. above.
marquezj commented 2 years ago

Hi @chapmancw: are there features of the S(Q,w) map not being represented by the grid? Otherwise, you can just interpolate the output of NCrystal to whatever grid you need.

tkittel commented 2 years ago

@chapmancw , just a quick note concerning your reply to item 1: vdoslus=5 implies Emax=12eV by default, which certainly "wastes" points in the 1eV..12eV region. But you can already now easily override this to e.g. Emax=1eV, and thus get all of your grid points in this lower energy region. Just add the line egrid 1.0 inside the relevant @DYNINFO section of your NCMAT file (see https://github.com/mctools/ncrystal/wiki/NCMAT-format). That will btw. automatically reduce the number of orders needed so it should also speed up initialisation. And don't worry, only orders not contributing to scattering of a 1eV neutron will be left out (NCrystal knows what part of the kernel is populated by which expansion order, so will decide the required maximum order exactly based on Emax).

ramic-k commented 1 year ago

Hi Thomas, tanks for implementing this in v3.1.0 release. Where can we find some instructions on how to actually use the new feature? ;)

Thanks for implementing this!

tkittel commented 1 year ago

Hi Kemal,

What?? :smiley: I didn't implement anything, as I was waiting to hear your response to the comments from me and @marquezj concerning possible workarounds (also, I had other work on my plate).

Cheers, Thomas

ramic-k commented 1 year ago

:D Did i read this wrong?

"Add internal utilities for evaluating S(alpha,beta) kernels at any (alpha,beta) point, and add utilities for what is essentially exact sampling or integration for cross sections. These utilities are useful as they allow us to validate scattering kernel models better, and are used directly in the new UCN-production models."

tkittel commented 1 year ago

I think you did partially, but it is certainly related. That note talks about a C++-only helper class, which is not exposed to python. And it only deals with interpolating in an existing kernel, not adding points before constructing the kernel.

But of course, it all depends what exactly you guys need :-)