1 /*============================================================================
2  * Physical properties management for groundwater flow module.
3  *============================================================================*/
4 
5 /*
6   This file is part of Code_Saturne, a general-purpose CFD tool.
7 
8   Copyright (C) 1998-2021 EDF S.A.
9 
10   This program is free software; you can redistribute it and/or modify it under
11   the terms of the GNU General Public License as published by the Free Software
12   Foundation; either version 2 of the License, or (at your option) any later
13   version.
14 
15   This program is distributed in the hope that it will be useful, but WITHOUT
16   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
18   details.
19 
20   You should have received a copy of the GNU General Public License along with
21   this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22   Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24 
25 /*----------------------------------------------------------------------------*/
26 
27 #include "cs_defs.h"
28 
29 
30 /*----------------------------------------------------------------------------*/
31 
32 /*----------------------------------------------------------------------------
33  * Standard C library headers
34  *----------------------------------------------------------------------------*/
35 
36 #include <assert.h>
37 #include <stdio.h>
38 #include <stdarg.h>
39 #include <stdlib.h>
40 #include <string.h>
41 #include <math.h>
42 
43 /*----------------------------------------------------------------------------
44  * Local headers
45  *----------------------------------------------------------------------------*/
46 
47 #include "bft_mem.h"
48 #include "bft_error.h"
49 #include "bft_printf.h"
50 
51 #include "cs_gwf_parameters.h"
52 #include "cs_math.h"
53 #include "cs_field.h"
54 #include "cs_field_pointer.h"
55 #include "cs_mesh.h"
56 #include "cs_mesh_quantities.h"
57 
58 /*----------------------------------------------------------------------------
59  * Header for the current file
60  *----------------------------------------------------------------------------*/
61 
62 #include "cs_gwf_physical_properties.h"
63 
64 /*----------------------------------------------------------------------------*/
65 
66 BEGIN_C_DECLS
67 
68 /*=============================================================================
69  * Public function definitions
70  *============================================================================*/
71 
72 /*----------------------------------------------------------------------------*/
73 /*!
74  * \brief  Update delay of all transported species (user scalars)
75  *
76  * Species transport is delayed by retention in solid phase.
77  * This delay is computed as follows:
78  * \f[
79  *    R = 1 + \rho \dfrac{K_d}{\theta} ;
80  * \f]
81  * where \f$ R \f$ is the delay factor, \f$ \rho \f$ the soil (bulk) density,
82  * \f$ K_d \f$ the contaminant distribution coefficient and \f$ \theta \f$
83  * the moisture content.
84  *
85  * Please refer to the
86  * <a href="../../theory.pdf#gwf_sp_transp"><b>dedicated section</b></a>
87  * of the theory guide for more informations.
88  *
89  */
90  /*---------------------------------------------------------------------------*/
91 
cs_gwf_delay_update(void)92 void cs_gwf_delay_update(void)
93 {
94   const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
95 
96   cs_field_t *sca, *delay, *sat, *rosoil, *kd;
97   cs_gwf_soilwater_partition_t sorption_scal;
98 
99   int key_part = cs_field_key_id("gwf_soilwater_partition");
100 
101   /* Get soil property fields */
102   sat = cs_field_by_name("saturation");
103   rosoil = cs_field_by_name("soil_density");
104 
105   /* Loop on user scalars */
106   for (int f_id = 0; f_id < cs_field_n_fields(); f_id++) {
107     sca = cs_field_by_id(f_id);
108     if ((sca->type & CS_FIELD_VARIABLE) && (sca->type & CS_FIELD_USER)) {
109       /* Sorption properties are needed to update delay */
110       cs_field_get_key_struct(sca, key_part, &sorption_scal);
111 
112       /* Update of delay */
113       kd = cs_field_by_id(sorption_scal.ikd);
114       delay = cs_field_by_id(sorption_scal.idel);
115 
116       for (int c_id = 0; c_id < n_cells; c_id++) {
117         delay->val[c_id] = 1. +  rosoil->val[c_id]
118                                * kd->val[c_id] / sat->val[c_id];
119       }
120     }
121   }
122 }
123 
124 /*----------------------------------------------------------------------------*/
125 /*!
126  * \brief  Add first-order decay to implicit part of source term array
127  *
128  * Corresponding EDO for decay phenomenon is:
129  * \f[
130  *    \dfrac{dc}{dt} = - decay_rate  c
131  * \f]
132  *
133  * \param[in]     f_id      field index of scalar on which decay is set
134  * \param[in,out] ts_imp    source term implicit part for scalar of index f_id
135  */
136  /*---------------------------------------------------------------------------*/
137 
cs_gwf_decay_rate(const int f_id,cs_real_t * ts_imp)138 void cs_gwf_decay_rate(const int        f_id,
139                        cs_real_t       *ts_imp)
140 {
141   const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
142   const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
143   const cs_real_t *vol = mesh_quantities->cell_vol;
144 
145   cs_field_t *sca = cs_field_by_id(f_id);
146   const int dr_id = cs_field_key_id("fo_decay_rate");
147   cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
148 
149   /* First-order decay_rate in implict term */
150   for (int c_id = 0; c_id < n_cells; c_id++)
151     ts_imp[c_id] -= decay_rate * vol[c_id];
152 }
153 
154 /*----------------------------------------------------------------------------*/
155 /*!
156  * \brief Update sorbed concentration for scalars with kinetic sorption.
157  *
158  * It is estimated by the following analytical expression :
159  * \f[
160  * S^{n+1} = S^n \exp(- (k^{-} + decay_rate) * \Delta t)
161  *           - C^n \dfrac{k^{+}}{k^{-}}
162  *             \left(\exp(- (k^{-} + decay_rate) * \Delta t) - 1 \right)
163  * \f]
164  *
165  * \param[in]   f_id   field index of scalar which properties are updated
166  */
167  /*---------------------------------------------------------------------------*/
168 
cs_gwf_sorbed_concentration_update(const int f_id)169 void cs_gwf_sorbed_concentration_update(const int f_id)
170 {
171   const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
172   cs_real_t *dt = CS_F_(dt)->val;
173 
174   /* Get scalar, sorbed concentration fields */
175   cs_field_t *sca = cs_field_by_id(f_id);
176   int key_sorb = cs_field_key_id("gwf_sorbed_concentration_id");
177   cs_field_t *sorb = cs_field_by_id(cs_field_get_key_int(sca, key_sorb));
178 
179   /* Get first-order decay rate */
180   const int dr_id = cs_field_key_id("fo_decay_rate");
181   cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
182 
183   /* Get soil-water partition structure */
184   cs_gwf_soilwater_partition_t sorption_scal;
185   const int key_part = cs_field_key_id("gwf_soilwater_partition");
186   cs_field_get_key_struct(sca, key_part, &sorption_scal);
187 
188   /* Get k+ and k- fields */
189   cs_field_t *kp = cs_field_by_id(sorption_scal.ikp);
190   cs_field_t *km = cs_field_by_id(sorption_scal.ikm);
191 
192   /* Update sorbed concentration */
193 
194   /* First choice : analytical resolution */
195   if (sorption_scal.anai){
196     for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++)
197       /* Case of reversible sorption or decay rate term */
198       if (km->val[c_id] + decay_rate > cs_math_epzero) {
199         cs_real_t expkdt = exp(-(km->val[c_id] + decay_rate) * dt[c_id]);
200         cs_real_t kpskm = kp->val[c_id] / (km->val[c_id] + decay_rate);
201         sorb->val[c_id] =  expkdt * sorb->val[c_id]
202                          - (expkdt-1.) * kpskm * sca->val[c_id];
203       }
204       else { /* Irreversible sorption and no decay rate */
205         sorb->val[c_id] =  sorb->val[c_id]
206                          + dt[c_id]*kp->val[c_id]*sca->val[c_id];
207       }
208   }
209   else { /* Second choice : explicit */
210     for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
211       sorb->val[c_id] +=  dt[c_id] * (kp->val[c_id] * sca->val[c_id]
212                         - (km->val[c_id] + decay_rate) * sorb->val[c_id]);
213     }
214   }
215 }
216 
217 /*----------------------------------------------------------------------------*/
218 /*!
219  * \brief Update liquid concentration according to precipitation phenomenon.
220  *
221  * If liquid concentration exceeds solubility index, then it is clipped and
222  * transferred in precipitated concentration.
223  *
224  * Note that the decay rate is applied before the clipping.
225  *
226  * \param[in]   f_id    field index of scalar which properties are updated
227  */
228  /*---------------------------------------------------------------------------*/
229 
cs_gwf_precipitation(const int f_id)230 void cs_gwf_precipitation(const int f_id)
231 {
232   const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
233 
234   cs_real_t *dt = CS_F_(dt)->val;
235   cs_real_t ctot_tmp = 0.;
236   cs_field_t *sca, *precip, *solub;
237   cs_gwf_soilwater_partition_t sorption_scal;
238 
239   int key_part = cs_field_key_id("gwf_soilwater_partition");
240   int key_pre = cs_field_key_id("gwf_precip_concentration_id");
241 
242   /* Get scalar, precipitation and solubility index fields */
243   sca = cs_field_by_id(f_id);
244   cs_field_get_key_struct(sca, key_part, &sorption_scal);
245   precip = cs_field_by_id(cs_field_get_key_int(sca, key_pre));
246   solub = cs_field_by_id(sorption_scal.imxsol);
247 
248   /* Get first-order decay rate */
249   const int dr_id =  cs_field_key_id("fo_decay_rate");
250   cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
251 
252   for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
253     /* Apply the radioactive decay rate to precipitation concentration before
254        clipping. To be consistent with decay treatment of liquid concentration,
255        an implicit scheme is applied. */
256     precip->val[c_id] *= 1. / (1. + decay_rate * dt[c_id]);
257 
258     /* Clipping of solute concentration, update of precip. concentration */
259     ctot_tmp          = sca->val[c_id] + precip->val[c_id];
260     sca->val[c_id]    = CS_MIN(ctot_tmp, solub->val[c_id]);
261     precip->val[c_id] = CS_MAX(0., ctot_tmp - solub->val[c_id]);
262   }
263 }
264 
265 /*----------------------------------------------------------------------------*/
266 /*!
267  * \brief Take into account kinetic chemical reaction in evolution equation
268  *        of total liquid concentration.
269  *
270  * \f[
271  * S^{n+1} = S^n \exp(- (k^{-} + decay_rate) * \Delta t)
272  *           - C^n \dfrac{k^{+}}{k^{-}}
273  *             \left(\exp(- (k^{-} + decay_rate) * \Delta t) - 1 \right)
274  * \f]
275  *
276  * \param[in]   f_id   field index of scalar which properties are updated
277  */
278  /*---------------------------------------------------------------------------*/
279 
cs_gwf_kinetic_reaction(const int f_id,cs_real_t * ts_imp,cs_real_t * ts_exp)280 void cs_gwf_kinetic_reaction(const int  f_id,
281                              cs_real_t *ts_imp,
282                              cs_real_t *ts_exp)
283 {
284   const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
285   const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
286   const cs_real_t *vol = mesh_quantities->cell_vol;
287 
288   cs_real_t *dt = CS_F_(dt)->val;
289 
290   /* Get soil density */
291   cs_field_t *rosoil = cs_field_by_name("soil_density");
292 
293   /* Get scalar and sorbed concentration fields */
294   cs_field_t *sca = cs_field_by_id(f_id);
295   int key_sorb = cs_field_key_id("gwf_sorbed_concentration_id");
296   cs_field_t *sorb = cs_field_by_id(cs_field_get_key_int(sca, key_sorb));
297 
298   /* Get first-order decay rate */
299   const int dr_id =  cs_field_key_id("fo_decay_rate");
300   cs_real_t decay_rate = cs_field_get_key_double(sca, dr_id);
301 
302   /* Soil-water partition structure */
303   int key_part = cs_field_key_id("gwf_soilwater_partition");
304   cs_gwf_soilwater_partition_t sorption_scal;
305   cs_field_get_key_struct(sca, key_part, &sorption_scal);
306 
307   /* Get k+ and k- fields */
308   cs_field_t *kp = cs_field_by_id(sorption_scal.ikp);
309   cs_field_t *km = cs_field_by_id(sorption_scal.ikm);
310 
311   /* First choice : analytical resolution */
312   if (sorption_scal.anai) {
313     for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++)
314       /* General case (reversible sorption or presence of decay rate) */
315       if (km->val[c_id] + decay_rate > cs_math_epzero) {
316         cs_real_t expkdt = exp(-(km->val[c_id]+decay_rate) * dt[c_id]);
317         cs_real_t kpskm = kp->val[c_id] / (km->val[c_id] + decay_rate);
318         ts_exp[c_id] += - vol[c_id]
319                          *(decay_rate * rosoil->val[c_id] * sorb->val[c_id]
320                            + rosoil->val[c_id]/dt[c_id] * (1-expkdt)
321                             *(kpskm*sca->val[c_id] - sorb->val[c_id]));
322         ts_imp[c_id] +=  vol[c_id] * rosoil->val[c_id] / dt[c_id]
323                         *(1-expkdt)*kpskm;
324       }
325       else { /* case of irreversible sorption without decay rate */
326         cs_real_t rokpl = rosoil->val[c_id] * kp->val[c_id];
327         ts_exp[c_id] += - vol[c_id] * rokpl * sca->val[c_id];
328         ts_imp[c_id] += + vol[c_id] * rokpl;
329       }
330   }
331   else { /* Second choice : explicit */
332     for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
333       ts_exp[c_id] += vol[c_id] * rosoil->val[c_id]
334                      *(  km->val[c_id] * sorb->val[c_id]
335                        - kp->val[c_id] * sca->val[c_id]);
336       ts_imp[c_id] += vol[c_id] * rosoil->val[c_id] * kp->val[c_id];
337     }
338   }
339 }
340 
341 /*----------------------------------------------------------------------------*/
342 
343 END_C_DECLS
344