CEMeNT-PSAAP / MCDC

MC/DC: Monte Carlo Dynamic Code
https://mcdc.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
20 stars 20 forks source link

added azurv1 convergence tests #205

Closed murrayaidanj closed 1 month ago

murrayaidanj commented 2 months ago

This PR does not change any of the functionality, I am merely adding an example of convergence demonstration, specifically to the benchmarks computed by William Bennett (http://dx.doi.org/10.48550/arXiv.2205.15783).

main.py contains the functions and an example input (at the bottom) that runs simulations for given initial conditions at given particle numbers, and asserts that the convergence is O($1/\sqrt{N}$) as expected (or a slope of -1/2 in log log form) to a given tolerance. main.py can be called with the same --mode=numba functionality as other inputs, but if using mpi (which is great for creating the data), making sure loud=False, as the plotting gets messy.

I figured examples was the best place for this, as it is not fast enough to be a test (especially since the data for each source is best computed separately given the limitations of numba compilation), and doesn't have much of an answer besides -1/2.

ilhamv commented 2 months ago

Thanks, @murrayaidanj !

I just changed the base into dev. @jpmorgan98 , can you please take a look if the base changing worked OK.

Aidan, I think this would make a great verification test. Please check MCDC/test/verification/analytic.

murrayaidanj commented 1 month ago

I've updated the format and placement - it should match the rest of the verification tests now

ilhamv commented 1 month ago

Thanks, @murrayaidanj !

There is an error with the black styling. You can check https://black.readthedocs.io/en/stable/ about it if you want; but I can fix that no problem.

@jpmorgan98 , there is a slight change in the pyproject.toml, I'm wondering why that is.

I'll review the rest of the PR now.

murrayaidanj commented 1 month ago

Of course - thank you so much for all your help!

azurv1_gaussianIC_flux azurv1_gaussiansource_flux azurv1_planeIC_flux azurv1_squareIC_flux azurv1_squaresource_flux

jpmorgan98 commented 1 month ago

I can't commit to this PR so I am attaching an example of the README.md description of AZURV1. Put this file in /tests/verification/analytic. Then just add some notes and equations for whatever makes scene for the tests you are adding README.md

murrayaidanj commented 1 month ago

Thank you for the feedback! I just added comments (and remembered to use black this time). I'm working on the README.md right now. Should I include the benchmark hdf5 that I use? I didn't originally only because it is a large (1.08GB) file.

jpmorgan98 commented 1 month ago

No please don't add anything north of like 50kb. If it's necessary to have something like that 1.08GB file for verification we can look at other solutions to host that

murrayaidanj commented 1 month ago

The README.md attached kept downloading as an empty file, so I created a basic version and added it where specified. To include the equations would entail using about half of an existing paper, so I referenced it where necessary, but I could also put it all in the README if that's better. If there's something in the file you sent that's drastically different, I can go ahead and match it once I'm able to see it.

jpmorgan98 commented 1 month ago

Great when the checks pass again I will merge this! Thanks again!!!!

jpmorgan98 commented 1 month ago

K all we need is to be back in black then we will be all good

ilhamv commented 1 month ago

Thanks @murrayaidanj !

It looks like some of the results do not converge well; they flatten when N is sufficiently high, which indicates that the MC solution leads to a different solution from the reference. Have you made sure that

  1. the reference solutions are not pointwise solutions but mesh average solutions, and
  2. sufficiently fine discretizations are used to specify the initial condition and source.

Another thing is that I found that other parts of the code are changed in this PR. I'm wondering why... Is it because the PR (and the origin branch) were initially made toward (and from) main not dev? What do you think @jpmorgan98 ?

murrayaidanj commented 1 month ago

The reference solutions are indeed pointwise, but I'd run many times more points and averaged them across time - I'm working on averaging them across x as well now. It'll take some time to run those solutions, but as soon as they do I'll let you know if that's the problem.

As for the changed files, I did originally fork from the main branch (I've never used github before ... I'm sorry about that) - I can switch my stuff into a branch off of dev and resubmit a request if that's better practice

ilhamv commented 1 month ago

@murrayaidanj, I see. Note that you can also use the Scipy numerical integration function to do the double integrals. You can check out MCDC/test/verification/analytic/azurv1/reference.py, which is how I do it for the original AZURV1.

As for the changed files, I did originally fork from the main branch (I've never used github before ... I'm sorry about that) - I can switch my stuff into a branch off of dev and resubmit a request if that's better practice

No worries! Let's try reapplying your changes to a new branch off of dev to make the integration go smoothly. Hopefully, it will not be too much of a work, especially since the changes are primarily adding new files/folders rather than modifying existing ones.

You can go ahead and "Close with comment" this PR and then make a new PR using the new branch off dev.

Thanks, Aidan!

murrayaidanj commented 1 month ago

Thank you for the suggestion - I'm trying it out now as I get results with a more refined x mesh. I'll open a new request once I can get the plots to converge better.