Test application of conductionQ to Moon terrain (extension of lsp_spice.py)

[steo85it]

Now that the wrapper is tested, it would be great to understand how to integrate it in our applications! :)

To get an initial sandbox, I extended lsp_spice.py with a (very preliminary) suggestion according to my basic understanding: based on the computed E at each time step, we first compute T with k*dT/dz = Fsurf = 0 (I added it as argument to flux/model.py) to get a preliminary Qn. We use Qn to set up the model, then enter a loop over E(t), compute Qnp1 and extrapolate Fsurf(t+1) to be used in the computation of T(t+1), and so on... we get "corrected" T that we can save and plot as usual.

It's not performing at all and I'm not even sure if it's used correctly (probably not), so please have a look and correct me where I got things wrong! I'll then improve it according to your suggestions, if any. :)

p.s.: accepting PR #21 should "clean up" this one, which in principle only contains the last 2 commits.

[nschorgh]

The sum of absorbed radiative fluxes is Q_absorbed = (1-albedo)*(Qdirect+Qrefl)+emiss*QIR. where Qrefl and QIR are obtained with two separate calls to solve_radiosity. Note that compute_steady_state_temp can't be used because it's the equilibrium temperature (no Fsurf). So, after two calls of solve_radiosity, the absorbed flux for the next time step is Qnp1 = (1-rho)*(E + Brefl) + emiss*BIR) Then, one call to conductionQ provides Fsurf. And the temperature is obtained with Qnp1 + Fsurf = eps*sigma*T^4

[nschorgh]

Alternatively, one doesn't have to use solve_radiosity at all and use eqs. (6), (7), and (1) in our paper draft, once at every time step.

If using eqs. (8) and (9) instead, it's should be something like: ` def compute_surface_temp(FF, E, rho, emiss, Fsurf, clamp=True):

# short-wavelength irradiance
B = solve_radiosity(FF, E, rho)[0]

# thermal irradiance
IR = FF@((1 - rho)*B + Fsurf)
Q = solve_radiosity(FF, IR)[0]

# surface energy balance
Qabs = (1 - rho)*B + emiss*Q + Fsurf
if clamp:
    Qabs = np.maximum(0, Qabs)

T = ( Qabs / (emiss*sigma) )**0.25

