pingu_98/cara
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'develop/temperature_dependency' into 'master'

Browse source 
Implementation of outside temperature time dependency

Closes #15

See merge request cara/cara!19
This commit is contained in:
Philip James Elson 2020-11-05 11:43:36 +00:00
commit aadbb8c1dd
3 changed files with 211 additions and 10 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
cara/apps/expert.py
View file

@ -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'],

91
cara/models.py
View file

@ -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):

125
cara/tests/test_known_quantities.py
View file

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