1 /*============================================================================
2 * Main functions dedicated to soil management in groundwater flows
3 *============================================================================*/
4
5 /* VERS */
6
7 /*
8 This file is part of Code_Saturne, a general-purpose CFD tool.
9
10 Copyright (C) 1998-2021 EDF S.A.
11
12 This program is free software; you can redistribute it and/or modify it under
13 the terms of the GNU General Public License as published by the Free Software
14 Foundation; either version 2 of the License, or (at your option) any later
15 version.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
20 details.
21
22 You should have received a copy of the GNU General Public License along with
23 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
24 Street, Fifth Floor, Boston, MA 02110-1301, USA.
25 */
26
27 /*----------------------------------------------------------------------------*/
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <ctype.h>
35 #include <float.h>
36 #include <math.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include <bft_mem.h>
45 #include <bft_printf.h>
46
47 #include "cs_field.h"
48 #include "cs_hodge.h"
49 #include "cs_log.h"
50 #include "cs_math.h"
51 #include "cs_mesh_location.h"
52 #include "cs_parall.h"
53 #include "cs_param_types.h"
54 #include "cs_post.h"
55 #include "cs_prototypes.h"
56 #include "cs_reco.h"
57 #include "cs_volume_zone.h"
58
59 /*----------------------------------------------------------------------------
60 * Header for the current file
61 *----------------------------------------------------------------------------*/
62
63 #include "cs_gwf_soil.h"
64
65 /*----------------------------------------------------------------------------*/
66
67 BEGIN_C_DECLS
68
69 /*!
70 \file cs_gwf_soil.c
71
72 \brief Main functions dedicated to soil management in groundwater flows when
73 using CDO schemes
74
75 */
76
77 /*============================================================================
78 * Local macro definitions
79 *============================================================================*/
80
81 #define CS_GWF_SOIL_DBG 0
82
83 /*============================================================================
84 * Structure definitions
85 *============================================================================*/
86
87 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
88
89 /*============================================================================
90 * Static global variables
91 *============================================================================*/
92
93 static const char _err_empty_soil[] =
94 " Stop execution. The structure related to a soil is empty.\n"
95 " Please check your settings.\n";
96
97 static int _n_soils = 0;
98 static cs_gwf_soil_t **_soils = NULL;
99
100 /* The following array enables to get the soil id related to each cell.
101 The array size is equal to n_cells */
102 static short int *_cell2soil_ids = NULL;
103
104 /*============================================================================
105 * Private function prototypes
106 *============================================================================*/
107
108 /*----------------------------------------------------------------------------*/
109 /*!
110 * \brief Function that compute the new values of the properties related to
111 * a soil with a Van Genuchten-Mualen.
112 * Case of an isotropic permeability and an unsteady Richards eq.
113 *
114 * \param[in] t_eval time at which one performs the evaluation
115 * \param[in] mesh pointer to a cs_mesh_t structure
116 * \param[in] connect pointer to a cs_cdo_connect_t structure
117 * \param[in] quant pointer to a cs_cdo_quantities_t structure
118 * \param[in] head_values array of values for head used in law
119 * \param[in] zone pointer to a cs_zone_t
120 * \param[in, out] soil_context pointer to a structure cast on-the-fly
121 */
122 /*----------------------------------------------------------------------------*/
123
124 static inline void
_update_soil_genuchten_iso(const cs_real_t t_eval,const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_zone_t * zone,void * soil_context)125 _update_soil_genuchten_iso(const cs_real_t t_eval,
126 const cs_mesh_t *mesh,
127 const cs_cdo_connect_t *connect,
128 const cs_cdo_quantities_t *quant,
129 const cs_zone_t *zone,
130 void *soil_context)
131 {
132 CS_UNUSED(t_eval);
133 CS_UNUSED(mesh);
134 CS_UNUSED(connect);
135 CS_UNUSED(quant);
136
137 cs_gwf_soil_context_genuchten_t *sc = soil_context;
138
139 /* Only isotropic values are considered in this case */
140
141 const double iso_satval = sc->saturated_permeability[0][0];
142 const double delta_moisture = sc->saturated_moisture - sc->residual_moisture;
143 const cs_real_t *head = sc->head_values;
144
145 /* Retrieve field values associated to properties to update */
146
147 cs_real_t *permeability = sc->permeability_values;
148 cs_real_t *moisture = sc->moisture_values;
149 cs_real_t *capacity = sc->capacity_values;
150
151 assert(capacity != NULL && permeability != NULL && moisture != NULL);
152
153 /* Main loop on cells belonging to this soil */
154
155 # pragma omp parallel for if (zone->n_elts > CS_THR_MIN) \
156 shared(head, zone, sc, permeability, moisture, capacity) \
157 firstprivate(iso_satval, delta_moisture)
158 for (cs_lnum_t i = 0; i < zone->n_elts; i++) {
159
160 const cs_lnum_t c_id = zone->elt_ids[i];
161 const cs_real_t h = head[c_id];
162
163 if (h < 0) { /* S_e(h) = [1 + |alpha*h|^n]^(-m) */
164
165 const double coef = pow(fabs(sc->scale * h), sc->n);
166 const double se = pow(1 + coef, -sc->m);
167 const double se_pow_overm = pow(se, 1/sc->m);
168 const double coef_base = 1 - pow(1 - se_pow_overm, sc->m);
169
170 /* Set the permeability value */
171
172 permeability[c_id] =
173 iso_satval* pow(se, sc->tortuosity) * coef_base*coef_base;
174
175 /* Set the moisture content */
176
177 moisture[c_id] = se*delta_moisture + sc->residual_moisture;
178
179 /* Set the soil capacity */
180
181 const double ccoef = -sc->n * sc->m * delta_moisture;
182 const double se_m1 = se/(1. + coef);
183
184 capacity[c_id] = ccoef * coef/h * se_m1;
185
186 }
187 else {
188
189 /* Set the permeability value to the saturated values */
190
191 permeability[c_id] = iso_satval;
192
193 /* Set the moisture content (Se = 1 in this case)*/
194
195 moisture[c_id] = delta_moisture + sc->residual_moisture;
196
197 /* Set the soil capacity */
198
199 capacity[c_id] = 0.;
200
201 }
202
203 } /* Loop on selected cells */
204
205 }
206
207 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
208
209 /*============================================================================
210 * Public function prototypes
211 *============================================================================*/
212
213 /*----------------------------------------------------------------------------*/
214 /*!
215 * \brief Get the number of allocated soils
216 *
217 * \return the number of allocated soils
218 */
219 /*----------------------------------------------------------------------------*/
220
221 int
cs_gwf_get_n_soils(void)222 cs_gwf_get_n_soils(void)
223 {
224 return _n_soils;
225 }
226
227 /*----------------------------------------------------------------------------*/
228 /*!
229 * \brief Retrieve a soil structure from its id
230 *
231 * \param[in] id id to look for
232 *
233 * \return a pointer to a cs_gwf_soil_t structure
234 */
235 /*----------------------------------------------------------------------------*/
236
237 cs_gwf_soil_t *
cs_gwf_soil_by_id(int id)238 cs_gwf_soil_by_id(int id)
239 {
240 if (id > -1 && id < _n_soils)
241 return _soils[id];
242 else
243 return NULL;
244 }
245
246 /*----------------------------------------------------------------------------*/
247 /*!
248 * \brief Retrieve a soil structure from its name
249 *
250 * \param[in] name name to look for
251 *
252 * \return a pointer to a cs_gwf_soil_t structure
253 */
254 /*----------------------------------------------------------------------------*/
255
256 cs_gwf_soil_t *
cs_gwf_soil_by_name(const char * name)257 cs_gwf_soil_by_name(const char *name)
258 {
259 if (name == NULL)
260 return NULL;
261
262 for (int i = 0; i < _n_soils; i++) {
263
264 cs_gwf_soil_t *s = _soils[i];
265 const cs_zone_t *zone = cs_volume_zone_by_id(s->zone_id);
266
267 if (strcmp(zone->name, name) == 0)
268 return s;
269 }
270
271 /* Not found among the list */
272 return NULL;
273 }
274
275 /*----------------------------------------------------------------------------*/
276 /*!
277 * \brief Get the saturated moisture for the given soil id
278 *
279 * \param[in] soil_id id of the requested soil
280 *
281 * \return the value of the saturated moisture
282 */
283 /*----------------------------------------------------------------------------*/
284
285 cs_real_t
cs_gwf_soil_get_saturated_moisture(int soil_id)286 cs_gwf_soil_get_saturated_moisture(int soil_id)
287 {
288 cs_gwf_soil_t *soil = cs_gwf_soil_by_id(soil_id);
289
290 if (soil == NULL)
291 bft_error(__FILE__, __LINE__, 0, "%s: Empty soil.\n", __func__);
292
293 return soil->saturated_moisture;
294 }
295
296 /*----------------------------------------------------------------------------*/
297 /*!
298 * \brief Check if all soils have been set as CS_GWF_SOIL_SATURATED
299 *
300 * \return true or false
301 */
302 /*----------------------------------------------------------------------------*/
303
304 bool
cs_gwf_soil_all_saturated(void)305 cs_gwf_soil_all_saturated(void)
306 {
307 for (int soil_id = 0; soil_id < _n_soils; soil_id++) {
308
309 const cs_gwf_soil_t *soil = _soils[soil_id];
310 if (soil->model != CS_GWF_SOIL_SATURATED)
311 return false;
312
313 }
314
315 return true;
316 }
317
318 /*----------------------------------------------------------------------------*/
319 /*!
320 * \brief Check that at least one soil has been defined and the model of soil
321 * exists.
322 * Raise an error if a problem is encoutered.
323 */
324 /*----------------------------------------------------------------------------*/
325
326 void
cs_gwf_soil_check(void)327 cs_gwf_soil_check(void)
328 {
329 if (_n_soils < 1)
330 bft_error(__FILE__, __LINE__, 0,
331 "%s: Groundwater module is activated but no soil is defined.",
332 __func__);
333 if (_soils == NULL)
334 bft_error(__FILE__, __LINE__, 0,
335 "%s: The soil structure is not allocated whereas %d soils"
336 " have been added.\n", __func__, _n_soils);
337
338 for (int i = 0; i < _n_soils; i++) {
339
340 if (_soils[i]->model == CS_GWF_SOIL_N_HYDRAULIC_MODELS) {
341 const cs_zone_t *z = cs_volume_zone_by_id(_soils[i]->zone_id);
342 bft_error(__FILE__, __LINE__, 0,
343 "%s: Invalid model of soil attached to zone %s\n",
344 __func__, z->name);
345 }
346
347 } /* Loop on soils */
348 }
349
350 /*----------------------------------------------------------------------------*/
351 /*!
352 * \brief Create a new cs_gwf_soil_t structure and add it to the array of
353 * soils. An initialization by default of all members is performed.
354 *
355 * \param[in] zone pointer to a volume zone structure
356 * \param[in] model type of modelling for the hydraulic behavior
357 * \param[in] perm_type type of permeability (iso/anisotropic)
358 * \param[in] sat_moisture value of the saturated moisture content
359 * \param[in] bulk_density value of the mass density
360 *
361 * \return a pointer to the new allocated structure
362 */
363 /*----------------------------------------------------------------------------*/
364
365 cs_gwf_soil_t *
cs_gwf_soil_create(const cs_zone_t * zone,cs_gwf_soil_hydraulic_model_t model,cs_property_type_t perm_type,double sat_moisture,double bulk_density)366 cs_gwf_soil_create(const cs_zone_t *zone,
367 cs_gwf_soil_hydraulic_model_t model,
368 cs_property_type_t perm_type,
369 double sat_moisture,
370 double bulk_density)
371 {
372 cs_gwf_soil_t *soil = NULL;
373
374 BFT_MALLOC(soil, 1, cs_gwf_soil_t);
375
376 soil->id = _n_soils;
377 soil->bulk_density = bulk_density;
378 soil->saturated_moisture = sat_moisture;
379 soil->model = model;
380
381 /* Attached a volume zone to the current soil */
382
383 assert(zone != NULL);
384 soil->zone_id = zone->id;
385 soil->context = NULL;
386
387 soil->update_properties = NULL;
388 soil->free_context = NULL;
389
390 switch (model) {
391
392 case CS_GWF_SOIL_SATURATED:
393 {
394 cs_gwf_soil_context_saturated_t *sc = NULL;
395
396 BFT_MALLOC(sc, 1, cs_gwf_soil_context_saturated_t);
397
398 /* Default initialization */
399
400 sc->saturated_moisture = sat_moisture;
401
402 for (int ki = 0; ki < 3; ki++)
403 for (int kj = 0; kj < 3; kj++)
404 sc->saturated_permeability[ki][kj] = 0.0;
405
406 sc->saturated_permeability[0][0] = 1.0;
407 sc->saturated_permeability[1][1] = 1.0;
408 sc->saturated_permeability[2][2] = 1.0;
409
410 soil->context = sc;
411 }
412 break;
413
414 case CS_GWF_SOIL_GENUCHTEN:
415 {
416 cs_gwf_soil_context_genuchten_t *sc = NULL;
417
418 BFT_MALLOC(sc, 1, cs_gwf_soil_context_genuchten_t);
419
420 sc->residual_moisture = 0.;
421 sc->saturated_moisture = sat_moisture;
422
423 for (int ki = 0; ki < 3; ki++)
424 for (int kj = 0; kj < 3; kj++)
425 sc->saturated_permeability[ki][kj] = 0.0;
426
427 sc->saturated_permeability[0][0] = 1.0;
428 sc->saturated_permeability[1][1] = 1.0;
429 sc->saturated_permeability[2][2] = 1.0;
430
431 sc->n = 1.25;
432 sc->m = 1 - 1./sc->n;
433 sc->scale = 1.;
434 sc->tortuosity = 1.;
435
436 /* Pointer to property values will be set after */
437
438 sc->permeability_values = NULL;
439 sc->head_values = NULL;
440 sc->moisture_values = NULL;
441 sc->capacity_values = NULL;
442
443 soil->context = sc;
444
445 if (perm_type & CS_PROPERTY_ISO)
446 soil->update_properties = _update_soil_genuchten_iso;
447 else
448 bft_error(__FILE__, __LINE__, 0,
449 "%s: Invalid type of property for the permeability.\n"
450 " Please check your settings.", __func__);
451 }
452 break;
453
454 case CS_GWF_SOIL_USER:
455 break;
456
457 default:
458 bft_error(__FILE__, __LINE__, 0,
459 "%s: Invalid type of soil model\n", __func__);
460 break; /* Nothing to do */
461
462 } /* Switch on soil modeling */
463
464 /* Store the new soils in the soil array */
465
466 _n_soils++;
467 BFT_REALLOC(_soils, _n_soils, cs_gwf_soil_t *);
468 _soils[soil->id] = soil;
469
470 return soil;
471 }
472
473 /*----------------------------------------------------------------------------*/
474 /*!
475 * \brief Build an array storing the associated soil for each cell
476 *
477 * \param[in] n_cells number of cells
478 */
479 /*----------------------------------------------------------------------------*/
480
481 void
cs_gwf_build_cell2soil(cs_lnum_t n_cells)482 cs_gwf_build_cell2soil(cs_lnum_t n_cells)
483 {
484 BFT_MALLOC(_cell2soil_ids, n_cells, short int);
485
486 if (_n_soils == 1)
487 memset(_cell2soil_ids, 0, sizeof(short int)*n_cells);
488
489 else {
490
491 assert(_n_soils > 1);
492 # pragma omp parallel for if (n_cells > CS_THR_MIN)
493 for (cs_lnum_t j = 0; j < n_cells; j++)
494 _cell2soil_ids[j] = -1; /* unset by default */
495
496 for (int soil_id = 0; soil_id < _n_soils; soil_id++) {
497
498 const cs_gwf_soil_t *soil = _soils[soil_id];
499 const cs_zone_t *z = cs_volume_zone_by_id(soil->zone_id);
500
501 assert(z != NULL);
502
503 # pragma omp parallel for if (z->n_elts > CS_THR_MIN)
504 for (cs_lnum_t j = 0; j < z->n_elts; j++)
505 _cell2soil_ids[z->elt_ids[j]] = soil_id;
506
507 } /* Loop on soils */
508
509 /* Chcek if every cells is associated to a soil */
510
511 for (cs_lnum_t j = 0; j < n_cells; j++)
512 if (_cell2soil_ids[j] == -1)
513 bft_error(__FILE__, __LINE__, 0,
514 " %s: At least cell %ld has no related soil.\n",
515 __func__, (long)j);
516
517 } /* n_soils > 1 */
518
519 }
520
521 /*----------------------------------------------------------------------------*/
522 /*!
523 * \brief Get the array storing the associated soil for each cell
524 */
525 /*----------------------------------------------------------------------------*/
526
527 const short int *
cs_gwf_get_cell2soil(void)528 cs_gwf_get_cell2soil(void)
529 {
530 return _cell2soil_ids;
531 }
532
533 /*----------------------------------------------------------------------------*/
534 /*!
535 * \brief Free all cs_gwf_soil_t structures
536 */
537 /*----------------------------------------------------------------------------*/
538
539 void
cs_gwf_soil_free_all(void)540 cs_gwf_soil_free_all(void)
541 {
542 if (_n_soils < 1)
543 return;
544 assert(_soils != NULL);
545
546 for (int i = 0; i < _n_soils; i++) {
547
548 cs_gwf_soil_t *soil = _soils[i];
549
550 if (soil->free_context != NULL)
551 soil->free_context(&(soil->context));
552
553 if (soil->context != NULL) {
554
555 switch (soil->model) {
556
557 case CS_GWF_SOIL_SATURATED:
558 {
559 cs_gwf_soil_context_saturated_t *sc = soil->context;
560
561 BFT_FREE(sc);
562 sc = NULL;
563 }
564 break;
565
566 case CS_GWF_SOIL_GENUCHTEN:
567 {
568 cs_gwf_soil_context_genuchten_t *sc = soil->context;
569
570 BFT_FREE(sc);
571 sc = NULL;
572 }
573 break;
574
575 default:
576 cs_base_warn(__FILE__, __LINE__);
577 bft_printf("%s: The context structure of a soil may not be freed.\n",
578 __func__);
579 break;
580
581 } /* Switch on predefined soil context */
582
583 }
584 BFT_FREE(soil);
585
586 } /* Loop on soils */
587
588 BFT_FREE(_soils);
589 BFT_FREE(_cell2soil_ids);
590 }
591
592 /*----------------------------------------------------------------------------*/
593 /*!
594 * \brief Summary of the settings related to all cs_gwf_soil_t structures
595 */
596 /*----------------------------------------------------------------------------*/
597
598 void
cs_gwf_soil_log_setup(void)599 cs_gwf_soil_log_setup(void)
600 {
601 cs_log_printf(CS_LOG_SETUP, " * GWF | Number of soils: %d\n", _n_soils);
602
603 char meta[64];
604 for (int i = 0; i < _n_soils; i++) {
605
606 const cs_gwf_soil_t *soil = _soils[i];
607 const cs_zone_t *z = cs_volume_zone_by_id(soil->zone_id);
608
609 cs_log_printf(CS_LOG_SETUP, "\n Soil.%d | Zone: %s\n",
610 soil->id, z->name);
611 cs_log_printf(CS_LOG_SETUP, "\n Soil.%d | Bulk.density: %6.3e\n",
612 soil->id, soil->bulk_density);
613
614 sprintf(meta, " Soil.%d |", soil->id);
615
616 switch (soil->model) {
617
618 case CS_GWF_SOIL_GENUCHTEN:
619 {
620 const cs_gwf_soil_context_genuchten_t *sc = soil->context;
621
622 cs_log_printf(CS_LOG_SETUP, "%s Model: VanGenuchten-Mualen\n", meta);
623 cs_log_printf(CS_LOG_SETUP, "%s Parameters:", meta);
624 cs_log_printf(CS_LOG_SETUP,
625 " residual_moisture %5.3e", sc->residual_moisture);
626 cs_log_printf(CS_LOG_SETUP,
627 " saturated_moisture %5.3e\n", sc->saturated_moisture);
628 cs_log_printf(CS_LOG_SETUP, "%s Parameters:", meta);
629 cs_log_printf(CS_LOG_SETUP, " n= %f, scale= %f, tortuosity= %f\n",
630 sc->n, sc->scale, sc->tortuosity);
631 cs_log_printf(CS_LOG_SETUP, "%s Saturated permeability\n", meta);
632 cs_log_printf(CS_LOG_SETUP, "%s [%-4.2e %4.2e %4.2e;\n", meta,
633 sc->saturated_permeability[0][0],
634 sc->saturated_permeability[0][1],
635 sc->saturated_permeability[0][2]);
636 cs_log_printf(CS_LOG_SETUP, "%s %-4.2e %4.2e %4.2e;\n", meta,
637 sc->saturated_permeability[1][0],
638 sc->saturated_permeability[1][1],
639 sc->saturated_permeability[1][2]);
640 cs_log_printf(CS_LOG_SETUP, "%s %-4.2e %4.2e %4.2e]\n", meta,
641 sc->saturated_permeability[2][0],
642 sc->saturated_permeability[2][1],
643 sc->saturated_permeability[2][2]);
644 }
645 break;
646
647 case CS_GWF_SOIL_SATURATED:
648 {
649 const cs_gwf_soil_context_saturated_t *sc = soil->context;
650
651 cs_log_printf(CS_LOG_SETUP, "%s Model: Saturated\n", meta);
652 cs_log_printf(CS_LOG_SETUP, "%s Parameters", meta);
653 cs_log_printf(CS_LOG_SETUP,
654 " saturated_moisture %5.3e\n", sc->saturated_moisture);
655 cs_log_printf(CS_LOG_SETUP, "%s Saturated permeability\n", meta);
656 cs_log_printf(CS_LOG_SETUP, "%s [%-4.2e %4.2e %4.2e;\n", meta,
657 sc->saturated_permeability[0][0],
658 sc->saturated_permeability[0][1],
659 sc->saturated_permeability[0][2]);
660 cs_log_printf(CS_LOG_SETUP, "%s %-4.2e %4.2e %4.2e;\n", meta,
661 sc->saturated_permeability[1][0],
662 sc->saturated_permeability[1][1],
663 sc->saturated_permeability[1][2]);
664 cs_log_printf(CS_LOG_SETUP, "%s %-4.2e %4.2e %4.2e]\n", meta,
665 sc->saturated_permeability[2][0],
666 sc->saturated_permeability[2][1],
667 sc->saturated_permeability[2][2]);
668 }
669 break;
670
671 case CS_GWF_SOIL_USER:
672 cs_log_printf(CS_LOG_SETUP, "%s Model: User-defined\n", meta);
673 break;
674
675 default:
676 bft_error(__FILE__, __LINE__, 0,
677 " Invalid model for groundwater module.\n"
678 " Please check your settings.");
679
680 } /* Switch model */
681
682 } /* Loop on soils */
683
684 cs_log_printf(CS_LOG_SETUP, "\n");
685 }
686
687 /*----------------------------------------------------------------------------*/
688 /*!
689 * \brief Set a soil defined by a saturated hydraulic model and attached to
690 * an isotropic permeability
691 *
692 * \param[in, out] soil pointer to a cs_gwf_soil_t structure
693 * \param[in] k_s value of the saturated permeability
694 */
695 /*----------------------------------------------------------------------------*/
696
697 void
cs_gwf_soil_set_iso_saturated(cs_gwf_soil_t * soil,double k_s)698 cs_gwf_soil_set_iso_saturated(cs_gwf_soil_t *soil,
699 double k_s)
700 {
701 if (soil == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_soil));
702
703 cs_gwf_soil_context_saturated_t *sc = soil->context;
704
705 if (soil->model != CS_GWF_SOIL_SATURATED)
706 bft_error(__FILE__, __LINE__, 0,
707 "%s: soil model is not saturated\n", __func__);
708 if (sc == NULL)
709 bft_error(__FILE__, __LINE__, 0,
710 "%s: soil context not allocated\n", __func__);
711
712 /* Default initialization is the identity matrix */
713
714 for (int i = 0; i < 3; i++)
715 sc->saturated_permeability[i][i] = k_s;
716 }
717
718 /*----------------------------------------------------------------------------*/
719 /*!
720 * \brief Set a soil defined by a saturated hydraulic model and attached to
721 * an anisotropic permeability
722 *
723 * \param[in, out] soil pointer to a cs_gwf_soil_t structure
724 * \param[in] k_s value of the anisotropic saturated permeability
725 */
726 /*----------------------------------------------------------------------------*/
727
728 void
cs_gwf_soil_set_aniso_saturated(cs_gwf_soil_t * soil,double k_s[3][3])729 cs_gwf_soil_set_aniso_saturated(cs_gwf_soil_t *soil,
730 double k_s[3][3])
731 {
732 if (soil == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_soil));
733
734 cs_gwf_soil_context_saturated_t *sc = soil->context;
735
736 if (soil->model != CS_GWF_SOIL_SATURATED)
737 bft_error(__FILE__, __LINE__, 0,
738 "%s : soil model is not saturated\n", __func__);
739 if (sc == NULL)
740 bft_error(__FILE__, __LINE__, 0,
741 "%s: soil context not allocated\n", __func__);
742
743 for (int i = 0; i < 3; i++)
744 for (int j = 0; j < 3; j++)
745 sc->saturated_permeability[i][j] = k_s[i][j];
746 }
747
748 /*----------------------------------------------------------------------------*/
749 /*!
750 * \brief Set a soil defined by a Van Genuchten-Mualen hydraulic model and
751 * attached to an isotropic saturated permeability
752 *
753 * The (effective) liquid saturation (also called moisture content)
754 * follows the identity
755 * S_l,eff = (S_l - theta_r)/(theta_s - theta_r)
756 * = (1 + |alpha . h|^n)^(-m)
757 *
758 * The isotropic relative permeability is defined as:
759 * k_r = S_l,eff^L * (1 - (1 - S_l,eff^(1/m))^m))^2
760 * where m = 1 - 1/n
761 *
762 * \param[in, out] soil pointer to a cs_gwf_soil_t structure
763 * \param[in] k_s value of the isotropic saturated permeability
764 * \param[in] theta_r residual moisture
765 * \param[in] alpha scale parameter (in m^-1)
766 * \param[in] n shape parameter
767 * \param[in] L turtuosity parameter
768 */
769 /*----------------------------------------------------------------------------*/
770
771 void
cs_gwf_soil_set_iso_genuchten(cs_gwf_soil_t * soil,double k_s,double theta_r,double alpha,double n,double L)772 cs_gwf_soil_set_iso_genuchten(cs_gwf_soil_t *soil,
773 double k_s,
774 double theta_r,
775 double alpha,
776 double n,
777 double L)
778 {
779 if (soil == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_soil));
780
781 cs_gwf_soil_context_genuchten_t *sc = soil->context;
782
783 if (soil->model != CS_GWF_SOIL_GENUCHTEN)
784 bft_error(__FILE__, __LINE__, 0,
785 "%s: soil model is not Van Genuchten\n", __func__);
786 if (sc == NULL)
787 bft_error(__FILE__, __LINE__, 0,
788 "%s: soil context not allocated\n", __func__);
789 if (n <= FLT_MIN)
790 bft_error(__FILE__, __LINE__, 0,
791 "%s: Invalid value for n = %6.4e (the shape parameter).\n"
792 "This value should be > 0.\n", __func__, n);
793
794 sc->residual_moisture = theta_r;
795
796 /* Default initialization is the identity matrix */
797
798 for (int i = 0; i < 3; i++)
799 sc->saturated_permeability[i][i] = k_s;
800
801 /* Additional advanced settings */
802
803 sc->n = n;
804 sc->m = 1 - 1/sc->n;
805 sc->scale = alpha;
806 sc->tortuosity = L;
807 }
808
809 /*----------------------------------------------------------------------------*/
810 /*!
811 * \brief Set a soil defined by a Van Genuchten-Mualen hydraulic model and
812 * attached to an anisotropic saturated permeability.
813 *
814 * The (effective) liquid saturation (also called moisture content)
815 * follows the identity
816 * S_l,eff = (S_l - theta_r)/(theta_s - theta_r)
817 * = (1 + |alpha . h|^n)^(-m)
818 *
819 * The isotropic relative permeability is defined as:
820 * k_r = S_l,eff^L * (1 - (1 - S_l,eff^(1/m))^m))^2
821 * where m = 1 - 1/n
822 *
823 * \param[in, out] soil pointer to a cs_gwf_soil_t structure
824 * \param[in] k_s value of the isotropic saturated permeability
825 * \param[in] theta_r residual moisture/liquid saturation
826 * \param[in] alpha scale parameter (in m^-1)
827 * \param[in] n shape parameter
828 * \param[in] L turtuosity parameter
829 */
830 /*----------------------------------------------------------------------------*/
831
832 void
cs_gwf_soil_set_aniso_genuchten(cs_gwf_soil_t * soil,double k_s[3][3],double theta_r,double alpha,double n,double L)833 cs_gwf_soil_set_aniso_genuchten(cs_gwf_soil_t *soil,
834 double k_s[3][3],
835 double theta_r,
836 double alpha,
837 double n,
838 double L)
839 {
840 if (soil == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_soil));
841
842 cs_gwf_soil_context_genuchten_t *sc = soil->context;
843
844 if (soil->model != CS_GWF_SOIL_GENUCHTEN)
845 bft_error(__FILE__, __LINE__, 0,
846 "%s: soil model is not Van Genuchten\n", __func__);
847 if (sc == NULL)
848 bft_error(__FILE__, __LINE__, 0,
849 "%s: soil context not allocated\n", __func__);
850 if (n <= FLT_MIN)
851 bft_error(__FILE__, __LINE__, 0,
852 "%s: Invalid value for n = %6.4e (the shape parameter).\n"
853 "This value should be > 0.\n", __func__, n);
854
855 sc->residual_moisture = theta_r;
856
857 for (int i = 0; i < 3; i++)
858 for (int j = 0; j < 3; j++)
859 sc->saturated_permeability[i][j] = k_s[i][j];
860
861 /* Additional advanced settings */
862
863 sc->n = n;
864 sc->m = 1 - 1/sc->n;
865 sc->scale = alpha;
866 sc->tortuosity = L;
867 }
868
869 /*----------------------------------------------------------------------------*/
870 /*!
871 * \brief Set a soil defined by a user-defined hydraulic model
872 *
873 * \param[in, out] soil pointer to a cs_gwf_soil_t structure
874 * \param[in] context pointer to a structure cast on-the-fly
875 * \param[in] update_func function pointer to update propoerties
876 * \param[in] free_context_func function pointer to free the context
877 */
878 /*----------------------------------------------------------------------------*/
879
880 void
cs_gwf_soil_set_user(cs_gwf_soil_t * soil,void * context,cs_gwf_soil_update_t * update_func,cs_gwf_soil_free_context_t * free_context_func)881 cs_gwf_soil_set_user(cs_gwf_soil_t *soil,
882 void *context,
883 cs_gwf_soil_update_t *update_func,
884 cs_gwf_soil_free_context_t *free_context_func)
885 {
886 if (soil == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_soil));
887
888 if (soil->model != CS_GWF_SOIL_USER)
889 bft_error(__FILE__, __LINE__, 0,
890 " %s: soil model is not user-defined.\n", __func__);
891
892 /* Set function pointers */
893 soil->context = context;
894 soil->update_properties = update_func;
895 soil->free_context = free_context_func;
896 }
897
898 /*----------------------------------------------------------------------------*/
899 /*!
900 * \brief Set the properties of the groundwater flow module in the case where
901 * all soils are considered as saturated.
902 *
903 * \param[in, out] permeability pointer to a cs_property_t structure
904 * \param[in, out] moisture_content pointer to a cs_property_t structure
905 */
906 /*----------------------------------------------------------------------------*/
907
908 void
cs_gwf_soil_saturated_set_properties(cs_property_t * permeability,cs_property_t * moisture_content)909 cs_gwf_soil_saturated_set_properties(cs_property_t *permeability,
910 cs_property_t *moisture_content)
911 {
912 assert(permeability != NULL && moisture_content != NULL);
913
914 for (int i = 0; i < _n_soils; i++) {
915
916 cs_gwf_soil_t *soil = _soils[i];
917
918 if (soil->model != CS_GWF_SOIL_SATURATED)
919 bft_error(__FILE__, __LINE__, 0,
920 " %s: Invalid way of setting soil parameter.\n"
921 " All soils are not considered as saturated.", __func__);
922
923 const cs_zone_t *z = cs_volume_zone_by_id(soil->zone_id);
924
925 cs_gwf_soil_context_saturated_t *sc = soil->context;
926
927 /* Set the permeability */
928
929 if (permeability->type & CS_PROPERTY_ISO)
930 cs_property_def_iso_by_value(permeability,
931 z->name,
932 sc->saturated_permeability[0][0]);
933
934 else if (permeability->type & CS_PROPERTY_ORTHO) {
935
936 cs_real_3_t val = {sc->saturated_permeability[0][0],
937 sc->saturated_permeability[1][1],
938 sc->saturated_permeability[2][2]};
939
940 cs_property_def_ortho_by_value(permeability, z->name, val);
941
942 }
943 else if (permeability->type & CS_PROPERTY_ANISO) {
944
945 cs_property_def_aniso_by_value(permeability,
946 z->name,
947 (double (*)[3])sc->saturated_permeability);
948
949 }
950 else
951 bft_error(__FILE__, __LINE__, 0,
952 " %s: Invalid type of property.\n", __func__);
953
954 /* Set the moisture content */
955
956 assert(fabs(sc->saturated_moisture - soil->saturated_moisture) < FLT_MIN);
957
958 cs_property_def_iso_by_value(moisture_content,
959 z->name,
960 sc->saturated_moisture);
961
962 } /* Loop on soils */
963
964 }
965
966 /*----------------------------------------------------------------------------*/
967 /*!
968 * \brief Set the different arrays used in soil context for a GWF model set
969 * to unsaturated single-phase flows in a porous media.
970 *
971 * \param[in] head pointer to the current head values in cells
972 * \param[in] permeability pointer to the current permeability values
973 * \param[in] moisture_content pointer to the current moisture content values
974 * \param[in] capacity pointer to the current soil capacity values
975 */
976 /*----------------------------------------------------------------------------*/
977
978 void
cs_gwf_soil_uspf_set_arrays(cs_real_t head[],cs_real_t permeability[],cs_real_t moisture_content[],cs_real_t capacity[])979 cs_gwf_soil_uspf_set_arrays(cs_real_t head[],
980 cs_real_t permeability[],
981 cs_real_t moisture_content[],
982 cs_real_t capacity[])
983 {
984 for (int i = 0; i < _n_soils; i++) {
985
986 cs_gwf_soil_t *soil = _soils[i];
987
988 switch (soil->model) {
989
990 case CS_GWF_SOIL_GENUCHTEN:
991 {
992 cs_gwf_soil_context_genuchten_t *sc = soil->context;
993
994 sc->permeability_values = permeability;
995 sc->head_values = head;
996 sc->moisture_values = moisture_content;
997 sc->capacity_values = capacity;
998 }
999 break;
1000
1001 default:
1002 break; /* Do nothing. For user-defined soils, one has to do similar
1003 things in cs_user_gwf.c */
1004
1005 }
1006
1007 } /* Loop on soils */
1008 }
1009
1010 /*----------------------------------------------------------------------------*/
1011 /*!
1012 * \brief Update the soil properties
1013 *
1014 * \param[in] time_eval time at which one evaluates properties
1015 * \param[in] mesh pointer to the mesh structure
1016 * \param[in] connect pointer to the cdo connectivity
1017 * \param[in] quant pointer to the cdo quantities
1018 */
1019 /*----------------------------------------------------------------------------*/
1020
1021 void
cs_gwf_soil_update(cs_real_t time_eval,const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)1022 cs_gwf_soil_update(cs_real_t time_eval,
1023 const cs_mesh_t *mesh,
1024 const cs_cdo_connect_t *connect,
1025 const cs_cdo_quantities_t *quant)
1026 {
1027 for (int i = 0; i < _n_soils; i++) {
1028
1029 const cs_gwf_soil_t *soil = _soils[i];
1030 const cs_zone_t *zone = cs_volume_zone_by_id(soil->zone_id);
1031
1032 switch (soil->model) {
1033
1034 case CS_GWF_SOIL_GENUCHTEN:
1035 case CS_GWF_SOIL_USER:
1036 soil->update_properties(time_eval,
1037 mesh, connect, quant,
1038 zone,
1039 soil->context);
1040 break;
1041
1042 default:
1043 break; /* Do nothing (for instance in the case of a saturated soil which
1044 is constant (steady and uniform) */
1045
1046 }
1047
1048 } /* Loop on soils */
1049
1050 }
1051
1052 /*----------------------------------------------------------------------------*/
1053
1054 END_C_DECLS
1055