This will require Embree AND CGAL to be installed.
Conditional import or try/except would solve this.

https://github.com/sampotter/python-flux/blob/700965ad79a8f387869459da433f63ad2629fe32/src/flux/shape.py#L1
https://github.com/sampotter/python-flux/blob/700965ad79a8f387869459da433f63ad2629fe32/src/flux/shape.py#L8

The issue seems to be with this row. If this is changed to indptr = [0], then the code doesn't crash, but it runs in what seems to be an infinite loop.

Just not to lose track of this potential issue @nschorgh mentioned here #21 (comment) (easy enough to solve, but since we are in the middle of several PR already...)

We should pick an open source license for this package. I have no preference as to which license we choose. It should be a license which makes it easy and straightforward for this library to be used in the communities in which it is most likely to gain traction.

There's also this: https://code.nasa.gov/#/guide

Thoughts?

@emazarico @nschorgh @steo85it

While making experiments with your libraries at the Moon's south pole, I noticed a (possible) issue when using lsp_spice.
This should be the interior of a crater (see elevation map of the region below), with a smaller crater within it even projecting a shadow of its own (which is clearly not the expected behavior...).

I think that this could be due to the fact of being in the southern hemisphere and the transformation to cartesian, which the code ending up "confusing" a crater for a hill.

(above is the elevation from the Moon mean radius, while below is Z ~ signed distance from the Moon center)

You can check it by setting

i0, i1 = 2425, 2663  # 600, 700
j0, j1 = 1912, 2172  # 650, 725

as indexes to select a region from the whole south_pole grd "LDEM_80S_150M_adjusted.grd" in lsp_make_shape_model and running the lsp_shape.py example for the 6-Sep-2021.
Since the region is a bit large, you can maybe subsample the map a bit by adding

IJ = IJ[::3]
V = np.array([x.ravel(), y.ravel(), z.ravel()]).T
V = V[::3]

to lsp_make_shape_model.py.

I tested several solutions/options:

• checked what happens if lat0 is set to 90 (North Pole) -> the illumination changes a bit (seasons?), but the craters persists behaving as hills
• using a (stereo-) projected map (instead of converting it to a cartesian point cloud) works in preserving the "crater behavior", but I don't know then how to get the actual illumination at the South Pole at the Moon (instead of at its center). At least in some tests, seasons seem to be opposite than expected (Sep 21 should be winter and Mar 21 Summer at the Moon South Pole, but I get a Sun with a higher elevation at these craters).
• following this test, I thought that in the end what counts is the sun direction at a given location, so that sun_dirs in lsp_spice.py could be tweaked to include the cartesian position of the crater at the South Pole. This doesn't seem to impact the result (the angular difference given by the radius of the Moon is too small? Or is there something I'm missing in how that sun_dir is then used by get_direct_irradiance?
• I also thought that the N(ormal) vector should be responsible of telling the code where the "outside" of the planet is, but that looks fine (~ 0,0,-1) in the non-working tests. Still, in some tests I noticed that introducing/removing the test

N[(N*P).sum(1) < 0] *= -1

can reproduce the same behavior of "hill vs crater".

Any hint here would be appreciated! :)

in ingersoll.py line 104, it should be if args.blocks and args.tol:, otherwise FF._root is not defined.

in ingersoll.py line 121, it should be '- wrote E.png'

@steo85it @nschorgh

I think we should try to improve the interface to the time-dependent thermal models a bit.

