API Reference

Crop modules

diffwofost.physical_models.crop.phenology.DVS_Phenology

Bases: SimulationObject

Implements the algorithms for phenologic development in WOFOST.

Phenologic development in WOFOST is expresses using a unitless scale which takes the values 0 at emergence, 1 at Anthesis (flowering) and 2 at maturity. This type of phenological development is mainly representative for cereal crops. All other crops that are simulated with WOFOST are forced into this scheme as well, although this may not be appropriate for all crops. For example, for potatoes development stage 1 represents the start of tuber formation rather than flowering.

Phenological development is mainly governed by temperature and can be modified by the effects of day length and vernalization during the period before Anthesis. After Anthesis, only temperature influences the development rate.

Simulation parameters

Name	Description	Type	Unit
TSUMEM	Temperature sum from sowing to emergence	SCr
TBASEM	Base temperature for emergence	SCr
TEFFMX	Maximum effective temperature for emergence	SCr
TSUM1	Temperature sum from emergence to anthesis	SCr
TSUM2	Temperature sum from anthesis to maturity	SCr
IDSL	Switch for development options: temp only (0), +daylength	SCr	-
	(1), +vernalization (>=2)
DLO	Optimal daylength for phenological development	SCr	hr
DLC	Critical daylength for phenological development	SCr	hr
DVSI	Initial development stage at emergence (may be >0 for	SCr	-
	transplanted crops)
DVSEND	Final development stage	SCr	-
DTSMTB	Daily increase in temperature sum as a function of daily	TCr
	mean temperature

State variables

Name	Description	Pbl	Unit
DVS	Development stage	Y	-
TSUM	Temperature sum	N
TSUME	Temperature sum for emergence	N
DOS	Day of sowing	N	-
DOE	Day of emergence	N	-
DOA	Day of Anthesis	N	-
DOM	Day of maturity	N	-
DOH	Day of harvest	N	-
STAGE	Current stage (`emerging\|vegetative\|reproductive\|mature`)	N	-

Rate variables

Name	Description	Pbl
DTSUME	Increase in temperature sum for emergence	N
DTSUM	Increase in temperature sum for anthesis or maturity	N
DVR	Development rate	Y

External dependencies:

None

Signals sent or handled

DVS_Phenology sends the crop_finish signal when maturity is reached and the end_type is 'maturity' or 'earliest'.

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
DVS	...
TSUM	...

[!NOTE] Notice that the gradient ∂DVS/∂TEFFMX is zero.

[!NOTE] The parameter IDSL it is not differentiable since it is a switch.

Methods:

initialize –

:param day: start date of the simulation
calc_rates –

Compute daily phenological development rates.
integrate –

Integrate phenology states and manage stage transitions.
get_variable –

Return the value of the specified state or rate variable.

Attributes:

device –

Get device from ComputeConfig.
dtype –

Get dtype from ComputeConfig.

device `property`

device

Get device from ComputeConfig.

dtype `property`

dtype

Get dtype from ComputeConfig.

initialize

initialize(day, kiosk, parvalues)

:param day: start date of the simulation

:param kiosk: variable kiosk of this PCSE instance :param parvalues: ParameterProvider object providing parameters as key/value pairs

Source code in src/diffwofost/physical_models/crop/phenology.py

def initialize(self, day, kiosk, parvalues):
    """:param day: start date of the simulation

    :param kiosk: variable kiosk of this PCSE  instance
    :param parvalues: `ParameterProvider` object providing parameters as
            key/value pairs
    """
    self.params = self.Parameters(parvalues)
    self.params_shape = _get_params_shape(self.params)

    # Initialize vernalisation for IDSL>=2
    # It has to be done in advance to get the correct params_shape
    IDSL = _broadcast_to(
        self.params.IDSL, self.params_shape, dtype=self.dtype, device=self.device
    )
    self.params.IDSL = IDSL
    if torch.any(IDSL >= 2):
        if self.params_shape != ():
            self.vernalisation = Vernalisation(
                day, kiosk, parvalues, dvs_shape=self.params_shape
            )
        else:
            self.vernalisation = Vernalisation(day, kiosk, parvalues)
        if self.vernalisation.params_shape != self.params_shape:
            self.params_shape = self.vernalisation.params_shape
    else:
        self.vernalisation = None

    # After Vernalisation initialization the final params_shape may have changed.
    self._cast_and_broadcast_params()

    # Create scalar constants once at the beginning to avoid recreating them
    self._ones = torch.ones(self.params_shape, dtype=self.dtype, device=self.device)
    self._zeros = torch.zeros(self.params_shape, dtype=self.dtype, device=self.device)
    self._epsilon = torch.tensor(1e-8, dtype=self.dtype, device=self.device)

    # Initialize rates and kiosk
    self.rates = self.RateVariables(kiosk)
    self.kiosk = kiosk

    self._connect_signal(self._on_CROP_FINISH, signal=signals.crop_finish)

    # Define initial states
    DVS, DOS, DOE, STAGE = self._get_initial_stage(day)
    DVS = _broadcast_to(DVS, self.params_shape, dtype=self.dtype, device=self.device)

    # Initialize all date tensors with -1 (not yet occurred)
    DOS = _broadcast_to(DOS, self.params_shape, dtype=self.dtype, device=self.device)
    DOE = _broadcast_to(DOE, self.params_shape, dtype=self.dtype, device=self.device)
    DOA = torch.full(self.params_shape, -1.0, dtype=self.dtype, device=self.device)
    DOM = torch.full(self.params_shape, -1.0, dtype=self.dtype, device=self.device)
    DOH = torch.full(self.params_shape, -1.0, dtype=self.dtype, device=self.device)
    STAGE = _broadcast_to(STAGE, self.params_shape, dtype=self.dtype, device=self.device)

    # Also ensure TSUM and TSUME are properly shaped
    TSUM = torch.zeros(
        self.params_shape, dtype=self.dtype, device=self.device, requires_grad=True
    )
    TSUME = torch.zeros(
        self.params_shape, dtype=self.dtype, device=self.device, requires_grad=True
    )

    self.states = self.StateVariables(
        kiosk,
        publish="DVS",
        TSUM=TSUM,
        TSUME=TSUME,
        DVS=DVS,
        DOS=DOS,
        DOE=DOE,
        DOA=DOA,
        DOM=DOM,
        DOH=DOH,
        STAGE=STAGE,
    )

calc_rates

calc_rates(day, drv)

Compute daily phenological development rates.

Parameters:

day (date) –

Current simulation date.
drv –

Meteorological driver object with at least TEMP and LAT.

Logic

Photoperiod reduction (DVRED) if IDSL >= 1 using daylength.
Vernalisation factor (VERNFAC) if IDSL >= 2 and in vegetative stage.
Stage-specific:
emerging: temperature sum for emergence (DTSUME), DVR via TSUMEM.
vegetative: temperature sum (DTSUM) scaled by VERNFAC and DVRED.
reproductive: temperature sum (DTSUM) only temperature-driven.
mature: all rates zero.

Sets

r.DTSUME, r.DTSUM, r.DVR.

Raises:

PCSEError –

If STAGE unrecognized.

Source code in src/diffwofost/physical_models/crop/phenology.py

