[ChasingNeutrons]

Hi everyone,

I have developed a new solver on OpenFOAM - it's used for neutron transport and, at least on the surface, isn't equivalent to finite volume methods.

A bit of a preface on how this solver works: it ray traces across a geometry from a number of different starting points and at a number of different angles. It then iterates on boundary flux values streaming along these rays to calculate cell values of the field I'm interested in, given a source. This field is then used to calculate an eigenvalue and a new source and the process is repeated until the eigenvalue and the field have converged. It's a method of characteristics neutron transport solver (or MOC).

My problem is that I would very much like to recast this program in a Foam-y way where I can put the field I'm interested in, the ray tracing information, the eigenvalue and other source terms in a form like:

solve
(
MOC::dds(boundaryFlux,rays) + (sigmaT, boundaryFlux) 
== MOC::Source(eigenvalue, chiNuSigmaEff, flux) + MOC::Source(SigmaS,flux)
)

My problem is I'm not entirely sure where to begin - where can I create my MOC::dds() for example? Where can I modify the solve function? Any help would be very much appreciated and I would happily share my code if it would help explain the situation.

Cheers.

Edit: To clarify, the above isn't necessarily how I want my solver to be written but somewhere along those lines would be ideal, subject to whatever constraints may exist.