robertsj / libdetran

deterministic transport utilities
MIT License
2 stars 2 forks source link

NestedGeoMesh #57

Open robertsj opened 1 year ago

robertsj commented 1 year ago

Imagine the following: given a typical IJK mesh, sweep the cells in the appropriate order for an octant and angle (within octant). Of course, that's what the various SN Sweeper's do (with analogs for MOC, were they implemented. :/ )

Now, imagine each of those "cells" is actually another IJK (or ijk?) mesh who happens to fit exactly within the first. Then, the fluxes leaving IJK's neighbors, say, from the right face of (I-1)JK, the forward face of I(J-1)K, and the top face of IJ(K-1), for a first-octant angle (mu, eta, xi > 0) provide the values for the incident fluxes psi(1jk), psi(i1k), and psi(ij1), for i, j, and k from 1 to their respective limits n_x, n_y, and n_z. As long as there is a sweep defined for the ijk mesh, the outgoing fluxes psi(n_x,jk), psi(i,n_y, k), psi(ij,n_k) can be passed to (I+1)JK, I(J+1)K, and IJ(K+1) or stored in appropriate temporaries.

The mesh (I-1)JK could have its own embedded mesh (i-1)jk. In such cases, the only requirement is a mapping between the fluxes passed from one embedded mesh to another. For the IJK=ijk case (no embedded mesh), psi(ijk) = psi(IJK) for all ijk. For all other cases, some form of averaging or interpolation is needed. For example, if two cells go into four, the two cell values could define a line to interpolate into 4. In the reverse, linear interpolation could be used, i.e., each of the two values would be weighted average of the closest pair of the four. Alternatively, linear regression would fit the four values to a line on which the 2 would be evaluated. Of course, these operations are not unrelated: going from many to one is o=(1/N) * [1 ... 1]' m, and from one to many is m=[1 ... 1] o. Call these operators M and O (for o = Mm). Then M = (1/N) O'. Or, one could be symmetric and say M=O', where O'=sqrt(1/N)[1...1]. More generically M goes from a some number o to m, i.e., M.shape = (m, o). Then, O.shape = (o, m) and O=M.T. The symmetric case is nice in theory, but internally, we might need the sweep to use "absolute" units so that various scaling is not required.

outgoing angular fluxes from IJK's coarse neighbors can be used to initialize the many larger number of boundary fluxes in ijk