@prepare_rates
def calc_rates(self, day, drv):
    """Compute daily phenological development rates.

    Args:
        day (datetime.date): Current simulation date.
        drv: Meteorological driver object with at least TEMP and LAT.

    Logic:
        1. Photoperiod reduction (DVRED) if IDSL >= 1 using daylength.
        2. Vernalisation factor (VERNFAC) if IDSL >= 2 and in vegetative stage.
        3. Stage-specific:
           - emerging: temperature sum for emergence (DTSUME), DVR via TSUMEM.
           - vegetative: temperature sum (DTSUM) scaled by VERNFAC and DVRED.
           - reproductive: temperature sum (DTSUM) only temperature-driven.
           - mature: all rates zero.

    Sets:
        r.DTSUME, r.DTSUM, r.DVR.

    Raises:
        PCSEError: If STAGE unrecognized.

    """
    p = self.params
    r = self.rates
    s = self.states
    shape = self.params_shape

    # Day length sensitivity
    DAYLP = daylength(day, drv.LAT)
    DAYLP_t = _broadcast_to(DAYLP, shape, dtype=self.dtype, device=self.device)
    # Compute DVRED conditionally based on IDSL >= 1
    safe_den = p.DLO - p.DLC
    safe_den = safe_den.sign() * torch.maximum(torch.abs(safe_den), self._epsilon)
    dvred_active = torch.clamp((DAYLP_t - p.DLC) / safe_den, 0.0, 1.0)
    DVRED = torch.where(p.IDSL >= 1, dvred_active, self._ones)

    # Vernalisation factor - always compute if module exists
    VERNFAC = self._ones
    if hasattr(self, "vernalisation") and self.vernalisation is not None:
        # Always call calc_rates (it handles stage internally now)
        self.vernalisation.calc_rates(day, drv)
        # Apply vernalisation only where IDSL >= 2 AND in vegetative stage
        is_vegetative = s.STAGE == 1
        VERNFAC = torch.where(
            (p.IDSL >= 2) & is_vegetative,
            self.kiosk["VERNFAC"],
            self._ones,
        )

    TEMP = _get_drv(drv.TEMP, shape, self.dtype, self.device)

    # Initialize all rate variables
    r.DTSUME = self._zeros
    r.DTSUM = self._zeros
    r.DVR = self._zeros

    # Compute rates for emerging stage (STAGE == 0)
    is_emerging = s.STAGE == 0
    if torch.any(is_emerging):
        temp_diff = TEMP - p.TBASEM
        # Ensure the maximum effective temperature difference is non-negative
        max_diff = torch.clamp(p.TEFFMX - p.TBASEM, min=0.0)
        dtsume_emerging = torch.clamp(temp_diff, min=0.0)
        dtsume_emerging = torch.minimum(dtsume_emerging, max_diff)
        safe_den = p.TSUMEM
        safe_den = safe_den.sign() * torch.maximum(torch.abs(safe_den), self._epsilon)
        dvr_emerging = 0.1 * dtsume_emerging / safe_den

        r.DTSUME = torch.where(is_emerging, dtsume_emerging, r.DTSUME)
        r.DVR = torch.where(is_emerging, dvr_emerging, r.DVR)

    # Compute rates for vegetative stage (STAGE == 1)
    is_vegetative = s.STAGE == 1
    if torch.any(is_vegetative):
        dtsum_vegetative = p.DTSMTB(TEMP) * VERNFAC * DVRED
        safe_den = p.TSUM1
        safe_den = safe_den.sign() * torch.maximum(torch.abs(safe_den), self._epsilon)
        dvr_vegetative = dtsum_vegetative / safe_den

        r.DTSUM = torch.where(is_vegetative, dtsum_vegetative, r.DTSUM)
        r.DVR = torch.where(is_vegetative, dvr_vegetative, r.DVR)

    # Compute rates for reproductive stage (STAGE == 2)
    is_reproductive = s.STAGE == 2
    if torch.any(is_reproductive):
        dtsum_reproductive = p.DTSMTB(TEMP)
        safe_den = p.TSUM2
        safe_den = safe_den.sign() * torch.maximum(torch.abs(safe_den), self._epsilon)
        dvr_reproductive = dtsum_reproductive / safe_den

        r.DTSUM = torch.where(is_reproductive, dtsum_reproductive, r.DTSUM)
        r.DVR = torch.where(is_reproductive, dvr_reproductive, r.DVR)

    # Mature stage (STAGE == 3) keeps zeros (already initialized)

    msg = "Finished rate calculation for %s"
    self.logger.debug(msg % day)

integrate

integrate(day, delt=1.0)

Integrate phenology states and manage stage transitions.

Parameters:

day (date) –

Current simulation day.
delt (float, default: 1.0 ) –

Timestep length in days (default 1.0).

Sequence

Integrates vernalisation module if active and in vegetative stage.
Accumulates TSUME, TSUM, advances DVS by DVR.
Checks threshold crossings to move through stages: emerging -> vegetative (DVS >= 0) vegetative -> reproductive (DVS >= 1) reproductive -> mature (DVS >= DVSEND)

Side Effects

Emits crop_emerged signal on emergence.
Emits crop_finish signal at maturity if end type matches.

Notes

Caps DVS at stage boundary values.

Raises:

PCSEError –

If STAGE undefined.

Source code in src/diffwofost/physical_models/crop/phenology.py

@prepare_states
def integrate(self, day, delt=1.0):
    """Integrate phenology states and manage stage transitions.

    Args:
        day (datetime.date): Current simulation day.
        delt (float, optional): Timestep length in days (default 1.0).

    Sequence:
        - Integrates vernalisation module if active and in vegetative stage.
        - Accumulates TSUME, TSUM, advances DVS by DVR.
        - Checks threshold crossings to move through stages:
            emerging -> vegetative (DVS >= 0)
            vegetative -> reproductive (DVS >= 1)
            reproductive -> mature (DVS >= DVSEND)

    Side Effects:
        - Emits crop_emerged signal on emergence.
        - Emits crop_finish signal at maturity if end type matches.

    Notes:
        Caps DVS at stage boundary values.

    Raises:
        PCSEError: If STAGE undefined.

    """
    p = self.params
    r = self.rates
    s = self.states
    shape = self.params_shape

    # Integrate vernalisation module
    if self.vernalisation:
        # Save a copy of state
        state_copy = _snapshot_state(self.vernalisation.states)
        mask_IDSL = p.IDSL >= 2

        # Check if any element is in vegetative stage i.e. stage 1
        mask_STAGE = mask_IDSL & (s.STAGE == 1)
        self.vernalisation.integrate(day, delt)
        state_integrated = _snapshot_state(self.vernalisation.states)

        # Restore original state
        _restore_state(self.vernalisation.states, state_copy)
        self.vernalisation.touch()
        state_touched = _snapshot_state(self.vernalisation.states)

        # Apply the masks
        for name in state_copy:
            # results of vernalisation module
            vernalisation_states = torch.where(
                mask_STAGE, state_integrated[name], state_touched[name]
            )
            setattr(
                self.vernalisation.states,
                name,
                torch.where(mask_IDSL, vernalisation_states, state_copy[name]),
            )

    # Integrate phenologic states
    s.TSUME = s.TSUME + r.DTSUME
    s.DVS = s.DVS + r.DVR
    s.TSUM = s.TSUM + r.DTSUM

    day_ordinal = torch.tensor(day.toordinal(), dtype=self.dtype, device=self.device)

    # Check transitions for emerging -> vegetative (STAGE 0 -> 1)
    is_emerging = s.STAGE == 0
    should_emerge = is_emerging & (s.DVS >= 0.0)
    s.STAGE = torch.where(
        should_emerge, torch.ones(shape, dtype=torch.long, device=self.device), s.STAGE
    )
    s.DOE = torch.where(
        should_emerge,
        torch.full(shape, day_ordinal, dtype=self.dtype, device=self.device),
        s.DOE,
    )
    s.DVS = torch.where(should_emerge, torch.clamp(s.DVS, max=0.0), s.DVS)

    # Send signal if any crop emerged (only once per day)
    if torch.any(should_emerge):
        self._send_signal(signals.crop_emerged)

    # Check transitions for vegetative -> reproductive (STAGE 1 -> 2)
    is_vegetative = s.STAGE == 1
    should_flower = is_vegetative & (s.DVS >= 1.0)
    s.STAGE = torch.where(
        should_flower, torch.full(shape, 2, dtype=torch.long, device=self.device), s.STAGE
    )
    s.DOA = torch.where(
        should_flower,
        torch.full(shape, day_ordinal, dtype=self.dtype, device=self.device),
        s.DOA,
    )
    s.DVS = torch.where(should_flower, torch.clamp(s.DVS, max=1.0), s.DVS)

    # Check transitions for reproductive -> mature (STAGE 2 -> 3)
    is_reproductive = s.STAGE == 2
    should_mature = is_reproductive & (s.DVS >= p.DVSEND)
    s.STAGE = torch.where(
        should_mature, torch.full(shape, 3, dtype=torch.long, device=self.device), s.STAGE
    )
    s.DOM = torch.where(
        should_mature,
        torch.full(shape, day_ordinal, dtype=self.dtype, device=self.device),
        s.DOM,
    )
    s.DVS = torch.where(should_mature, torch.minimum(s.DVS, p.DVSEND), s.DVS)

    # Send crop_finish signal if maturity reached for one.
    # assumption is that all elements mature simultaneously
    # TODO: revisit this when fixing engine for agromanager
    if torch.any(should_mature) and p.CROP_END_TYPE in ["maturity", "earliest"]:
        self._send_signal(
            signal=signals.crop_finish,
            day=day,
            finish_type="maturity",
            crop_delete=True,
        )

    msg = "Finished state integration for %s"
    self.logger.debug(msg % day)

get_variable

get_variable(varname)

Return the value of the specified state or rate variable.

:param varname: Name of the variable.

Note that the get_variable() will searches for varname exactly as specified (case sensitive).

Source code in src/diffwofost/physical_models/crop/phenology.py

def get_variable(self, varname):
    # TODO: should be removed while fixing #49. this is needed because
    # conditions are applied on STAGE in pcse.crop.wofost72.py
    """Return the value of the specified state or rate variable.

    :param varname: Name of the variable.

    Note that the `get_variable()` will searches for `varname` exactly
    as specified (case sensitive).
    """
    if varname == "STAGE":
        # Return string representation of current stage
        stage_map = {
            0: "emerging",
            1: "vegetative",
            2: "reproductive",
            3: "mature",
        }
        stage_value = self.states.STAGE
        if stage_value.dim() != 0:
            stage_id = stage_value.flatten()[0].item()
        else:
            stage_id = stage_value.item()
        return stage_map[stage_id]

    # Search for variable in the current object, then traverse the hierarchy
    value = None
    if hasattr(self.states, varname):
        value = getattr(self.states, varname)
    elif hasattr(self.rates, varname):
        value = getattr(self.rates, varname)
    # Query individual sub-SimObject for existence of variable v
    else:
        for simobj in self.subSimObjects:
            value = simobj.get_variable(varname)
            if value is not None:
                break
    return value

diffwofost.physical_models.crop.partitioning.DVS_Partitioning

Bases: _BaseDVSPartitioning

Class for assimilate partitioning based on development stage (DVS).

DVS_Partitioning calculates the partitioning of the assimilates to roots, stems, leaves and storage organs using fixed partitioning tables as a function of crop development stage. The available assimilates are first split into below-ground and aboveground using the values in FRTB. In a second stage they are split into leaves (FLTB), stems (FSTB) and storage organs (FOTB).

Since the partitioning fractions are derived from the state variable DVS they are regarded state variables as well.

Simulation parameters (To be provided in cropdata dictionary):

Name	Description	Type	Unit
FRTB	Partitioning to roots as a function of development stage	TCr	-
FSTB	Partitioning to stems as a function of development stage	TCr	-
FLTB	Partitioning to leaves as a function of development stage	TCr	-
FOTB	Partitioning to storage organs as a function of development stage	TCr	-

State variables

Name	Description	Pbl	Unit
FR	Fraction partitioned to roots	Y	-
FS	Fraction partitioned to stems	Y	-
FL	Fraction partitioned to leaves	Y	-
FO	Fraction partitioned to storage organs	Y	-
PF	Partitioning factors packed in tuple	N	-

Rate variables

None

External dependencies:

Name	Description	Provided by	Unit
DVS	Crop development stage	DVS_Phenology	-

Outputs

Name	Description	Pbl	Unit
FR	Fraction partitioned to roots	Y	-
FL	Fraction partitioned to leaves	Y	-
FS	Fraction partitioned to stems	Y	-
FO	Fraction partitioned to storage organs	Y	-

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
FR	FRTB, DVS
FL	FLTB, DVS
FS	FSTB, DVS
FO	FOTB, DVS

Exceptions raised

A PartitioningError is raised if the partitioning coefficients to leaves, stems and storage organs on a given day do not add up to 1.

Methods:

initialize –

Initialize the DVS_Partitioning simulation object.
integrate –

Update partitioning factors based on development stage (DVS).
calc_rates –

Return partitioning factors based on current DVS.

initialize

initialize(day, kiosk, parvalues)

Initialize the DVS_Partitioning simulation object.

Parameters:

day –

Start date of the simulation.
kiosk (VariableKiosk) –

Variable kiosk of this PCSE instance.
parvalues (ParameterProvider) –

Object providing parameters as key/value pairs.

Source code in src/diffwofost/physical_models/crop/partitioning.py

def initialize(self, day, kiosk, parvalues):
    """Initialize the DVS_Partitioning simulation object.

    Args:
        day: Start date of the simulation.
        kiosk (VariableKiosk): Variable kiosk of this PCSE instance.
        parvalues (ParameterProvider): Object providing parameters as
            key/value pairs.
    """
    self._initialize_from_tables(kiosk, parvalues)

integrate

integrate(day, delt=1.0)

Update partitioning factors based on development stage (DVS).

Source code in src/diffwofost/physical_models/crop/partitioning.py

@prepare_states
def integrate(self, day, delt=1.0):
    """Update partitioning factors based on development stage (DVS)."""
    self._update_from_tables()

calc_rates

calc_rates(day, drv)

Return partitioning factors based on current DVS.

Rate calculation does nothing for partitioning as it is a derived state.

Source code in src/diffwofost/physical_models/crop/partitioning.py

def calc_rates(self, day, drv):
    """Return partitioning factors based on current DVS.

    Rate calculation does nothing for partitioning as it is a derived state.
    """
    return self.states.PF

diffwofost.physical_models.crop.assimilation.WOFOST72_Assimilation

Bases: SimulationObject

Class implementing a WOFOST/SUCROS style assimilation routine.

WOFOST calculates the daily gross CO2 assimilation rate of a crop from the absorbed radiation and the photosynthesis-light response curve of individual leaves. This response is dependent on temperature and leaf age. The absorbed radiation is calculated from the total incoming radiation and the leaf area. Daily gross CO2 assimilation is obtained by integrating the assimilation rates over the leaf layers and over the day.

Simulation parameters (provide in cropdata dictionary)

Name	Description	Type	Unit
AMAXTB	Max. leaf CO2 assimilation rate as function of DVS	TCr	kg CO2 ha⁻¹ leaf h⁻¹
EFFTB	Light use effic. single leaf as a function of daily mean temperature	TCr	kg CO2 ha⁻¹ h⁻¹ /(J m⁻² s⁻¹)
KDIFTB	Extinction coefficient for diffuse visible light as function of DVS	TCr	-
TMPFTB	Reduction factor on AMAX as function of daily mean temperature	TCr	-
TMNFTB	Reduction factor on AMAX as function of daily minimum temperature	TCr	-

Rate variables This class returns the potential gross assimilation rate 'PGASS' directly from the __call__() method, but also includes it as a rate variable.

Name	Description	Pbl	Unit
PGASS	Potential gross assimilation	Y	kg CH2O ha⁻¹ d⁻¹

External dependencies

Name	Description	Provided by	Unit
DVS	Crop development stage	DVS_Phenology	-
LAI	Leaf area index	Leaf_dynamics	-

Weather inputs used

Name	Description	Unit
IRRAD	Daily shortwave radiation	J m⁻² d⁻¹
DTEMP	Daily mean temperature	°C
TMIN	Daily minimum temperature	°C
LAT	Latitude	degrees

Outputs

Name	Description	Pbl	Unit
PGASS	Potential gross assimilation	Y	kg CH2O ha⁻¹ d⁻¹

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
PGASS	AMAXTB, EFFTB, KDIFTB, TMPFTB, TMNFTB

Methods:

initialize –

Initialize the assimilation module.
calc_rates –

Compute the potential gross assimilation rate (PGASS).
__call__ –

Calculate and return the potential gross assimilation rate (PGASS).
integrate –

No state variables to integrate for this module.

Attributes:

device –

Get device from ComputeConfig.
dtype –

Get dtype from ComputeConfig.

device `property`

device

Get device from ComputeConfig.

dtype `property`

dtype

Get dtype from ComputeConfig.

initialize

initialize(day: date, kiosk: VariableKiosk, parvalues: ParameterProvider) -> None

Initialize the assimilation module.

Source code in src/diffwofost/physical_models/crop/assimilation.py

def initialize(
    self, day: datetime.date, kiosk: VariableKiosk, parvalues: ParameterProvider
) -> None:
    """Initialize the assimilation module."""
    self.kiosk = kiosk
    self.params = self.Parameters(parvalues)
    self.params_shape = _get_params_shape(self.params)
    self.rates = self.RateVariables(kiosk, publish=["PGASS"])

    # 7-day running average buffer for TMIN (stored as tensors).
    self._tmn_window = deque(maxlen=7)
    self._tmn_window_mask = deque(maxlen=7)
    # Reused scalar constants
    self._epsilon = torch.tensor(1e-12, dtype=self.dtype, device=self.device)

calc_rates

calc_rates(day: date = None, drv: WeatherDataContainer = None) -> None

Compute the potential gross assimilation rate (PGASS).

Source code in src/diffwofost/physical_models/crop/assimilation.py

@prepare_rates
def calc_rates(self, day: datetime.date = None, drv: WeatherDataContainer = None) -> None:
    """Compute the potential gross assimilation rate (PGASS)."""
    p = self.params
    r = self.rates
    k = self.kiosk

    _exist_required_external_variables(k)

    # External states
    dvs = _broadcast_to(k["DVS"], self.params_shape, dtype=self.dtype, device=self.device)
    lai = _broadcast_to(k["LAI"], self.params_shape, dtype=self.dtype, device=self.device)

    # Weather drivers
    irrad = _get_drv(drv.IRRAD, self.params_shape, dtype=self.dtype, device=self.device)
    dtemp = _get_drv(drv.DTEMP, self.params_shape, dtype=self.dtype, device=self.device)
    tmin = _get_drv(drv.TMIN, self.params_shape, dtype=self.dtype, device=self.device)

    # Assimilation is zero before crop emergence (DVS < 0)
    dvs_mask = (dvs >= 0).to(dtype=self.dtype)
    # 7-day running average of TMIN
    self._tmn_window.appendleft(tmin * dvs_mask)
    self._tmn_window_mask.appendleft(dvs_mask)
    tmin_stack = torch.stack(list(self._tmn_window), dim=0)
    mask_stack = torch.stack(list(self._tmn_window_mask), dim=0)
    tminra = tmin_stack.sum(dim=0) / (mask_stack.sum(dim=0) + 1e-8)

    # Astronomical variables (computed with PCSE util; then broadcast to tensors)
    lat = _as_python_float(drv.LAT)
    irrad_for_astro = _as_python_float(drv.IRRAD)
    dayl, _daylp, sinld, cosld, difpp, _atmtr, dsinbe, _angot = astro(day, lat, irrad_for_astro)

    dayl_t = _broadcast_to(dayl, self.params_shape, dtype=self.dtype, device=self.device)
    sinld_t = _broadcast_to(sinld, self.params_shape, dtype=self.dtype, device=self.device)
    cosld_t = _broadcast_to(cosld, self.params_shape, dtype=self.dtype, device=self.device)
    difpp_t = _broadcast_to(difpp, self.params_shape, dtype=self.dtype, device=self.device)
    dsinbe_t = _broadcast_to(dsinbe, self.params_shape, dtype=self.dtype, device=self.device)

    # Parameter tables
    amax = p.AMAXTB(dvs)
    amax = amax * p.TMPFTB(dtemp)
    kdif = p.KDIFTB(dvs)
    eff = p.EFFTB(dtemp)

    dtga = totass7(
        dayl_t,
        amax,
        eff,
        lai,
        kdif,
        irrad,
        difpp_t,
        dsinbe_t,
        sinld_t,
        cosld_t,
        epsilon=self._epsilon,
    )

    # Correction for low minimum temperature potential
    dtga = dtga * p.TMNFTB(tminra)

    # Convert kg CO2 -> kg CH2O
    pgass = dtga * (30.0 / 44.0)

    # Assimilation is zero before crop emergence (DVS < 0)
    r.PGASS = pgass * dvs_mask
    return r.PGASS

