Crop modules

diffwofost.physical_models.crop.wofost72.Wofost72

Bases: SimulationObject

Top level object organizing the different components of WOFOST.

The CropSimulation object organizes the different processes of the crop simulation. Moreover, it contains the parameters, rate and state variables which are relevant at the level of the entire crop. The processes that are implemented as embedded simulation objects consist of:

1. Phenology (self.pheno)
2. Partitioning (self.part)
3. Assimilation (self.assim)
4. Maintenance respiration (self.mres)
5. Evapotranspiration (self.evtra)
6. Leaf dynamics (self.lv_dynamics)
7. Stem dynamics (self.st_dynamics)
8. Root dynamics (self.ro_dynamics)
9. Storage organ dynamics (self.so_dynamics)

Simulation parameters:

Name	Description	Type	Unit
CVL	Conversion factor for assimilates to leaves	SCr	-
CVO	Conversion factor for assimilates to storage organs	SCr	-
CVR	Conversion factor for assimilates to roots	SCr	-
CVS	Conversion factor for assimilates to stems	SCr	-

State variables:

Name	Description	Type	Unit
TAGP	Total above-ground Production	N	kg ha-1
GASST	Total gross assimilation	N	kg CH2O ha-1
MREST	Total gross maintenance respiration	N	kg CH2O ha-1
CTRAT	Total crop transpiration accumulated over the crop cycle	N	cm
CEVST	Total soil evaporation accumulated over the crop cycle	N	cm
HI	Harvest Index (only calculated during finalize())	N	-
DOF	Date representing the day of finish of the crop simulation	N	-
FINISH_TYPE	String representing the reason for finishing the simulation:	N	-
	maturity, harvest, leave death, etc.

Rate variables:

Name	Description	Type	Unit
GASS	Assimilation rate corrected for water stress	N	kg CH2O ha-1 d-1
MRES	Actual maintenance respiration rate, ...	N	kg CH2O ha-1 d-1
	... taking into account that MRES <= GASS
ASRC	Net available assimilates (GASS - MRES)	N	kg CH2O ha-1 d-1
DMI	Total dry matter increase, ...	Y	kg ha-1 d-1
	... calculated as ASRC times a weighted conversion efficiency
ADMI	Aboveground dry matter increase	Y	kg ha-1 d-1

Methods:

• initialize –

  Initialize the crop simulation and its embedded components.

• calc_rates –

  Calculate the rates of change of the state variables.

• integrate –

  Integrate the state variables using the rates of change.

• finalize –

  Finalize the crop simulation by computing the Harvest Index.

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, shape: tuple | Size | None = None, component_overrides: dict | None = None) -> None

Initialize the crop simulation and its embedded components.

Parameters:

• day (date) –

  Start date of the simulation.

• kiosk (VariableKiosk) –

  Variable kiosk used to read and publish crop state.

• parvalues (ParameterProvider) –

  Parameter provider containing the physical-model parameters for the crop.

• shape (tuple | Size | None, default: None ) –

  Target tensor shape for state and rate variables.

• component_overrides (dict | None, default: None ) –

  Optional mapping used to replace one or more internal WOFOST components at construction time. The component_overrides is a dictionary containing: - "component_class": The class to use for the component. - "model": The model to use for the component, if specified in the override. - "kwargs": Any additional keyword arguments to pass to the component constructor.

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

def initialize(
    self,
    day: datetime.date,
    kiosk: VariableKiosk,
    parvalues: ParameterProvider,
    shape: tuple | torch.Size | None = None,
    component_overrides: dict | None = None,
) -> None:
    """Initialize the crop simulation and its embedded components.

    Args:
        day: Start date of the simulation.
        kiosk: Variable kiosk used to read and publish crop state.
        parvalues: Parameter provider containing the physical-model
            parameters for the crop.
        shape: Target tensor shape for state and rate variables.
        component_overrides: Optional mapping used to replace one or more
            internal WOFOST components at construction time.
            The ``component_overrides`` is a dictionary containing:
            - "component_class": The class to use for the component.
            - "model": The model to use for the component, if specified in the override.
            - "kwargs": Any additional keyword arguments to pass to the component
            constructor.
    """
    self.params = self.Parameters(parvalues, shape=shape)
    self.rates = self.RateVariables(
        kiosk, publish=["DMI", "ADMI", "REALLOC_LV", "REALLOC_ST", "REALLOC_SO"], shape=shape
    )
    self.kiosk = kiosk

    # Initialize components of the crop
    # This will add attributes to self for each component, e.g. self.pheno, self.part, etc.
    initialize_all_components(
        self,
        day,
        kiosk,
        parvalues,
        shape=shape,
        component_overrides=component_overrides,
    )

    # Initial total (living+dead) above-ground biomass of the crop
    TAGP = self.kiosk.TWLV + self.kiosk.TWST + self.kiosk.TWSO
    self.states = self.StateVariables(
        kiosk,
        publish=["TAGP", "GASST", "MREST", "HI"],
        TAGP=TAGP,
        GASST=0.0,
        MREST=0.0,
        CTRAT=0.0,
        CEVST=0.0,
        HI=0.0,
        DOF=None,
        FINISH_TYPE=None,
        shape=shape,
    )

    # Check partitioning of TDWI over plant organs
    checksum = parvalues["TDWI"] - self.states.TAGP - self.kiosk["TWRT"]
    if torch.any(torch.abs(checksum) > 0.0001):
        msg = "Error in partitioning of initial biomass (TDWI)!"
        raise exc.PartitioningError(msg)

    # assign handler for CROP_FINISH signal
    self._connect_signal(self._on_CROP_FINISH, signal=signals.crop_finish)

calc_rates

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

