Numerical integration of the soil module #191

jacobcook1995 · 2023-03-08T15:21:52Z

Description

This pull request adds functionality to integrate the soil model (or at least the simple dummy model that currently exists). This involves a function to set up and run the integration process, and a function that combines the existing functions into a form suitable for integration using scipy's solve_ivp function.

One thing that definitely needs thought is where DataArrays get used and where plain numpy arrays are used. I basically just adopted the structure that worked that was the smallest change from what I had previously. This might be an area where it's worth being more systematic as I am guessing this has performance implications (at least if data is being copied/altered).

This pull request adds dotmap as a dependency. I added it because I needed a container for mocking integration output for one of the tests. I considered using a dataclass rather than a DotMap for but it seemed like it would end up being significantly more longwinded. I am happy to try and implement this if adding dotmap as a dependancy is a problem.

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

…unction

…find it

dalonsoa

This is a good starting point, and you have definitely got the idea of wrapping the functionality to be used by the solver. I've just made a suggestion to make the code more generalisable down the line.

virtual_rainforest/core/base_model.py

virtual_rainforest/models/soil/carbon.py

virtual_rainforest/models/soil/soil_model.py

davidorme

Looking good - and I think this is the only way to go. Some minor comments. I don't think we need DotMap but it isn't the most expensive thing in the world.

Oh wait. I've just got why it is being used - because you're mocking, you need to pass something like the real output of solve_ivp forward to integrate_soil_model. I have to admit that for a simple integration test, I'd just run solve_ivp rather than mocking it. If we do mock it, then we can import the actual expected return class, which is basically just scipys own implementation of a dict with dot attribute access.

from scipy.integrate._ivp.ivp import OdeResult
# or avoiding the private special subclass
from scipy.optimize import OptimizeResult

and then we just define the expected output as:

OptimizeResult(
                success=True,
                y=np.array(
                    [
                        [5.000e-02, 3.210e-02],
                        [2.000e-02, 1.921e-02],
                        [4.500e00, 4.509e00],
                        [5.000e-01, 5.000e-01],
                        [3.210e-02, 5.000e-02],
                        [1.921e-02, 2.000e-02],
                        [4.509e00, 4.500e00],
                        [5.000e-01, 5.000e-01],
                    ]))

tests/models/soil/test_soil_model.py

virtual_rainforest/models/soil/soil_model.py

alexdewar

I've made some minor suggestions in places.

virtual_rainforest/models/soil/soil_model.py

alexdewar · 2023-03-13T10:47:02Z

virtual_rainforest/models/soil/soil_model.py


        # Update carbon pools (attributes and data object)
        # n.b. this also updates the data object automatically
-        self.carbon.update_soil_carbon_pools(self.data, carbon_pool_updates)
+        self.replace_soil_pools(updated_carbon_pools)

        # Finally increment timing
        self.next_update += self.update_interval

    def cleanup(self) -> None:
        """Placeholder function for soil model cleanup."""


(I know this isn't part of this PR, but I've only just noticed it.)

Do all models have a cleanup() method? What's it needed for? You don't normally have to worry about manually cleaning up resources with Python.

Not 100% sure tbh, we decided on that set of core functions for the standard model api a while back. We settled on this as one of the key functions but I'm not sure what its roll was envisioned as being. It was originally intended as being the opposite of spinup (I think we initially called it tear_down) but not sure what that entails.

In a similar vein, the setup function is a bit redundant now that from_config is defined

virtual_rainforest/models/soil/soil_model.py

… doesn't matter

dalonsoa

I've made a few suggestions - formatting will be wrong - to make the code more concise. I'm not an enthusiast of using locals, as you do in line 72, but I cannot think on a better, scalable method at the moment.

virtual_rainforest/models/soil/soil_model.py

jacobcook1995 · 2023-03-20T13:44:12Z

I do agree that using locals() is a bit weird, but I couldn't come up with an alternative (apart from vars() which is either worse or much the same)

jacobcook1995 · 2023-03-22T11:57:54Z

@dalonsoa I've removed the use of locals by passing the pool order as a dictionary rather than a list (so that values can be assigned to it).

I'm not 100% whether this is a reasonable thing to do in python, as don't really know how memory allocation works under the hood. Relatedly I'm also unsure whether

delta_pools_ordered = {
    str(name): np.array([])
    for name in self.data.data.keys()
    if str(name).startswith("soil_c_pool_")
}

or

delta_pools_ordered = {
    str(name): np.zeros(no_cells)
    for name in self.data.data.keys()
    if str(name).startswith("soil_c_pool_")
}

will work out being more efficient when memory is reallocated, e.g. when this step is carried out

delta_pools_ordered["soil_c_pool_lmwc"] = -lmwc_to_maom
delta_pools_ordered["soil_c_pool_maom"] = lmwc_to_maom