call

__call__(day: date = None, drv: WeatherDataContainer = None) -> Tensor

Calculate and return the potential gross assimilation rate (PGASS).

Source code in src/diffwofost/physical_models/crop/assimilation.py

def __call__(self, day: datetime.date = None, drv: WeatherDataContainer = None) -> torch.Tensor:
    """Calculate and return the potential gross assimilation rate (PGASS)."""
    return self.calc_rates(day, drv)

integrate

integrate(day: date = None, delt=1.0) -> None

No state variables to integrate for this module.

Source code in src/diffwofost/physical_models/crop/assimilation.py

@prepare_states
def integrate(self, day: datetime.date = None, delt=1.0) -> None:
    """No state variables to integrate for this module."""
    return

diffwofost.physical_models.crop.leaf_dynamics.WOFOST_Leaf_Dynamics

Bases: SimulationObject

Leaf dynamics for the WOFOST crop model.

Implementation of biomass partitioning to leaves, growth and senenscence of leaves. WOFOST keeps track of the biomass that has been partitioned to the leaves for each day (variable LV), which is called a leaf class). For each leaf class the leaf age (variable 'LVAGE') and specific leaf area (variable SLA) are also registered. Total living leaf biomass is calculated by summing the biomass values for all leaf classes. Similarly, leaf area is calculated by summing leaf biomass times specific leaf area (LV * SLA).

Senescense of the leaves can occur as a result of physiological age, drought stress or self-shading.

Simulation parameters (provide in cropdata dictionary)

Name	Description	Type	Unit
RGRLAI	Maximum relative increase in LAI.	SCr	ha ha⁻¹ d⁻¹
SPAN	Life span of leaves growing at 35 Celsius	SCr	d
TBASE	Lower threshold temp. for ageing of leaves	SCr	C
PERDL	Max. relative death rate of leaves due to water stress	SCr
TDWI	Initial total crop dry weight	SCr	kg ha⁻¹
KDIFTB	Extinction coefficient for diffuse visible light as function of DVS	TCr
SLATB	Specific leaf area as a function of DVS	TCr	ha kg⁻¹

State variables

Name	Description	Pbl	Unit
LV	Leaf biomass per leaf class	N	kg ha⁻¹
SLA	Specific leaf area per leaf class	N	ha kg⁻¹
LVAGE	Leaf age per leaf class	N	d
LVSUM	Sum of LV	N	kg ha⁻¹
LAIEM	LAI at emergence	N	-
LASUM	Total leaf area as sum of LV*SLA, not including stem and pod area	N	-
LAIEXP	LAI value under theoretical exponential growth	N	-
LAIMAX	Maximum LAI reached during growth cycle	N	-
LAI	Leaf area index, including stem and pod area	Y	-
WLV	Dry weight of living leaves	Y	kg ha⁻¹
DWLV	Dry weight of dead leaves	N	kg ha⁻¹
TWLV	Dry weight of total leaves (living + dead)	Y	kg ha⁻¹

Rate variables

Name	Description	Pbl	Unit
GRLV	Growth rate leaves	N	kg ha⁻¹ d⁻¹
DSLV1	Death rate leaves due to water stress	N	kg ha⁻¹ d⁻¹
DSLV2	Death rate leaves due to self-shading	N	kg ha⁻¹ d⁻¹
DSLV3	Death rate leaves due to frost kill	N	kg ha⁻¹ d⁻¹
DSLV	Maximum of DSLV1, DSLV2, DSLV3	N	kg ha⁻¹ d⁻¹
DALV	Death rate leaves due to aging	N	kg ha⁻¹ d⁻¹
DRLV	Death rate leaves as a combination of DSLV and DALV	N	kg ha⁻¹ d⁻¹
SLAT	Specific leaf area for current time step, adjusted for source/sink limited leaf expansion rate	N	ha kg⁻¹
FYSAGE	Increase in physiological leaf age	N	-
GLAIEX	Sink-limited leaf expansion rate (exponential curve)	N	ha ha⁻¹ d⁻¹
GLASOL	Source-limited leaf expansion rate (biomass increase)	N	ha ha⁻¹ d⁻¹

External dependencies

Name	Description	Provided by	Unit
DVS	Crop development stage	DVS_Phenology	-
FL	Fraction biomass to leaves	DVS_Partitioning	-
FR	Fraction biomass to roots	DVS_Partitioning	-
SAI	Stem area index	WOFOST_Stem_Dynamics	-
PAI	Pod area index	WOFOST_Storage_Organ_Dynamics	-
TRA	Transpiration rate	Evapotranspiration	cm day⁻¹ ?
TRAMX	Maximum transpiration rate	Evapotranspiration	cm day⁻¹ ?
ADMI	Above-ground dry matter increase	CropSimulation	kg ha⁻¹ d⁻¹
RFTRA	Reduction factor for transpiration (water & oxygen)	Y	-
RF_FROST	Reduction factor frost kill	FROSTOL (optional)	-

Outputs

Name	Description	Pbl	Unit
LAI	Leaf area index, including stem and pod area	Y	-
TWLV	Dry weight of total leaves (living + dead)	Y	kg ha⁻¹

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
LAI	TDWI, SPAN, RGRLAI, TBASE, KDIFTB, SLATB
TWLV	TDWI, PERDL

[!NOTE] Notice that the following gradients are zero: - ∂SPAN/∂TWLV - ∂PERDL/∂TWLV - ∂KDIFTB/∂LAI

Methods:

initialize –

Initialize the WOFOST_Leaf_Dynamics simulation object.
calc_rates –

Calculate the rates of change for the leaf dynamics.
integrate –

Integrate the leaf dynamics state variables.

Attributes:

device –

Get device from ComputeConfig.
dtype –

Get dtype from ComputeConfig.

device `property`

device

Get device from ComputeConfig.

dtype `property`

dtype

Get dtype from ComputeConfig.

initialize

initialize(day: date, kiosk: VariableKiosk, parvalues: ParameterProvider) -> None

Initialize the WOFOST_Leaf_Dynamics simulation object.

Parameters:

day (date) –

The starting date of the simulation.
kiosk (VariableKiosk) –

A container for registering and publishing (internal and external) state variables. See PCSE documentation for details.
parvalues (ParameterProvider) –

A dictionary-like container holding all parameter sets (crop, soil, site) as key/value. The values are arrays or scalars. See PCSE documentation for details.

Source code in src/diffwofost/physical_models/crop/leaf_dynamics.py

def initialize(
    self, day: datetime.date, kiosk: VariableKiosk, parvalues: ParameterProvider
) -> None:
    """Initialize the WOFOST_Leaf_Dynamics simulation object.

    Args:
        day (datetime.date): The starting date of the simulation.
        kiosk (VariableKiosk): A container for registering and publishing
            (internal and external) state variables. See PCSE documentation for
            details.
        parvalues (ParameterProvider): A dictionary-like container holding
            all parameter sets (crop, soil, site) as key/value. The values are
            arrays or scalars. See PCSE documentation for details.
    """
    self.START_DATE = day
    self.kiosk = kiosk
    # TODO check if parvalues are already torch.nn.Parameters
    self.params = self.Parameters(parvalues)
    self.rates = self.RateVariables(kiosk)

    # Create scalar constants once at the beginning to avoid recreating them
    self._zero = torch.tensor(0.0, dtype=self.dtype, device=self.device)
    self._epsilon = torch.tensor(1e-12, dtype=self.dtype, device=self.device)
    self._sigmoid_sharpness = torch.tensor(1000, dtype=self.dtype, device=self.device)
    self._sigmoid_epsilon = torch.tensor(1e-14, dtype=self.dtype, device=self.device)

    # CALCULATE INITIAL STATE VARIABLES
    # check for required external variables
    _exist_required_external_variables(self.kiosk)
    # TODO check if external variables are already torch tensors

    # Get kiosk values and ensure they are on the correct device
    FL = torch.as_tensor(self.kiosk["FL"], dtype=self.dtype, device=self.device)
    FR = torch.as_tensor(self.kiosk["FR"], dtype=self.dtype, device=self.device)
    DVS = torch.as_tensor(self.kiosk["DVS"], dtype=self.dtype, device=self.device)

    params = self.params
    self.params_shape = _get_params_shape(params)

    # Initial leaf biomass
    TDWI = _broadcast_to(params.TDWI, self.params_shape, dtype=self.dtype, device=self.device)
    WLV = (TDWI * (1 - FR)) * FL
    DWLV = torch.zeros(self.params_shape, dtype=self.dtype, device=self.device)
    TWLV = WLV + DWLV

    # Initialize leaf classes (SLA, age and weight)
    SLA = torch.zeros((*self.params_shape, self.MAX_DAYS), dtype=self.dtype, device=self.device)
    LVAGE = torch.zeros(
        (*self.params_shape, self.MAX_DAYS), dtype=self.dtype, device=self.device
    )
    LV = torch.zeros((*self.params_shape, self.MAX_DAYS), dtype=self.dtype, device=self.device)
    SLA[..., 0] = params.SLATB(DVS).to(dtype=self.dtype, device=self.device)
    LV[..., 0] = WLV

    # Initial values for leaf area
    LAIEM = LV[..., 0] * SLA[..., 0]
    LASUM = LAIEM
    LAIEXP = LAIEM
    LAIMAX = LAIEM
    SAI = torch.as_tensor(self.kiosk["SAI"], dtype=self.dtype, device=self.device)
    PAI = torch.as_tensor(self.kiosk["PAI"], dtype=self.dtype, device=self.device)
    LAI = LASUM + SAI + PAI

    # Initialize StateVariables object
    self.states = self.StateVariables(
        kiosk,
        publish=["LAI", "TWLV", "WLV"],
        LV=LV,
        SLA=SLA,
        LVAGE=LVAGE,
        LAIEM=LAIEM,
        LASUM=LASUM,
        LAIEXP=LAIEXP,
        LAIMAX=LAIMAX,
        LAI=LAI,
        WLV=WLV,
        DWLV=DWLV,
        TWLV=TWLV,
    )