Calculate the rates of change of the state variables.

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/wofost72.py

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

    Args:
        day (datetime.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.
    """
    p = self.params
    r = self.rates
    k = self.kiosk

    # Phenology
    self.pheno.calc_rates(day, drv)

    # if before emergence there is no need to continue
    # because only the phenology is running.
    # TODO: revisit this when fixing #60
    if torch.all(self.pheno.states.STAGE == 0):
        self.evtra(day, drv)
        return

    # Potential assimilation
    PGASS = self.assim(day, drv)

    # (evapo)transpiration rates
    self.evtra(day, drv)

    # water stress reduction
    r.GASS = PGASS * k.RFTRA

    # Respiration
    PMRES = self.mres(day, drv)
    r.MRES = torch.minimum(r.GASS, PMRES)

    # Net available assimilates
    r.ASRC = r.GASS - r.MRES

    # DM partitioning factors (pf), conversion factor (CVF),
    # dry matter increase (DMI) and check on carbon balance
    pf = self.part.calc_rates(day, drv)
    CVF = 1.0 / (
        (pf.FL / p.CVL + pf.FS / p.CVS + pf.FO / p.CVO) * (1.0 - pf.FR) + pf.FR / p.CVR
    )
    r.DMI = CVF * r.ASRC
    self._check_carbon_balance(day, r.DMI, r.GASS, r.MRES, CVF, pf)

    # distribution over plant organ

    # Reallocation from stems/leaves not applicable in WOFOST72
    r.REALLOC_LV = torch.zeros_like(r.DMI)
    r.REALLOC_ST = torch.zeros_like(r.DMI)
    r.REALLOC_SO = torch.zeros_like(r.DMI)

    # Below-ground dry matter increase and root dynamics
    self.ro_dynamics.calc_rates(day, drv)
    # Aboveground dry matter increase and distribution over stems,
    # leaves, organs
    r.ADMI = (1.0 - pf.FR) * r.DMI
    self.st_dynamics.calc_rates(day, drv)
    self.so_dynamics.calc_rates(day, drv)
    self.lv_dynamics.calc_rates(day, drv)

integrate

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

Integrate the state variables using the rates of change.

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/wofost72.py

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

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

    # Capture stage *before* integration (phenology will advance it).
    # STAGE == 0 means "emerging"; read directly from the tensor to avoid
    # a GPU→CPU sync from get_variable() / .item().
    crop_stage_before = self.pheno.states.STAGE.clone()

    # Phenology
    self.pheno.integrate(day, delt)

    # if before emergence there is no need to continue
    # because only the phenology is running.
    # Just run a touch() to to ensure that all state variables are available
    # in the kiosk
    # TODO: revisit this when fixing #60
    if torch.all(crop_stage_before == 0):
        self.touch()
        return

    # Partitioning
    self.part.integrate(day, delt)

    # Integrate states on leaves, storage organs, stems and roots
    self.ro_dynamics.integrate(day, delt)
    self.so_dynamics.integrate(day, delt)
    self.st_dynamics.integrate(day, delt)
    self.lv_dynamics.integrate(day, delt)

    # Integrate total (living+dead) above-ground biomass of the crop
    states.TAGP = self.kiosk.TWLV + self.kiosk.TWST + self.kiosk.TWSO

    # total gross assimilation and maintenance respiration
    states.GASST = states.GASST + rates.GASS
    states.MREST = states.MREST + rates.MRES

    # total crop transpiration and soil evaporation
    states.CTRAT = states.CTRAT + self.kiosk.TRA
    states.CEVST = states.CEVST + self.kiosk.EVS

finalize

finalize(day: date) -> None

Finalize the crop simulation by computing the Harvest Index.

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

def finalize(self, day: datetime.date) -> None:
    """Finalize the crop simulation by computing the Harvest Index."""
    # Calculate Harvest Index
    TAGP = self.states.TAGP
    if torch.any(TAGP <= 0):
        msg = "Cannot calculate Harvest Index because TAGP=0"
        self.logger.warning(msg)
    self.states.HI = torch.where(
        TAGP > 0, self.kiosk.TWSO / torch.clamp(TAGP, min=1e-10), torch.full_like(TAGP, -1.0)
    )

    SimulationObject.finalize(self, day)

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	Unit
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	TSUMEM, TBASEM, TEFFMX,TSUM1, TSUM2, DLO, DLC, DVSI,
	DVSEND, DTSMTB, VERNSAT, VERNBASE, VERNDVS
TSUM	DVSI, DVSEND, DTSMTB, VERNSAT, VERNBASE, VERNDVS
TSUME	TBASEM, TEFFMX

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

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

Methods:

• initialize –

  Initialize the DVS_Phenology module.

• calc_rates –

  Compute daily phenological development rates.

• integrate –

  Integrate phenology states and manage stage transitions.

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, shape=None)

Initialize the DVS_Phenology module.

Parameters:

• day –

  start date of the simulation

• kiosk –

  variable kiosk of this PCSE instance

• parvalues –

  ParameterProvider object providing parameters as key/value pairs

• shape –

  optional shape for state and rate tensors (default None for scalar)

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

def initialize(self, day, kiosk, parvalues, shape=None):
    """Initialize the DVS_Phenology module.

    Args:
        day: start date of the simulation
        kiosk: variable kiosk of this PCSE instance
        parvalues: `ParameterProvider` object providing parameters as
            key/value pairs
        shape: optional shape for state and rate tensors (default None for scalar)
    """
    self._device = ComputeConfig.get_device()
    self._dtype = ComputeConfig.get_dtype()

    self.params = self.Parameters(parvalues, shape=shape)

    # Initialize vernalisation for IDSL>=2
    if torch.any(self.params.IDSL >= 2):
        self.vernalisation = Vernalisation(day, kiosk, parvalues, shape=shape)
    else:
        self.vernalisation = None

    # 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, shape=shape)
    self.kiosk = kiosk

    # see issue #60
    self._connect_signal(self._on_CROP_FINISH, signal=signals.crop_finish)

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

    self.states = self.StateVariables(
        kiosk,
        publish="DVS",
        TSUM=0.0,
        TSUME=0.0,
        DVS=DVS,
        DOS=DOS,
        DOE=DOE,
        DOA=-1.0,  # not yet occurred
        DOM=-1.0,  # not yet occurred
        DOH=-1.0,  # not yet occurred
        STAGE=STAGE,
        shape=shape,
    )

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

1. Photoperiod reduction (DVRED) if IDSL >= 1 using daylength.
2. Vernalisation factor (VERNFAC) if IDSL >= 2 and in vegetative stage.
3. Stage-specific:
4. emerging: temperature sum for emergence (DTSUME), DVR via TSUMEM.
5. vegetative: temperature sum (DTSUM) scaled by VERNFAC and DVRED.
6. reproductive: temperature sum (DTSUM) only temperature-driven.
7. 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

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

    # Day length sensitivity
    # daylength returns a Tensor directly; broadcast to parameter shape.
    DAYLP = daylength(day, drv.LAT, dtype=self.dtype, device=self.device)
    DAYLP_t = _broadcast_to(DAYLP, p.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 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, p.shape, self.dtype, self.device)

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

    # Emerging stage (STAGE == 0)
    is_emerging = s.STAGE == 0
    temp_diff = TEMP - p.TBASEM
    max_diff = torch.clamp(p.TEFFMX - p.TBASEM, min=0.0)
    dtsume_emerging = torch.minimum(torch.clamp(temp_diff, min=0.0), max_diff)
    safe_den_em = p.TSUMEM.sign() * torch.maximum(torch.abs(p.TSUMEM), self._epsilon)
    dvr_emerging = 0.1 * dtsume_emerging / safe_den_em
    r.DTSUME = torch.where(is_emerging, dtsume_emerging, r.DTSUME)
    r.DVR = torch.where(is_emerging, dvr_emerging, r.DVR)

    # Vegetative stage (STAGE == 1)
    is_vegetative = s.STAGE == 1
    dtsum_vegetative = p.DTSMTB(TEMP) * VERNFAC * DVRED
    safe_den_v1 = p.TSUM1.sign() * torch.maximum(torch.abs(p.TSUM1), self._epsilon)
    dvr_vegetative = dtsum_vegetative / safe_den_v1
    r.DTSUM = torch.where(is_vegetative, dtsum_vegetative, r.DTSUM)
    r.DVR = torch.where(is_vegetative, dvr_vegetative, r.DVR)

    # Reproductive stage (STAGE == 2)
    is_reproductive = s.STAGE == 2
    dtsum_reproductive = p.DTSMTB(TEMP)
    safe_den_v2 = p.TSUM2.sign() * torch.maximum(torch.abs(p.TSUM2), self._epsilon)
    dvr_reproductive = dtsum_reproductive / safe_den_v2
    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 initialised)

    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

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

    # 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, 1.0, s.STAGE)
    s.DOE = torch.where(should_emerge, day_ordinal, 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, 2.0, s.STAGE)
    s.DOA = torch.where(should_flower, day_ordinal, 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, 3.0, s.STAGE)
    s.DOM = torch.where(should_mature, day_ordinal, 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, see issue #60
    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)

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, shape=None)

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.

• shape (tuple | Size | None, default: None ) –

  Target shape for the state and rate variables.

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

def initialize(self, day, kiosk, parvalues, shape=None):
    """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.
        shape (tuple | torch.Size | None): Target shape for the state and rate variables.
    """
    self._initialize_from_tables(kiosk, parvalues, shape=shape)

integrate

integrate(day, delt=1.0)

Update partitioning factors based on development stage (DVS).

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

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, shape: tuple | Size | None = None) -> 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,
    shape: tuple | torch.Size | None = None,
) -> None:
    """Initialize the assimilation module."""
    self._device = ComputeConfig.get_device()
    self._dtype = ComputeConfig.get_dtype()

    self.kiosk = kiosk
    self.params = self.Parameters(parvalues, shape=shape)
    self.rates = self.RateVariables(kiosk, publish=["PGASS"], shape=shape)

    # 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)
    # Cache for astro() results keyed by (day, lat).  astro() only depends
    # on day and latitude so the same result can be reused across batch
    # elements (which share the same weather driver).
    self._astro_cache: dict = {}

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

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
    # 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 via vectorized torch astro routine.
    # latitude and radiation are passed directly – they may be scalars or
    # tensors; the function returns torch.Tensor results in all cases.
    dayl, _daylp, sinld, cosld, difpp, _atmtr, dsinbe, _angot = astro(
        day, drv.LAT, drv.IRRAD, dtype=self.dtype, device=self.device
    )

    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,
        dtype=self.dtype,
        device=self.device,
    )

    # 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

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, shape: tuple | Size | None = None) -> 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.

• shape (tuple | Size | None, default: None ) –

  Target shape for the state and rate variables.

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

def initialize(
    self,
    day: datetime.date,
    kiosk: VariableKiosk,
    parvalues: ParameterProvider,
    shape: tuple | torch.Size | None = None,
) -> 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.
        shape (tuple | torch.Size | None): Target shape for the state and rate variables.
    """
    self.START_DATE = day
    self.kiosk = kiosk

    # Get defaults from ComputeConfig if not already set
    self._device = ComputeConfig.get_device()
    self._dtype = ComputeConfig.get_dtype()

    # TODO check if parvalues are already torch.nn.Parameters
    self.params = self.Parameters(parvalues, shape=shape)
    self.rates = self.RateVariables(kiosk, shape=shape)

    # 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)

    params = self.params

    # Initial leaf biomass
    WLV = (params.TDWI * (1 - self.kiosk["FR"])) * self.kiosk["FL"]
    DWLV = 0.0
    TWLV = WLV + DWLV

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

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

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

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

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_mask = k["DVS"] >= self._zero

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

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

    # death due to self shading cause by high LAI
    LAICR = 3.2 / p.KDIFTB(k["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)

    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.

    # 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 - p.SPAN - self._sigmoid_epsilon) / self._sigmoid_sharpness
        )

        # originial hard mask
        hard_mask = s.LVAGE > p.SPAN

        # 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 > p.SPAN

    r.DALV = torch.sum(span_mask * s.LV, dim=0)
    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, p.shape, self.dtype, self.device)

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

    # specific leaf area of leaves per time step
    r.SLAT = dvs_mask * p.SLATB(k["DVS"])

    # leaf area not to exceed exponential growth curve
    is_lai_exp = s.LAIEXP < 6.0
    DTEFF = torch.clamp(TEMP - p.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.

    r.GLAIEX = torch.where(
        dvs_mask,
        torch.where(is_lai_exp, s.LAIEXP * p.RGRLAI * DTEFF, r.GLAIEX),
        self._zero,
    )

    # source-limited increase in leaf area
    r.GLASOL = torch.where(
        dvs_mask,
        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,
        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

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()

    # 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=0) - rates.DRLV
    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=0, keepdim=True)
    new_biomass = torch.take_along_dim(weight_cumsum, indices=idx_oldest, dim=0)
    tLV = torch.scatter(tLV, dim=0, 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
    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=0)
    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=0)
    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, shape: tuple | Size | None = None) -> 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.

• shape (tuple | Size | None, default: None ) –

  Target shape for the state and rate variables.

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

def initialize(
    self,
    day: datetime.date,
    kiosk: VariableKiosk,
    parvalues: ParameterProvider,
    shape: tuple | torch.Size | None = None,
) -> 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.
        shape (tuple | torch.Size | None): Target shape for the state and rate variables.
    """
    self._device = ComputeConfig.get_device()
    self._dtype = ComputeConfig.get_dtype()

    self.kiosk = kiosk
    self.params = self.Parameters(parvalues, shape=shape)
    self.rates = self.RateVariables(kiosk, publish=["DRRT", "GRRT"], shape=shape)

    # INITIAL STATES
    params = self.params

    # Initial root depth states
    RDM = torch.maximum(params.RDI, torch.minimum(params.RDMCR, params.RDMSOL))
    RD = params.RDI

    # Initial root biomass states
    WRT = params.TDWI * self.kiosk["FR"]
    DWRT = 0.0
    TWRT = WRT + DWRT

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

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

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_mask = k["DVS"] >= 0

    # Increase in root biomass
    r.GRRT = dvs_mask * k["FR"] * k["DMI"]
    r.DRRT = dvs_mask * s.WRT * p.RDRRTB(k["DVS"])
    r.GWRT = r.GRRT - r.DRRT

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

    # Do not let the roots growth if partioning to the roots
    # (variable FR) is zero.
    mask = k["FR"] > 0.0
    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

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

diffwofost.physical_models.crop.storage_organ_dynamics.WOFOST_Storage_Organ_Dynamics

Bases: SimulationObject

Implementation of storage organ dynamics.

Storage organs are the most simple component of the plant in WOFOST and consist of a static pool of biomass. Growth of the storage organs is the result of assimilate partitioning. Death of storage organs is not implemented and the corresponding rate variable (DRSO) is always set to zero.

Pods are green elements of the plant canopy and can as such contribute to the total photosynthetic active area. This is expressed as the Pod Area Index which is obtained by multiplying pod biomass with a fixed Specific Pod Area (SPA).

Simulation parameters

| Name | Description | Type | Unit | |------|===============================================|========|=============| | TDWI | Initial total crop dry weight | SCr | kg ha⁻¹ | | SPA | Specific Pod Area | SCr | ha kg⁻¹ |

State variables

| Name | Description | Pbl | Unit | |------|==================================================|======|=============| | PAI | Pod Area Index | Y | - | | WSO | Weight of living storage organs | Y | kg ha⁻¹ | | DWSO | Weight of dead storage organs | N | kg ha⁻¹ | | TWSO | Total weight of storage organs | Y | kg ha⁻¹ |

Rate variables

| Name | Description | Pbl | Unit | |------|==================================================|======|=============| | GRSO | Growth rate storage organs | N | kg ha⁻¹ d⁻¹ | | DRSO | Death rate storage organs | N | kg ha⁻¹ d⁻¹ | | GWSO | Net change in storage organ biomass | N | kg ha⁻¹ d⁻¹ |

Signals send or handled

None

External dependencies

| Name | Description | Provided by | Unit | |------|====================================|=====================|=============| | ADMI | Above-ground dry matter increase | CropSimulation | kg ha⁻¹ d⁻¹ | | FO | Fraction biomass to storage organs | DVS_Partitioning | - | | FR | Fraction biomass to roots | DVS_Partitioning | - |

Outputs:

Name	Description	Provided by	Unit
PAI	Pod Area Index	Y	-
TWSO	Total weight storage organs	Y	kg ha⁻¹
WSO	Weight living storage organs	Y	kg ha⁻¹

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
PAI	SPA
TWSO	TDWI
WSO	TDWI

Methods:

• initialize –

  Initialize the storage organ dynamics model.

• calc_rates –

  Calculate the rates of change of the state variables.

• integrate –

  Integrate the 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, shape: tuple | Size | None = None) -> None

Initialize the storage organ dynamics 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.

• shape (tuple | Size | None, default: None ) –

  Target shape for the state and rate variables.

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

def initialize(
    self,
    day: datetime.date,
    kiosk: VariableKiosk,
    parvalues: ParameterProvider,
    shape: tuple | torch.Size | None = None,
) -> None:
    """Initialize the storage organ dynamics 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.
        shape (tuple | torch.Size | None): Target shape for the state and rate variables.
    """
    self.kiosk = kiosk
    self.params = self.Parameters(parvalues, shape=shape)
    self.rates = self.RateVariables(kiosk, publish=["GRSO"])

    self._drso_zeros = torch.zeros(self.params.shape, dtype=self.dtype, device=self.device)

    # Initial storage organ biomass
    TDWI = self.params.TDWI
    SPA = self.params.SPA
    FO = self.kiosk["FO"]
    FR = self.kiosk["FR"]

    WSO = (TDWI * (1 - FR)) * FO
    DWSO = self._drso_zeros
    TWSO = WSO + DWSO
    # Initial Pod Area Index
    PAI = WSO * SPA

    self.states = self.StateVariables(
        kiosk, publish=["TWSO", "WSO", "PAI"], WSO=WSO, DWSO=DWSO, TWSO=TWSO, PAI=PAI
    )

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.

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

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.
    """
    rates = self.rates
    k = self.kiosk

    FO = k["FO"]
    ADMI = k["ADMI"]
    REALLOC_SO = k.get("REALLOC_SO", self._drso_zeros)

    # Growth/death rate organs
    rates.GRSO = ADMI * FO
    rates.DRSO = self._drso_zeros
    rates.GWSO = rates.GRSO - rates.DRSO + REALLOC_SO

integrate

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

Integrate the state variables.

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/storage_organ_dynamics.py

def integrate(self, day: datetime.date = None, delt=1.0) -> None:
    """Integrate the 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.
    """
    params = self.params
    rates = self.rates
    states = self.states

    # Stem biomass (living, dead, total)
    states.WSO = states.WSO + rates.GWSO
    states.DWSO = states.DWSO + rates.DRSO
    states.TWSO = states.WSO + states.DWSO

    # Calculate Pod Area Index (SAI)
    states.PAI = states.WSO * params.SPA

diffwofost.physical_models.crop.respiration.WOFOST_Maintenance_Respiration

Bases: SimulationObject

Maintenance respiration in WOFOST.

WOFOST calculates the maintenance respiration as proportional to the dry weights of the plant organs to be maintained, where each plant organ can be assigned a different maintenance coefficient. Multiplying organ weight with the maintenance coeffients yields the relative maintenance respiration (RMRES) which is than corrected for senescence (parameter RFSETB). Finally, the actual maintenance respiration rate is calculated using the daily mean temperature, assuming a relative increase for each 10 degrees increase in temperature as defined by Q10.

Simulation parameters (provide in cropdata dictionary)

Name	Description	Type	Unit
Q10	Relative increase in maintenance respiration rate with	SCr	-
	each 10 degrees increase in temperature		-
RMR	Relative maintenance respiration rate for roots	SCr	kg CH₂O kg⁻¹ d⁻¹
RMS	Relative maintenance respiration rate for stems	SCr	kg CH₂O kg⁻¹ d⁻¹
RML	Relative maintenance respiration rate for leaves	SCr	kg CH₂O kg⁻¹ d⁻¹
RMO	Relative maintenance respiration rate for storage organs	SCr	kg CH₂O kg⁻¹ d⁻¹
RFSETB	Reduction factor for senescence	TCr	-

Rate variables

Name	Description	Pbl	Unit
PMRES	Potential maintenance respiration rate	N	kg CH₂O ha⁻¹ d⁻¹

Signals send or handled

None

External dependencies

Name	Description	Provided by	Unit
DVS	Crop development stage	DVS_Phenology	-
WRT	Dry weight of living roots	WOFOST_Root_Dynamics	kg ha⁻¹
WST	Dry weight of living stems	WOFOST_Stem_Dynamics	kg ha⁻¹
WLV	Dry weight of living leaves	WOFOST_Leaf_Dynamics	kg ha⁻¹
WSO	Dry weight of living storage organs	WOFOST_Storage_Organ_Dynamics	kg ha⁻¹

Outputs

Name	Description	Pbl	Unit
PMRES	Potential maintenance respiration rate	N	kg CH₂O ha⁻¹ d⁻¹

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
PMRES	Q10, RMR, RML, RMS, RMO, RFSETB

Methods:

• initialize –

  Initialize the maintenance respiration module.

• calc_rates –

  Calculate maintenance respiration rates.

• __call__ –

  Calculate and return maintenance respiration (PMRES).

• 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, shape: tuple | None = None)

Initialize the maintenance respiration module.

Parameters:

• day (date) –

  Start date of the simulation

• kiosk (VariableKiosk) –

  Variable kiosk of this PCSE instance

• parvalues (ParameterProvider) –

  ParameterProvider object providing parameters as key/value pairs

• shape (tuple | None, default: None ) –

  Shape of the parameters tensors (optional)

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