(if this isn't premature optimisation).

dalonsoa

I've made a small suggestion and a few questions, but otherwise it looks good!

dalonsoa · 2023-03-22T13:10:30Z

virtual_rainforest/models/soil/soil_model.py

@@ -47,8 +54,8 @@ class SoilModel(BaseModel):
    model_name = "soil"
    """An internal name used to register the model and schema"""
    required_init_vars = (


A probably out of place question, but just worth the check here, as it has come in a conversation with @vgro:

Are these variables the ones that need to be present in the data object to be able to initialise this class, i.e. to run the __init__ method?

If so, are these variables (these arrays) loaded from files when creating the data object?

As the soil_schema.json do not indicate that those files should be present in the input config, how is the user informed that they should be present - keeping aside trying to run the thing and getting an error when trying to initialise these model?

How should they be included in the input file

In my case these don't have to be present in the data object to initialise the class as my class doesn't have any attributes. I've just put all the things that need to be in data for the functions used by the SoilModel to work in here. This does cause initialisation to fail when the data object technically contains everything needed to initialise the class, which might not be ideal behaviour.

I don't have any data for the soil variables as of yet, so they are not at present loaded, but they should be loaded when the data object is created.

On the latter 2 questions, I haven't really given them much thought, and as far as I know there isn't any functionality in place to handle them as of yet.

dalonsoa · 2023-03-22T13:13:31Z

virtual_rainforest/models/soil/soil_model.py

+        self.data["soil_c_pool_lmwc"] = new_pools["soil_c_pool_lmwc"]
+        self.data["soil_c_pool_maom"] = new_pools["soil_c_pool_maom"]


For what I understand here, and based on what required_init_vars says, the value of these pools over time is not kept in the data object - or kept at all - right? In other words, you are no interested into tracking how these evolve over time.

I would definitely be interested in tracking the value of these pools over time. However, I couldn't see an easy way to add a time dimension to the existing data object so decided to just overwrite it for this first attempt at integrating the model

@jacobcook1995 Let me know if you want a chat about adding that time dimension! But I agree that it is fine to park that for now.

virtual_rainforest/models/soil/carbon.py

dalonsoa · 2023-03-22T13:22:51Z

virtual_rainforest/models/soil/carbon.py

+    delta_pools_ordered["soil_c_pool_lmwc"] = -lmwc_to_maom
+    delta_pools_ordered["soil_c_pool_maom"] = lmwc_to_maom


Assigning these to the dictionary keys seems a bit overkilling with just 2 pools, but I think it will make sense once there are more pools involved. In that sense, using a dictionary to ensure the order is the same than the one used in the inputs seems like the right move.

As you are totally replacing the existing, dummy array assigned to each dictionary entry, using during creation np.array([]) makes sense as you are just using as place holder.

alexdewar

I've got a couple of small queries/suggestions, but otherwise LGTM!

tests/models/soil/test_soil_model.py

virtual_rainforest/models/soil/soil_model.py

alexdewar · 2023-03-22T14:18:47Z

tests/models/soil/test_soil_model.py

@@ -6,7 +6,7 @@

 import numpy as np
 import pytest
-from dotmap import DotMap  # type: ignore
+from scipy.optimize import OptimizeResult  # type: ignore


Why are you ignoring the type here?

As far as I can tell scipy.optimize is not typed even though other scipy modules are typed. I could have missed them but also couldn't track down type stubs

alexdewar · 2023-03-22T14:21:48Z

tests/models/soil/test_soil_model.py


    log_check(caplog, expected_log)


+def test_order_independance(dummy_carbon_data, soil_model_fixture):


This is fine as is, but it could be part of a parametrized test, with an order argument that determines what order the data is in.

virtual_rainforest/models/soil/soil_model.py

alexdewar

Sorry, I meant to approve before...

…or gets documented

davidorme

LGTM

jacobcook1995 added 22 commits March 3, 2023 15:21

Moved bounds checking into soil model

213fbf4

Removed soil carbon class

5136be8

Moved pool update function to soil_model

3b99540

Added additional soil model tests to catch bad data objects

d24bf6a

Made soil pool update function a part of the SoilModel class

c142b4d

Added skelton function to run the integration

fb82d76

Added basic structure to the integrate function

916592d

Removed dependance on dt from calculate_soil_carbon_updates function

9bc5e07

Removed direct dependance of soil carbon functions on data object

4c68cf1

Setup extraction of the variables from a single numpy array

b9765a2

Switched to returning a 1D numpy data array rather than a named dataset

c576eaa

Added step to generate initial numpy array

e10e081

Switched to extracting number of cells from the Grid object

88b2399

Documented data unpacking order

962ac7e

Added error handling for failed numerical integration

2e7f3ec

Added test for construct_full_soil_model function

d643431

Starting returning a populated data array from the soil integration f…

037f91e

…unction

Added dotmap as a dependancy

0e5b2d3

Added test of soil model integration function

a76b48c

Changed pool update function to replace rather than increment pools

8802c5d

Moved where Integration error is defined so that sphinx can actually …

3ad710a

…find it

Switched to using integrate function in soil model

c4c6d0b

jacobcook1995 requested review from dalonsoa, davidorme and alexdewar March 8, 2023 15:22

dalonsoa requested changes Mar 10, 2023

View reviewed changes

davidorme requested changes Mar 13, 2023

View reviewed changes

tests/models/soil/test_soil_model.py Outdated Show resolved Hide resolved

virtual_rainforest/models/soil/soil_model.py Outdated Show resolved Hide resolved

alexdewar requested changes Mar 13, 2023

View reviewed changes

jacobcook1995 added 2 commits March 14, 2023 14:56

Removed unnecessary temporary variables

4b2d1fe

Removed another unnecessary temporary variable

89707da

jacobcook1995 added 7 commits March 17, 2023 09:59

Integration result is now unpacked based on order within data object

78393b9

Updated naming of soil carbon pools

605165a

Made slice making list comp a function in its own right

934a290

Removed explicit ordering from construct_full_soil_model function

0ec4370

Added explict test that order that pools are added to the data object…

bfdef41

… doesn't matter

Added ordering to function to calculate the soil carbon pool updates

62d7b7c

Moved pool_order creation higher in function flow

8ead71b

jacobcook1995 requested review from alexdewar, davidorme and dalonsoa March 20, 2023 13:11

dalonsoa requested changes Mar 20, 2023

View reviewed changes

virtual_rainforest/models/soil/soil_model.py Outdated Show resolved Hide resolved

virtual_rainforest/models/soil/soil_model.py Outdated Show resolved Hide resolved

virtual_rainforest/models/soil/soil_model.py Outdated Show resolved Hide resolved

jacobcook1995 added 4 commits March 20, 2023 13:45

Made Diego's suggested edits to make ordering more consise

7a9a88b

Moved IntegrationError to more sensible location

2dcc35e

Improved docstring for SoilModel class

de21f63

Removed use of locals() by passing dictonary rather than list

f2f2f17

jacobcook1995 requested a review from dalonsoa March 22, 2023 11:58

dalonsoa approved these changes Mar 22, 2023

View reviewed changes

Switch to Diego's suggested method of condensing the cocatenation

58bb2e5

alexdewar requested changes Mar 22, 2023

View reviewed changes

alexdewar self-requested a review March 22, 2023 14:24

alexdewar approved these changes Mar 22, 2023

View reviewed changes

jacobcook1995 added 4 commits March 22, 2023 15:12

Removed unnecessary temporary variable

324a7ea

Simplified mocking by using patch.object

ffb224d

Merge branch 'develop' into feature/integrate_soil_model

90877bd

Changed soil_model autodocumentation to make sure that IntegrationErr…

f91ab52

…or gets documented

davidorme approved these changes Mar 28, 2023

View reviewed changes

jacobcook1995 merged commit 7ce9b32 into develop Mar 28, 2023

jacobcook1995 deleted the feature/integrate_soil_model branch March 28, 2023 09:35

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Numerical integration of the soil module #191

Numerical integration of the soil module #191

jacobcook1995 commented Mar 8, 2023

dalonsoa left a comment

davidorme left a comment

alexdewar left a comment

alexdewar Mar 13, 2023

jacobcook1995 Mar 14, 2023

dalonsoa left a comment

jacobcook1995 commented Mar 20, 2023

jacobcook1995 commented Mar 22, 2023 •

edited

Loading

dalonsoa left a comment

dalonsoa Mar 22, 2023

jacobcook1995 Mar 22, 2023 •

edited

Loading

dalonsoa Mar 22, 2023

jacobcook1995 Mar 22, 2023

davidorme Mar 28, 2023

dalonsoa Mar 22, 2023

alexdewar left a comment

alexdewar Mar 22, 2023

jacobcook1995 Mar 22, 2023

alexdewar Mar 22, 2023

alexdewar left a comment

davidorme left a comment

		self.data["soil_c_pool_lmwc"] = new_pools["soil_c_pool_lmwc"]
		self.data["soil_c_pool_maom"] = new_pools["soil_c_pool_maom"]

		delta_pools_ordered["soil_c_pool_lmwc"] = -lmwc_to_maom
		delta_pools_ordered["soil_c_pool_maom"] = lmwc_to_maom


		log_check(caplog, expected_log)


		def test_order_independance(dummy_carbon_data, soil_model_fixture):

Numerical integration of the soil module #191

Numerical integration of the soil module #191

Conversation

jacobcook1995 commented Mar 8, 2023

Description

Type of change

Key checklist

Further checks

dalonsoa left a comment

Choose a reason for hiding this comment

davidorme left a comment

Choose a reason for hiding this comment

alexdewar left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dalonsoa left a comment

Choose a reason for hiding this comment

jacobcook1995 commented Mar 20, 2023

jacobcook1995 commented Mar 22, 2023 • edited Loading

dalonsoa left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jacobcook1995 Mar 22, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alexdewar left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alexdewar left a comment

Choose a reason for hiding this comment

davidorme left a comment

Choose a reason for hiding this comment

jacobcook1995 commented Mar 22, 2023 •

edited

Loading

jacobcook1995 Mar 22, 2023 •

edited

Loading