diff --git a/docs/source/usage/how_to_run.rst b/docs/source/usage/how_to_run.rst index e2a9e31..1519dfe 100644 --- a/docs/source/usage/how_to_run.rst +++ b/docs/source/usage/how_to_run.rst @@ -73,6 +73,16 @@ The following are inputs for the overall simulation: * ``agent.air_travel_int`` (`integer`, default ``-1``) The number of time steps between air travel events. Set to -1 to disable all air travel events. Currently this is implemented only for ``ic_type = census``. +* ``agent.weather_int`` (`integer`, default ``-1``) + The number of time steps between weather data updates. Set to -1 to disable the weather-dependent transmission model. + When enabled, ``agent.weather_filename`` and ``agent.startdate`` must also be provided. +* ``agent.weather_filename`` (`string`) + Path to the CSV file containing weather data (temperature and humidity by county and week). + Required when ``agent.weather_int > 0``. Example data files are provided in ``ExaEpi/data/``. +* ``agent.startdate`` (`string`) + Start date of the simulation in ``YYYY-MM-DD`` format (e.g. ``2020-01-01``). + Used to align the weather data with the simulation timeline. + Required when ``agent.weather_int > 0``. * ``agent.aggregated_diag_int`` (`integer`, default ``-1``) The number of time steps between writing aggregated data, for example wastewater data. Set to -1 to disable writing. * ``agent.aggregated_diag_prefix`` (`string`, default ``cases``) @@ -254,6 +264,29 @@ can be specified as follows: where ``[disease name]`` is any of the names specified in ``agent.disease_names`` (or the default value), and ``[key]`` is any of the parameters listed above. +The following inputs configure the weather-dependent transmission model. These parameters are only +used when ``agent.weather_int > 0``. The model computes a weather-based scale factor for the +transmission probability as a function of near-surface air temperature :math:`T` (°C) and absolute +humidity :math:`AH` (g/m³): + +.. math:: + + p(T,\,AH) = p_{\max} \exp\!\bigl(-\beta_{AH}\,AH - \alpha_T\,\max(0,\,T - T_0)\bigr) + +* ``weather_transmission.p_max`` (`float`, default ``1.0``) + Maximum transmission scale factor (dimensionless). At low temperature and low humidity the + scale factor approaches this value. It is multiplied by the baseline transmission probability + ``disease.p_trans`` to give the effective transmission probability. +* ``weather_transmission.beta_AH`` (`float`, default ``0.18``) + Absolute-humidity sensitivity coefficient in units of m³/g. Larger values cause transmission + to decrease more rapidly as humidity increases. +* ``weather_transmission.T0`` (`float`, default ``5.0``) + Temperature threshold in °C above which the temperature penalty is applied. Below this + temperature the model applies no temperature-dependent reduction. +* ``weather_transmission.alpha_T`` (`float`, default ``0.015``) + Temperature sensitivity coefficient in units of °C⁻¹. Larger values cause transmission + to decrease more rapidly above ``T0``. + In addition to the ExaEpi inputs, there are also a number of runtime options that can be configured for AMReX itself. Please see `__ for more information on these options. diff --git a/src/AgentContainer.H b/src/AgentContainer.H index 6483dda..21cbf6d 100644 --- a/src/AgentContainer.H +++ b/src/AgentContainer.H @@ -26,6 +26,7 @@ #include "InteractionModelLibrary.H" #include "Utils.H" #include "WeatherData.H" +#include "WeatherTransmissionModel.H" /*! \brief Derived class from ParticleContainer that defines agents and their functions */ class AgentContainer : public amrex::ParticleContainer<0, 0, RealIdx::nattribs, IntIdx::nattribs> { @@ -131,6 +132,11 @@ class AgentContainer : public amrex::ParticleContainer<0, 0, RealIdx::nattribs, int activeWeatherWeek; ActiveWeather* awd; + /*! Weather-dependent transmission probability model */ + WeatherTransmissionModel m_weather_model; + /*! Controls weather coupling: <= 0 disables weather scaling, > 0 enables it */ + int m_weather_int = -1; + protected: amrex::Real m_shelter_compliance = 0.95_rt; /*!< Shelter-in-place compliance rate */ diff --git a/src/AgentContainer.cpp b/src/AgentContainer.cpp index 3cd240d..5d4e4e6 100644 --- a/src/AgentContainer.cpp +++ b/src/AgentContainer.cpp @@ -101,6 +101,13 @@ AgentContainer::AgentContainer (const amrex::Geometry& a_geom, addAttributes(); + awd = nullptr; + { + amrex::ParmParse pp("agent"); + pp.query("weather_int", m_weather_int); + } + m_weather_model.readInputs("weather_transmission"); + amrex::ParmParse pp("agent"); pp.query("shelter_compliance", m_shelter_compliance); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 355409d..d865a6b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,7 +32,9 @@ set(_sources UrbanPopData.H UrbanPopData.cpp Utils.H - Utils.cpp) + Utils.cpp + WeatherTransmissionModel.H + WeatherTransmissionModel.cpp) # List of input files set(_input_files) diff --git a/src/InteractionModHome.H b/src/InteractionModHome.H index bc18d28..5226d51 100644 --- a/src/InteractionModHome.H +++ b/src/InteractionModHome.H @@ -177,7 +177,10 @@ void InteractionModHome::fastInteractHome (PCType& agent auto prob_ptr = soa.GetRealData(RealIdx::nattribs + r0(d) + RealIdxDisease::prob).data(); auto lparm = agents.getDiseaseParameters_d(d); auto lparm_h = agents.getDiseaseParameters_h(d); - Real scale = 1.0_rt; // TODO this should vary based on cell + auto weatherIdxPtr = soa.GetIntData(IntIdx::weatherLookup).data(); + auto varVec = agents.awd != nullptr ? agents.awd->varVec.data() : nullptr; + auto weather_model = agents.m_weather_model; + int a_weather_int = agents.m_weather_int; Real infect = (1.0_rt - lparm_h->vac_eff); // loop to count infectious agents in each group auto d_ptr = getCommunityIndex.comm_to_local_index_d.data(); @@ -235,8 +238,14 @@ void InteractionModHome::fastInteractHome (PCType& agent int family_i = community * max_family + family_ptr[i]; int num_infected_family = infected_family_d_ptr[family_i]; int num_infected_family_asymp = infected_family_asymp_d_ptr[family_i]; - Real family_prob = 1.0_rt - infect * xmit_family_prob * scale; - Real family_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * xmit_family_prob * scale; + Real weather_scale = 1.0_rt; + if (a_weather_int > 0 && varVec != nullptr && weatherIdxPtr[i] >= 0) { + weatherVars wvar = varVec[weatherIdxPtr[i]]; + weather_scale = weather_model(wvar); + } + Real family_prob = 1.0_rt - infect * xmit_family_prob * weather_scale; + Real family_prob_asymp = + 1.0_rt - lparm->asymp_relative_inf * infect * xmit_family_prob * weather_scale; prob_ptr[i] *= static_cast(std::pow(family_prob, num_infected_family)); prob_ptr[i] *= static_cast(std::pow(family_prob_asymp, num_infected_family_asymp)); if (!ptd.m_idata[IntIdx::withdrawn][i]) { @@ -245,7 +254,7 @@ void InteractionModHome::fastInteractHome (PCType& agent int cluster_i = community * max_hcs + hh_cluster_ptr[i]; int num_infected_nc = infected_nc_d_ptr[cluster_i] - num_infected_family_not_withdrawn; AMREX_ALWAYS_ASSERT(num_infected_nc >= 0); - Real nc_prob = 1.0_rt - infect * xmit_nc_prob * scale; + Real nc_prob = 1.0_rt - infect * xmit_nc_prob * weather_scale; prob_ptr[i] *= static_cast(std::pow(nc_prob, num_infected_nc)); // same for asymptomatic int num_infected_family_not_withdrawn_asymp = infected_family_not_withdrawn_asymp_d_ptr[family_i]; @@ -253,7 +262,7 @@ void InteractionModHome::fastInteractHome (PCType& agent int num_infected_nc_asymp = infected_nc_asymp_d_ptr[cluster_i] - num_infected_family_not_withdrawn_asymp; AMREX_ALWAYS_ASSERT(num_infected_nc_asymp >= 0); - Real nc_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * xmit_nc_prob * scale; + Real nc_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * xmit_nc_prob * weather_scale; prob_ptr[i] *= static_cast(std::pow(nc_prob_asymp, num_infected_nc_asymp)); } } diff --git a/src/InteractionModHomeNborhood.H b/src/InteractionModHomeNborhood.H index 0d64f96..b6240ba 100644 --- a/src/InteractionModHomeNborhood.H +++ b/src/InteractionModHomeNborhood.H @@ -138,7 +138,10 @@ void InteractionModHomeNborhood::fastInteractHomeNborhoo auto prob_ptr = soa.GetRealData(RealIdx::nattribs + r0(d) + RealIdxDisease::prob).data(); auto lparm = agents.getDiseaseParameters_d(d); auto lparm_h = agents.getDiseaseParameters_h(d); - Real scale = 1.0_rt; // TODO this should vary based on cell + auto weatherIdxPtr = soa.GetIntData(IntIdx::weatherLookup).data(); + auto varVec = agents.awd != nullptr ? agents.awd->varVec.data() : nullptr; + auto weather_model = agents.m_weather_model; + int a_weather_int = agents.m_weather_int; Real infect = 1.0_rt - lparm_h->vac_eff; auto d_ptr = getCommunityIndex.comm_to_local_index_d.data(); @@ -176,15 +179,23 @@ void InteractionModHomeNborhood::fastInteractHomeNborhoo infected_community_asymp_d_ptr[community] - num_infected_nborhood_asymp; AMREX_ALWAYS_ASSERT(num_infected_community >= 0); - Real comm_prob = 1.0_prt - infect * lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * scale; + Real weather_scale = 1.0_rt; + if (a_weather_int > 0 && varVec != nullptr && weatherIdxPtr[i] >= 0) { + weatherVars wvar = varVec[weatherIdxPtr[i]]; + weather_scale = weather_model(wvar); + } + Real comm_prob = 1.0_prt - infect * lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * weather_scale; Real comm_prob_asymp = 1.0_prt - lparm->asymp_relative_inf * infect * - lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * scale; + lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * + weather_scale; prob_ptr[i] *= static_cast(std::pow(comm_prob, num_infected_community)); prob_ptr[i] *= static_cast(std::pow(comm_prob_asymp, num_infected_community_asymp)); - Real nborhood_prob = 1.0_prt - infect * lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * scale; + Real nborhood_prob = + 1.0_prt - infect * lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * weather_scale; Real nborhood_prob_asymp = 1.0_prt - lparm->asymp_relative_inf * infect * - lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * scale; + lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * + weather_scale; prob_ptr[i] *= static_cast(std::pow(nborhood_prob, num_infected_nborhood)); prob_ptr[i] *= static_cast(std::pow(nborhood_prob_asymp, num_infected_nborhood_asymp)); } diff --git a/src/InteractionModSchool.H b/src/InteractionModSchool.H index 6da658c..9dee903 100644 --- a/src/InteractionModSchool.H +++ b/src/InteractionModSchool.H @@ -170,7 +170,10 @@ void InteractionModSchool::fastInteractSchool (PCType& a auto prob_ptr = soa.GetRealData(RealIdx::nattribs + r0(d) + RealIdxDisease::prob).data(); auto lparm = agents.getDiseaseParameters_d(d); auto lparm_h = agents.getDiseaseParameters_h(d); - Real scale = 1.0_rt; // TODO this should vary based on cell + auto weatherIdxPtr = soa.GetIntData(IntIdx::weatherLookup).data(); + auto varVec = agents.awd != nullptr ? agents.awd->varVec.data() : nullptr; + auto weather_model = agents.m_weather_model; + int a_weather_int = agents.m_weather_int; Real infect = (1.0_rt - lparm_h->vac_eff); auto d_ptr = getCommunityIndex.comm_to_local_index_d.data(); @@ -205,17 +208,23 @@ void InteractionModSchool::fastInteractSchool (PCType& a auto tidx = getTileIndex(iv, valid_box, true, bin_size, tbx); auto community = d_ptr[tidx]; int pos = (community * max_school_id + school_id_ptr[i]) * max_school_grade + school_grade_ptr[i]; + Real weather_scale = 1.0_rt; + if (a_weather_int > 0 && varVec != nullptr && weatherIdxPtr[i] >= 0) { + weatherVars wvar = varVec[weatherIdxPtr[i]]; + weather_scale = weather_model(wvar); + } if (getSchoolType(school_grade_ptr[i]) == SchoolType::daycare) { if (adults && age_group_ptr[i] > AgeGroups::a5to17) { // adults-adults interaction accounted for // in InteractionModWork.H // do nothing! } else { int num_infected_daycare = infected_daycare_d_ptr[pos]; - Real daycare_prob = 1.0_rt - infect * lparm->xmit_school[SchoolType::daycare] * scale; + Real daycare_prob = 1.0_rt - infect * lparm->xmit_school[SchoolType::daycare] * weather_scale; prob_ptr[i] *= static_cast(std::pow(daycare_prob, num_infected_daycare)); int num_infected_daycare_asymp = infected_daycare_asymp_d_ptr[pos]; Real daycare_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * - lparm->xmit_school[SchoolType::daycare] * scale; + lparm->xmit_school[SchoolType::daycare] * + weather_scale; prob_ptr[i] *= static_cast(std::pow(daycare_prob_asymp, num_infected_daycare_asymp)); } @@ -238,9 +247,9 @@ void InteractionModSchool::fastInteractSchool (PCType& a xmit = lparm->xmit_school_c2a[getSchoolType(school_grade_ptr[i])]; } } - Real school_prob = 1.0_rt - infect * xmit * scale; + Real school_prob = 1.0_rt - infect * xmit * weather_scale; prob_ptr[i] *= static_cast(std::pow(school_prob, num_infected_school)); - Real school_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * xmit * scale; + Real school_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * xmit * weather_scale; prob_ptr[i] *= static_cast(std::pow(school_prob_asymp, num_infected_school_asymp)); } } diff --git a/src/InteractionModWork.H b/src/InteractionModWork.H index 166e48c..db4ead6 100644 --- a/src/InteractionModWork.H +++ b/src/InteractionModWork.H @@ -137,7 +137,10 @@ void InteractionModWork::fastInteractWork (PCType& agent auto prob_ptr = soa.GetRealData(RealIdx::nattribs + r0(d) + RealIdxDisease::prob).data(); auto lparm = agents.getDiseaseParameters_d(d); auto lparm_h = agents.getDiseaseParameters_h(d); - Real scale = 1.0_rt; // TODO this should vary based on cell + auto weatherIdxPtr = soa.GetIntData(IntIdx::weatherLookup).data(); + auto varVec = agents.awd != nullptr ? agents.awd->varVec.data() : nullptr; + auto weather_model = agents.m_weather_model; + int a_weather_int = agents.m_weather_int; Real infect = 1.0_rt - lparm_h->vac_eff; auto d_ptr = getCommunityIndex.comm_to_local_index_d.data(); @@ -165,10 +168,16 @@ void InteractionModWork::fastInteractWork (PCType& agent auto community = d_ptr[tidx]; int wgroup_i = (community * max_workgroup + workgroup_ptr[i]) * max_naics + naics_ptr[i]; int num_infected_workgroup = infected_workgroup_d_ptr[wgroup_i]; - Real workgroup_prob = 1.0_prt - infect * lparm->xmit_work * scale; + Real weather_scale = 1.0_rt; + if (a_weather_int > 0 && varVec != nullptr && weatherIdxPtr[i] >= 0) { + weatherVars wvar = varVec[weatherIdxPtr[i]]; + weather_scale = weather_model(wvar); + } + Real workgroup_prob = 1.0_prt - infect * lparm->xmit_work * weather_scale; prob_ptr[i] *= static_cast(std::pow(workgroup_prob, num_infected_workgroup)); int num_infected_workgroup_asymp = infected_workgroup_asymp_d_ptr[wgroup_i]; - Real workgroup_prob_asymp = 1.0_prt - lparm->asymp_relative_inf * infect * lparm->xmit_work * scale; + Real workgroup_prob_asymp = + 1.0_prt - lparm->asymp_relative_inf * infect * lparm->xmit_work * weather_scale; prob_ptr[i] *= static_cast(std::pow(workgroup_prob_asymp, num_infected_workgroup_asymp)); } }); diff --git a/src/InteractionModWorkNborhood.H b/src/InteractionModWorkNborhood.H index ded19a0..f94accd 100644 --- a/src/InteractionModWorkNborhood.H +++ b/src/InteractionModWorkNborhood.H @@ -144,7 +144,10 @@ void InteractionModWorkNborhood::fastInteractWorkNborhoo auto prob_ptr = soa.GetRealData(RealIdx::nattribs + r0(d) + RealIdxDisease::prob).data(); auto lparm = agents.getDiseaseParameters_d(d); auto lparm_h = agents.getDiseaseParameters_h(d); - Real scale = 1.0_rt; // TODO this should vary based on cell + auto weatherIdxPtr = soa.GetIntData(IntIdx::weatherLookup).data(); + auto varVec = agents.awd != nullptr ? agents.awd->varVec.data() : nullptr; + auto weather_model = agents.m_weather_model; + int a_weather_int = agents.m_weather_int; Real infect = (1.0_rt - lparm_h->vac_eff); auto d_ptr = getCommunityIndex.comm_to_local_index_d.data(); @@ -180,17 +183,25 @@ void InteractionModWorkNborhood::fastInteractWorkNborhoo int num_infected_nborhood_asymp = infected_nborhood_asymp_d_ptr[community * max_nborhood + nborhood]; int num_infected_community_asymp = infected_community_asymp_d_ptr[community]; AMREX_ALWAYS_ASSERT(num_infected_community >= num_infected_nborhood); - Real comm_prob = 1.0_rt - infect * lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * scale; + Real weather_scale = 1.0_rt; + if (a_weather_int > 0 && varVec != nullptr && weatherIdxPtr[i] >= 0) { + weatherVars wvar = varVec[weatherIdxPtr[i]]; + weather_scale = weather_model(wvar); + } + Real comm_prob = 1.0_rt - infect * lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * weather_scale; prob_ptr[i] *= static_cast(std::pow(comm_prob, num_infected_community - num_infected_nborhood)); - Real nborhood_prob = 1.0_rt - infect * lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * scale; + Real nborhood_prob = + 1.0_rt - infect * lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * weather_scale; prob_ptr[i] *= static_cast(std::pow(nborhood_prob, num_infected_nborhood)); Real comm_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * - lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * scale; + lparm->xmit_comm[ptd.m_idata[IntIdx::age_group][i]] * + weather_scale; prob_ptr[i] *= static_cast( std::pow(comm_prob_asymp, num_infected_community_asymp - num_infected_nborhood_asymp)); Real nborhood_prob_asymp = 1.0_rt - lparm->asymp_relative_inf * infect * - lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * scale; + lparm->xmit_hood[ptd.m_idata[IntIdx::age_group][i]] * + weather_scale; prob_ptr[i] *= static_cast(std::pow(nborhood_prob_asymp, num_infected_nborhood_asymp)); } }); diff --git a/src/InteractionModel.H b/src/InteractionModel.H index 7b766b2..82b1d60 100644 --- a/src/InteractionModel.H +++ b/src/InteractionModel.H @@ -8,6 +8,7 @@ #include #include "AgentDefinitions.H" +#include "WeatherData.H" #include #include #include @@ -148,6 +149,11 @@ void interactAgentsImpl (IModel& interaction_model, AgentContainer& agents, int auto& soa = ptile.GetStructOfArrays(); auto num_tiles = numTilesInBox(mfi.tilebox(), true, bin_size); + auto weatherIdxPtr = soa.GetIntData(IntIdx::weatherLookup).data(); + auto varVec = agents.awd != nullptr ? agents.awd->varVec.data() : nullptr; + auto weather_model = agents.m_weather_model; + int a_weather_int = agents.m_weather_int; + int max_group = 0; if (bin_idx != -1) { max_group = agents.getMaxGroup(bin_idx); @@ -182,7 +188,6 @@ void interactAgentsImpl (IModel& interaction_model, AgentContainer& agents, int auto prob_ptr = soa.GetRealData(RealIdx::nattribs + r0(d) + RealIdxDisease::prob).data(); auto lparm = agents.getDiseaseParameters_d(d); auto lparm_h = agents.getDiseaseParameters_h(d); - Real scale = 1.0_prt; // TODO this should vary based on cell Real infect = (1.0_rt - lparm_h->vac_eff); Real asymp_relative_inf = lparm_h->asymp_relative_inf; @@ -205,9 +210,14 @@ void interactAgentsImpl (IModel& interaction_model, AgentContainer& agents, int AMREX_ALWAYS_ASSERT((Long)susceptible_i < np); if (infectious_i != susceptible_i && isSusceptible(susceptible_i, ptd, d) && isCandidate(susceptible_i, ptd)) { - ParticleReal prob = - 1.0_prt - - relative_inf * infect * binaryInteraction(infectious_i, susceptible_i, ptd, lparm, scale); + Real weather_scale = 1.0_prt; + if (a_weather_int > 0 && varVec != nullptr && weatherIdxPtr[susceptible_i] >= 0) { + weatherVars wvar = varVec[weatherIdxPtr[susceptible_i]]; + weather_scale = weather_model(wvar); + } + ParticleReal prob = 1.0_prt - relative_inf * infect * + binaryInteraction(infectious_i, susceptible_i, ptd, lparm, + weather_scale); // The atomic operation is needed because we find all the susceptible for each infectious in turn. // It can be eliminated by switching the order of infectious and susceptible. Gpu::Atomic::Multiply(&prob_ptr[susceptible_i], prob); diff --git a/src/WeatherData.H b/src/WeatherData.H index 5048591..f4a2741 100644 --- a/src/WeatherData.H +++ b/src/WeatherData.H @@ -75,8 +75,9 @@ struct ActiveWeather { return simWeek * numUnitsWithDataOnThisProc + lunit; } weatherVars* lookupWeatherVars (int unit, int simWeek) { - if (computeIndex(unit, simWeek) == -1) { return NULL; } - return &varVec[computeIndex(unit, simWeek)]; + int idx = computeIndex(unit, simWeek); + if (idx == -1) { return NULL; } + return &varVec[idx]; } void print (int unit, int simWeek) { // gpu lookupWeatherVars(unit, simWeek)->print(); diff --git a/src/WeatherData.cpp b/src/WeatherData.cpp index 7e0454e..9acf4df 100644 --- a/src/WeatherData.cpp +++ b/src/WeatherData.cpp @@ -20,7 +20,7 @@ void WeatherData::readDataFromFile (const std::string& fname) { std::istringstream lis(line); std::string temp[17]; int i = 0; - while (std::getline(lis, temp[i], ',')) { + while (i < 17 && std::getline(lis, temp[i], ',')) { i++; } weatherVars vars; @@ -138,6 +138,9 @@ bool WeatherData::lookupWeatherVars (int stateFP, int countyFP, date d, int& wee } void WeatherData::extractActiveData (DemographicData& demo, int startWeek, int numSimWeeks) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(startWeek + numSimWeeks <= numWeeks, + "Weather data does not cover the full simulation period: " + "reduce nsteps or provide a weather file with more data."); int numUnitsOnThisProc = 0; for (int i = 0; i < demo.Nunit; i++) { if (demo.Unit_on_proc[i]) { numUnitsOnThisProc++; } diff --git a/src/WeatherTransmissionModel.H b/src/WeatherTransmissionModel.H new file mode 100644 index 0000000..fd20f83 --- /dev/null +++ b/src/WeatherTransmissionModel.H @@ -0,0 +1,104 @@ +/*! @file WeatherTransmissionModel.H + \brief #WeatherTransmissionModel struct definition +*/ + +#ifndef WEATHER_TRANSMISSION_MODEL_H_ +#define WEATHER_TRANSMISSION_MODEL_H_ + +#include +#include +#include + +#include "WeatherData.H" + +#include +#include + +/*! \brief Weather-dependent transmission probability model. + * + * Computes the transmission probability as a function of near-surface air + * temperature and absolute humidity: + * + * \f[ + * p(T, AH) = p_{\max} \exp\!\bigl(-\beta_{AH}\,AH + * -\alpha_T\,\max(0,\,T - T_0)\bigr) + * \f] + * + * where \f$T\f$ is temperature in degrees Celsius and \f$AH\f$ is absolute + * humidity in g/m^3. The compute function (#operator()) is decorated with + * AMREX_GPU_HOST_DEVICE so it may be called from within amrex::ParallelFor + * kernels on both CPU and GPU backends. + * + * Parameters are parsed from the inputs file under the prefix supplied to + * #readInputs (typically \c "weather_transmission"): + * | Key | Default | Units | + * |---------------------------------|---------|-------------| + * | \c weather_transmission.p_max | 1.00 | — | + * | \c weather_transmission.beta_AH | 0.18 | m^3/g | + * | \c weather_transmission.T0 | 5.0 | °C | + * | \c weather_transmission.alpha_T | 0.015 | °C^{-1} | + */ +struct WeatherTransmissionModel { + + /*! Maximum transmission probability (dimensionless). */ + /*! This will be multiplied by p_trans later */ + amrex::Real p_max = amrex::Real(1.0); + + /*! Absolute-humidity sensitivity coefficient (m^3/g). + * Higher values make transmission drop faster with increasing humidity. */ + amrex::Real beta_AH = amrex::Real(0.18); + + /*! Temperature threshold above which the temperature penalty applies (°C). */ + amrex::Real T0 = amrex::Real(5.0); + + /*! Temperature sensitivity coefficient (°C^{-1}). + * Higher values make transmission drop faster above T0. */ + amrex::Real alpha_T = amrex::Real(0.015); + + /*! \brief Parse model parameters from the AMReX inputs file. + * + * Call this on the host during initialisation, passing the ParmParse + * prefix (e.g. \c "weather_transmission"). All four parameters are + * optional; unspecified keys retain their default values. + * + * \param pp_str ParmParse prefix string. + */ + void readInputs(const std::string& pp_str); + + /*! \brief Evaluate the transmission probability. + * + * This function is callable from device code and may therefore be used + * inside \c amrex::ParallelFor lambdas. + * + * \param T_celsius Near-surface air temperature (°C). + * \param AH_g_per_m3 Absolute humidity (g/m^3). + * \return Transmission probability \f$p(T,AH)\f$. + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real operator()(amrex::Real T_celsius, amrex::Real AH_g_per_m3) const noexcept { + amrex::Real T_excess = amrex::max(amrex::Real(0), T_celsius - T0); + return p_max * std::exp(-beta_AH * AH_g_per_m3 - alpha_T * T_excess); + } + + /*! \brief Evaluate the transmission probability directly from a #weatherVars record. + * + * Converts \c wvar.tas (surface temperature, K) to Celsius and + * \c wvar.huss (specific humidity, kg/kg) to absolute humidity (g/m^3) + * using a standard sea-level air density of 1.225 kg/m^3, then + * delegates to the two-argument overload. + * + * \param wvar Weather variables for the agent's location. + * \return Transmission probability \f$p(T,AH)\f$. + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + amrex::Real operator()(const weatherVars& wvar) const noexcept { + // tas is in Kelvin; convert to Celsius + amrex::Real T_celsius = static_cast(wvar.tas) - amrex::Real(273.15); + // huss is specific humidity [kg vapor / kg moist air]; multiply by standard + // air density (1.225 kg/m^3) and convert to g/m^3 to get absolute humidity + amrex::Real AH_g_per_m3 = static_cast(wvar.huss) * amrex::Real(1225.0); + return (*this)(T_celsius, AH_g_per_m3); + } +}; + +#endif diff --git a/src/WeatherTransmissionModel.cpp b/src/WeatherTransmissionModel.cpp new file mode 100644 index 0000000..550bc6a --- /dev/null +++ b/src/WeatherTransmissionModel.cpp @@ -0,0 +1,16 @@ +/*! @file WeatherTransmissionModel.cpp + \brief Function implementations for #WeatherTransmissionModel +*/ + +#include "WeatherTransmissionModel.H" + +#include + +/*! \brief Read weather-transmission model parameters from the inputs file. */ +void WeatherTransmissionModel::readInputs (const std::string& pp_str) { + amrex::ParmParse pp(pp_str); + pp.query("p_max", p_max); + pp.query("beta_AH", beta_AH); + pp.query("T0", T0); + pp.query("alpha_T", alpha_T); +}