the get_predicted() (aka tendency) logic for Moist environments

[slayoo]

CC: @tluettm

The Moist environments feature a logic which defines two sets of thermodynamic state variables: "current" and "predicted".
These are swapped at the end of each timestep from Particulator._notify_observers() in which it is ensured that the environment.notify() is the last one called among all notify() methods:

PySDM/PySDM/particulator.py

Lines 49 to 60 in 4ff7b64

	def run(self, steps):
	for _ in range(steps):
	for key, dynamic in self.dynamics.items():
	with self.timers[key]:
	dynamic()
	self.n_steps += 1
	self._notify_observers()
	def _notify_observers(self):
	reversed_order_so_that_environment_is_last = reversed(self.observers)
	for observer in reversed_order_so_that_environment_is_last:
	observer.notify()

All environments inheriting from Moist (like Parcel, Kinematic_1D and Kinematic_2D) require that the AmbientThermodynamics "dynamic" is added to the builder upon Particulator instantiation. The AmbientThermodynamics.__call__() method is where the environment.sync() is called which is expected to update the thermodynamic state with new "forcings" (be it due to advection in multi-dimensional envs, or due to adiabatic ascent in the parcel env):

PySDM/PySDM/dynamics/ambient_thermodynamics.py

Lines 9 to 17 in 4ff7b64

	class AmbientThermodynamics:
	def __init__(self):
	self.particulator = None
	def register(self, builder):
	self.particulator = builder.particulator
	def __call__(self):
	self.particulator.environment.sync()

This updated state goes into the "predicted" dataset:

PySDM/PySDM/environments/impl/moist.py

Line 100 in 4ff7b64

self._values["predicted"] = target

So, in each timestep, the microphysics dynamics have access to both "current" and "predicted" thermodynamic states. The "current" state should not be modified, while modifying the "predicted" state effectively constitutes the definition of the thermodynamic tendencies:

PySDM/PySDM/particulator.py

Lines 119 to 120 in 4ff7b64

	This function should only change environment's predicted `thd` and
	`water_vapour_mixing_ratio` (and not `rhod`).

Multiple dynamics may modify the "predicted" state. Access to the "predicted" vars is provided with Moist.get_predicted(), while access to the "current" state is provided by Moist.__getitem__(). The swapping of these two states, which happens in Moist.notify() effectively leads to "application" of the accumulated tendencies:

PySDM/PySDM/environments/impl/moist.py

Lines 110 to 116 in 4ff7b64

	def notify(self):
	if self._values["predicted"] is None:
	return
	self._tmp = self._values["current"]
	self._values["current"] = self._values["predicted"]
	self._values["predicted"] = None

[tluettm]

Thanks a lot @slayoo

What I don't understand is this part in the implicit step within the condensation methods backend. Here 'water_vapour_mixing_ratio' is modified directly. It might be that the naming of 'water_vapour_mixing_ratio' in this case is misleading and is in fact the predicted variable or a temporary variable for sub stepping. I just want to confirm that the method is applied correctly here since its hard to follow function calls in this routine.

PySDM/PySDM/backends/impl_numba/methods/condensation_methods.py