return T`

The first option does one order of reflection at each time step. The second option does infinitely many reflections at each time step.

[steo85it]

Hi Norbert, thanks a lot for the clarification: I had a lot of confusion up there and needed precisely your answers to straighten things up!

I was indeed following the last option you presented (and eqs. 1, 6-9 of the paper draft) when setting up my modifications to model.py: despite the name, when introducing a parameter Fsurf != 0 (which is not appropriate, I agree) and then replacing tot with Qabs in the case E.ndim == 1, the algorithm should be the same as the one you sketched (just checking).

Ok, one should probably rename the function to compute_surface_temp, and/or make a more specific compute_steady_state one which would be the same but with Fsurf=0.

I'll update things accordingly tomorrow (I'm still not 100% sure about the initialization step, but it seems ok), then we can think about "computational efficiency" (if I understand correctly, one needs to loop over this function for quite a long time, and preferably lots of time steps, to reach the steady state...).

[steo85it]

So, would this be correct (assuming I'm using the compute_surface_temp function you sketched above)? My issue is that I need 2 values of Q (at t and t+1) to call conductionQ and compute the first Fsurf. Either I use Qprev=0, or I call compute_surface_temp twice (with E(t) and E(t+1)), both with Fsurf=0, before getting the first Fsurf, correcting T, etc...

# compute provisional T and Q from IR only (Fsurf=kdT/dz=0)
T = compute_surface_temp(FF, E[:, 0], rho, emiss, Fsurf=0.) # this contains the 2 radiosity calls, I still don't have any Fsurf
Qn = sigma * emiss * (T ** 4.) # it would be nice to get Q from compute_surface_temp, w/o reverting T every time

# set up model at t=0 (with Qprev = Q from provisional model) ==> this is what I'm unsure about, what value to use as initial Q
model = PccThermalModel1D(nfaces=E.shape[0], z=z, T0=210, ti=120., rhoc=960000.,# we have rho, but what value for c?
                              emissivity=emiss, Fgeotherm=0., Qprev=Qn.astype(np.double), bcond='Q')

# loop over time-steps/sun angles
for i in range(len(sun_dirs) - 1): # len - 1 since we always take E(t_{i+1})=E(np1)
    # get Q(np1), including contribution Fsurf=k dT/dz (except at first iter); convert to Qnp1
    T = compute_surface_temp(FF, E[:, i + 1], rho, emiss, Fsurf=model.Fsurf) => here we 'advance time', using Fsurf from the previous time step
    Qnp1 = sigma * emiss * (Tloc ** 4.)
    # extrapolate model at t_{i+1}
    model.step(stepet, Qnp1) => updates the value of model.Fsurf to be used in the next `compute_surface_temp`, sets model.Qprev = Qnp1

and basically stay in the loop until T(t) - T(t-1) < eps (meaning that I'd need to have stored enough sun_dirs, but that's another issue)?

[nschorgh]

Qn=Qnp1 is a better initialization than Qn=0. I use Qn=(1-rho)*E(n) and Qnp1=(1-rho)*E(n+1) for initialization. If one uses solve_radiosity then one needs an initial value for Fsurf, which can be zero.

Without solve_radiosity, no initialization of Fsurf is required. A Fortran version of this second approach is here, but I'm not sure that makes it any clearer: https://github.com/nschorgh/Planetary-Code-Collection/blob/master/Topo3D/cratersQ_moon.f90

[nschorgh]

I spelled out the equations in the Overleaf manuscript, (13)-(19). With this approach Fsurf isn't needed at all, because conductionQ already includes the surface energy balance equation and outputs Tsurf^(n+1).

[steo85it]

Nice, thanks for the suggestions/clarifications! Now it seems more "linear" (that loop will be quite slow in python though, and no way to vectorize on the temporal axis if each step depends on the previous one... mmm...)

# compute provisional Tsurf and Q from direct illumination only
Qabs0 = (1-rho)*E[:,0] + Fgeoth
Tsurf0 = (Qabs0/(sigma * emiss))**0.25
# set up model (with Qprev = Qabs0)
model = PccThermalModel1D(nfaces=E.shape[0], z=z, T0=210, ti=120., rhoc=960000.,# we have rho, but how do we get c?
                              emissivity=emiss, Fgeotherm=Fgeoth, Qprev=Qabs0.astype(np.double), bcond='Q')
# loop over time-steps/sun angles
for i in range(len(sun_dirs) - 1):
    # get Q(np1) as in Eq. 17-19
    if i == 0:
        Qnp1 = compute_absorbed_flux(FF, E[:, i + 1], rho, emiss, Tsurf=Tsurf0)
    else:
        Qnp1 = compute_absorbed_flux(FF, E[:, i + 1], rho, emiss, Tsurf=model.T[:,0])
    # extrapolate model at t_{i+1}, get model.T(i+1)
    model.step(stepet, Qnp1)
    # add check for steady-state (model.T(i) - model.T(i-1) < eps) to break loop

[nschorgh]

That loop might not be slow, because in Python 3 range is the same as xrange.

[sampotter]

I'm going to hold off on accepting this until accepting #21 first.

However, one quick thing: let's avoid making too many new files. Ideally, fewer scripts in the examples directory is better. So instead of adding lsp_spice_condQ.py, it would be better to modify lsp_spice.py to use PccThermalModel1d for a thermal model. Then it will be an end-to-end example and more interesting for someone trying to figure out how to use the library.

[steo85it]

I perfectly agree + the only difference with lsp_spice.py should be the call to conductionQ (which can be depending on some parameter in the example). I set it up in this way as this was more for my understanding at first, and I didn't want to "contaminate" other functions. I can amend that.

[steo85it]

I removed lsp_spice_condQ.py and updated lsp_spice.py to include the option to consider subsurface heat conduction. The commit now also includes the update_incoming_radiances function in flux/model.py to compute Eq. 17-18 of the radiosity manuscript, as discussed. It runs smoothly, but the output plots are still a bit useless right now (the surface temperature during spin-up of the model).

[steo85it]

superseeded by PR #30