CAMP 1.0.0
Chemistry Across Multiple Phases
rxn_SIMPOL_phase_transfer.c
Go to the documentation of this file.
1/* Copyright (C) 2021 Barcelona Supercomputing Center and University of
2 * Illinois at Urbana-Champaign
3 * SPDX-License-Identifier: MIT
4 *
5 * Phase Transfer reaction solver functions
6 *
7 */
8/** \file
9 * \brief Phase Transfer reaction solver functions
10 */
11#include <math.h>
12#include <stdio.h>
13#include <stdlib.h>
14#include "../aero_rep_solver.h"
15#include "../rxns.h"
16#include "../sub_model_solver.h"
17#include "../util.h"
18
19// TODO Lookup environmental indices during initialization
20#define TEMPERATURE_K_ env_data[0]
21#define PRESSURE_PA_ env_data[1]
22
23// Jacobian set indices
24#define JAC_GAS 0
25#define JAC_AERO 1
26
27// Aerosol mass concentration types
28#define PER_PARTICLE_MASS 0
29#define TOTAL_PARTICLE_MASS 1
30
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]
76
77/** \brief Flag Jacobian elements used by this reaction
78 *
79 * \param model_data Pointer to the model data
80 * \param rxn_int_data Pointer to the reaction integer data
81 * \param rxn_float_data Pointer to the reaction floating-point data
82 * \param jac Jacobian
83 */
85 int *rxn_int_data,
86 double *rxn_float_data,
87 Jacobian *jac) {
88 int *int_data = rxn_int_data;
89 double *float_data = rxn_float_data;
90
91 bool *aero_jac_elem =
92 (bool *)malloc(sizeof(bool) * model_data->n_per_cell_state_var);
93 if (aero_jac_elem == NULL) {
94 printf(
95 "\n\nERROR allocating space for 1D Jacobian structure array for "
96 "SIMPOL phase transfer reaction\n\n");
97 exit(1);
98 }
99
101 for (int i_aero_phase = 0; i_aero_phase < NUM_AERO_PHASE_; i_aero_phase++) {
104 jacobian_register_element(jac, AERO_SPEC_(i_aero_phase),
105 AERO_SPEC_(i_aero_phase));
106
107 if (AERO_ACT_ID_(i_aero_phase) > 0) {
109 jacobian_register_element(jac, AERO_SPEC_(i_aero_phase),
110 AERO_ACT_ID_(i_aero_phase));
111 }
112
113 for (int i_elem = 0; i_elem < model_data->n_per_cell_state_var; ++i_elem)
114 aero_jac_elem[i_elem] = false;
115
116 int n_jac_elem =
117 aero_rep_get_used_jac_elem(model_data, AERO_REP_ID_(i_aero_phase),
118 AERO_PHASE_ID_(i_aero_phase), aero_jac_elem);
119 if (n_jac_elem > NUM_AERO_PHASE_JAC_ELEM_(i_aero_phase)) {
120 printf(
121 "\n\nERROR Received more Jacobian elements than expected for SIMPOL "
122 "partitioning reaction. Got %d, expected <= %d",
123 n_jac_elem, NUM_AERO_PHASE_JAC_ELEM_(i_aero_phase));
124 exit(1);
125 }
126 int i_used_elem = 0;
127 for (int i_elem = 0; i_elem < model_data->n_per_cell_state_var; ++i_elem) {
128 if (aero_jac_elem[i_elem] == true) {
130 jacobian_register_element(jac, AERO_SPEC_(i_aero_phase), i_elem);
131 PHASE_JAC_ID_(i_aero_phase, JAC_GAS, i_used_elem) = i_elem;
132 PHASE_JAC_ID_(i_aero_phase, JAC_AERO, i_used_elem) = i_elem;
133 ++i_used_elem;
134 }
135 }
136 for (; i_used_elem < NUM_AERO_PHASE_JAC_ELEM_(i_aero_phase);
137 ++i_used_elem) {
138 PHASE_JAC_ID_(i_aero_phase, JAC_GAS, i_used_elem) = -1;
139 PHASE_JAC_ID_(i_aero_phase, JAC_AERO, i_used_elem) = -1;
140 }
141 if (i_used_elem != n_jac_elem) {
142 printf(
143 "\n\nERROR setting used Jacobian elements in SIMPOL phase "
144 "transfer reaction %d %d\n\n",
145 i_used_elem, n_jac_elem);
146 rxn_SIMPOL_phase_transfer_print(rxn_int_data, rxn_float_data);
147 exit(1);
148 }
149 }
150
151 free(aero_jac_elem);
152 return;
153}
154
155/** \brief Update the time derivative and Jacbobian array indices
156 *
157 * \param model_data Pointer to the model data for finding sub model ids
158 * \param deriv_ids Id of each state variable in the derivative array
159 * \param jac Jacobian
160 * \param rxn_int_data Pointer to the reaction integer data
161 * \param rxn_float_data Pointer to the reaction floating-point data
162 */
163void rxn_SIMPOL_phase_transfer_update_ids(ModelData *model_data, int *deriv_ids,
164 Jacobian jac, int *rxn_int_data,
165 double *rxn_float_data) {
166 int *int_data = rxn_int_data;
167 double *float_data = rxn_float_data;
168
169 // Update the time derivative ids
170 DERIV_ID_(0) = deriv_ids[GAS_SPEC_];
171 for (int i = 0; i < NUM_AERO_PHASE_; i++)
172 DERIV_ID_(i + 1) = deriv_ids[AERO_SPEC_(i)];
173
174 // Update the Jacobian ids
175 int i_jac = 0;
177 for (int i_aero_phase = 0; i_aero_phase < NUM_AERO_PHASE_; i_aero_phase++) {
178 JAC_ID_(i_jac++) =
179 jacobian_get_element_id(jac, AERO_SPEC_(i_aero_phase), GAS_SPEC_);
180 JAC_ID_(i_jac++) =
181 jacobian_get_element_id(jac, GAS_SPEC_, AERO_SPEC_(i_aero_phase));
182 JAC_ID_(i_jac++) = jacobian_get_element_id(jac, AERO_SPEC_(i_aero_phase),
183 AERO_SPEC_(i_aero_phase));
184 if (AERO_ACT_ID_(i_aero_phase) > 0) {
185 GAS_ACT_JAC_ID_(i_aero_phase) =
188 jac, AERO_SPEC_(i_aero_phase), AERO_ACT_ID_(i_aero_phase));
189 } else {
190 GAS_ACT_JAC_ID_(i_aero_phase) = -1;
191 AERO_ACT_JAC_ID_(i_aero_phase) = -1;
192 }
193 for (int i_elem = 0; i_elem < NUM_AERO_PHASE_JAC_ELEM_(i_aero_phase);
194 ++i_elem) {
195 if (PHASE_JAC_ID_(i_aero_phase, JAC_GAS, i_elem) > 0) {
196 PHASE_JAC_ID_(i_aero_phase, JAC_GAS, i_elem) = jacobian_get_element_id(
197 jac, GAS_SPEC_, PHASE_JAC_ID_(i_aero_phase, JAC_GAS, i_elem));
198 }
199 if (PHASE_JAC_ID_(i_aero_phase, JAC_AERO, i_elem) > 0) {
200 PHASE_JAC_ID_(i_aero_phase, JAC_AERO, i_elem) = jacobian_get_element_id(
201 jac, AERO_SPEC_(i_aero_phase),
202 PHASE_JAC_ID_(i_aero_phase, JAC_AERO, i_elem));
203 }
204 }
205 }
206
207 return;
208}
209
210/** \brief Update reaction data for new environmental conditions
211 *
212 * For Phase Transfer reaction this only involves recalculating the rate
213 * constant.
214 *
215 * \param model_data Pointer to the model data
216 * \param rxn_int_data Pointer to the reaction integer data
217 * \param rxn_float_data Pointer to the reaction floating-point data
218 * \param rxn_env_data Pointer to the environment-dependent parameters
219 */
221 int *rxn_int_data,
222 double *rxn_float_data,
223 double *rxn_env_data) {
224 int *int_data = rxn_int_data;
225 double *float_data = rxn_float_data;
226 double *env_data = model_data->grid_cell_env;
227
228 // Calculate the mass accomodation coefficient if the N* parameter
229 // was provided, otherwise set it to 0.1 (per Zaveri 2008)
230 ALPHA_ = 0.1;
231 if (DELTA_H_ != 0.0 || DELTA_S_ != 0.0) {
232 double del_G = DELTA_H_ - TEMPERATURE_K_ * DELTA_S_;
233 ALPHA_ = exp(-del_G / (UNIV_GAS_CONST_ * TEMPERATURE_K_));
234 ALPHA_ = ALPHA_ / (1.0 + ALPHA_);
235 }
236
237 // replaced by transition-regime rate equation
238#if 0
239 // Save c_rms * mass_acc for use in mass transfer rate calc
240 MFP_M_ = PRE_C_AVG_ * sqrt(TEMPERATURE_K_) * mass_acc; // [m/s]
241#endif
242
243 /// save the mean free path [m] for calculating condensation rates
245
246 // SIMPOL.1 vapor pressure [Pa]
247 double vp = B1_ / TEMPERATURE_K_ + B2_ + B3_ * TEMPERATURE_K_ +
248 B4_ * log(TEMPERATURE_K_);
249 vp = 101325.0 * pow(10, vp);
250
251 // Calculate the conversion from kg_x/m^3 -> ppm_x
253
254 // Calculate the partitioning coefficient K_eq (ppm_x/kg_x*kg_tot)
255 // such that for partitioning species X at equilibrium:
256 // [X]_gas = [X]_aero * activity_coeff_X * K_eq * MW_tot_aero / [tot]_aero
257 // where 'tot' indicates all species within an aerosol phase combined
258 // with []_gas in (ppm) and []_aero in (kg/m^3)
259 EQUIL_CONST_ = vp // Pa_x / (mol_x_aero/mol_tot_aero)
260 / PRESSURE_PA_ // 1/Pa_air
261 / MW_ // mol_x_aero/kg_x_aero
262 * 1.0e6; // ppm_x / (Pa_x/P_air)
263 // = ppm_x * mol_tot_aero / kg_x_aero
264
265 return;
266}
267
268/** \brief Calculate contributions to the time derivative \f$f(t,y)\f$ from
269 * this reaction.
270 *
271 * \param model_data Pointer to the model data, including the state array
272 * \param time_deriv TimeDerivative object
273 * \param rxn_int_data Pointer to the reaction integer data
274 * \param rxn_float_data Pointer to the reaction floating-point data
275 * \param rxn_env_data Pointer to the environment-dependent parameters
276 * \param time_step Current time step being computed (s)
277 */
278#ifdef CAMP_USE_SUNDIALS
280 ModelData *model_data, TimeDerivative time_deriv, int *rxn_int_data,
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;
284 double *state = model_data->grid_cell_state;
285 double *env_data = model_data->grid_cell_env;
286
287 // Calculate derivative contributions for each aerosol phase
288 for (int i_phase = 0; i_phase < NUM_AERO_PHASE_; i_phase++) {
289 // Get the particle effective radius (m)
290 realtype radius;
292 model_data, // model data
293 AERO_REP_ID_(i_phase), // aerosol representation index
294 AERO_PHASE_ID_(i_phase), // aerosol phase index
295 &radius, // particle effective radius (m)
296 NULL); // partial derivative
297
298 // Check the aerosol concentration type (per-particle or total per-phase
299 // mass)
300 int aero_conc_type = aero_rep_get_aero_conc_type(
301 model_data, // model data
302 AERO_REP_ID_(i_phase), // aerosol representation index
303 AERO_PHASE_ID_(i_phase)); // aerosol phase index
304
305 // Get the particle number concentration (#/m3)
306 realtype number_conc;
308 model_data, // model data
309 AERO_REP_ID_(i_phase), // aerosol representation index
310 AERO_PHASE_ID_(i_phase), // aerosol phase index
311 &number_conc, // particle number conc (#/m3)
312 NULL); // partial derivative
313
314 // Get the total mass of the aerosol phase (kg/m3)
315 realtype aero_phase_mass;
317 model_data, // model data
318 AERO_REP_ID_(i_phase), // aerosol representation index
319 AERO_PHASE_ID_(i_phase), // aerosol phase index
320 &aero_phase_mass, // total aerosol-phase mass (kg/m3)
321 NULL); // partial derivatives
322
323 // Get the total mass of the aerosol phase (kg/mol)
324 realtype aero_phase_avg_MW;
326 model_data, // model data
327 AERO_REP_ID_(i_phase), // aerosol representation index
328 AERO_PHASE_ID_(i_phase), // aerosol phase index
329 &aero_phase_avg_MW, // avg MW in the aerosol phase (kg/mol)
330 NULL); // partial derivatives
331
332 // This was replaced with the transition-regime condensation rate
333 // equations
334#if 0
335 long double cond_rate =
336 ((long double)1.0) / (radius * radius / (3.0 * DIFF_COEFF_) +
337 4.0 * radius / (3.0 * MFP_M_));
338#endif
339
340 // Calculate the rate constant for diffusion limited mass transfer to the
341 // aerosol phase (m3/#/s)
342 long double cond_rate =
344
345 // Calculate the evaporation rate constant (ppm_x*m^3/kg_x/s)
346 long double evap_rate =
347 cond_rate * (EQUIL_CONST_ * aero_phase_avg_MW / aero_phase_mass);
348
349 // Get the activity coefficient (if one exists)
350 long double act_coeff = 1.0;
351 if (AERO_ACT_ID_(i_phase) > -1) {
352 act_coeff = state[AERO_ACT_ID_(i_phase)];
353 }
354
355 // Calculate aerosol-phase evaporation rate (ppm/s)
356 evap_rate *= act_coeff;
357
358 // Calculate the evaporation and condensation rates
359 cond_rate *= state[GAS_SPEC_];
360 evap_rate *= state[AERO_SPEC_(i_phase)];
361
362 // per-particle mass concentration rates
363 if (aero_conc_type == PER_PARTICLE_MASS) {
364 // Change in the gas-phase is evaporation - condensation (ppm/s)
365 if (DERIV_ID_(0) >= 0) {
367 number_conc * evap_rate);
369 -number_conc * cond_rate);
370 }
371
372 // Change in the aerosol-phase species is condensation - evaporation
373 // (kg/m^3/s)
374 if (DERIV_ID_(1 + i_phase) >= 0) {
375 time_derivative_add_value(time_deriv, DERIV_ID_(1 + i_phase),
376 -evap_rate / KGM3_TO_PPM_);
377 time_derivative_add_value(time_deriv, DERIV_ID_(1 + i_phase),
378 cond_rate / KGM3_TO_PPM_);
379 }
380
381 // total-aerosol mass concentration rates
382 } else {
383 // Change in the gas-phase is evaporation - condensation (ppm/s)
384 if (DERIV_ID_(0) >= 0) {
386 number_conc * evap_rate);
388 -number_conc * cond_rate);
389 }
390
391 // Change in the aerosol-phase species is condensation - evaporation
392 // (kg/m^3/s)
393 if (DERIV_ID_(1 + i_phase) >= 0) {
394 time_derivative_add_value(time_deriv, DERIV_ID_(1 + i_phase),
395 -number_conc * evap_rate / KGM3_TO_PPM_);
396 time_derivative_add_value(time_deriv, DERIV_ID_(1 + i_phase),
397 number_conc * cond_rate / KGM3_TO_PPM_);
398 }
399 }
400 }
401 return;
402}
403#endif
404
405/** \brief Calculate contributions to the Jacobian from this reaction
406 *
407 * \param model_data Pointer to the model data
408 * \param jac Reaction Jacobian
409 * \param rxn_int_data Pointer to the reaction integer data
410 * \param rxn_float_data Pointer to the reaction floating-point data
411 * \param rxn_env_data Pointer to the environment-dependent parameters
412 * \param time_step Current time step being calculated (s)
413 */
414#ifdef CAMP_USE_SUNDIALS
416 Jacobian jac, int *rxn_int_data,
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;
422 double *state = model_data->grid_cell_state;
423 double *env_data = model_data->grid_cell_env;
424
425 // Calculate derivative contributions for each aerosol phase
426 for (int i_phase = 0; i_phase < NUM_AERO_PHASE_; i_phase++) {
427 // Get the particle effective radius (m)
428 realtype radius;
430 model_data, // model data
431 AERO_REP_ID_(i_phase), // aerosol representation index
432 AERO_PHASE_ID_(i_phase), // aerosol phase index
433 &radius, // particle effective radius (m)
434 &(EFF_RAD_JAC_ELEM_(i_phase, 0))); // partial derivative
435
436 // Check the aerosol concentration type (per-particle or total per-phase
437 // mass)
438 int aero_conc_type = aero_rep_get_aero_conc_type(
439 model_data, // model data
440 AERO_REP_ID_(i_phase), // aerosol representation index
441 AERO_PHASE_ID_(i_phase)); // aerosol phase index
442
443 // Get the particle number concentration (#/m3)
444 realtype number_conc;
446 model_data, // model data
447 AERO_REP_ID_(i_phase), // aerosol representation index
448 AERO_PHASE_ID_(i_phase), // aerosol phase index
449 &number_conc, // particle number conc (#/m3)
450 &(NUM_CONC_JAC_ELEM_(i_phase, 0))); // partial derivative
451
452 // Get the total mass of the aerosol phase (kg/m3)
453 realtype aero_phase_mass;
455 model_data, // model data
456 AERO_REP_ID_(i_phase), // aerosol representation index
457 AERO_PHASE_ID_(i_phase), // aerosol phase index
458 &aero_phase_mass, // total aerosol-phase mass (kg/m3)
459 &(MASS_JAC_ELEM_(i_phase, 0))); // partial derivatives
460
461 // Get the total average MW of the aerosol phase (kg/mol)
462 realtype aero_phase_avg_MW;
464 model_data, // model data
465 AERO_REP_ID_(i_phase), // aerosol representation index
466 AERO_PHASE_ID_(i_phase), // aerosol phase index
467 &aero_phase_avg_MW, // avg MW in the aerosol phase (kg/mol)
468 &(MW_JAC_ELEM_(i_phase, 0))); // partial derivatives
469
470 // This was replaced with the transition-regime condensation rate
471 // equations
472#if 0
473 long double cond_rate =
474 ((long double)1.0) / (radius * radius / (3.0 * DIFF_COEFF_) +
475 4.0 * radius / (3.0 * MFP_M_));
476#endif
477
478 // Calculate the rate constant for diffusion limited mass transfer to the
479 // aerosol phase (m3/#/s)
480 long double cond_rate =
482
483 // Calculate the evaporation rate constant (ppm_x*m^3/kg_x/s)
484 long double evap_rate =
485 cond_rate * (EQUIL_CONST_ * aero_phase_avg_MW / aero_phase_mass);
486
487 // Get the activity coefficient (if one exists)
488 long double act_coeff = 1.0;
489 if (AERO_ACT_ID_(i_phase) > -1) {
490 act_coeff = state[AERO_ACT_ID_(i_phase)];
491 }
492
493 // per-particle mass concentrations
494 if (aero_conc_type == PER_PARTICLE_MASS) {
495 // Change in the gas-phase is evaporation - condensation (ppm/s)
496 if (JAC_ID_(1 + i_phase * 3 + 1) >= 0) {
497 jacobian_add_value(jac, (unsigned int)JAC_ID_(1 + i_phase * 3 + 1),
499 number_conc * evap_rate * act_coeff);
500 }
501 if (JAC_ID_(0) >= 0) {
502 jacobian_add_value(jac, (unsigned int)JAC_ID_(0), JACOBIAN_LOSS,
503 number_conc * cond_rate);
504 }
505
506 // Change in the aerosol-phase species is condensation - evaporation
507 // (kg/m^3/s)
508 if (JAC_ID_(1 + i_phase * 3) >= 0) {
509 jacobian_add_value(jac, (unsigned int)JAC_ID_(1 + i_phase * 3),
510 JACOBIAN_PRODUCTION, cond_rate / KGM3_TO_PPM_);
511 }
512 if (JAC_ID_(1 + i_phase * 3 + 2) >= 0) {
513 jacobian_add_value(jac, (unsigned int)JAC_ID_(1 + i_phase * 3 + 2),
514 JACOBIAN_LOSS, evap_rate * act_coeff / KGM3_TO_PPM_);
515 }
516
517 // Activity coefficient contributions
518 if (GAS_ACT_JAC_ID_(i_phase) > 0) {
520 jac, (unsigned int)GAS_ACT_JAC_ID_(i_phase), JACOBIAN_PRODUCTION,
521 number_conc * evap_rate * state[AERO_SPEC_(i_phase)]);
522 }
523 if (AERO_ACT_JAC_ID_(i_phase) > 0) {
525 jac, (unsigned int)AERO_ACT_JAC_ID_(i_phase), JACOBIAN_LOSS,
526 evap_rate / KGM3_TO_PPM_ * state[AERO_SPEC_(i_phase)]);
527 }
528
529 // total-particle mass concentrations
530 } else {
531 // Change in the gas-phase is evaporation - condensation (ppm/s)
532 if (JAC_ID_(1 + i_phase * 3 + 1) >= 0) {
533 jacobian_add_value(jac, (unsigned int)JAC_ID_(1 + i_phase * 3 + 1),
535 number_conc * evap_rate * act_coeff);
536 }
537 if (JAC_ID_(0) >= 0) {
538 jacobian_add_value(jac, (unsigned int)JAC_ID_(0), JACOBIAN_LOSS,
539 number_conc * cond_rate);
540 }
541
542 // Change in the aerosol-phase species is condensation - evaporation
543 // (kg/m^3/s)
544 if (JAC_ID_(1 + i_phase * 3) >= 0) {
545 jacobian_add_value(jac, (unsigned int)JAC_ID_(1 + i_phase * 3),
547 number_conc * cond_rate / KGM3_TO_PPM_);
548 }
549 if (JAC_ID_(1 + i_phase * 3 + 2) >= 0) {
550 jacobian_add_value(jac, (unsigned int)JAC_ID_(1 + i_phase * 3 + 2),
552 number_conc * evap_rate * act_coeff / KGM3_TO_PPM_);
553 }
554
555 // Activity coefficient contributions
556 if (GAS_ACT_JAC_ID_(i_phase) > 0) {
558 jac, (unsigned int)GAS_ACT_JAC_ID_(i_phase), JACOBIAN_PRODUCTION,
559 number_conc * evap_rate * state[AERO_SPEC_(i_phase)]);
560 }
561 if (AERO_ACT_JAC_ID_(i_phase) > 0) {
562 jacobian_add_value(jac, (unsigned int)AERO_ACT_JAC_ID_(i_phase),
564 number_conc * evap_rate / KGM3_TO_PPM_ *
565 state[AERO_SPEC_(i_phase)]);
566 }
567 }
568
569 // Get the overall rates
570 evap_rate *= act_coeff;
571 cond_rate *= state[GAS_SPEC_];
572 evap_rate *= state[AERO_SPEC_(i_phase)];
573
574 // Calculate partial derivatives
575
576 // this was replaced with the transition regime rate equations
577#if 0
578 realtype d_cond_d_radius =
579 -(2.0 * radius / (3.0 * DIFF_COEFF_) + 4.0 / (3.0 * MFP_M_)) *
580 cond_rate * cond_rate / state[GAS_SPEC_];
581#endif
583 DIFF_COEFF_, MFP_M_, radius, ALPHA_) *
584 state[GAS_SPEC_];
585 realtype d_evap_d_radius = d_cond_d_radius / state[GAS_SPEC_] *
586 EQUIL_CONST_ * aero_phase_avg_MW /
587 aero_phase_mass * state[AERO_SPEC_(i_phase)];
588 realtype d_evap_d_mass = -evap_rate / aero_phase_mass;
589 realtype d_evap_d_MW = evap_rate / aero_phase_avg_MW;
590
591 // per-particle mass concentrations
592 if (aero_conc_type == PER_PARTICLE_MASS) {
593 // Loop through Jac elements and update
594 for (int i_elem = 0; i_elem < NUM_AERO_PHASE_JAC_ELEM_(i_phase);
595 ++i_elem) {
596 // Gas-phase species dependencies
597 if (PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem) > 0) {
598 // species involved in effective radius calculations
600 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
602 number_conc * d_evap_d_radius *
603 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
605 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
607 number_conc * d_cond_d_radius *
608 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
609
610 // species involved in number concentration
612 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
614 evap_rate * NUM_CONC_JAC_ELEM_(i_phase, i_elem));
616 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
617 JACOBIAN_LOSS, cond_rate * NUM_CONC_JAC_ELEM_(i_phase, i_elem));
618
619 // species involved in mass calculations
621 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
623 number_conc * d_evap_d_mass * MASS_JAC_ELEM_(i_phase, i_elem));
624
625 // species involved in average MW calculations
627 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
629 number_conc * d_evap_d_MW * MW_JAC_ELEM_(i_phase, i_elem));
630 }
631
632 // Aerosol-phase species dependencies
633 if (PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem) > 0) {
634 // species involved in effective radius calculations
636 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
638 d_evap_d_radius / KGM3_TO_PPM_ *
639 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
641 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
643 d_cond_d_radius / KGM3_TO_PPM_ *
644 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
645
646 // species involved in mass calculations
648 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
650 d_evap_d_mass / KGM3_TO_PPM_ * MASS_JAC_ELEM_(i_phase, i_elem));
651
652 // species involved in average MW calculations
654 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
656 d_evap_d_MW / KGM3_TO_PPM_ * MW_JAC_ELEM_(i_phase, i_elem));
657 }
658 }
659
660 // total-particle mass concentrations
661 } else {
662 // Loop through Jac elements and update
663 for (int i_elem = 0; i_elem < NUM_AERO_PHASE_JAC_ELEM_(i_phase);
664 ++i_elem) {
665 // Gas-phase species dependencies
666 if (PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem) > 0) {
667 // species involved in effective radius calculations
669 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
671 number_conc * d_evap_d_radius *
672 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
674 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
676 number_conc * d_cond_d_radius *
677 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
678
679 // species involved in number concentration
681 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
683 evap_rate * NUM_CONC_JAC_ELEM_(i_phase, i_elem));
685 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
686 JACOBIAN_LOSS, cond_rate * NUM_CONC_JAC_ELEM_(i_phase, i_elem));
687
688 // species involved in mass calculations
690 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
692 number_conc * d_evap_d_mass * MASS_JAC_ELEM_(i_phase, i_elem));
693
694 // species involved in average MW calculations
696 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_GAS, i_elem),
698 number_conc * d_evap_d_MW * MW_JAC_ELEM_(i_phase, i_elem));
699 }
700
701 // Aerosol-phase species dependencies
702 if (PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem) > 0) {
703 // species involved in effective radius calculations
705 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
707 number_conc * d_evap_d_radius / KGM3_TO_PPM_ *
708 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
710 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
712 number_conc * d_cond_d_radius / KGM3_TO_PPM_ *
713 EFF_RAD_JAC_ELEM_(i_phase, i_elem));
714
715 // species involved in number concentration
717 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
719 evap_rate / KGM3_TO_PPM_ * NUM_CONC_JAC_ELEM_(i_phase, i_elem));
721 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
723 cond_rate / KGM3_TO_PPM_ * NUM_CONC_JAC_ELEM_(i_phase, i_elem));
724
725 // species involved in mass calculations
727 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
729 number_conc * d_evap_d_mass / KGM3_TO_PPM_ *
730 MASS_JAC_ELEM_(i_phase, i_elem));
731
732 // species involved in average MW calculations
734 jac, (unsigned int)PHASE_JAC_ID_(i_phase, JAC_AERO, i_elem),
736 number_conc * d_evap_d_MW / KGM3_TO_PPM_ *
737 MW_JAC_ELEM_(i_phase, i_elem));
738 }
739 }
740 }
741 }
742 return;
743}
744#endif
745
746/** \brief Print the Phase Transfer reaction parameters
747 *
748 * \param rxn_int_data Pointer to the reaction integer data
749 * \param rxn_float_data Pointer to the reaction floating-point data
750 */
751void rxn_SIMPOL_phase_transfer_print(int *rxn_int_data,
752 double *rxn_float_data) {
753 int *int_data = rxn_int_data;
754 double *float_data = rxn_float_data;
755
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_);
761 printf("\nNumber of aerosol phases: %d", NUM_AERO_PHASE_);
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 ***");
766 for (int i = 0; i < NUM_AERO_PHASE_; ++i) {
767 printf("\n Aerosol species id: %d", AERO_SPEC_(i));
768 printf("\n Activity coefficient id: %d", AERO_ACT_ID_(i));
769 printf("\n Aerosol phase id: %d", AERO_PHASE_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));
772 printf("\n dGas/dAct coeff Jac id: %d", GAS_ACT_JAC_ID_(i));
773 printf("\n dAero/dAct coeff Jac id: %d", AERO_ACT_JAC_ID_(i));
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:");
780 for (int j = 0; j < NUM_AERO_PHASE_JAC_ELEM_(i); ++j)
781 printf(" %d", PHASE_JAC_ID_(i, JAC_GAS, j));
782 printf("\n dAero/dx ids:");
783 for (int j = 0; j < NUM_AERO_PHASE_JAC_ELEM_(i); ++j)
784 printf(" %d", PHASE_JAC_ID_(i, JAC_AERO, j));
785 printf("\n Effective radius Jac elem:");
786 for (int j = 0; j < NUM_AERO_PHASE_JAC_ELEM_(i); ++j)
787 printf(" %le", EFF_RAD_JAC_ELEM_(i, j));
788 printf("\n Number concentration Jac elem:");
789 for (int j = 0; j < NUM_AERO_PHASE_JAC_ELEM_(i); ++j)
790 printf(" %le", NUM_CONC_JAC_ELEM_(i, j));
791 printf("\n Aerosol mass Jac elem:");
792 for (int j = 0; j < NUM_AERO_PHASE_JAC_ELEM_(i); ++j)
793 printf(" %le", MASS_JAC_ELEM_(i, j));
794 printf("\n Average MW Jac elem:");
795 for (int j = 0; j < NUM_AERO_PHASE_JAC_ELEM_(i); ++j)
796 printf(" %le", MW_JAC_ELEM_(i, j));
797 }
798
799 return;
800}
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.
Definition Jacobian.c:200
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.
Definition Jacobian.c:234
void jacobian_register_element(Jacobian *jac, unsigned int dep_id, unsigned int ind_id)
Adds an element to the sparse matrix.
Definition Jacobian.c:105
#define JACOBIAN_LOSS
Definition Jacobian.h:19
#define JACOBIAN_PRODUCTION
Definition Jacobian.h:18
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 AERO_ACT_ID_(x)
#define JAC_AERO
#define MW_JAC_ELEM_(x, e)
#define DELTA_H_
#define MW_
#define JAC_GAS
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 PRESSURE_PA_
#define AERO_PHASE_ID_(x)
#define B2_
#define B1_
#define EQUIL_CONST_
void rxn_SIMPOL_phase_transfer_print(int *rxn_int_data, double *rxn_float_data)
Print the Phase Transfer reaction parameters.
#define PRE_C_AVG_
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 MFP_M_
#define KGM3_TO_PPM_
#define CONV_
#define ALPHA_
#define GAS_ACT_JAC_ID_(x)
#define GAS_SPEC_
#define TEMPERATURE_K_
#define PER_PARTICLE_MASS
#define B4_
#define JAC_ID_(x)
#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 B3_
#define NUM_AERO_PHASE_JAC_ELEM_(x)
#define DIFF_COEFF_
#define NUM_AERO_PHASE_
#define AERO_REP_ID_(x)
#define EFF_RAD_JAC_ELEM_(x, e)
#define AERO_SPEC_(x)
#define DELTA_S_
#define NUM_CONC_JAC_ELEM_(x, e)
#define DERIV_ID_(x)
int n_per_cell_state_var
Definition camp_common.h:63
double * grid_cell_env
double * grid_cell_state
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)
Definition util.h:154
static double gas_aerosol_transition_rxn_rate_constant(double diffusion_coeff__m2_s, double mean_free_path__m, double radius__m, double alpha)
Definition util.h:129
static double mean_free_path__m(double diffusion_coeff__m2_s, double temperature__K, double mw__kg_mol)
Definition util.h:50
#define UNIV_GAS_CONST_
Definition util.h:15