diff --git a/cara/apps/expert.py b/cara/apps/expert.py index c6ee010e..c02c56df 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -288,8 +288,9 @@ baseline_model = models.Model( room=models.Room(volume=75), 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, + inside_temp=models.PiecewiseConstant((0,24),(293,)), + outside_temp=models.PiecewiseConstant((0,24),(283,)), + cd_b=0.6, window_height=1.6, opening_length=0.6, ), infected=models.InfectedPerson( virus=models.Virus.types['SARS_CoV_2'], diff --git a/cara/models.py b/cara/models.py index 95aabe08..805993ac 100644 --- a/cara/models.py +++ b/cara/models.py @@ -6,6 +6,34 @@ from abc import abstractmethod from dataclasses import dataclass +# average temperature of each month, hour per hour (from midnight to 11 pm) +Geneva_hourly_temperatures_celsius_per_hour = { + 'Jan': [0.2, -0.3, -0.5, -0.9, -1.1, -1.4, -1.5, -1.5, -1.1, 0.1, 1.5, + 2.8, 3.8, 4.4, 4.5, 4.4, 4.4, 3.9, 3.1, 2.7, 2.2, 1.7, 1.5, 1.1], + 'Feb': [0.9, 0.3, 0.0, -0.5, -0.7, -1.1, -1.2, -1.1, -0.7, 0.8, 2.5, + 4.2, 5.4, 6.2, 6.3, 6.2, 6.1, 5.5, 4.5, 4.1, 3.5, 2.8, 2.5, 2.0], + 'Mar': [4.2, 3.5, 3.1, 2.5, 2.1, 1.6, 1.5, 1.6, 2.2, 4.0, 6.3, 8.4, + 10.0, 11.1, 11.2, 11.1, 11.0, 10.2, 8.9, 8.3, 7.5, 6.7, 6.3, 5.6], + 'Apr': [7.4, 6.7, 6.2, 5.5, 5.2, 4.7, 4.5, 4.6, 5.3, 7.2, 9.6, 11.9, + 13.7, 14.8, 14.9, 14.8, 14.7, 13.8, 12.4, 11.8, 10.9, 10.1, 9.6, 8.9], + 'May': [11.8, 11.1, 10.6, 9.9, 9.5, 8.9, 8.8, 8.9, 9.6, 11.6, 14.2, 16.6, + 18.4, 19.6, 19.7, 19.6, 19.4, 18.6, 17.1, 16.5, 15.6, 14.6, 14.2, 13.4], + 'Jun': [15.2, 14.4, 13.9, 13.2, 12.7, 12.2, 12.0, 12.1, 12.8, 15.0, 17.7, + 20.2, 22.1, 23.3, 23.5, 23.4, 23.2, 22.3, 20.8, 20.1, 19.1, 18.2, 17.7, 16.9], + 'Jul': [17.6, 16.7, 16.1, 15.3, 14.9, 14.3, 14.1, 14.2, 15.0, 17.3, 20.2, + 23.0, 25.0, 26.3, 26.5, 26.4, 26.2, 25.2, 23.6, 22.8, 21.8, 20.8, 20.2, 19.4], + 'Aug': [17.1, 16.2, 15.7, 14.9, 14.5, 13.9, 13.7, 13.8, 14.6, 16.9, 19.7, + 22.4, 24.4, 25.6, 25.8, 25.7, 25.5, 24.5, 22.9, 22.2, 21.2, 20.2, 19.7, 18.9], + 'Sep': [13.4, 12.7, 12.2, 11.5, 11.2, 10.7, 10.5, 10.6, 11.3, 13.2, 15.6, + 17.9, 19.6, 20.8, 20.9, 20.8, 20.7, 19.8, 18.4, 17.8, 16.9, 16.1, 15.6, 14.9], + 'Oct': [9.4, 8.8, 8.5, 7.9, 7.6, 7.2, 7.1, 7.2, 7.7, 9.3, 11.2, 13.0, + 14.4, 15.3, 15.4, 15.3, 15.2, 14.5, 13.4, 12.9, 12.2, 11.6, 11.2, 10.6], + 'Nov': [4.0, 3.6, 3.3, 2.9, 2.6, 2.3, 2.2, 2.2, 2.7, 3.9, 5.5, 6.9, 8.0, + 8.7, 8.8, 8.7, 8.7, 8.1, 7.2, 6.8, 6.3, 5.7, 5.5, 5.0], + 'Dec': [1.4, 1.0, 0.8, 0.4, 0.2, -0.0, -0.1, -0.1, 0.3, 1.3, 2.6, 3.8, + 4.7, 5.2, 5.3, 5.2, 5.2, 4.7, 4.0, 3.7, 3.2, 2.8, 2.6, 2.2] + } + @dataclass(frozen=True) class Room: @@ -72,6 +100,50 @@ class PeriodicInterval(Interval): return tuple(result) +@dataclass(frozen=True) +class PiecewiseConstant: + + #: transition times at which the function changes value (hours). + transition_times: typing.Tuple[float, ...] + + #: values of the function between transitions + values: typing.Tuple[float, ...] + + def __post_init__(self): + if len(self.transition_times) != len(self.values)+1: + raise ValueError("transition_times should contain one more element than values") + if tuple(sorted(set(self.transition_times))) != self.transition_times: + raise ValueError("transition_times should not contain duplicated elements and should be sorted") + + def value(self,time) -> float: + if time <= self.transition_times[0]: + return self.values[0] + if time > self.transition_times[-1]: + return self.values[-1] + + for t1,t2,value in zip(self.transition_times[:-1], + self.transition_times[1:],self.values): + if time > t1 and time <= t2: + return value + + def interval(self) -> Interval: + # build an Interval object + present_times = [] + for t1,t2,value in zip(self.transition_times[:-1], + self.transition_times[1:],self.values): + if value: + present_times.append((t1,t2)) + return SpecificInterval(present_times=present_times) + + +# Geneva hourly temperatures as piecewise constant function (in Kelvin) +GenevaTemperatures = { + month: PiecewiseConstant(tuple(np.arange(25.)), + tuple(273.15+np.array(temperatures))) + for month,temperatures in Geneva_hourly_temperatures_celsius_per_hour.items() +} + + @dataclass(frozen=True) class Ventilation: """ @@ -87,6 +159,9 @@ class Ventilation: #: The times at which the air exchange is taking place. active: Interval + def transition_times(self): + return self.active.transition_times() + @abstractmethod def air_exchange(self, room: Room, time: float) -> float: """ @@ -106,8 +181,8 @@ 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) + inside_temp: PiecewiseConstant #: The temperature inside the room (Kelvin) + outside_temp: PiecewiseConstant #: The temperature outside of the window (Kelvin) window_height: float #: The height of the window @@ -115,6 +190,12 @@ class WindowOpening(Ventilation): cd_b: float = 0.6 #: Discharge coefficient: what portion effective area is used to exchange air (0 <= cd_b <= 1) + def transition_times(self): + transitions = super().transition_times() + transitions.update(self.inside_temp.interval().transition_times()) + transitions.update(self.outside_temp.interval().transition_times()) + return transitions + def air_exchange(self, room: Room, time: float) -> float: # If the window is shut, no air is being exchanged. if not self.active.triggered(time): @@ -122,7 +203,8 @@ class WindowOpening(Ventilation): # Reminder, no dependence on time in the resulting calculation. - temp_delta = abs(self.inside_temp - self.outside_temp) / self.outside_temp + temp_delta = abs(self.inside_temp.value(time) - + self.outside_temp.value(time)) / self.outside_temp.value(time) 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 @@ -321,7 +403,8 @@ class Model: """ state_change_times = set() state_change_times.update(self.infected.presence.transition_times()) - state_change_times.update(self.ventilation.active.transition_times()) + state_change_times.update(self.ventilation.transition_times()) + return sorted(state_change_times) def last_state_change(self, time: float): diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 82c3ac71..4069bacf 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -29,8 +29,9 @@ def baseline_model(): room=models.Room(volume=75), 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, + inside_temp=models.PiecewiseConstant((0,24),(293,)), + outside_temp=models.PiecewiseConstant((0,24),(283,)), + cd_b=0.6, window_height=1.6, opening_length=0.6, ), infected=models.InfectedPerson( virus=models.Virus.types['SARS_CoV_2'], @@ -50,8 +51,9 @@ def baseline_model(): def baseline_periodic_window(): 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, + inside_temp=models.PiecewiseConstant((0,24),(293,)), + outside_temp=models.PiecewiseConstant((0,24),(283,)), + cd_b=0.6, window_height=1.6, opening_length=0.6, ) @@ -148,3 +150,118 @@ def test_expiration_aerosols(): exp2 = models.Expiration((0.059, 0.0139, 0.751, 0.139), particle_sizes = (5.5e-4, 3.5e-4, 0.8e-4, 1.8e-4)) npt.assert_allclose(exp1.aerosols(mask), exp2.aerosols(mask), rtol=1e-5) + + +def test_piecewiseconstantfunction_wrongarguments(): + # number of values should be 1+number of transition times + pytest.raises(ValueError,models.PiecewiseConstant,(0,1),(0,0)) + pytest.raises(ValueError,models.PiecewiseConstant,(0,),(0,0)) + # two transition times cannot be equal + pytest.raises(ValueError,models.PiecewiseConstant,(0,2,2),(0,0)) + # unsorted transition times are not allowed + pytest.raises(ValueError,models.PiecewiseConstant,(2,0),(0,0)) + + +def test_piecewiseconstant(): + transition_times = (0,8,16,24) + values = (2,5,8) + fun = models.PiecewiseConstant(transition_times,values) + assert (fun.value(10) == 5) and (fun.value(20.5) == 8) and \ + (fun.value(8) == 2) and (fun.value(0) == 2) and \ + (fun.value(24) == 8) and (fun.value(-1) == 2) and (fun.value(25) == 8) + + +def test_constantfunction(): + transition_times = (0,24) + values = (20,) + fun = models.PiecewiseConstant(transition_times,values) + for t in [0,1,8,10,16,20.1,24]: + assert (fun.value(t) == 20) + + +def test_piecewiseconstant_vs_interval(): + transition_times = (0,8,16,24) + values = (0,1,0) + fun = models.PiecewiseConstant(transition_times,values) + interval = models.SpecificInterval(present_times=[(8,16)]) + assert interval.transition_times() == fun.interval().transition_times() + for t in [0,1,8,10,16,20.1,24]: + assert fun.interval().triggered(t) == interval.triggered(t) + + +def test_windowopening(): + tempOutside = models.PiecewiseConstant((0,10,24),(273.15,283.15)) + tempInside = models.PiecewiseConstant((0,24),(293.15,)) + w = models.WindowOpening(active=models.SpecificInterval([(0,24)]), + inside_temp=tempInside,outside_temp=tempOutside, + window_height=1.,opening_length=0.6) + npt.assert_allclose(w.air_exchange(models.Room(volume=68),16.), + 3.7393925,rtol=1e-5) + npt.assert_allclose(w.air_exchange(models.Room(volume=68),8.), + 5.3842316,rtol=1e-5) + + +def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),)): + model = models.Model( + room=models.Room(volume=75), + ventilation=models.WindowOpening( + active=models.SpecificInterval(intervals_open), + inside_temp=models.PiecewiseConstant((0,24),(293,)), + outside_temp=models.GenevaTemperatures[month], + cd_b=0.6, window_height=1.6, opening_length=0.6, + ), + infected=models.InfectedPerson( + virus=models.Virus.types['SARS_CoV_2'], + presence=models.SpecificInterval(((0, 4), (5, 7.5))), + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Light exercise'], + expiration=models.Expiration.types['Unmodulated Vocalization'], + ), + infected_occupants=1, + exposed_occupants=10, + exposed_activity=models.Activity.types['Light exercise'], + ) + return model + + +def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)): + model = models.Model( + room=models.Room(volume=75), + ventilation=models.WindowOpening( + active=models.SpecificInterval(intervals_open), + inside_temp=models.PiecewiseConstant((0,24),(293,)), + outside_temp=models.PiecewiseConstant((0,24),(outside_temp,)), + cd_b=0.6, window_height=1.6, opening_length=0.6, + ), + infected=models.InfectedPerson( + virus=models.Virus.types['SARS_CoV_2'], + presence=models.SpecificInterval(((0, 4), (5, 7.5))), + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Light exercise'], + expiration=models.Expiration.types['Unmodulated Vocalization'], + ), + infected_occupants=1, + exposed_occupants=10, + exposed_activity=models.Activity.types['Light exercise'], + ) + return model + + +def test_concentrations_hourly_dep_startup(): + # The concentrations should be the same up to 8 AM (time when the + # temperature changes DURING the window opening). + for month,temperatures in models.Geneva_hourly_temperatures_celsius_per_hour.items(): + m1 = build_hourly_dependent_model(month) + m2 = build_constant_temp_model(temperatures[7]+273.15) + for t in [0.5, 1.2, 2., 3.5, 5., 6.5, 7.5, 7.9, 8.]: + npt.assert_allclose(m1.concentration(t), m2.concentration(t), rtol=1e-5) + + +def test_concentrations_hourly_dep_adding_artificial_transitions(): + # Adding a second opening inside the first one should not change anything + for month,temperatures in models.Geneva_hourly_temperatures_celsius_per_hour.items(): + m1 = build_hourly_dependent_model(month,intervals_open=((7.5, 8.5),)) + m2 = build_hourly_dependent_model(month,intervals_open=((7.5, 8.5),(8.,8.1))) + for t in [0.5, 1.2, 2., 3.5, 5., 6.5, 7.5, 7.9, 8., 8.5, 9., 12.]: + npt.assert_allclose(m1.concentration(t), m2.concentration(t), rtol=1e-5) +