calc_rates

calc_rates(day: date, drv: WeatherDataContainer) -> None

Calculate the rates of change for the leaf dynamics.

Parameters:

day (date) –

The current date of the simulation.
drv (WeatherDataContainer) –

A dictionary-like container holding weather data elements as key/value. The values are arrays or scalars. See PCSE documentation for details.

Source code in src/diffwofost/physical_models/crop/leaf_dynamics.py

@prepare_rates
def calc_rates(self, day: datetime.date, drv: WeatherDataContainer) -> None:
    """Calculate the rates of change for the leaf dynamics.

    Args:
        day (datetime.date, optional): The current date of the simulation.
        drv (WeatherDataContainer, optional): A dictionary-like container holding
            weather data elements as key/value. The values are
            arrays or scalars. See PCSE documentation for details.
    """
    r = self.rates
    s = self.states
    p = self.params
    k = self.kiosk

    # If DVS < 0, the crop has not yet emerged, so we zerofy the rates using mask
    # A mask (0 if DVS < 0, 1 if DVS >= 0)
    DVS = torch.as_tensor(k["DVS"], dtype=self.dtype, device=self.device)
    dvs_mask = (DVS >= 0).to(dtype=self.dtype).to(device=self.device)

    # Growth rate leaves
    # weight of new leaves
    r.GRLV = dvs_mask * k.ADMI * k.FL

    # death of leaves due to water/oxygen stress
    RFTRA = _broadcast_to(k.RFTRA, self.params_shape, dtype=self.dtype, device=self.device)
    PERDL = _broadcast_to(p.PERDL, self.params_shape, dtype=self.dtype, device=self.device)
    r.DSLV1 = dvs_mask * s.WLV * (1.0 - RFTRA) * PERDL

    # death due to self shading cause by high LAI
    DVS = _broadcast_to(
        self.kiosk["DVS"], self.params_shape, dtype=self.dtype, device=self.device
    )
    KDIFTB = p.KDIFTB.to(device=self.device, dtype=self.dtype)
    LAICR = 3.2 / KDIFTB(DVS)
    r.DSLV2 = dvs_mask * s.WLV * torch.clamp(0.03 * (s.LAI - LAICR) / LAICR, 0.0, 0.03)

    # Death of leaves due to frost damage as determined by
    # Reduction Factor Frost "RF_FROST"
    if "RF_FROST" in self.kiosk:
        r.DSLV3 = s.WLV * k.RF_FROST
    else:
        r.DSLV3 = torch.zeros_like(s.WLV, dtype=self.dtype)

    r.DSLV3 = dvs_mask * r.DSLV3

    # leaf death equals maximum of water stress, shading and frost
    r.DSLV = torch.maximum(torch.maximum(r.DSLV1, r.DSLV2), r.DSLV3)
    r.DSLV = dvs_mask * r.DSLV

    # Determine how much leaf biomass classes have to die in states.LV,
    # given the a life span > SPAN, these classes will be accumulated
    # in DALV.
    # Note that the actual leaf death is imposed on the array LV during the
    # state integration step.
    tSPAN = _broadcast_to(
        p.SPAN, s.LVAGE.shape, dtype=self.dtype, device=self.device
    )  # Broadcast to same shape

    # Using a sigmoid here instead of a conditional statement on the value of
    # SPAN because the latter would not allow for the gradient to be tracked.
    # the if statement `p.SPAN.requires_grad` to avoid unnecessary
    # approximation when SPAN is not a learnable parameter.
    # here we use STE (straight through estimator) method.
    # TODO: sharpness can be exposed as a parameter
    if p.SPAN.requires_grad:
        # soft mask using sigmoid
        soft_mask = torch.sigmoid(
            (s.LVAGE - tSPAN - self._sigmoid_epsilon) / self._sigmoid_sharpness
        ).to(dtype=self.dtype)

        # originial hard mask
        hard_mask = (s.LVAGE > tSPAN).to(dtype=self.dtype)

        # STE method. Here detach is used to stop the gradient flow. This
        # way, during backpropagation, the gradient is computed only through
        # the `soft_mask``, while during the forward pass, the `hard_mask``
        # is used.
        span_mask = hard_mask.detach() + soft_mask - soft_mask.detach()
    else:
        span_mask = (s.LVAGE > tSPAN).to(dtype=self.dtype)

    r.DALV = torch.sum(span_mask * s.LV, dim=-1)
    r.DALV = dvs_mask * r.DALV

    # Total death rate leaves
    r.DRLV = torch.maximum(r.DSLV, r.DALV)

    # Get the temperature from the drv
    TEMP = _get_drv(drv.TEMP, self.params_shape, self.dtype, self.device)

    # physiologic ageing of leaves per time step
    TBASE = _broadcast_to(p.TBASE, self.params_shape, dtype=self.dtype, device=self.device)
    FYSAGE = (TEMP - TBASE) / (35.0 - TBASE)
    r.FYSAGE = dvs_mask * torch.clamp(FYSAGE, 0.0)

    # specific leaf area of leaves per time step
    SLATB = p.SLATB.to(device=self.device, dtype=self.dtype)
    r.SLAT = dvs_mask * SLATB(DVS)

    # leaf area not to exceed exponential growth curve
    is_lai_exp = s.LAIEXP < 6.0
    DTEFF = torch.clamp(TEMP - TBASE, 0.0)

    # NOTE: conditional statements do not allow for the gradient to be
    # tracked through the condition. Thus, the gradient with respect to
    # parameters that contribute to `is_lai_exp` (e.g. RGRLAI and TBASE)
    # are expected to be incorrect.
    RGRLAI = _broadcast_to(p.RGRLAI, self.params_shape, dtype=self.dtype, device=self.device)
    r.GLAIEX = torch.where(
        dvs_mask.bool(),
        torch.where(is_lai_exp, s.LAIEXP * RGRLAI * DTEFF, r.GLAIEX),
        self._zero,
    )

    # source-limited increase in leaf area
    r.GLASOL = torch.where(
        dvs_mask.bool(),
        torch.where(is_lai_exp, r.GRLV * r.SLAT, r.GLASOL),
        self._zero,
    )

    # sink-limited increase in leaf area
    GLA = torch.minimum(r.GLAIEX, r.GLASOL)

    # adjustment of specific leaf area of youngest leaf class
    r.SLAT = torch.where(
        dvs_mask.bool(),
        torch.where(
            is_lai_exp & (r.GRLV > self._epsilon), GLA / (r.GRLV + self._epsilon), r.SLAT
        ),
        self._zero,
    )

integrate

integrate(day: date, delt=1.0) -> None

Integrate the leaf dynamics state variables.

Parameters:

day (date) –

The current date of the simulation.
delt (float, default: 1.0 ) –

The time step for integration. Defaults to 1.0.

Source code in src/diffwofost/physical_models/crop/leaf_dynamics.py

@prepare_states
def integrate(self, day: datetime.date, delt=1.0) -> None:
    """Integrate the leaf dynamics state variables.

    Args:
        day (datetime.date, optional): The current date of the simulation.
        delt (float, optional): The time step for integration. Defaults to 1.0.
    """
    # TODO check if DVS < 0 and skip integration needed
    rates = self.rates
    states = self.states

    # --------- leave death ---------
    tLV = states.LV.clone()
    tSLA = states.SLA.clone()
    tLVAGE = states.LVAGE.clone()
    tDRLV = _broadcast_to(rates.DRLV, tLV.shape, dtype=self.dtype, device=self.device)

    # Leaf death is imposed on leaves from the oldest ones.
    # Calculate the cumulative sum of weights after leaf death, and
    # find out which leaf classes are dead (negative weights)
    weight_cumsum = tLV.cumsum(dim=-1) - tDRLV
    is_alive = weight_cumsum >= 0

    # Adjust value of oldest leaf class, i.e. the first non-zero
    # weight along the time axis (the last dimension).
    # Cast argument to int because torch.argmax requires it to be numeric
    idx_oldest = torch.argmax(is_alive.type(torch.int), dim=-1, keepdim=True).to(
        device=self.device
    )
    new_biomass = torch.take_along_dim(weight_cumsum, indices=idx_oldest, dim=-1)
    tLV = torch.scatter(tLV, dim=-1, index=idx_oldest, src=new_biomass)

    # Integration of physiological age
    # Zero out all dead leaf classes
    # NOTE: conditional statements do not allow for the gradient to be
    # tracked through the condition. Thus, the gradient with respect to
    # parameters that contribute to `is_alive` are expected to be incorrect.
    tLV = torch.where(is_alive, tLV, 0.0)
    tLVAGE = tLVAGE + rates.FYSAGE.unsqueeze(-1)
    tLVAGE = torch.where(is_alive, tLVAGE, 0.0)
    tSLA = torch.where(is_alive, tSLA, 0.0)

    # --------- leave growth ---------
    idx = int((day - self.START_DATE).days / delt)
    tLV[..., idx] = rates.GRLV
    tSLA[..., idx] = rates.SLAT
    tLVAGE[..., idx] = 0.0

    # calculation of new leaf area
    states.LASUM = torch.sum(tLV * tSLA, dim=-1)
    states.LAI = self._calc_LAI()
    states.LAIMAX = torch.maximum(states.LAI, states.LAIMAX)

    # exponential growth curve
    states.LAIEXP = states.LAIEXP + rates.GLAIEX

    # Update leaf biomass states
    states.WLV = torch.sum(tLV, dim=-1)
    states.DWLV = states.DWLV + rates.DRLV
    states.TWLV = states.WLV + states.DWLV

    # Store final leaf biomass deques
    self.states.LV = tLV
    self.states.SLA = tSLA
    self.states.LVAGE = tLVAGE

