Extract the concept of an interval for the Ventilation schemes, and correctly handle the time dependence of Ventilation in the concentration calculations.
This commit is contained in:
Phil Elson
parent
943a65b063
commit
a70f899b52
3 changed files with 134 additions and 58 deletions
27
cara/apps.py
27
cara/apps.py
|
|
@ -160,14 +160,14 @@ class WidgetView:
|
|||
return widget
|
||||
|
||||
def _build_window(self, node):
|
||||
period = widgets.IntSlider(value=node.period, min=0, max=240)
|
||||
interval = widgets.IntSlider(value=node.duration, min=0, max=240)
|
||||
period = widgets.IntSlider(value=node.active.period, min=0, max=240)
|
||||
interval = widgets.IntSlider(value=node.active.duration, min=0, max=240)
|
||||
|
||||
def on_period_change(change):
|
||||
node.period = change['new']
|
||||
node.active.period = change['new']
|
||||
|
||||
def on_interval_change(change):
|
||||
node.duration = change['new']
|
||||
node.active.duration = change['new']
|
||||
|
||||
# TODO: Link the state back to the widget, not just the other way around.
|
||||
period.observe(on_period_change, names=['value'])
|
||||
|
|
@ -181,14 +181,14 @@ class WidgetView:
|
|||
)
|
||||
|
||||
def _build_hepa(self, node):
|
||||
period = widgets.IntSlider(value=node.period, min=0, max=240)
|
||||
interval = widgets.IntSlider(value=node.duration, min=0, max=240)
|
||||
period = widgets.IntSlider(value=node.active.period, min=0, max=240, step=5)
|
||||
interval = widgets.IntSlider(value=node.active.duration, min=0, max=240, step=5)
|
||||
|
||||
def on_period_change(change):
|
||||
node.period = change['new']
|
||||
node.active.period = change['new']
|
||||
|
||||
def on_interval_change(change):
|
||||
node.duration = change['new']
|
||||
node.active.duration = change['new']
|
||||
|
||||
# TODO: Link the state back to the widget, not just the other way around.
|
||||
period.observe(on_period_change, names=['value'])
|
||||
|
|
@ -286,8 +286,9 @@ class WidgetView:
|
|||
|
||||
baseline_model = models.Model(
|
||||
room=models.Room(volume=75),
|
||||
ventilation=models.PeriodicWindow(
|
||||
period=120, duration=120, inside_temp=293, outside_temp=283, cd_b=0.6,
|
||||
ventilation=models.WindowOpening(
|
||||
active=models.PeriodicInterval(period=120, duration=120),
|
||||
inside_temp=293, outside_temp=283, cd_b=0.6,
|
||||
window_height=1.6, opening_length=0.6,
|
||||
),
|
||||
infected=models.InfectedPerson(
|
||||
|
|
@ -313,14 +314,14 @@ class CARAStateBuilder(state.StateBuilder):
|
|||
def build_type_Ventilation(self, _: dataclasses.Field):
|
||||
s = state.DataclassStateNamed(
|
||||
states={
|
||||
'Natural': self.build_generic(models.PeriodicWindow),
|
||||
'HEPA': self.build_generic(models.PeriodicHEPA),
|
||||
'Natural': self.build_generic(models.WindowOpening),
|
||||
'HEPA': self.build_generic(models.HEPAFilter),
|
||||
},
|
||||
state_builder=self,
|
||||
)
|
||||
# Initialise the HEPA state
|
||||
s._states['HEPA'].dcs_update_from(
|
||||
models.PeriodicHEPA(120, 120, 500.)
|
||||
models.HEPAFilter(models.PeriodicInterval(120, 120), 500.)
|
||||
)
|
||||
return s
|
||||
|
||||
|
|
|
|||
144
cara/models.py
144
cara/models.py
|
|
@ -13,11 +13,68 @@ class Room:
|
|||
volume: int
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Interval:
|
||||
"""
|
||||
Represents a collection of times in which a "thing" happens.
|
||||
|
||||
The "thing" may be when an action is taken, such as opening a window, or
|
||||
entering a room.
|
||||
|
||||
"""
|
||||
def triggered(self, time: float) -> bool:
|
||||
"""Whether the given time falls inside this interval."""
|
||||
return False
|
||||
|
||||
def boundaries(self) -> typing.Set[float]:
|
||||
"""Returns the edges of this interval."""
|
||||
return set()
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class PeriodicInterval(Interval):
|
||||
#: How often does the interval occur.
|
||||
period: int
|
||||
|
||||
#: How long does the interval occur for.
|
||||
#: A value greater than :data:`period` signifies the event is permanently
|
||||
#: occurring, a value of 0 signifies that the event never happens.
|
||||
duration: int
|
||||
|
||||
def triggered(self, time: float) -> bool:
|
||||
period = self.period / 60.
|
||||
duration = self.duration / 60.
|
||||
return (time % period) >= (period - duration)
|
||||
|
||||
def boundaries(self) -> typing.Set[float]:
|
||||
state_changes = set()
|
||||
period_h = self.period / 60
|
||||
duration_h = self.duration / 60
|
||||
# The current implementation starts at the end of the first period, not at 0.
|
||||
start = period_h
|
||||
# Take as many steps as we need to get to a full day.
|
||||
for i in np.arange(start, 24, period_h):
|
||||
state_changes.add(i)
|
||||
if duration_h < period_h:
|
||||
state_changes.add(i - duration_h)
|
||||
|
||||
return state_changes
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Ventilation:
|
||||
"""
|
||||
An abstract class for ventilation schemes
|
||||
Represents a mechanism by which air can be exchanged (replaced/filtered)
|
||||
in a time dependent manner.
|
||||
|
||||
The nature of the various air exchange schemes means that it is expected
|
||||
for subclasses of Ventilation to exist. Known subclasses include
|
||||
WindowOpening for window based air exchange, and HEPAFilter, for
|
||||
mechanical air exchange through a filter.
|
||||
|
||||
"""
|
||||
#: The times at which the air exchange is taking place.
|
||||
active: Interval
|
||||
|
||||
@abstractmethod
|
||||
def air_exchange(self, room: Room, time: float) -> float:
|
||||
|
|
@ -26,15 +83,20 @@ class Ventilation:
|
|||
cubic meter at a given time (in hours).
|
||||
|
||||
"""
|
||||
pass
|
||||
return 0.
|
||||
|
||||
def times_of_state_change(self) -> typing.Set[float]:
|
||||
"""
|
||||
Returns the times at which a change in ventilation occurs.
|
||||
|
||||
"""
|
||||
return self.active.boundaries()
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class PeriodicWindow(Ventilation):
|
||||
|
||||
# The window is opened for <duration> minutes every <period> minutes
|
||||
period: int #: How often the window is opened (minutes)
|
||||
duration: int #: How long the window remains opened for (minutes)
|
||||
class WindowOpening(Ventilation):
|
||||
#: The interval in which the window is open.
|
||||
active: Interval
|
||||
|
||||
inside_temp: float #: The temperature inside the room (Kelvin)
|
||||
outside_temp: float #: The temperature outside of the window (Kelvin)
|
||||
|
|
@ -46,36 +108,27 @@ class PeriodicWindow(Ventilation):
|
|||
cd_b: float = 0.6 #: Discharge coefficient: what portion effective area is used to exchange air (0 <= cd_b <= 1)
|
||||
|
||||
def air_exchange(self, room: Room, time: float) -> float:
|
||||
period = self.period / 60.
|
||||
duration = self.duration / 60.
|
||||
# Returns the rate at which air is being exchanged in the given room per cubic meter at a given time
|
||||
# If the window is closed, no air is being exchanged
|
||||
if (time % period) < (period - duration):
|
||||
return 0
|
||||
# If the window is shut, no air is being exchanged.
|
||||
if not self.active.triggered(time):
|
||||
return 0.
|
||||
|
||||
root = np.sqrt(9.81 * self.window_height * (abs(self.inside_temp - self.outside_temp)) / self.outside_temp)
|
||||
temp_delta = abs(self.inside_temp - self.outside_temp) / self.outside_temp
|
||||
root = np.sqrt(9.81 * self.window_height * temp_delta)
|
||||
|
||||
return (3600 / (3 * room.volume)) * self.cd_b * self.window_height * self.opening_length * root
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class PeriodicHEPA(Ventilation):
|
||||
|
||||
# The HEPA is switched on for <duration> minutes every <period> minutes
|
||||
period: int #: How often the HEPA is switched on (minutes)
|
||||
duration: int #: How long the HEPA remains switched on for (minutes)
|
||||
class HEPAFilter(Ventilation):
|
||||
#: The interval in which the HEPA filter is operating.
|
||||
active: Interval
|
||||
|
||||
q_air_mech: float #: The rate at which the HEPA exchanges air (when switched on)
|
||||
|
||||
def air_exchange(self, room: Room, time: float) -> float:
|
||||
# Returns the rate at which air is being exchanged in the given room per cubic meter at a given time
|
||||
|
||||
period = self.period / 60.
|
||||
duration = self.duration / 60.
|
||||
|
||||
# If the HEPA is off, no air is being exchanged
|
||||
if (time % period) < (period - duration):
|
||||
return 0
|
||||
# If the HEPA is off, no air is being exchanged.
|
||||
if not self.active.triggered(time):
|
||||
return 0.
|
||||
|
||||
return self.q_air_mech / room.volume
|
||||
|
||||
|
|
@ -266,12 +319,29 @@ class Model:
|
|||
|
||||
return k + self.virus.decay_constant + self.ventilation.air_exchange(self.room, time)
|
||||
|
||||
@functools.lru_cache()
|
||||
def collect_time_state_changes(self):
|
||||
"""
|
||||
All time dependent entities on this model must provide information about
|
||||
the times at which their state changes.
|
||||
|
||||
"""
|
||||
state_change_times = set()
|
||||
for start, end in self.infected.present_times:
|
||||
state_change_times.add(start)
|
||||
state_change_times.add(end)
|
||||
state_change_times.update(self.ventilation.times_of_state_change())
|
||||
return sorted(state_change_times)
|
||||
|
||||
def last_state_change(self, time: float):
|
||||
"""
|
||||
Find the most recent state change.
|
||||
|
||||
"""
|
||||
for change_time in self.collect_time_state_changes()[::-1]:
|
||||
if time >= change_time:
|
||||
return change_time
|
||||
return 0
|
||||
|
||||
@functools.lru_cache()
|
||||
def concentration(self, time: float) -> float:
|
||||
|
|
@ -280,20 +350,17 @@ class Model:
|
|||
V = self.room.volume
|
||||
Ni = self.infected_occupants
|
||||
ER = self.infected.emission_rate(time)
|
||||
t0, t1 = self.infected.start_end_of_presence(time)
|
||||
|
||||
if t == 0:
|
||||
return 0.0
|
||||
elif t0 < t <= t1:
|
||||
# Concentration while infected present.
|
||||
init_concentration = self.concentration(t0)
|
||||
fac = np.exp(IVRR * (t0 - t))
|
||||
return ((ER * Ni) / (IVRR * V)) * (1 - fac) + init_concentration * fac
|
||||
t_last_state_change = self.last_state_change(time)
|
||||
if t != t_last_state_change:
|
||||
concentration_at_last_state_change = self.concentration(t_last_state_change)
|
||||
else:
|
||||
# Concentration while infected not present.
|
||||
end_concentration = self.concentration(t1)
|
||||
fac = np.exp(IVRR * (t1 - t))
|
||||
return end_concentration * fac
|
||||
concentration_at_last_state_change = self.concentration(t - 0.0001)
|
||||
|
||||
fac = np.exp(IVRR * (t_last_state_change - t))
|
||||
concentration_limit = (ER * Ni) / (IVRR * V)
|
||||
return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac
|
||||
|
||||
def infection_probability(self):
|
||||
# Infection probability
|
||||
|
|
@ -317,4 +384,3 @@ class Model:
|
|||
|
||||
# Probability of infection.
|
||||
return (1 - np.exp(-inf_aero)) * 100
|
||||
|
||||
|
|
|
|||
|
|
@ -26,8 +26,11 @@ def test_no_mask_emission_rate(baseline_model):
|
|||
def baseline_model():
|
||||
model = models.Model(
|
||||
room=models.Room(volume=75),
|
||||
ventilation=models.PeriodicWindow(period=120, duration=120, inside_temp=293, outside_temp=283, cd_b=0.6,
|
||||
window_height=1.6, opening_length=0.6),
|
||||
ventilation=models.WindowOpening(
|
||||
active=models.PeriodicInterval(period=120, duration=120),
|
||||
inside_temp=293, outside_temp=283, cd_b=0.6,
|
||||
window_height=1.6, opening_length=0.6,
|
||||
),
|
||||
infected=models.InfectedPerson(
|
||||
virus=models.Virus.types['SARS_CoV_2'],
|
||||
present_times=((0, 4), (5, 8)),
|
||||
|
|
@ -44,8 +47,11 @@ def baseline_model():
|
|||
|
||||
@pytest.fixture
|
||||
def baseline_periodic_window():
|
||||
return models.PeriodicWindow(period=120, duration=15, inside_temp=293, outside_temp=283, cd_b=0.6,
|
||||
window_height=1.6, opening_length=0.6)
|
||||
return models.WindowOpening(
|
||||
active=models.PeriodicInterval(period=120, duration=15),
|
||||
inside_temp=293, outside_temp=283, cd_b=0.6,
|
||||
window_height=1.6, opening_length=0.6,
|
||||
)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
|
|
@ -55,7 +61,10 @@ def baseline_room():
|
|||
|
||||
@pytest.fixture
|
||||
def baseline_periodic_hepa():
|
||||
return models.PeriodicHEPA(period=120, duration=15, q_air_mech=514.74)
|
||||
return models.HEPAFilter(
|
||||
active=models.PeriodicInterval(period=120, duration=15),
|
||||
q_air_mech=514.74,
|
||||
)
|
||||
|
||||
|
||||
def test_concentrations(baseline_model):
|
||||
|
|
@ -63,7 +72,7 @@ def test_concentrations(baseline_model):
|
|||
concentrations = [baseline_model.concentration(t) for t in ts]
|
||||
npt.assert_allclose(
|
||||
concentrations,
|
||||
[0.000000e+00, 2.891970e-01, 1.266287e-04, 2.891969e-01, 5.544607e-08],
|
||||
[0.000000e+00, 2.891970e-01, 1.267267e-04, 2.891969e-01, 5.548897e-08],
|
||||
rtol=1e-5
|
||||
)
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue