14#include "../aero_rep_solver.h"
16#include "../sub_model_solver.h"
20#define TEMPERATURE_K_ env_data[0]
21#define PRESSURE_PA_ env_data[1]
28#define PER_PARTICLE_MASS 0
29#define TOTAL_PARTICLE_MASS 1
31#define DELTA_H_ float_data[0]
32#define DELTA_S_ float_data[1]
33#define DIFF_COEFF_ float_data[2]
34#define PRE_C_AVG_ float_data[3]
35#define B1_ float_data[4]
36#define B2_ float_data[5]
37#define B3_ float_data[6]
38#define B4_ float_data[7]
39#define CONV_ float_data[8]
40#define MW_ float_data[9]
41#define NUM_AERO_PHASE_ int_data[0]
42#define GAS_SPEC_ (int_data[1] - 1)
43#define MFP_M_ rxn_env_data[0]
44#define ALPHA_ rxn_env_data[1]
45#define EQUIL_CONST_ rxn_env_data[2]
46#define KGM3_TO_PPM_ rxn_env_data[3]
47#define NUM_INT_PROP_ 2
48#define NUM_FLOAT_PROP_ 10
49#define NUM_ENV_PARAM_ 4
50#define AERO_SPEC_(x) (int_data[NUM_INT_PROP_ + x] - 1)
51#define AERO_ACT_ID_(x) (int_data[NUM_INT_PROP_ + NUM_AERO_PHASE_ + x] - 1)
52#define AERO_PHASE_ID_(x) \
53 (int_data[NUM_INT_PROP_ + 2 * (NUM_AERO_PHASE_) + x] - 1)
54#define AERO_REP_ID_(x) \
55 (int_data[NUM_INT_PROP_ + 3 * (NUM_AERO_PHASE_) + x] - 1)
56#define DERIV_ID_(x) (int_data[NUM_INT_PROP_ + 4 * (NUM_AERO_PHASE_) + x])
57#define GAS_ACT_JAC_ID_(x) \
58 int_data[NUM_INT_PROP_ + 1 + 5 * (NUM_AERO_PHASE_) + x]
59#define AERO_ACT_JAC_ID_(x) \
60 int_data[NUM_INT_PROP_ + 1 + 6 * (NUM_AERO_PHASE_) + x]
61#define JAC_ID_(x) (int_data[NUM_INT_PROP_ + 1 + 7 * (NUM_AERO_PHASE_) + x])
62#define PHASE_INT_LOC_(x) \
63 (int_data[NUM_INT_PROP_ + 2 + 10 * (NUM_AERO_PHASE_) + x] - 1)
64#define PHASE_FLOAT_LOC_(x) \
65 (int_data[NUM_INT_PROP_ + 2 + 11 * (NUM_AERO_PHASE_) + x] - 1)
66#define NUM_AERO_PHASE_JAC_ELEM_(x) (int_data[PHASE_INT_LOC_(x)])
67#define PHASE_JAC_ID_(x, s, e) \
68 int_data[PHASE_INT_LOC_(x) + 1 + (s) * NUM_AERO_PHASE_JAC_ELEM_(x) + e]
69#define EFF_RAD_JAC_ELEM_(x, e) float_data[PHASE_FLOAT_LOC_(x) + e]
70#define NUM_CONC_JAC_ELEM_(x, e) \
71 float_data[PHASE_FLOAT_LOC_(x) + NUM_AERO_PHASE_JAC_ELEM_(x) + e]
72#define MASS_JAC_ELEM_(x, e) \
73 float_data[PHASE_FLOAT_LOC_(x) + 2 * NUM_AERO_PHASE_JAC_ELEM_(x) + e]
74#define MW_JAC_ELEM_(x, e) \
75 float_data[PHASE_FLOAT_LOC_(x) + 3 * NUM_AERO_PHASE_JAC_ELEM_(x) + e]
86 double *rxn_float_data,
88 int *int_data = rxn_int_data;
89 double *float_data = rxn_float_data;
93 if (aero_jac_elem == NULL) {
95 "\n\nERROR allocating space for 1D Jacobian structure array for "
96 "SIMPOL phase transfer reaction\n\n");
101 for (
int i_aero_phase = 0; i_aero_phase <
NUM_AERO_PHASE_; i_aero_phase++) {
114 aero_jac_elem[i_elem] =
false;
121 "\n\nERROR Received more Jacobian elements than expected for SIMPOL "
122 "partitioning reaction. Got %d, expected <= %d",
128 if (aero_jac_elem[i_elem] ==
true) {
141 if (i_used_elem != n_jac_elem) {
143 "\n\nERROR setting used Jacobian elements in SIMPOL phase "
144 "transfer reaction %d %d\n\n",
145 i_used_elem, n_jac_elem);
165 double *rxn_float_data) {
166 int *int_data = rxn_int_data;
167 double *float_data = rxn_float_data;
177 for (
int i_aero_phase = 0; i_aero_phase <
NUM_AERO_PHASE_; i_aero_phase++) {
222 double *rxn_float_data,
223 double *rxn_env_data) {
224 int *int_data = rxn_int_data;
225 double *float_data = rxn_float_data;
249 vp = 101325.0 * pow(10, vp);
278#ifdef CAMP_USE_SUNDIALS
281 double *rxn_float_data,
double *rxn_env_data, realtype time_step) {
282 int *int_data = rxn_int_data;
283 double *float_data = rxn_float_data;
306 realtype number_conc;
315 realtype aero_phase_mass;
324 realtype aero_phase_avg_MW;
335 long double cond_rate =
336 ((
long double)1.0) / (radius * radius / (3.0 *
DIFF_COEFF_) +
337 4.0 * radius / (3.0 *
MFP_M_));
342 long double cond_rate =
346 long double evap_rate =
347 cond_rate * (
EQUIL_CONST_ * aero_phase_avg_MW / aero_phase_mass);
350 long double act_coeff = 1.0;
356 evap_rate *= act_coeff;
367 number_conc * evap_rate);
369 -number_conc * cond_rate);
386 number_conc * evap_rate);
388 -number_conc * cond_rate);
414#ifdef CAMP_USE_SUNDIALS
417 double *rxn_float_data,
418 double *rxn_env_data,
419 realtype time_step) {
420 int *int_data = rxn_int_data;
421 double *float_data = rxn_float_data;
444 realtype number_conc;
453 realtype aero_phase_mass;
462 realtype aero_phase_avg_MW;
473 long double cond_rate =
474 ((
long double)1.0) / (radius * radius / (3.0 *
DIFF_COEFF_) +
475 4.0 * radius / (3.0 *
MFP_M_));
480 long double cond_rate =
484 long double evap_rate =
485 cond_rate * (
EQUIL_CONST_ * aero_phase_avg_MW / aero_phase_mass);
488 long double act_coeff = 1.0;
496 if (
JAC_ID_(1 + i_phase * 3 + 1) >= 0) {
499 number_conc * evap_rate * act_coeff);
503 number_conc * cond_rate);
508 if (
JAC_ID_(1 + i_phase * 3) >= 0) {
512 if (
JAC_ID_(1 + i_phase * 3 + 2) >= 0) {
521 number_conc * evap_rate * state[
AERO_SPEC_(i_phase)]);
532 if (
JAC_ID_(1 + i_phase * 3 + 1) >= 0) {
535 number_conc * evap_rate * act_coeff);
539 number_conc * cond_rate);
544 if (
JAC_ID_(1 + i_phase * 3) >= 0) {
549 if (
JAC_ID_(1 + i_phase * 3 + 2) >= 0) {
559 number_conc * evap_rate * state[
AERO_SPEC_(i_phase)]);
570 evap_rate *= act_coeff;
578 realtype d_cond_d_radius =
580 cond_rate * cond_rate / state[
GAS_SPEC_];
585 realtype d_evap_d_radius = d_cond_d_radius / state[
GAS_SPEC_] *
588 realtype d_evap_d_mass = -evap_rate / aero_phase_mass;
589 realtype d_evap_d_MW = evap_rate / aero_phase_avg_MW;
602 number_conc * d_evap_d_radius *
607 number_conc * d_cond_d_radius *
629 number_conc * d_evap_d_MW *
MW_JAC_ELEM_(i_phase, i_elem));
671 number_conc * d_evap_d_radius *
676 number_conc * d_cond_d_radius *
698 number_conc * d_evap_d_MW *
MW_JAC_ELEM_(i_phase, i_elem));
752 double *rxn_float_data) {
753 int *int_data = rxn_int_data;
754 double *float_data = rxn_float_data;
756 printf(
"\n\nSIMPOL.1 Phase Transfer reaction\n");
757 printf(
"\ndelta H: %le delta S: %le diffusion coeff: %le Pre C_avg: %le",
759 printf(
"\nB1: %le B2: %le B3:%le B4: %le",
B1_,
B2_,
B3_,
B4_);
760 printf(
"\nconv: %le MW: %le",
CONV_,
MW_);
762 printf(
"\nGas-phase species id: %d",
GAS_SPEC_);
763 printf(
"\nGas-phase derivative id: %d",
DERIV_ID_(0));
764 printf(
"\ndGas/dGas Jac id: %d",
JAC_ID_(0));
765 printf(
"\n*** Aerosol phase data ***");
767 printf(
"\n Aerosol species id: %d",
AERO_SPEC_(i));
768 printf(
"\n Activity coefficient id: %d",
AERO_ACT_ID_(i));
770 printf(
"\n Aerosol representation id: %d",
AERO_REP_ID_(i));
771 printf(
"\n Aerosol species derivative id: %d",
DERIV_ID_(i + 1));
774 printf(
"\n dAero/dGas Jac id: %d",
JAC_ID_(1 + i * 3));
775 printf(
"\n dGas/dAero Jac id: %d",
JAC_ID_(2 + i * 3));
776 printf(
"\n dAero/dAero Jac id: %d",
JAC_ID_(3 + i * 3));
777 printf(
"\n Number of aerosol-phase species Jac elements: %d",
779 printf(
"\n dGas/dx ids:");
782 printf(
"\n dAero/dx ids:");
785 printf(
"\n Effective radius Jac elem:");
788 printf(
"\n Number concentration Jac elem:");
791 printf(
"\n Aerosol mass Jac elem:");
794 printf(
"\n Average MW Jac elem:");
unsigned int jacobian_get_element_id(Jacobian jac, unsigned int dep_id, unsigned int ind_id)
Get an element id in the Jacobian data arrays.
void jacobian_add_value(Jacobian jac, unsigned int elem_id, unsigned int prod_or_loss, long double jac_contribution)
Add a contribution to the Jacobian.
void jacobian_register_element(Jacobian *jac, unsigned int dep_id, unsigned int ind_id)
Adds an element to the sparse matrix.
#define JACOBIAN_PRODUCTION
int aero_rep_get_aero_conc_type(ModelData *model_data, int aero_rep_idx, int aero_phase_idx)
Check whether aerosol concentrations are per-particle or total for each phase.
int aero_rep_get_used_jac_elem(ModelData *model_data, int aero_rep_idx, int aero_phase_idx, bool *jac_struct)
Flag Jacobian elements used to calculated mass, volume, etc.
void aero_rep_get_aero_phase_mass__kg_m3(ModelData *model_data, int aero_rep_idx, int aero_phase_idx, double *aero_phase_mass, double *partial_deriv)
Get the total mass of an aerosol phase in this representation ( )
void aero_rep_get_aero_phase_avg_MW__kg_mol(ModelData *model_data, int aero_rep_idx, int aero_phase_idx, double *aero_phase_avg_MW, double *partial_deriv)
Get the average molecular weight of an aerosol phase in this representation ( )
void aero_rep_get_number_conc__n_m3(ModelData *model_data, int aero_rep_idx, int aero_phase_idx, double *number_conc, double *partial_deriv)
Get the particle number concentration ( )
void aero_rep_get_effective_radius__m(ModelData *model_data, int aero_rep_idx, int aero_phase_idx, double *radius, double *partial_deriv)
Get the effective particle radius, (m)
#define MW_JAC_ELEM_(x, e)
void rxn_SIMPOL_phase_transfer_get_used_jac_elem(ModelData *model_data, int *rxn_int_data, double *rxn_float_data, Jacobian *jac)
Flag Jacobian elements used by this reaction.
void rxn_SIMPOL_phase_transfer_update_env_state(ModelData *model_data, int *rxn_int_data, double *rxn_float_data, double *rxn_env_data)
Update reaction data for new environmental conditions.
void rxn_SIMPOL_phase_transfer_calc_jac_contrib(ModelData *model_data, Jacobian jac, int *rxn_int_data, double *rxn_float_data, double *rxn_env_data, realtype time_step)
Calculate contributions to the Jacobian from this reaction.
#define AERO_PHASE_ID_(x)
void rxn_SIMPOL_phase_transfer_print(int *rxn_int_data, double *rxn_float_data)
Print the Phase Transfer reaction parameters.
void rxn_SIMPOL_phase_transfer_calc_deriv_contrib(ModelData *model_data, TimeDerivative time_deriv, int *rxn_int_data, double *rxn_float_data, double *rxn_env_data, realtype time_step)
Calculate contributions to the time derivative from this reaction.
#define GAS_ACT_JAC_ID_(x)
#define PER_PARTICLE_MASS
#define PHASE_JAC_ID_(x, s, e)
void rxn_SIMPOL_phase_transfer_update_ids(ModelData *model_data, int *deriv_ids, Jacobian jac, int *rxn_int_data, double *rxn_float_data)
Update the time derivative and Jacbobian array indices.
#define AERO_ACT_JAC_ID_(x)
#define MASS_JAC_ELEM_(x, e)
#define NUM_AERO_PHASE_JAC_ELEM_(x)
#define EFF_RAD_JAC_ELEM_(x, e)
#define NUM_CONC_JAC_ELEM_(x, e)
void time_derivative_add_value(TimeDerivative time_deriv, unsigned int spec_id, long double rate_contribution)
Add a contribution to the time derivative.
static double d_gas_aerosol_transition_rxn_rate_constant_d_radius(double diffusion_coeff__m2_s, double mean_free_path__m, double radius__m, double alpha)
static double gas_aerosol_transition_rxn_rate_constant(double diffusion_coeff__m2_s, double mean_free_path__m, double radius__m, double alpha)
static double mean_free_path__m(double diffusion_coeff__m2_s, double temperature__K, double mw__kg_mol)