def initialize(
    self,
    day: datetime.date,
    kiosk: VariableKiosk,
    parvalues: ParameterProvider,
    shape: tuple | None = None,
):
    """Initialize the maintenance respiration module.

    Args:
        day: Start date of the simulation
        kiosk: Variable kiosk of this PCSE instance
        parvalues: ParameterProvider object providing parameters as key/value pairs
        shape: Shape of the parameters tensors (optional)
    """
    self.params = self.Parameters(parvalues, shape=shape)
    self.rates = self.RateVariables(kiosk, shape=shape)
    self.kiosk = kiosk

calc_rates

calc_rates(day: date, drv: WeatherDataContainer)

Calculate maintenance respiration rates.

Parameters:

• day (date) –

  Current date

• drv (WeatherDataContainer) –

  Weather data for the current day

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

def calc_rates(self, day: datetime.date, drv: WeatherDataContainer):
    """Calculate maintenance respiration rates.

    Args:
        day: Current date
        drv: Weather data for the current day
    """
    p = self.params
    kk = self.kiosk
    r = self.rates

    Q10 = p.Q10
    RMR = p.RMR
    RML = p.RML
    RMS = p.RMS
    RMO = p.RMO

    WRT = kk["WRT"]
    WLV = kk["WLV"]
    WST = kk["WST"]
    WSO = kk["WSO"]
    # [!] DVS needs to be broadcasted explicetly because it is used
    # in torch.where and the kiosk does not format it correctly
    # TODO see #22
    DVS = _broadcast_to(kk["DVS"], p.shape, self.dtype, self.device)

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

    RMRES = RMR * WRT + RML * WLV + RMS * WST + RMO * WSO
    RMRES = RMRES * p.RFSETB(DVS)
    TEFF = Q10 ** ((TEMP - 25.0) / 10.0)
    PMRES = RMRES * TEFF

    # No maintenance respiration before emergence (DVS < 0).
    r.PMRES = torch.where(DVS < 0, torch.zeros_like(PMRES), PMRES)

__call__

__call__(day: date, drv: WeatherDataContainer)

Calculate and return maintenance respiration (PMRES).

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

def __call__(self, day: datetime.date, drv: WeatherDataContainer):
    """Calculate and return maintenance respiration (PMRES)."""
    self.calc_rates(day, drv)
    return self.rates.PMRES

integrate

integrate(day: date, delt: float = 1.0)

No state variables to integrate for this module.

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

def integrate(self, day: datetime.date, delt: float = 1.0):
    """No state variables to integrate for this module."""
    return

diffwofost.physical_models.crop.evapotranspiration.Evapotranspiration

Bases: _BaseEvapotranspirationNonLayered

Potential evaporation (water and soil) rates and crop transpiration rate.

Simulation parameters

Name	Description	Type	Unit
CFET	Correction factor for potential transpiration rate	SCr	-
DEPNR	Dependency number for crop sensitivity to soil moisture stress.	SCr	-
KDIFTB	Extinction coefficient for diffuse visible light vs DVS	TCr	-
IAIRDU	Switch airducts on (1) or off (0)	SCr	-
IOX	Switch oxygen stress on (1) or off (0)	SCr	-
CRAIRC	Critical air content for root aeration	SSo	-
SM0	Soil porosity	SSo	-
SMW	Volumetric soil moisture at wilting point	SSo	-
SMFCF	Volumetric soil moisture at field capacity	SSo	-

State variables

Name	Description	Pbl	Unit
IDWST	Number of days with water stress	N	-
IDOST	Number of days with oxygen stress	N	-

Rate variables

Name	Description	Pbl	Unit
EVWMX	Max evaporation rate from open water surface	Y	cm day⁻¹
EVSMX	Max evaporation rate from wet soil surface	Y	cm day⁻¹
TRAMX	Max transpiration rate from canopy	Y	cm day⁻¹
TRA	Actual transpiration rate from canopy	Y	cm day⁻¹
IDOS	Indicates oxygen stress on this day (True	False)	N
IDWS	Indicates water stress on this day (True	False)	N
RFWS	Reduction factor for water stress	N	-
RFOS	Reduction factor for oxygen stress	N	-
RFTRA	Combined reduction factor for transpiration	Y	-

External dependencies

Name	Description	Provided by	Unit
DVS	Crop development stage	Phenology	-
LAI	Leaf area index	Leaf dynamics	-
SM	Volumetric soil moisture content	Waterbalance	-

Outputs

Name	Description	Pbl	Unit
TRA	Actual transpiration rate from canopy	Y	cm day⁻¹
TRAMX	Max transpiration rate from canopy	Y	cm day⁻¹
EVWMX	Max evaporation rate from open water surface	Y	cm day⁻¹
EVSMX	Max evaporation rate from wet soil surface	Y	cm day⁻¹
RFTRA	Combined reduction factor for transpiration	Y	-

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
EVWMX	KDIFTB
EVSMX	KDIFTB
TRAMX	CFET, KDIFTB
TRA	CFET, KDIFTB, DEPNR, SMFCF, SMW, CRAIRC, SM0
RFTRA	CFET, DEPNR, SMFCF, SMW, CRAIRC, SM0

Methods:

• initialize –

  Initialize the standard evapotranspiration module (no CO2 effects).

initialize

initialize(day: date, kiosk: VariableKiosk, parvalues: ParameterProvider, shape: tuple | None = None) -> None

Initialize the standard evapotranspiration module (no CO2 effects).

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.

• shape (tuple | Size | None, default: None ) –

  Target shape for the state and rate variables.

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

def initialize(
    self,
    day: datetime.date,
    kiosk: VariableKiosk,
    parvalues: ParameterProvider,
    shape: tuple | None = None,
) -> None:
    """Initialize the standard evapotranspiration module (no CO2 effects).

    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.
        shape (tuple | torch.Size | None): Target shape for the state and rate variables.
    """
    self._initialize_base(
        day,
        kiosk,
        parvalues,
        publish_rates=["EVWMX", "EVSMX", "TRAMX", "TRA", "RFTRA"],
        shape=shape,
    )

Soil modules

diffwofost.physical_models.soil.classic_waterbalance.WaterbalancePP

Bases: SimulationObject

Fake waterbalance for simulation under potential production.

Keeps the soil moisture content at field capacity and only accumulates crop transpiration and soil evaporation rates through the course of the simulation.

Methods:

• initialize –

  Initialize the potential-production waterbalance.

• calc_rates –

  Calculate soil evaporation and transpiration rates.

• integrate –

  Integrate state variables over one time step.

initialize

initialize(day, kiosk, parvalues, shape: tuple | Size | None = None)

Initialize the potential-production waterbalance.

Parameters:

• day –

  start date of the simulation

• kiosk –

  variable kiosk of this PCSE instance

• parvalues –

  ParameterProvider object containing all parameters

• shape (tuple | Size | None, default: None ) –

  optional shape for state and rate tensors (default None for scalar)

This waterbalance keeps the soil moisture always at field capacity. Therefore WaterbalancePP has only one parameter (SMFCF: the field capacity of the soil) and one state variable (SM: the volumetric soil moisture).

Source code in src/diffwofost/physical_models/soil/classic_waterbalance.py

def initialize(self, day, kiosk, parvalues, shape: tuple | torch.Size | None = None):
    """Initialize the potential-production waterbalance.

    Args:
        day: start date of the simulation
        kiosk: variable kiosk of this PCSE instance
        parvalues: ParameterProvider object containing all parameters
        shape: optional shape for state and rate tensors (default None for scalar)

    This waterbalance keeps the soil moisture always at field capacity. Therefore
    `WaterbalancePP` has only one parameter (`SMFCF`: the field capacity of the
    soil) and one state variable (`SM`: the volumetric soil moisture).
    """
    self.params = self.Parameters(parvalues, shape=shape)
    self.rates = self.RateVariables(kiosk, publish="EVS", shape=shape)
    self.states = self.StateVariables(
        kiosk,
        SM=self.params.SMFCF,
        publish="SM",
        EVST=0.0,
        WTRAT=0.0,
        shape=shape,
    )

calc_rates

calc_rates(day, drv)

Calculate soil evaporation and transpiration rates.

Source code in src/diffwofost/physical_models/soil/classic_waterbalance.py

def calc_rates(self, day, drv):
    """Calculate soil evaporation and transpiration rates."""
    r = self.rates
    # Transpiration and maximum soil and surface water evaporation rates
    # are calculated by the crop Evapotranspiration module.
    # However, if the crop is not yet emerged then set TRA=0 and use
    # the potential soil/water evaporation rates directly because there is
    # no shading by the canopy.
    if "TRA" not in self.kiosk:
        r.WTRA = torch.zeros_like(torch.as_tensor(drv.ES0))
        EVSMX = torch.as_tensor(drv.ES0)
    else:
        r.WTRA = self.kiosk["TRA"]
        EVSMX = torch.as_tensor(self.kiosk["EVSMX"])

    # Actual evaporation rates

    rain_ge_1 = torch.as_tensor(self.RAINold) >= 1.0

    # If rainfall amount >= 1cm on previous day assume maximum soil evaporation.
    # Else soil evaporation is a function days-since-last-rain (DSLR).
    dslr_inc = torch.as_tensor(self.DSLR) + 1.0
    evsmxt = EVSMX * (torch.sqrt(dslr_inc) - torch.sqrt(dslr_inc - 1.0))
    evs_else = torch.minimum(EVSMX, evsmxt + torch.as_tensor(self.RAINold))

    r.EVS = torch.where(rain_ge_1, EVSMX, evs_else)
    self.DSLR = torch.where(rain_ge_1, torch.ones_like(dslr_inc), dslr_inc)

    # Hold rainfall amount to keep track of soil surface wetness and reset self.DSLR if needed
    self.RAINold = torch.as_tensor(drv.RAIN)

integrate

integrate(day, delt=1.0)

Integrate state variables over one time step.

Source code in src/diffwofost/physical_models/soil/classic_waterbalance.py

def integrate(self, day, delt=1.0):
    """Integrate state variables over one time step."""
    # Keep soil moisture on field capacity
    self.states.SM = self.params.SMFCF

    # Accumulated transpiration and soil evaporation amounts
    # Avoid in-place updates because TensorStatesTemplate may broadcast
    # state tensors using `expand()`, which forbids in-place writes.
    self.states.EVST = self.states.EVST + (self.rates.EVS * delt)
    self.states.WTRAT = self.states.WTRAT + (self.rates.WTRA * delt)

diffwofost.physical_models.soil.classic_waterbalance.WaterbalanceFD

Bases: SimulationObject

Waterbalance for freely draining soils under water-limited production.

The purpose of the soil water balance calculations is to estimate the daily value of the soil moisture content. The soil moisture content influences soil moisture uptake and crop transpiration.

The dynamic calculations are carried out in two sections, one for the calculation of rates of change per timestep (= 1 day) and one for the calculation of summation variables and state variables. The water balance is driven by rainfall, possibly buffered as surface storage, and evapotranspiration. The processes considered are infiltration, soil water retention, percolation (here conceived as downward water flow from rooted zone to second layer), and the loss of water beyond the maximum root zone.

The textural profile of the soil is conceived as homogeneous. Initially the soil profile consists of two layers, the actually rooted soil and the soil immediately below the rooted zone until the maximum rooting depth is reached by roots(soil and crop dependent). The extension of the root zone from the initial rooting depth to maximum rooting depth is described in Root_Dynamics class. From the moment that the maximum rooting depth is reached the soil profile may be described as a one layer system depending if the roots are able to penetrate the entire profile. If not a non-rooted part remains at the bottom of the profile.

The class WaterbalanceFD is derived from WATFD.FOR in WOFOST7.1 with the exception that the depth of the soil is now completely determined by the maximum soil depth (RDMSOL) and not by the minimum of soil depth and crop maximum rooting depth (RDMCR).

Simulation parameters:

Name	Description	Type	Unit
SMFCF	Field capacity of the soil	SSo	-
SM0	Porosity of the soil	SSo	-
SMW	Wilting point of the soil	SSo	-
CRAIRC	Soil critical air content (waterlogging)	SSo	-
SOPE	Maximum percolation rate root zone	SSo	cm day⁻¹
KSUB	Maximum percolation rate subsoil	SSo	cm day⁻¹
RDMSOL	Soil rootable depth	SSo	cm
IFUNRN	Indicates whether non-infiltrating fraction of rain is a function of storm size (1) or not (0)	SSi	-
SSMAX	Maximum surface storage	SSi	cm
SSI	Initial surface storage	SSi	cm
WAV	Initial amount of water in total soil profile	SSi	cm
NOTINF	Maximum fraction of rain not infiltrating into the soil	SSi	-
SMLIM	Initial maximum moisture content in initial rooting depth zone	SSi	-

State variables:

Name	Description	Pbl	Unit
SM	Volumetric moisture content in root zone	Y	-
SS	Surface storage (layer of water on surface)	N	cm
SSI	Initial surface storage	N	cm
W	Amount of water in root zone	N	cm
WI	Initial amount of water in the root zone	N	cm
WLOW	Amount of water in the subsoil between current rooting depth and maximum rootable depth	N	cm
WLOWI	Initial amount of water in the subsoil	N	cm
WWLOW	Total amount of water in the soil profile; WWLOW = WLOW + W	N	cm
WTRAT	Total water lost as transpiration from the water balance; can differ from CTRAT which only counts transpiration for a crop cycle	N	cm
EVST	Total evaporation from the soil surface	N	cm
EVWT	Total evaporation from a water surface	N	cm
TSR	Total surface runoff	N	cm
RAINT	Total amount of rainfall (effective + non-effective)	N	cm
WDRT	Amount of water added to root zone by increase of root growth	N	cm
TOTINF	Total amount of infiltration	N	cm
TOTIRR	Total amount of effective irrigation	N	cm
PERCT	Total amount of water percolating from rooted zone to subsoil	N	cm
LOSST	Total amount of water lost to deeper soil	N	cm
DSOS	Days since oxygen stress, accumulating consecutive days of oxygen stress	Y	-
WBALRT	Checksum for root zone water balance; computed in finalize(), abs(WBALRT) > 0.0001 raises a WaterBalanceError	N	cm
WBALTT	Checksum for total water balance; computed in finalize(), abs(WBALTT) > 0.0001 raises a WaterBalanceError	N	cm

Rate variables:

Name	Description	Pbl	Unit
EVS	Actual evaporation rate from soil	N	cm day⁻¹
EVW	Actual evaporation rate from water surface	N	cm day⁻¹
WTRA	Actual transpiration rate from plant canopy, directly derived from TRA in the evapotranspiration module	N	cm day⁻¹
RAIN_INF	Infiltrating rainfall rate for current day	N	cm day⁻¹
RAIN_NOTINF	Non-infiltrating rainfall rate for current day	N	cm day⁻¹
RIN	Infiltration rate for current day	N	cm day⁻¹
RIRR	Effective irrigation rate for current day, computed as irrigation amount times efficiency	N	cm day⁻¹
PERC	Percolation rate to non-rooted zone	N	cm day⁻¹
LOSS	Rate of water loss to deeper soil	N	cm day⁻¹
DW	Change in amount of water in rooted zone due to infiltration, transpiration, and evaporation	N	cm day⁻¹
DWLOW	Change in amount of water in subsoil	N	cm day⁻¹
DTSR	Change in surface runoff	N	cm day⁻¹
DSS	Change in surface storage	N	cm day⁻¹

External dependencies:

Name	Description	Provided by	Unit
TRA	Crop transpiration rate	Evapotranspiration	cm day⁻¹
EVSMX	Maximum evaporation rate from a soil surface below the crop canopy	Evapotranspiration	cm day⁻¹
EVWMX	Maximum evaporation rate from a water surface below the crop canopy	Evapotranspiration	cm day⁻¹
RD	Rooting depth	Root_dynamics	cm

Outputs

Name	Description	Pbl	Unit
SM	Volumetric soil moisture in root zone	Y	-
EVS	Actual evaporation rate from soil	Y	cm day⁻¹
DSOS	Days since oxygen stress	Y	-

Gradient mapping (which parameters have a gradient):

Output	Parameters influencing it
SM	SMFCF, SMW, SM0
EVS	SMFCF, SMW, SM0
DSOS	None; thresholded diagnostic output

Exceptions raised:

A WaterbalanceError is raised when the waterbalance is not closing at the end of the simulation cycle (e.g water has "leaked" away).

Methods:

• initialize –

  Initialize the freely-draining water balance.

• calc_rates –

  Calculate the rates of change for all water balance components.

• integrate –

  Integrate state variables over one time step.

• finalize –

  Calculate and check water balance checksums at end of simulation.

initialize

initialize(day, kiosk, parvalues, shape=None)

Initialize the freely-draining water balance.

Parameters:

• day –

  Start date of the simulation.

• kiosk –

  Variable kiosk used to read and publish crop state.

• parvalues –

  Parameter provider containing the physical-model parameters for the soil.

• shape –

  Target tensor shape for state and rate variables.

Source code in src/diffwofost/physical_models/soil/classic_waterbalance.py