We have a few different ideas at this point for how to do the time-dependent simulation (basically Norbert's "just do one MVP" approach vs solving the full radiosity system twice at each step). I think it's worth doing a direct comparison of these. But I'm getting lost in the details of how to set them up.

I would like to be able to write something like:

from flux.model import ThermalModel

thermal_model = ThermalModel(
    FF,
    times, # len(times) == N, times[i] = ith time point
    sun_positions, # sun_positions[0] = ith sun position corresponding to times[i],
    method='norbert', # should be == 'norbert' or 'sam', just for example :-)
    **physical_parameter_kwargs, # passed on to PccThermalModel1D
)

The variable thermal_model should then be a Python iterator or generator. That is, we should be able to do things like:

for T in thermal_model:
    # T is an nz x shape_model.num_faces ndarray of temperatures for each skin depth
    ... do things with T ...

I'm going to put together a rough draft of this hopefully today or tomorrow and then it would be helpful to get you guys' feedback.

Picking up from the conversation in issue #6.

Just as a reminder, the README needs to be updated to allow CGAL compile (and some parts clarified for embree, too).
(I can also do it, not assigning to anyone specific)

https://github.com/sampotter/python-flux/blob/8cb4e46368b1278975746899401872b72a29bf17/src/flux/model.py#L8

At some point we decided to rename rho to albedo here, because in the paper rho is left of the sum.

Right now the strategy for assembling the form factor matrix involves building the complete sparse form factor matrix for the coarsest blocks in the spatial partition. This takes O(N^2) memory, so assembly takes O(N^2) memory.

To build really large compressed form factor matrices, we need a form factor MVP which only uses O(N) memory. This is easy enough (just compute each MVP directly using a pair of for loops, redoing the raytracing and form factor computations for each pair of elements), but the goal will be to do this as fast as humanly possible so that assembly doesn't run 100 times slower than what we're doing now.

In the radiosity equation rho is outside of the sum, but the albedo is inside the sum, so our solver is only correct for spatially uniform albedo.

We will either have to solve for (1/A - F) B=E or store two form factor matrices F and A F. Some asteroids are quite dark, so 1/A might grow large, and the second approach might be better, but the second approach can't handle time-dependent albedos.

This is low priority, because we can do a lot of time-dependent calculations without radiosity solvers.

I think that while updating shape.py with the is_occluded function, the option for an extended light source (D.ndim == 2, with D the direction between the vector array to the "Sun") got broken.

The current _is_occluded doesn't work in either cgal or embree version. For the cgal version, the issue is visible at

https://github.com/sampotter/python-flux/blob/700965ad79a8f387869459da433f63ad2629fe32/src/flux/shape.py#L263-L267

where D[i] cycles on the number of shape-model's facets instead of the number of "Sun facets/vectors" (in my case, D.shape=(100,3)).

For reference, see the old function get_direct_irradiance in the attached file (for a working embree version).
shape.txt

Any help is welcome @sampotter

Please consider this addition to shape.py, which allows each pixel to have a different incidence angle.
This will interpret an array of sun directions of the same length as the number of topographic pixels as a lat/lon dependence.

       # Compute the direct irradiance
        if Dsun.ndim == 1:
            E = np.zeros(n, dtype=self.dtype)
            E[I] = F0*np.maximum(0, self.N[I]@Dsun)
        elif n==m: #### new option
            E = np.zeros(n, dtype=self.dtype)
            for i in range(n):
                if I[i]:
                    E[i] = F0*np.maximum(0, self.N[i,:]@Dsun[i,:])
        else:
            E = np.zeros((n, m), dtype=self.dtype)
            I = I.reshape(m, n)
            for i, d in enumerate(Dsun):
                E[I[i], i] = F0*np.maximum(0, self.N[I[i]]@d)
        return E

And the value of m will have to be assigned outside of the previous if construct.

I'm pretty sure a good solver for the radiosity system will be a "block Jacobi iteration". Going to try to sketch out what this in the issue before trying to implement it.

For a system Ax = b, the Jacobi iteration separates A = D + L + U, D = diag(A), splitting A into D and L + U terms, then relaxing. This works very reliably for this system (seems to converge to machine precision in at most 8 iterations). It's also pretty simple since D = rho*I for the radiosity system.

A "block Jacobi iteration" would separate A = D + L + U where now D = diag(D1, ..., Dn), where the Dis are diagonal blocks at a certain depth. The idea here would be to relax as before, i.e. x <- D\(b - (L + U)x), except now D\ can be done by recursively applying the block Jacobi iteration, since each Di is a smaller radiosity system. At a fine enough level we can just use a sparse direct solver and avoid iterating. This leads to something reminiscent of a multigrid V-cycle (maybe there's another name for this in this context...). It seems likely that we would only need one or two of these "V-cycles" to converge fully.

We should be able to extend this library to handle arbitrary BRDFs (although it makes sense to focus on important special cases). This is a thread to collect ideas related to this.

Here's a quick to-do list:

• Pick a shape model and a small region of interest (~1K-10K triangles---note that 10K tris => ~1T triples of points!).
• Pick a suitable BRDF.
• Try visualizing the irradiance & compare w/ Lambertian to see the difference.

Try setting vmin and vmax by passing them as kwargs to flux.plot.tripcolor_vector. Doesn't work.

One difference leading to differences in #19 is the initialization method: in the F90 code, both Qprev and Qnp1 are given as input, while here only Qnp1.
To initialize Qprev, I run model.step(0,Qprev), but this changes Tsurf at t_0. A more consistent way would be to allow for directly setting model.Qprev = X (which now gives an error).
I'm not sure if this issue really causes discrepancies in #19, but it's a difference...

@sampotter At some point, it might be useful to get a review and publish this code through the Journal of Open Source Software. The review process is quite interesting, well structured, and useful to improve code readability, reliability, get suggestions of best practices (for testing, benchmarking, examples, docu, etc...), ...
Just pinning an idea for the future.

@steo85it

What has been done with this so far? I would like to push this along

@nschorgh @steo85it

Callling shape_model.get_visibility(I,J) with I and J different than the full mesh resulted in getting the wrong orientations of the visible faces and thus in a wrong evaluation of the visibility matrix. Moreover, the resulting orientations depended on the arrays I and J.

I think debug-compressed-form-factor-matrix branch has an issue similar to #4 .
When running haworth.py from that branch, I get the left result in the figure below (E, above; T, below) - on the other hand, when I run the same experiment from a9f3dc1 , I get the result on the right. The latter makes more sense, considering the shape of the crater and shadowed areas therein.

I checked and the issue which solved #4 still looks "solved": now that I know that the problem doesn't come from my own code/meshes/changes, I will try to track it down checking differences among the commits.
Ah, this happens with both embree or CGAL: it has to be related to some of the clean-ups, but not to this "fundamental change". Also, I went back in history in the master branch only to avoid having to debug unrelated errors in the haworth.py test, not because the issue showed up (I'll check now at which point it appears in master, if any).

@sampotter

Hi Sam,
while testing my "outer mesh" stuff for ray tracing (as discussed last week), I found some weird artifacts in the T results.

The "2x2 grid" appears both in my tests and using your lsp_spice "base example" (before my mods, too), both with quad and octree. Did you ever notice this? It's also visible on T (obviously) but not on the E plot.
This particular output results from selecting
i0, i1 = 929, 1831
j0, j1 = 500, 1435
from lunar_south_pole_80mpp_curvature.grd and

utc0 = '2021 FEB 15 00:00:00.00'
utc1 = '2021 MAR 15 00:00:00.00'

et0 = spice.str2et(utc0)
et1 = spice.str2et(utc1)
stepet = 3*24.*3600
nbet = int(np.ceil((et1 - et0) / stepet))
et = np.linspace(et0, et1, nbet, endpoint=False)

I have several examples of this happening. Also, I checked if anything would change when removing the if not self._blocks[i, j].is_empty_leaf from compressed_form_factors.py (just in case, as it is a recent change), but it doesn't seem to be related to it.

Still not 100% sure that it's not related to something I did, but your base example seems to present the same kind of issue...

Right now we have duplicated versions of C files from @nschorgh's PlanetaryCodeCollection in this repository, but eventually the ideal setup will be one of the following:

Option 1:

1. made PCC a submodule of this repository
2. clone PCC & compile
3. tell setup.py the location of the compiled static or shared libraries

Option 2:

1. clone and compile PCC
2. install static or shared libraries + headers system wide
3. tell setup.py to look for these libraries and headers

Then we won't have any duplicated code between the libraries and we can pull in changes that Norbert makes directly in the PCC repo.

It would be helpful to add some instructions for how to compile and install this library using spack. @steo85it pinging you since you just had some luck with this. I'm going to use spack on the clusters at work once I start running some larger tests for the paper (most likely this week).

@steo85it

IIRC, tqdm is some internal utility library of yours. There's some stuff that's using tqdm that's checked in now. Can we fix this dependency? If you'd like to just copy the relevant parts of tqdm into python-flux so that all of the examples run, that's fine. We can clean it up later if necessary.

Somehow shape_model.N in the program and shape_model.N in FF.bin are not the same.
I haven't pinpointed the reason yet, but the normals in FF.bin are not normalized.
The other variables, V, F, P, and A, are identical.

FF2 = CompressedFormFactorMatrix.from_file('FF.bin')   
shape_model2 = FF2.shape_model
print(min(shape_model.N[:,2] - shape_model2.N[:,2]))

Title says it all. Pickling causes problems when the API changes.

@steo85it

When you have a chance, please make a PR to include your GUI. It's OK if it's in a rough draft state... I have some tests I'm doing for issue #6 and it would be very helpful. ;-)

(No rush! It's the weekend...)

L201 and L205 take care of the same basemesh is None, so one of them is presumably redundant.
https://github.com/sampotter/python-flux/blob/700965ad79a8f387869459da433f63ad2629fe32/src/flux/shape.py#L201