diffwofost.physical_models.crop.root_dynamics.WOFOST_Root_Dynamics

Bases: SimulationObject

Root biomass dynamics and rooting depth.

Root growth and root biomass dynamics in WOFOST are separate processes, with the only exception that root growth stops when no more biomass is sent to the root system.

Root biomass increase results from the assimilates partitioned to the root system. Root death is defined as the current root biomass multiplied by a relative death rate (RDRRTB). The latter as a function of the development stage (DVS).

Increase in root depth is a simple linear expansion over time until the maximum rooting depth (RDM) is reached.

Simulation parameters

Name	Description	Type	Unit
RDI	Initial rooting depth	SCr	cm
RRI	Daily increase in rooting depth	SCr	cm day⁻¹
RDMCR	Maximum rooting depth of the crop	SCR	cm
RDMSOL	Maximum rooting depth of the soil	SSo	cm
TDWI	Initial total crop dry weight	SCr	kg ha⁻¹
IAIRDU	Presence of air ducts in the root (1) or not (0)	SCr	-
RDRRTB	Relative death rate of roots as a function of development stage	TCr	-

State variables

Name	Description	Pbl	Unit
RD	Current rooting depth	Y	cm
RDM	Maximum attainable rooting depth at the minimum of the soil and crop maximum rooting depth	N	cm
WRT	Weight of living roots	Y	kg ha⁻¹
DWRT	Weight of dead roots	N	kg ha⁻¹
TWRT	Total weight of roots	Y	kg ha⁻¹

Rate variables

Name	Description	Pbl	Unit
RR	Growth rate root depth	N	cm
GRRT	Growth rate root biomass	N	kg ha⁻¹ d⁻¹
DRRT	Death rate root biomass	N	kg ha⁻¹ d⁻¹
GWRT	Net change in root biomass	N	kg ha⁻¹ d⁻¹

Signals send or handled

None

External dependencies:

Name	Description	Provided by	Unit
DVS	Crop development stage	DVS_Phenology	-
DMI	Total dry matter increase	CropSimulation	kg ha⁻¹ d⁻¹
FR	Fraction biomass to roots	DVS_Partitioning	-

Outputs:

Name	Description	Provided by	Unit
RD	Current rooting depth	Y	cm
TWRT	Total weight of roots	Y	kg ha⁻¹

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
RD	RDI, RRI, RDMCR, RDMSOL
TWRT	TDWI, RDRRTB

[!NOTE] Notice that the gradient ∂TWRT/∂RDRRTB is zero.

IMPORTANT NOTICE

Currently root development is linear and depends only on the fraction of assimilates send to the roots (FR) and not on the amount of assimilates itself. This means that roots also grow through the winter when there is no assimilation due to low temperatures. There has been a discussion to change this behaviour and make root growth dependent on the assimilates send to the roots: so root growth stops when there are no assimilates available for growth.

Finally, we decided not to change the root model and keep the original WOFOST approach because of the following reasons: - A dry top layer in the soil could create a large drought stress that reduces the assimilates to zero. In this situation the roots would not grow if dependent on the assimilates, while water is available in the zone just below the root zone. Therefore a dependency on the amount of assimilates could create model instability in dry conditions (e.g. Southern-Mediterranean, etc.). - Other solutions to alleviate the problem above were explored: only put this limitation after a certain development stage, putting a dependency on soil moisture levels in the unrooted soil compartment. All these solutions were found to introduce arbitrary parameters that have no clear explanation. Therefore all proposed solutions were discarded.

We conclude that our current knowledge on root development is insufficient to propose a better and more biophysical approach to root development in WOFOST.

Methods:

initialize –

Initialize the model.
calc_rates –

Calculate the rates of change of the state variables.
integrate –

Integrate the state variables using the rates of change.

Attributes:

device –

Get device from ComputeConfig.
dtype –

Get dtype from ComputeConfig.

device `property`

device

Get device from ComputeConfig.

dtype `property`

dtype

Get dtype from ComputeConfig.

initialize

initialize(day: date, kiosk: VariableKiosk, parvalues: ParameterProvider) -> None

Initialize the model.

Parameters:

day (date) –

The starting date of the simulation.
kiosk (VariableKiosk) –

A container for registering and publishing (internal and external) state variables. See PCSE documentation for details.
parvalues (ParameterProvider) –

A dictionary-like container holding all parameter sets (crop, soil, site) as key/value. The values are arrays or scalars. See PCSE documentation for details.

Source code in src/diffwofost/physical_models/crop/root_dynamics.py

def initialize(
    self, day: datetime.date, kiosk: VariableKiosk, parvalues: ParameterProvider
) -> None:
    """Initialize the model.

    Args:
        day (datetime.date): The starting date of the simulation.
        kiosk (VariableKiosk): A container for registering and publishing
            (internal and external) state variables. See PCSE documentation for
            details.
        parvalues (ParameterProvider): A dictionary-like container holding
            all parameter sets (crop, soil, site) as key/value. The values are
            arrays or scalars. See PCSE documentation for details.
    """
    self.kiosk = kiosk
    self.params = self.Parameters(parvalues)
    self.rates = self.RateVariables(kiosk, publish=["DRRT", "GRRT"])

    # INITIAL STATES
    params = self.params
    self.params_shape = _get_params_shape(params)
    shape = self.params_shape

    # Initial root depth states
    RDI = _broadcast_to(params.RDI, shape, dtype=self.dtype, device=self.device)
    RDMCR = _broadcast_to(params.RDMCR, shape, dtype=self.dtype, device=self.device)
    RDMSOL = _broadcast_to(params.RDMSOL, shape, dtype=self.dtype, device=self.device)

    rdmax = torch.maximum(RDI, torch.minimum(RDMCR, RDMSOL))
    RDM = rdmax
    RD = RDI

    # Initial root biomass states
    TDWI = _broadcast_to(params.TDWI, shape, dtype=self.dtype, device=self.device)
    FR = _broadcast_to(self.kiosk["FR"], shape, dtype=self.dtype, device=self.device)
    WRT = TDWI * FR
    DWRT = torch.zeros(shape, dtype=self.dtype, device=self.device)
    TWRT = WRT + DWRT

    self.states = self.StateVariables(
        kiosk, publish=["RD", "WRT", "TWRT"], RD=RD, RDM=RDM, WRT=WRT, DWRT=DWRT, TWRT=TWRT
    )

calc_rates

calc_rates(day: date = None, drv: WeatherDataContainer = None) -> None

Calculate the rates of change of the state variables.

Parameters:

day (date, default: None ) –

The current date of the simulation.
drv (WeatherDataContainer, default: None ) –

A dictionary-like container holding weather data elements as key/value. The values are arrays or scalars. See PCSE documentation for details.

Source code in src/diffwofost/physical_models/crop/root_dynamics.py

@prepare_rates
def calc_rates(self, day: datetime.date = None, drv: WeatherDataContainer = None) -> None:
    """Calculate the rates of change of the state variables.

    Args:
        day (datetime.date, optional): The current date of the simulation.
        drv (WeatherDataContainer, optional): A dictionary-like container holding
            weather data elements as key/value. The values are
            arrays or scalars. See PCSE documentation for details.
    """
    p = self.params
    r = self.rates
    s = self.states
    k = self.kiosk

    # If DVS < 0, the crop has not yet emerged, so we zerofy the rates using mask.
    # Make a mask (0 if DVS < 0, 1 if DVS >= 0)
    DVS = _broadcast_to(k["DVS"], self.params_shape, dtype=self.dtype, device=self.device)
    dvs_mask = (DVS >= 0).to(dtype=self.dtype)

    # Increase in root biomass
    FR = _broadcast_to(k["FR"], self.params_shape, dtype=self.dtype, device=self.device)
    DMI = _broadcast_to(k["DMI"], self.params_shape, dtype=self.dtype, device=self.device)
    RDRRTB = p.RDRRTB.to(device=self.device, dtype=self.dtype)

    r.GRRT = dvs_mask * FR * DMI
    r.DRRT = dvs_mask * s.WRT * RDRRTB(DVS)
    r.GWRT = r.GRRT - r.DRRT

    # Increase in root depth
    RRI = _broadcast_to(p.RRI, self.params_shape, dtype=self.dtype, device=self.device)
    r.RR = dvs_mask * torch.minimum((s.RDM - s.RD), RRI)

    # Do not let the roots growth if partioning to the roots
    # (variable FR) is zero.
    mask = (FR > 0.0).to(dtype=self.dtype)
    r.RR = r.RR * mask * dvs_mask