def initialize(self, day, kiosk, parvalues, shape=None):
    """Initialize the freely-draining water balance.

    Args:
        day: Start date of the simulation.
        kiosk: Variable kiosk used to read and publish crop state.
        parvalues: Parameter provider containing the physical-model
            parameters for the soil.
        shape: Target tensor shape for state and rate variables.
    """
    dtype = ComputeConfig.get_dtype()
    device = ComputeConfig.get_device()

    self.params = self.Parameters(parvalues, shape=shape)

    # Check validity of SMLIM (site parameter, not stored in self.params)
    SMLIM_raw = torch.as_tensor(parvalues["SMLIM"], dtype=dtype, device=device)
    SMLIM = torch.maximum(self.params.SMW, torch.minimum(self.params.SM0, SMLIM_raw))
    if torch.any(SMLIM != SMLIM_raw):
        self.logger.warning("SMLIM not in valid range, adjusted.")

    # Default rooting depth: 10 cm; derive maximum rootable depth
    self.RDold = torch.as_tensor(self.DEFAULT_RD, dtype=dtype, device=device)
    self.RDM = torch.maximum(self.RDold, self.params.RDMSOL)

    # Initial soil moisture content in rooted zone, clamped to [SMW, SMLIM]
    SM = torch.clamp(
        self.params.SMW + self.params.WAV / self.RDold,
        min=self.params.SMW,
        max=SMLIM,
    )
    W = SM * self.RDold

    # Initial water in subsoil (below rooted zone to maximum rootable depth)
    WLOW = torch.clamp(
        self.params.WAV + self.RDM * self.params.SMW - W,
        min=torch.zeros_like(W),
        max=self.params.SM0 * (self.RDM - self.RDold),
    )

    # Days since last rain: 1 if SM is above halfway between SMW and SMFCF, else 5
    halfway = self.params.SMW + 0.5 * (self.params.SMFCF - self.params.SMW)
    self.DSLR = torch.where(SM >= halfway, torch.ones_like(SM), torch.full_like(SM, 5.0))

    # Initialize helper variables
    self.RINold = torch.zeros_like(SM)
    self.in_crop_cycle = False
    self.rooted_layer_needs_reset = False
    self.NINFTB = Afgen([0.0, 0.0, 0.5, 0.0, 1.5, 1.0])
    self._increments_w = []

    # Initialize state and rate variables
    self.states = self.StateVariables(
        kiosk,
        publish=["SM", "DSOS"],
        SM=SM,
        SS=self.params.SSI,
        SSI=self.params.SSI,
        W=W,
        WI=W,
        WLOW=WLOW,
        WLOWI=WLOW,
        WWLOW=(W + WLOW),  # Total water in soil column
        WTRAT=0.0,
        EVST=0.0,
        EVWT=0.0,
        TSR=0.0,
        RAINT=0.0,
        WDRT=0.0,
        TOTINF=0.0,
        TOTIRR=0.0,
        DSOS=0.0,
        PERCT=0.0,
        LOSST=0.0,
        WBALRT=-999.0,
        WBALTT=-999.0,
        shape=shape,
    )
    self.rates = self.RateVariables(kiosk, publish="EVS", shape=shape)

    # Connect signals
    self._connect_signal(self._on_CROP_START, signals.crop_start)
    self._connect_signal(self._on_CROP_FINISH, signals.crop_finish)
    self._connect_signal(self._on_IRRIGATE, signals.irrigate)

calc_rates

calc_rates(day, drv)

Calculate the rates of change for all water balance components.

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/soil/classic_waterbalance.py

def calc_rates(self, day, drv):
    """Calculate the rates of change for all water balance components.

    Args:
        day (datetime.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.

    """
    s = self.states
    p = self.params
    r = self.rates
    k = self.kiosk
    dtype = ComputeConfig.get_dtype()
    device = ComputeConfig.get_device()

    # Rate of irrigation (RIRR)
    r.RIRR = self._RIRR
    self._RIRR = torch.zeros_like(self._RIRR)

    # Transpiration and maximum evaporation rates from crop module.
    # Before emergence there is no canopy shading yet, so the water balance
    # must fall back to the weather-driven soil and surface evaporation.
    weather_wtra = torch.zeros_like(torch.as_tensor(drv.ES0, dtype=dtype, device=device))
    weather_evwmx = torch.as_tensor(drv.E0, dtype=dtype, device=device)
    weather_evsmx = torch.as_tensor(drv.ES0, dtype=dtype, device=device)

    if "TRA" not in k:
        r.WTRA = weather_wtra
        EVWMX = weather_evwmx
        EVSMX = weather_evsmx
    else:
        crop_wtra = torch.as_tensor(k["TRA"], dtype=dtype, device=device)
        crop_evwmx = torch.as_tensor(k["EVWMX"], dtype=dtype, device=device)
        crop_evsmx = torch.as_tensor(k["EVSMX"], dtype=dtype, device=device)

        if "DVS" in k:
            emerged_mask = torch.as_tensor(k["DVS"], dtype=dtype, device=device) >= 0.0
            r.WTRA = torch.where(emerged_mask, crop_wtra, weather_wtra)
            EVWMX = torch.where(emerged_mask, crop_evwmx, weather_evwmx)
            EVSMX = torch.where(emerged_mask, crop_evsmx, weather_evsmx)
        else:
            r.WTRA = crop_wtra
            EVWMX = crop_evwmx
            EVSMX = crop_evsmx

    # Actual evaporation rates
    # If SS > 1 cm: evaporate from surface water; else from soil
    ss_gt_1 = s.SS > 1.0
    rin_ge_1 = self.RINold >= 1.0
    DSLR_t = self.DSLR
    dslr_inc = DSLR_t + 1.0
    # Soil evaporation as function of days-since-last-rain (DSLR)
    EVSMXT = EVSMX * (torch.sqrt(dslr_inc) - torch.sqrt(DSLR_t))
    evs_dslr = torch.minimum(EVSMX, EVSMXT + self.RINold)
    EVS_no_ss = torch.where(rin_ge_1, EVSMX, evs_dslr)

    r.EVW = torch.where(ss_gt_1, EVWMX, torch.zeros_like(EVWMX))
    r.EVS = torch.where(ss_gt_1, torch.zeros_like(EVS_no_ss), EVS_no_ss)

    # Update DSLR:
    #   - SS > 1: unchanged (no soil evaporation)
    #   - SS <= 1 and RINold >= 1: reset to 1
    #   - SS <= 1 and RINold < 1: increment by 1
    self.DSLR = torch.where(
        ss_gt_1,
        DSLR_t,
        torch.where(rin_ge_1, torch.ones_like(DSLR_t), dslr_inc),
    )

    # Potentially infiltrating rainfall
    RAIN_t = torch.as_tensor(drv.RAIN, dtype=dtype, device=device)
    RINPRE_fixed = (1.0 - p.NOTINF) * RAIN_t
    RINPRE_storm = (1.0 - p.NOTINF * self.NINFTB(RAIN_t)) * RAIN_t
    # IFUNRN: 0 = fixed non-infiltrating fraction, 1 = function of storm size
    RINPRE = torch.where(p.IFUNRN == 0, RINPRE_fixed, RINPRE_storm)

    # Second stage: add surface storage and irrigation
    RINPRE = RINPRE + r.RIRR + s.SS
    # With surface storage, infiltration is limited by SOPE
    AVAIL = RINPRE - r.EVW
    RINPRE = torch.where(s.SS > 0.1, torch.minimum(p.SOPE, AVAIL), RINPRE)

    RD = self._determine_rooting_depth()

    # Equilibrium water in rooted zone
    WE = p.SMFCF * RD
    # Percolation: excess moisture above field capacity, capped at SOPE
    PERC1 = torch.maximum(
        torch.zeros_like(s.W),
        torch.minimum(p.SOPE, (s.W - WE) - r.WTRA - r.EVS),
    )

    # Loss from lower zone
    WELOW = p.SMFCF * (self.RDM - RD)
    r.LOSS = torch.maximum(
        torch.zeros_like(s.WLOW),
        torch.minimum(p.KSUB, s.WLOW - WELOW + PERC1),
    )

    # Percolation capped by subsoil uptake capacity
    PERC2 = (self.RDM - RD) * p.SM0 - s.WLOW + r.LOSS
    r.PERC = torch.minimum(PERC1, PERC2)

    # Adjusted infiltration rate
    r.RIN = torch.minimum(RINPRE, (p.SM0 - s.SM) * RD + r.WTRA + r.EVS + r.PERC)
    self.RINold = r.RIN

    # Rates of change in water amounts
    r.DW = r.RIN - r.WTRA - r.EVS - r.PERC
    r.DWLOW = r.PERC - r.LOSS

    # Prevent W from going negative by reducing EVS
    Wtmp = s.W + r.DW
    under_zero = Wtmp < 0.0
    r.EVS = torch.where(under_zero, torch.maximum(torch.zeros_like(r.EVS), r.EVS + Wtmp), r.EVS)
    r.DW = torch.where(under_zero, -s.W, r.DW)

    # Surface storage and runoff
    SStmp = RAIN_t + r.RIRR - r.EVW - r.RIN
    r.DSS = torch.minimum(SStmp, p.SSMAX - s.SS)
    r.DTSR = SStmp - r.DSS
    # incoming rainfall rate
    r.DRAINT = RAIN_t

integrate

integrate(day, delt=1.0)

Integrate state variables over one time step.

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/soil/classic_waterbalance.py

def integrate(self, day, delt=1.0):
    """Integrate state variables over one time step.

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

    # INTEGRALS OF THE WATERBALANCE: SUMMATIONS AND STATE VARIABLES

    # total transpiration
    s.WTRAT = s.WTRAT + r.WTRA * delt

    # total evaporation from surface water layer and/or soil
    s.EVWT = s.EVWT + r.EVW * delt
    s.EVST = s.EVST + r.EVS * delt

    # totals for rainfall, irrigation and infiltration
    s.RAINT = s.RAINT + r.DRAINT * delt
    s.TOTINF = s.TOTINF + r.RIN * delt
    s.TOTIRR = s.TOTIRR + r.RIRR * delt

    # Update surface storage and total surface runoff (TSR)
    s.SS = s.SS + r.DSS * delt
    s.TSR = s.TSR + r.DTSR * delt

    # amount of water in rooted zone
    s.W = s.W + r.DW * delt

    # total percolation and loss of water by deep leaching
    s.PERCT = s.PERCT + r.PERC * delt
    s.LOSST = s.LOSST + r.LOSS * delt

    # amount of water in unrooted, lower part of rootable zone
    s.WLOW = s.WLOW + r.DWLOW * delt
    # total amount of water in the whole rootable zone
    s.WWLOW = s.W + s.WLOW

    # CHANGE OF ROOTZONE SUBSYSTEM BOUNDARY

    RD = self._determine_rooting_depth()
    RDchange = RD - self.RDold
    self._redistribute_water(RDchange)

    # mean soil moisture content in rooted zone
    s.SM = s.W / RD

    # Accumulate days since oxygen stress
    stress_mask = s.SM >= (p.SM0 - p.CRAIRC)
    s.DSOS = torch.where(stress_mask, s.DSOS + 1.0, torch.zeros_like(s.DSOS))

    # save rooting depth
    self.RDold = RD

finalize

finalize(day)

Calculate and check water balance checksums at end of simulation.

Source code in src/diffwofost/physical_models/soil/classic_waterbalance.py

def finalize(self, day):
    """Calculate and check water balance checksums at end of simulation."""
    s = self.states

    wbal_increments = sum(self._increments_w) if self._increments_w else 0.0

    # Checksums for rootzone (RT) and whole system (TT)
    s.WBALRT = s.TOTINF + s.WI + s.WDRT - s.EVST - s.WTRAT - s.PERCT - s.W + wbal_increments
    s.WBALTT = (
        s.SSI
        + s.RAINT
        + s.TOTIRR
        + s.WI
        - s.W
        + wbal_increments
        + s.WLOWI
        - s.WLOW
        - s.WTRAT
        - s.EVWT
        - s.EVST
        - s.TSR
        - s.LOSST
        - s.SS
    )

    if torch.any(torch.abs(s.WBALRT) > 0.0001):
        msg = "Water balance for root zone does not close."
        raise exc.WaterBalanceError(msg)

    if torch.any(torch.abs(s.WBALTT) > 0.0001):
        total_in = (s.WI + s.WLOWI + s.SSI + s.TOTIRR + s.RAINT).tolist()
        total_out = (s.W + s.WLOW + s.SS + s.EVWT + s.EVST + s.WTRAT + s.TSR + s.LOSST).tolist()
        msg = "Water balance for complete soil profile does not close.\n"
        msg += f"Total INIT + IN:   {total_in}\n"
        msg += f"Total FINAL + OUT: {total_out}"
        raise exc.WaterBalanceError(msg)

    SimulationObject.finalize(self, day)

Utility (under development)

diffwofost.physical_models.config.Configuration dataclass

Configuration(CROP: type[SimulationObject], CROP_COMPONENTS: dict | None = None, CROP_NN_MODEL: type[Module] | None = None, 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.

Used by:

• Utility (under development)
  • Engine
  • EngineTestHelper

Methods:

• __post_init__ –

  Validate config data based on CROP.initialize signature.

• 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.

__post_init__

__post_init__()

Validate config data based on CROP.initialize signature.

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

def __post_init__(self):
    """Validate config data based on CROP.initialize signature."""
    sig_arguments = inspect.signature(self.CROP.initialize).parameters

    # Nullify CROP_NN_MODEL and CROP_COMPONENTS, if not compatible with CROP.initialize
    for field_value, sig_key, attr_name in [
        (self.CROP_NN_MODEL, "nn_model", "CROP_NN_MODEL"),
        (self.CROP_COMPONENTS, "component_overrides", "CROP_COMPONENTS"),
    ]:
        if field_value is not None and sig_key not in sig_arguments:
            setattr(self, attr_name, None)

    # Validate component overrides have "class" key with non-None value
    for component_name, override in (self.CROP_COMPONENTS or {}).items():
        self._validate_component_override(component_name, override)

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 acts as a factory for default configuration settings that are captured by simulation objects upon initialization. This enables precise control over where (device) and how (dtype) each model computation occurs, allowing for multiple models with different configurations to coexist.

Key Concept: Configuration Capture

When a simulation object (e.g., WOFOST_Leaf_Dynamics) is initialized, it queries ComputeConfig for the current device and dtype. The model captures and stores these settings for its lifetime. Subsequent changes to ComputeConfig will only affect newly created objects, leaving existing ones unchanged.

Default Behavior:

• Device: Defaults to torch.get_default_device()
• Dtype: Defaults to torch.get_default_dtype()

Basic Usage:

>>> from diffwofost.physical_models.config import ComputeConfig
>>> import torch
>>>
>>> # Configure defaults for new models
>>> ComputeConfig.set_device('cuda')
>>> ComputeConfig.set_dtype(torch.float32)
>>>
>>> # Get current defaults
>>> device = ComputeConfig.get_device()
>>> dtype = ComputeConfig.get_dtype()

Creating Models with Different Settings:

Because models capture the configuration at initialization, you can create instances with different settings in the same process:

>>> from diffwofost.physical_models.crop.leaf_dynamics import WOFOST_Leaf_Dynamics
>>>
>>> # Create a model on GPU (float32)
>>> ComputeConfig.set_device('cuda')
>>> ComputeConfig.set_dtype(torch.float32)
>>> model_gpu = WOFOST_Leaf_Dynamics(...)
>>>
>>> # Create a model on CPU (float64)
>>> ComputeConfig.set_device('cpu')
>>> ComputeConfig.set_dtype(torch.float64)
>>> model_cpu = WOFOST_Leaf_Dynamics(...)
>>>
>>> # model_gpu remains on cuda, model_cpu stays on cpu.

Setting the model properties like model.device = torch.device("cpu") or model.dtype = torch.float64 returns AttributeError. Always use ComputeConfig.set_device(...) and ComputeConfig.set_dtype(...).

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(config: str | Path | Configuration | None = None)

Bases: Engine

PCSE engine wrapper with reusable runtime state.

The engine accepts either a loaded Configuration instance or a path to a PCSE configuration file. Calling setup(...) resets all run-specific state so the same Engine object can be reused safely for multiple simulations.

Parameters:

• config (str | Path | Configuration | None, default: None ) –

  Model configuration as a loaded Configuration instance or a path to a PCSE configuration file.

Raises:

• TypeError –

  If no configuration is provided.

Subclassed by:

• Utility (under development) EngineTestHelper

Methods:

• setup –

  Set up the engine for a new simulation run.

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

def __init__(
    self,
    config: str | Path | Configuration | None = None,
):
    """Initialize the engine with a model configuration.

    Args:
        config (str | pathlib.Path | Configuration | None): Model
            configuration as a loaded Configuration instance or a path to a
            PCSE configuration file.

    Raises:
        TypeError: If no configuration is provided.
    """
    BaseEngine.__init__(self)

    if config is None:
        msg = "A model configuration must be provided when initializing the engine."
        raise TypeError(msg)

    # 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

    # Get the arguments of the CROP class's initialize method
    self._components_overrides = None
    if self.mconf.CROP_COMPONENTS is not None:
        self._components_overrides = normalize_components(
            self.mconf.CROP_COMPONENTS,
            self.mconf.CROP.COMPONENT_SPECS,
        )

setup

setup(parameterprovider, weatherdataprovider, agromanagement, external_states=None)

Set up the engine for a new simulation run.

Parameters:

• parameterprovider –

  Provider with crop and soil parameter values.

• weatherdataprovider –

  Provider used to retrieve daily driving weather variables.

• agromanagement –

  AgroManagement definition passed to the configured agromanagement component.

• external_states (list[dict] | None, default: None ) –

  Optional list of day-keyed external states to inject through the VariableKiosk.

Returns:

• Engine –

  The configured engine instance.

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

def setup(
    self,
    parameterprovider,
    weatherdataprovider,
    agromanagement,
    external_states=None,
):
    """Set up the engine for a new simulation run.

    Args:
        parameterprovider: Provider with crop and soil parameter values.
        weatherdataprovider: Provider used to retrieve daily driving
            weather variables.
        agromanagement: AgroManagement definition passed to the configured
            agromanagement component.
        external_states (list[dict] | None): Optional list of day-keyed
            external states to inject through the VariableKiosk.

    Returns:
        Engine: The configured engine instance.
    """
    self._reset_runtime_state()

    self.parameterprovider = parameterprovider
    self._shape = _get_params_shape(self.parameterprovider)

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

    # 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()
    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)
    return self

diffwofost.physical_models.utils.EngineTestHelper

EngineTestHelper(config: str | Path | Configuration | None = None)

Bases: Engine

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

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

def __init__(
    self,
    config: str | Path | Configuration | None = None,
):
    """Initialize the engine with a model configuration.

    Args:
        config (str | pathlib.Path | Configuration | None): Model
            configuration as a loaded Configuration instance or a path to a
            PCSE configuration file.

    Raises:
        TypeError: If no configuration is provided.
    """
    BaseEngine.__init__(self)

    if config is None:
        msg = "A model configuration must be provided when initializing the engine."
        raise TypeError(msg)

    # 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

    # Get the arguments of the CROP class's initialize method
    self._components_overrides = None
    if self.mconf.CROP_COMPONENTS is not None:
        self._components_overrides = normalize_components(
            self.mconf.CROP_COMPONENTS,
            self.mconf.CROP.COMPONENT_SPECS,
        )

Other classes (for developers)

diffwofost.physical_models.base.states_rates.TensorStatesTemplate

TensorStatesTemplate(kiosk=None, publish=None, shape=None, do_not_broadcast=None, **kwargs)

Bases: TensorContainer, StatesTemplate

Template for storing state variable values as tensors.

It includes functionality to broadcast state variables to a common shape. See diffwofost.base.states_rates.TensorContainer and pcse.base.states_rates.StatesTemplate for details.

Source code in src/diffwofost/physical_models/base/states_rates.py

def __init__(self, kiosk=None, publish=None, shape=None, do_not_broadcast=None, **kwargs):
    self._shape = ()
    self._do_not_broadcast = [] if do_not_broadcast is None else do_not_broadcast
    StatesTemplate.__init__(self, kiosk=kiosk, publish=publish, **kwargs)
    self._broadcast(shape)

diffwofost.physical_models.base.states_rates.TensorRatesTemplate

TensorRatesTemplate(kiosk=None, publish=None, shape=None, do_not_broadcast=None)

Bases: TensorContainer, RatesTemplate

Template for storing rate variable values as tensors.

It includes functionality to broadcast rate variables to a common shape. See diffwofost.base.states_rates.TensorContainer and pcse.base.states_rates.RatesTemplate for details.

Source code in src/diffwofost/physical_models/base/states_rates.py

def __init__(self, kiosk=None, publish=None, shape=None, do_not_broadcast=None):
    self._shape = ()
    self._do_not_broadcast = [] if do_not_broadcast is None else do_not_broadcast
    RatesTemplate.__init__(self, kiosk=kiosk, publish=publish)
    self._broadcast(shape)

diffwofost.physical_models.base.states_rates.TensorParamTemplate

TensorParamTemplate(parvalues, shape=None, do_not_broadcast=None)

Bases: TensorContainer, ParamTemplate

Template for storing parameter values as tensors.

It includes functionality to broadcast parameters to a common shape. See diffwofost.base.states_rates.TensorContainer and pcse.base.states_rates.ParamTemplate for details.

Source code in src/diffwofost/physical_models/base/states_rates.py

def __init__(self, parvalues, shape=None, do_not_broadcast=None):
    self._shape = ()
    self._do_not_broadcast = [] if do_not_broadcast is None else do_not_broadcast
    ParamTemplate.__init__(self, parvalues=parvalues)
    self._broadcast(shape)

diffwofost.physical_models.base.states_rates.TensorContainer

TensorContainer(shape=None, do_not_broadcast=None, **variables)

Bases: HasTraits

It includes functionality to broadcast variables to a common shape. This common shape can be inferred from the container's tensor and AFGEN variables, or it can be set as an input argument.

Parameters:

• shape (tuple | Size, default: None ) –

  Shape to which the variables in the container are broadcasted. If given, it should match the shape of all the input variables that already have dimensions. Defaults to None.

• do_not_broadcast (list, default: None ) –

  Name of the variables that are not broadcasted to the container shape. Defaults to None, which means that all variables are broadcasted.

• variables (dict, default: {} ) –

  Collection of variables to initialize the container, as key-value pairs.

Subclassed by:

• Other classes (for developers)
  • TensorParamTemplate
  • TensorRatesTemplate
  • TensorStatesTemplate

Attributes:

• shape –

  Base shape of the variables in the container.

Source code in src/diffwofost/physical_models/base/states_rates.py

def __init__(self, shape=None, do_not_broadcast=None, **variables):
    """Container of tensor variables.

    It includes functionality to broadcast variables to a common shape. This common shape can
    be inferred from the container's tensor and AFGEN variables, or it can be set as an input
    argument.

    Args:
        shape (tuple | torch.Size, optional): Shape to which the variables in the container
            are broadcasted. If given, it should match the shape of all the input variables that
            already have dimensions. Defaults to None.
        do_not_broadcast (list, optional): Name of the variables that are not broadcasted
            to the container shape. Defaults to None, which means that all variables are
            broadcasted.
        variables (dict): Collection of variables to initialize the container, as key-value
            pairs.
    """
    self._shape = ()
    self._do_not_broadcast = [] if do_not_broadcast is None else do_not_broadcast
    HasTraits.__init__(self, **variables)
    self._broadcast(shape)

shape property writable

shape

Base shape of the variables in the container.

diffwofost.physical_models.traitlets.Tensor

Tensor(default_value=Undefined, allow_none=False, read_only=None, help=None, config=None, dtype=None, **kwargs)

Bases: TraitType

Methods:

• validate –

  Validate input object, recasting it into a tensor if possible.

• from_string –

  Casting tensor from string is not supported for now.

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

def __init__(
    self,
    default_value=Undefined,
    allow_none=False,
    read_only=None,
    help=None,
    config=None,
    dtype=None,
    **kwargs,
):
    super().__init__(
        default_value=default_value,
        allow_none=allow_none,
        read_only=read_only,
        help=help,
        config=config,
        **kwargs,
    )
    self.dtype = dtype

validate

validate(obj, value)

Validate input object, recasting it into a tensor if possible.

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

def validate(self, obj, value):
    """Validate input object, recasting it into a tensor if possible."""
    device = ComputeConfig.get_device()
    dtype = ComputeConfig.get_dtype() if self.dtype is None else self.dtype
    if isinstance(value, torch.Tensor):
        casted = value.to(dtype=dtype, device=device)
        return casted
    try:
        # Try casting value into a tensor, raise validation error if it fails
        return torch.tensor(value, dtype=dtype, device=device)
    except:  # noqa: E722
        self.error(obj, value)

from_string

from_string(s)

Casting tensor from string is not supported for now.

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

def from_string(self, s):
    """Casting tensor from string is not supported for now."""
    raise NotImplementedError

Crop modules (ML-Based)

diffwofost.ml_models.crop.partitioning.DVS_Partitioning_NN

Bases: _BaseDVSPartitioning

Drop-in wrapper for using a neural-network model as a partitioning component.

DVS_Partitioning_NN mirrors the interface of DVS_Partitioning so that a neural-network model can be plugged into Wofost72 anywhere the standard rule-based partitioning module is expected.

The wrapped model must accept DVS as input and return a :class:PartioningFactors namedtuple containing FR, FL, FS, and FO. The wrapper handles publishing these values to the kiosk and storing them in the same state layout as the physical partitioning implementation.

External dependencies:

Name	Description	Provided by	Unit
DVS	Crop development stage	DVS_Phenology	-

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	-

Methods:

• initialize –

  Initialize the DVS_Partitioning_NN simulation object.

• integrate –

  Update partitioning factors by running the network on the current DVS.

• calc_rates –

  Return the current partitioning factors.

initialize

initialize(day, kiosk, nn_model, shape=None)

Initialize the DVS_Partitioning_NN simulation object.

Parameters:

• day –

  Start date of the simulation.

• kiosk (VariableKiosk) –

  Variable kiosk of this PCSE instance.

• nn_model (Module) –

  Network that maps DVS → PartioningFactors.

• shape (tuple | Size | None, default: None ) –

  Target shape for batch simulations.

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

def initialize(self, day, kiosk, nn_model, shape=None):
    """Initialize the DVS_Partitioning_NN simulation object.

    Args:
        day: Start date of the simulation.
        kiosk (VariableKiosk): Variable kiosk of this PCSE instance.
        nn_model (torch.nn.Module): Network that maps DVS → PartioningFactors.
        shape (tuple | torch.Size | None): Target shape for batch simulations.
    """
    self._device = ComputeConfig.get_device()
    self._dtype = ComputeConfig.get_dtype()
    self.nn_model = nn_model
    self.kiosk = kiosk

    DVS = _broadcast_to(self.kiosk["DVS"], shape or ())
    pf = self.nn_model(DVS)
    self.states = self.StateVariables(
        kiosk,
        publish=["FR", "FL", "FS", "FO"],
        FR=pf.FR,
        FL=pf.FL,
        FS=pf.FS,
        FO=pf.FO,
        PF=PartioningFactors(pf.FR, pf.FL, pf.FS, pf.FO),
        shape=shape,
    )

integrate

integrate(day, delt=1.0)

Update partitioning factors by running the network on the current DVS.

Parameters:

• day –

  Current simulation day.

• delt (float, default: 1.0 ) –

  Integration step size. Included for interface compatibility and ignored in this implementation.

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

def integrate(self, day, delt=1.0):
    """Update partitioning factors by running the network on the current DVS.

    Args:
        day: Current simulation day.
        delt (float): Integration step size. Included for interface
            compatibility and ignored in this implementation.
    """
    DVS = _broadcast_to(self.kiosk["DVS"], self.states.shape)
    pf = self.nn_model(DVS)
    self._set_partitioning_states(pf.FR, pf.FL, pf.FS, pf.FO)

calc_rates

calc_rates(day, drv)

Return the current partitioning factors.

Parameters:

• day –

  Current simulation day.

• drv –

  Driving variables object. Included for interface compatibility.

Returns:

• PartioningFactors –

  Current partitioning factors stored in the state.

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

def calc_rates(self, day, drv):
    """Return the current partitioning factors.

    Args:
        day: Current simulation day.
        drv: Driving variables object. Included for interface compatibility.

    Returns:
        PartioningFactors: Current partitioning factors stored in the state.
    """
    return self.states.PF