Lines 257 to 355 in 4ff7b64

	def step_impl( # pylint: disable=too-many-arguments,too-many-locals
	attributes,
	cell_idx,
	thd,
	water_vapour_mixing_ratio,
	rhod,
	dthd_dt_pred,
	d_water_vapour_mixing_ratio__dt_predicted,
	drhod_dt,
	m_d,
	air_density,
	air_dynamic_viscosity,
	rtol_x,
	timestep,
	n_substeps,
	fake,
	):
	timestep /= n_substeps
	ml_old = calculate_ml_old(
	attributes.water_mass, attributes.multiplicity, cell_idx
	)
	count_activating, count_deactivating, count_ripening = 0, 0, 0
	RH_max = 0
	success = True
	for _ in range(n_substeps):
	# note: no example yet showing that the trapezoidal scheme brings any improvement
	thd += timestep * dthd_dt_pred / 2
	water_vapour_mixing_ratio += (
	timestep * d_water_vapour_mixing_ratio__dt_predicted / 2
	)
	rhod += timestep * drhod_dt / 2
	T = formulae.state_variable_triplet__T(rhod, thd)
	p = formulae.state_variable_triplet__p(
	rhod, T, water_vapour_mixing_ratio
	)
	pv = formulae.state_variable_triplet__pv(p, water_vapour_mixing_ratio)
	lv = formulae.latent_heat__lv(T)
	pvs = formulae.saturation_vapour_pressure__pvs_water(T)
	DTp = formulae.diffusion_thermics__D(T, p)
	RH = pv / pvs
	Sc = formulae.trivia__air_schmidt_number(
	dynamic_viscosity=air_dynamic_viscosity,
	diffusivity=DTp,
	density=air_density,
	)
	(
	ml_new,
	success_within_substep,
	n_activating,
	n_deactivating,
	n_ripening,
	) = calculate_ml_new(
	attributes,
	timestep,
	fake,
	T,
	p,
	RH,
	Sc,
	cell_idx,
	lv,
	pvs,
	DTp,
	formulae.diffusion_thermics__K(T, p),
	rtol_x,
	)
	dml_dt = (ml_new - ml_old) / timestep
	d_water_vapour_mixing_ratio__dt_corrected = -dml_dt / m_d
	dthd_dt_corr = formulae.state_variable_triplet__dthd_dt(
	rhod=rhod,
	thd=thd,
	T=T,
	d_water_vapour_mixing_ratio__dt=d_water_vapour_mixing_ratio__dt_corrected,
	lv=lv,
	)
	thd += timestep * (dthd_dt_pred / 2 + dthd_dt_corr)
	water_vapour_mixing_ratio += timestep * (
	d_water_vapour_mixing_ratio__dt_predicted / 2
	+ d_water_vapour_mixing_ratio__dt_corrected
	)
	rhod += timestep * drhod_dt / 2
	ml_old = ml_new
	count_activating += n_activating
	count_deactivating += n_deactivating
	count_ripening += n_ripening
	RH_max = max(RH_max, RH)
	success = success and success_within_substep
	return (
	water_vapour_mixing_ratio,
	thd,
	count_activating,
	count_deactivating,
	count_ripening,
	RH_max,
	success,
	)
	return step_impl

Where is the predicted vars set to the current vars? In Moist.notify() predicted and current are not swapped, instead predicted is set to None. Likely so that the function returns directly next call in case the predicted vars was not needed. But would it not make sense to make to set predicted here for the case that multiple dynamics want to accesses the predicted vars?

PySDM/PySDM/environments/impl/moist.py

Line 116 in 4ff7b64

self._values["predicted"] = None

[slayoo]

What I don't understand is this part in the implicit step within the condensation methods backend. Here 'water_vapour_mixing_ratio' is modified directly. It might be that the naming of 'water_vapour_mixing_ratio' in this case is misleading and is in fact the predicted variable or a temporary variable for sub stepping. I just want to confirm that the method is applied correctly here since its hard to follow function calls in this routine.

Modification of the predicted state by condensation happens here (the second value in the returned tuple gets written to the predicted state):

PySDM/PySDM/backends/impl_numba/methods/condensation_methods.py

Lines 139 to 148 in 4ff7b64

	(
	success[cell_id],
	cdt_predicted_water_vapour_mixing_ratio[cell_id],
	cdt_pthd[cell_id],
	cnt_n_substeps[cell_id],
	cnt_n_activating[cell_id],
	cnt_n_deactivating[cell_id],
	cnt_n_ripening[cell_id],
	RH_max[cell_id],
	) = solver(

Indeed, the variable naming is misleading, as what is returned is internally called "water_vapour_mixing_ratio" as you point out:

PySDM/PySDM/backends/impl_numba/methods/condensation_methods.py

Lines 695 to 697 in 4ff7b64

	return (
	success,
	water_vapour_mixing_ratio,

Where is the predicted vars set to the current vars? In Moist.notify() predicted and current are not swapped, instead predicted is set to None. Likely so that the function returns directly next call in case the predicted vars was not needed. But would it not make sense to make to set predicted here for the case that multiple dynamics want to accesses the predicted vars?

The point is that the notify logic ensures that environment.notiy() is called last, so no other dynamic sees the None. The first thing to do in the a timestep is to call AmbientThermodynamics.__call__() which populates the predicted values with data.