integrate

integrate(day: date = None, delt=1.0) -> None

Integrate the state variables using the rates of change.

Parameters:

day (date, default: None ) –

The current date of the simulation.
delt (float, default: 1.0 ) –

The time step for integration. Defaults to 1.0.

Source code in src/diffwofost/physical_models/crop/root_dynamics.py

@prepare_states
def integrate(self, day: datetime.date = None, delt=1.0) -> None:
    """Integrate the state variables using the rates of change.

    Args:
        day (datetime.date, optional): The current date of the simulation.
        delt (float, optional): The time step for integration. Defaults to 1.0.
    """
    rates = self.rates
    states = self.states

    # Dry weight of living roots
    states.WRT = states.WRT + rates.GWRT

    # Dry weight of dead roots
    states.DWRT = states.DWRT + rates.DRRT

    # Total weight dry + living roots
    states.TWRT = states.WRT + states.DWRT

    # New root depth
    states.RD = states.RD + rates.RR

Utility (under development)

diffwofost.physical_models.config.Configuration `dataclass`

Configuration(CROP: type[SimulationObject], SOIL: type[SimulationObject] | None = None, AGROMANAGEMENT: type[AncillaryObject] = AgroManager, OUTPUT_VARS: list = list(), SUMMARY_OUTPUT_VARS: list = list(), TERMINAL_OUTPUT_VARS: list = list(), OUTPUT_INTERVAL: str = 'daily', OUTPUT_INTERVAL_DAYS: int = 1, OUTPUT_WEEKDAY: int = 0, model_config_file: str | Path | None = None, description: str | None = None)

Class to store model configuration from a PCSE configuration files.

Methods:

from_pcse_config_file –

Load the model configuration from a PCSE configuration file.
update_output_variable_lists –

Updates the lists of output variables that are defined in the configuration file.

from_pcse_config_file `classmethod`

from_pcse_config_file(filename: str | Path) -> Self

Load the model configuration from a PCSE configuration file.

Parameters:

filename (str | Path) –

Path to the configuraiton file. The path is first interpreted with respect to the current working directory and, if not found, it will then be interpreted with respect to the conf folder in the PCSE package.

Returns:

Configuration ( Self ) –

Model configuration instance

Raises:

FileNotFoundError –

if the configuraiton file does not exist
RuntimeError –

if parsing the configuration file fails

Source code in src/diffwofost/physical_models/config.py

@classmethod
def from_pcse_config_file(cls, filename: str | Path) -> Self:
    """Load the model configuration from a PCSE configuration file.

    Args:
        filename (str | pathlib.Path): Path to the configuraiton file. The path is first
            interpreted with respect to the current working directory and, if not found, it will
            then be interpreted with respect to the `conf` folder in the PCSE package.

    Returns:
        Configuration: Model configuration instance

    Raises:
        FileNotFoundError: if the configuraiton file does not exist
        RuntimeError: if parsing the configuration file fails
    """
    config = {}

    path = Path(filename)
    if path.is_absolute() or path.is_file():
        model_config_file = path
    else:
        pcse_dir = Path(pcse.__path__[0])
        model_config_file = pcse_dir / "conf" / path
    model_config_file = model_config_file.resolve()

    # check that configuration file exists
    if not model_config_file.exists():
        msg = f"PCSE model configuration file does not exist: {model_config_file.name}"
        raise FileNotFoundError(msg)
    # store for later use
    config["model_config_file"] = model_config_file

    # Load file using execfile
    try:
        loc = {}
        bytecode = compile(open(model_config_file).read(), model_config_file, "exec")
        exec(bytecode, {}, loc)
    except Exception as e:
        msg = f"Failed to load configuration from file {model_config_file}"
        raise RuntimeError(msg) from e

    # Add the descriptive header for later use
    if "__doc__" in loc:
        desc = loc.pop("__doc__")
        if len(desc) > 0:
            description = desc
            if description[-1] != "\n":
                description += "\n"
        config["description"] = description

    # Loop through the attributes in the configuration file
    for key, value in loc.items():
        if key.isupper():
            config[key] = value
    return cls(**config)

update_output_variable_lists

update_output_variable_lists(output_vars: str | list | tuple | set | None = None, summary_vars: str | list | tuple | set | None = None, terminal_vars: str | list | tuple | set | None = None)

Updates the lists of output variables that are defined in the configuration file.

This is useful because sometimes you want the flexibility to get access to an additional model variable which is not in the standard list of variables defined in the model configuration file. The more elegant way is to define your own configuration file, but this adds some flexibility particularly for use in jupyter notebooks and exploratory analysis.

Note that there is a different behaviour given the type of the variable provided. List and string inputs will extend the list of variables, while set/tuple inputs will replace the current list.

Parameters:

output_vars (str | list | tuple | set | None, default: None ) –

the variable names to add/replace for the OUTPUT_VARS configuration variable
summary_vars (str | list | tuple | set | None, default: None ) –

the variable names to add/replace for the SUMMARY_OUTPUT_VARS configuration variable
terminal_vars (str | list | tuple | set | None, default: None ) –

the variable names to add/replace for the TERMINAL_OUTPUT_VARS configuration variable

Raises:

TypeError –

if the type of the input arguments is not recognized

Source code in src/diffwofost/physical_models/config.py

def update_output_variable_lists(
    self,
    output_vars: str | list | tuple | set | None = None,
    summary_vars: str | list | tuple | set | None = None,
    terminal_vars: str | list | tuple | set | None = None,
):
    """Updates the lists of output variables that are defined in the configuration file.

    This is useful because sometimes you want the flexibility to get access to an additional
    model variable which is not in the standard list of variables defined in the model
    configuration file. The more elegant way is to define your own configuration file, but this
    adds some flexibility particularly for use in jupyter notebooks and exploratory analysis.

    Note that there is a different behaviour given the type of the variable provided. List and
    string inputs will extend the list of variables, while set/tuple inputs will replace the
    current list.

    Args:
        output_vars: the variable names to add/replace for the OUTPUT_VARS configuration
            variable
        summary_vars: the variable names to add/replace for the SUMMARY_OUTPUT_VARS
            configuration variable
        terminal_vars: the variable names to add/replace for the TERMINAL_OUTPUT_VARS
            configuration variable

    Raises:
        TypeError: if the type of the input arguments is not recognized
    """
    config_varnames = ["OUTPUT_VARS", "SUMMARY_OUTPUT_VARS", "TERMINAL_OUTPUT_VARS"]
    for varitems, config_varname in zip(
        [output_vars, summary_vars, terminal_vars], config_varnames, strict=True
    ):
        if varitems is None:
            continue
        else:
            if isinstance(varitems, str):  # A string: we extend the current list
                getattr(self, config_varname).extend(varitems.split())
            elif isinstance(varitems, list):  # a list: we extend the current list
                getattr(self, config_varname).extend(varitems)
            elif isinstance(varitems, tuple | set):  # tuple/set we replace the current list
                attr = getattr(self, config_varname)
                attr.clear()
                attr.extend(list(varitems))
            else:
                msg = f"Unrecognized input for `output_vars` to engine(): {output_vars}"
                raise TypeError(msg)

diffwofost.physical_models.config.ComputeConfig

Central configuration for device and dtype settings.

This class provides a centralized way to control PyTorch device and dtype settings across all simulation objects in diffWOFOST. Instead of setting device and dtype individually for each class, use this central configuration to apply settings globally.

Default Behavior:

Device: Automatically defaults to 'cuda' if available, otherwise 'cpu'
Dtype: Defaults to torch.float64

Basic Usage:

>>> from diffwofost.physical_models.config import ComputeConfig
>>> import torch
>>>
>>> # Set device to CPU
>>> ComputeConfig.set_device('cpu')
>>>
>>> # Or use a torch.device object
>>> ComputeConfig.set_device(torch.device('cuda'))
>>>
>>> # Set dtype to float32
>>> ComputeConfig.set_dtype(torch.float32)
>>>
>>> # Get current settings
>>> device = ComputeConfig.get_device()  # Returns: torch.device('cpu')
>>> dtype = ComputeConfig.get_dtype()    # Returns: torch.float32

Using with Simulation Objects:

All simulation objects (e.g., WOFOST_Leaf_Dynamics, WOFOST_Phenology) automatically use the settings from ComputeConfig. No changes needed to instantiation code:

>>> from diffwofost.physical_models.config import ComputeConfig
>>> from diffwofost.physical_models.crop.leaf_dynamics import WOFOST_Leaf_Dynamics
>>>
>>> # Set global compute settings
>>> ComputeConfig.set_device('cuda')
>>> ComputeConfig.set_dtype(torch.float32)
>>>
>>> # Instantiate objects - they automatically use global settings
>>> leaf_dynamics = WOFOST_Leaf_Dynamics()

Switching Between Devices:

Useful for switching between GPU training and CPU evaluation:

>>> # Train on GPU
>>> ComputeConfig.set_device('cuda')
>>> ComputeConfig.set_dtype(torch.float32)
>>> # ... run training ...
>>>
>>> # Evaluate on CPU
>>> ComputeConfig.set_device('cpu')
>>> ComputeConfig.set_dtype(torch.float64)
>>> # ... run evaluation ...

Resetting to Defaults:

>>> ComputeConfig.reset_to_defaults()

Methods:

get_device –

Get the current device setting.
set_device –

Set the device to use for tensor operations.
get_dtype –

Get the current dtype setting.
set_dtype –

Set the dtype to use for tensor creation.
reset_to_defaults –

Reset device and dtype to their default values.

get_device `classmethod`

get_device() -> device

Get the current device setting.

Returns:

device –

torch.device: The current device (cuda or cpu)

Source code in src/diffwofost/physical_models/config.py

@classmethod
def get_device(cls) -> torch.device:
    """Get the current device setting.

    Returns:
        torch.device: The current device (cuda or cpu)
    """
    cls._initialize_defaults()
    return cls._device

set_device `classmethod`

set_device(device: str | device) -> None

Set the device to use for tensor operations.

Parameters:

device (str | device) –

Device to use ('cuda', 'cpu', or torch.device object)

Example

ComputeConfig.set_device('cuda') ComputeConfig.set_device(torch.device('cpu'))

Source code in src/diffwofost/physical_models/config.py

@classmethod
def set_device(cls, device: str | torch.device) -> None:
    """Set the device to use for tensor operations.

    Args:
        device (str | torch.device): Device to use ('cuda', 'cpu', or torch.device object)

    Example:
        >>> ComputeConfig.set_device('cuda')
        >>> ComputeConfig.set_device(torch.device('cpu'))
    """
    if isinstance(device, str):
        cls._device = torch.device(device)
    else:
        cls._device = device

get_dtype `classmethod`

get_dtype() -> dtype

Get the current dtype setting.

Returns:

dtype –

torch.dtype: The current dtype (e.g., torch.float32, torch.float64)

Source code in src/diffwofost/physical_models/config.py

@classmethod
def get_dtype(cls) -> torch.dtype:
    """Get the current dtype setting.

    Returns:
        torch.dtype: The current dtype (e.g., torch.float32, torch.float64)
    """
    cls._initialize_defaults()
    return cls._dtype

set_dtype `classmethod`

set_dtype(dtype: dtype) -> None

Set the dtype to use for tensor creation.

Parameters:

dtype (dtype) –

PyTorch dtype (torch.float32, torch.float64, etc.)

Example

ComputeConfig.set_dtype(torch.float32)

Source code in src/diffwofost/physical_models/config.py

@classmethod
def set_dtype(cls, dtype: torch.dtype) -> None:
    """Set the dtype to use for tensor creation.

    Args:
        dtype (torch.dtype): PyTorch dtype (torch.float32, torch.float64, etc.)

    Example:
        >>> ComputeConfig.set_dtype(torch.float32)
    """
    cls._dtype = dtype

reset_to_defaults `classmethod`

reset_to_defaults() -> None

Reset device and dtype to their default values.

Source code in src/diffwofost/physical_models/config.py

@classmethod
def reset_to_defaults(cls) -> None:
    """Reset device and dtype to their default values."""
    cls._device = None
    cls._dtype = None
    cls._initialize_defaults()

diffwofost.physical_models.engine.Engine

Engine(parameterprovider, weatherdataprovider, agromanagement, config: str | Path | Configuration)

Bases: Engine

Source code in src/diffwofost/physical_models/engine.py

def __init__(
    self,
    parameterprovider,
    weatherdataprovider,
    agromanagement,
    config: str | Path | Configuration,
):
    BaseEngine.__init__(self)

    # If a path is given, load the model configuration from a PCSE config file
    if isinstance(config, str | Path):
        self.mconf = Configuration.from_pcse_config_file(config)
    else:
        self.mconf = config

    self.parameterprovider = parameterprovider

    # Variable kiosk for registering and publishing variables
    self.kiosk = VariableKiosk()

    # Placeholder for variables to be saved during a model run
    self._saved_output = []
    self._saved_summary_output = []
    self._saved_terminal_output = {}

    # register handlers for starting/finishing the crop simulation, for
    # handling output and terminating the system
    self._connect_signal(self._on_CROP_START, signal=signals.crop_start)
    self._connect_signal(self._on_CROP_FINISH, signal=signals.crop_finish)
    self._connect_signal(self._on_OUTPUT, signal=signals.output)
    self._connect_signal(self._on_TERMINATE, signal=signals.terminate)

    # Component for agromanagement
    self.agromanager = self.mconf.AGROMANAGEMENT(self.kiosk, agromanagement)
    start_date = self.agromanager.start_date
    end_date = self.agromanager.end_date

    # Timer: starting day, final day and model output
    self.timer = Timer(self.kiosk, start_date, end_date, self.mconf)
    self.day, _ = self.timer()

    # Driving variables
    self.weatherdataprovider = weatherdataprovider
    self.drv = self._get_driving_variables(self.day)

    # Component for simulation of soil processes
    if self.mconf.SOIL is not None:
        self.soil = self.mconf.SOIL(self.day, self.kiosk, parameterprovider)

    # Call AgroManagement module for management actions at initialization
    self.agromanager(self.day, self.drv)

    # Calculate initial rates
    self.calc_rates(self.day, self.drv)

diffwofost.physical_models.utils.EngineTestHelper

EngineTestHelper(parameterprovider, weatherdataprovider, agromanagement, config, external_states=None, device=None, dtype=None)

Bases: Engine

An engine which is purely for running the YAML unit tests.

Source code in src/diffwofost/physical_models/utils.py

def __init__(
    self,
    parameterprovider,
    weatherdataprovider,
    agromanagement,
    config,
    external_states=None,
    device=None,
    dtype=None,
):
    BaseEngine.__init__(self)

    # If a path is given, load the model configuration from a PCSE config file
    if isinstance(config, str | Path):
        self.mconf = Configuration.from_pcse_config_file(config)
    else:
        self.mconf = config

    self.parameterprovider = parameterprovider

    # Configure device and dtype on crop module class if it supports them
    if hasattr(self.mconf.CROP, "device") and device is not None:
        self.mconf.CROP.device = device
    if hasattr(self.mconf.CROP, "dtype") and dtype is not None:
        self.mconf.CROP.dtype = dtype

    # Variable kiosk for registering and publishing variables
    self.kiosk = VariableKioskTestHelper(external_states)

    # Placeholder for variables to be saved during a model run
    self._saved_output = list()
    self._saved_summary_output = list()
    self._saved_terminal_output = dict()

    # register handlers for starting/finishing the crop simulation, for
    # handling output and terminating the system
    self._connect_signal(self._on_CROP_START, signal=signals.crop_start)
    self._connect_signal(self._on_CROP_FINISH, signal=signals.crop_finish)
    self._connect_signal(self._on_OUTPUT, signal=signals.output)
    self._connect_signal(self._on_TERMINATE, signal=signals.terminate)

    # Component for agromanagement
    self.agromanager = self.mconf.AGROMANAGEMENT(self.kiosk, agromanagement)
    start_date = self.agromanager.start_date
    end_date = self.agromanager.end_date

    # Timer: starting day, final day and model output
    self.timer = Timer(self.kiosk, start_date, end_date, self.mconf)
    self.day, delt = self.timer()
    # Update external states in the kiosk
    self.kiosk(self.day)

    # Driving variables
    self.weatherdataprovider = weatherdataprovider
    self.drv = self._get_driving_variables(self.day)

    # Component for simulation of soil processes
    if self.mconf.SOIL is not None:
        self.soil = self.mconf.SOIL(self.day, self.kiosk, parameterprovider)

    # Call AgroManagement module for management actions at initialization
    self.agromanager(self.day, self.drv)

    # Calculate initial rates
    self.calc_rates(self.day, self.drv)

Crop modules

diffwofost.physical_models.crop.phenology.DVS_Phenology

device property

dtype property

initialize

calc_rates

integrate

get_variable

diffwofost.physical_models.crop.partitioning.DVS_Partitioning

initialize

integrate

calc_rates

diffwofost.physical_models.crop.assimilation.WOFOST72_Assimilation

device property

dtype property

initialize

calc_rates

__call__

integrate

diffwofost.physical_models.crop.leaf_dynamics.WOFOST_Leaf_Dynamics

device property

dtype property

initialize

calc_rates

integrate

diffwofost.physical_models.crop.root_dynamics.WOFOST_Root_Dynamics

device property

dtype property

initialize

calc_rates

integrate

Utility (under development)

diffwofost.physical_models.config.Configuration dataclass

from_pcse_config_file classmethod

update_output_variable_lists

diffwofost.physical_models.config.ComputeConfig

get_device classmethod

set_device classmethod

get_dtype classmethod

set_dtype classmethod

reset_to_defaults classmethod

diffwofost.physical_models.engine.Engine

diffwofost.physical_models.utils.EngineTestHelper

device `property`

dtype `property`

device `property`

dtype `property`

call

device `property`

dtype `property`

device `property`

dtype `property`

diffwofost.physical_models.config.Configuration `dataclass`

from_pcse_config_file `classmethod`

get_device `classmethod`

set_device `classmethod`

get_dtype `classmethod`

set_dtype `classmethod`

reset_to_defaults `classmethod`