848 abiotic model add latent heat flux to derivative of temperature update

[vgro]

The PR started out as adding a missing term to the derivative of the energy balance to fix #848 , I also made some small changes to the structure to make it more readable and updated the documentation.

The docs are built here

Fixes #848 #781

Type of change

• New feature (non-breaking change which adds functionality)
• Optimization (back-end change that speeds up the code)
• Bug fix (non-breaking change which fixes an issue)

Key checklist

• Make sure you've run the pre-commit checks: $ pre-commit run -a
• All tests pass: $ poetry run pytest

Further checks

• Code is commented, particularly in hard-to-understand areas
• Tests added that prove fix is effective or that feature works
• Relevant documentation reviewed and updated

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 93.18182% with 6 lines in your changes missing coverage. Please review.

Project coverage is 93.90%. Comparing base (f47818d) to head (30ea662).

Files with missing lines	Patch %	Lines
virtual_ecosystem/models/abiotic/energy_balance.py	92.68%	3 Missing ⚠️
virtual_ecosystem/models/abiotic/microclimate.py	92.85%	3 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #897      +/-   ##
===========================================
- Coverage    93.96%   93.90%   -0.07%     
===========================================
  Files           77       77              
  Lines         5864     5902      +38     
===========================================
+ Hits          5510     5542      +32     
- Misses         354      360       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jacobcook1995]

Might be worth giving either a citation to or writing out the temperature function that is being differentiated here

[vgro]

Will need to look into this again, I noted this down in a separate issue #919

[jacobcook1995]

I had a bunch of comments which I'm adding so we can discuss them when we meet

[vgro]

Thanks @jacobcook1995 for your feedback! I made a few more changes:

• updated documentation to have the same variable names, including the suggested change from EB to dQ/dt. I hope I caugt everything, will have a separate go at this when the insides of the model are ok.
• simplified calculation of net radiation
• incorporated scipy.optimise.newton to solve canopy temperature. It is not converging so I currently write out the best estimate after maxiter. In the course of this deleted my own derivative function.
• changed how air temperature is updated

Variable values are not in the range of physically possible but there needs to be some calibration done. There are also still open TODOs regarding unit conversions and time intervals, will address this separately to not make this more complicated here...

[jacobcook1995]

LGTM! Just had one small query

[jacobcook1995]

Are these loops so you can solve the problem for each canopy layer for each cell individually?

[vgro]

yes, the scipy function only takes scalar values

[jacobcook1995]

oh didn't realise that (does make sense). Was going to suggest flagging this as a potential performance bottleneck that you could improve by using a vector rather than scalar solver but if the solver only takes scalars then that is moot really

[vgro]

I thinks it's really fast but I can add a flag, good point

[davidorme]

Late to this party, but I suspect that this might be a bottleneck - many definitions of a cell specific function. I think a gridded version of a secant solver (which is what is being used here) might be faster, but this is definitely something to circle back to rather than worry about now.

[davidorme]

Actually, having said that, I've got a similar gridded root finding in pyrealm.splash which has been intermittently bothering me, so I've looked back into this and it is possible to parallelise this. I think this is going to be more efficient because it can operate over all cells and layers simultaneously - and numpy is great at handling those mass calculations rather than iterating over cells and layers. I suspect I'd want @dalonsoa or @alexdewar to weigh in on the implementation, but I think this is doable.

It makes more sense to use canned versions like scipy.optimise.newton - unless they impose restrictions like only scalar inputs.

from typing import Callable

import matplotlib.pyplot as plt
import numpy as np
from numpy.typing import NDArray


def secant_method(
    f: Callable,
    x0: NDArray[np.floating],
    x1: NDArray[np.floating],
    tol: float,
    **kwargs,
) -> NDArray[np.floating]:
    """Secant root for array inputs"""

    # Starting difference from root
    unsolved = np.full_like(x0, True, dtype=np.bool_)
    n_iter = 0

    while np.any(unsolved):
        # Calculate fx0 and fx1
        fx0 = f(x0, **kwargs)
        fx1 = f(x1, **kwargs)

        # Update using secant method
        x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0)

        # Update unsolved values
        unsolved = np.abs(fx0) > tol
        x0, x1 = np.where(unsolved, x1, x0), np.where(unsolved, x2, x1)

        n_iter += 1
        print(f"{unsolved.sum()} unsolved at iteration {n_iter}")

    return x2

def func(x: NDArray[np.floating], b: NDArray[np.floating]):
    return x**2 - 10 + b

# Cell specific variation
b = np.array([[1, 2, 3, 4]])
# Starting points
x0 = np.full_like(b, 0)
x1 = np.full_like(b, 40)

# Solve
root = secant_method(func, x0=x0, x1=x1, tol=1e-8, b=b)

# Plot
x = np.broadcast_to(np.arange(0, 5, 0.01)[:, None], (500, 4))
plt.plot(x, func(x, b=b))
ax = plt.gca()
style = dict(linewidth=0.2, color="black")
ax.vlines(root, ymin=0, ymax=1, transform=ax.get_xaxis_transform(), **style)
ax.hlines(0, xmin=0, xmax=1, transform=ax.get_yaxis_transform(), **style)
plt.show()