1 /*============================================================================
2 * Main functions dedicated to 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 #include "cs_defs.h"
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <ctype.h>
37 #include <float.h>
38 #include <math.h>
39 #include <stdlib.h>
40 #include <string.h>
41
42
43 /*----------------------------------------------------------------------------
44 * Local headers
45 *----------------------------------------------------------------------------*/
46
47 #include <bft_mem.h>
48
49 #include "cs_boundary_zone.h"
50 #include "cs_cdovb_scaleq.h"
51 #include "cs_equation_bc.h"
52 #include "cs_evaluate.h"
53 #include "cs_field.h"
54 #include "cs_gwf_priv.h"
55 #include "cs_hodge.h"
56 #include "cs_log.h"
57 #include "cs_math.h"
58 #include "cs_mesh_location.h"
59 #include "cs_parall.h"
60 #include "cs_param_types.h"
61 #include "cs_physical_constants.h"
62 #include "cs_post.h"
63 #include "cs_reco.h"
64 #include "cs_zone.h"
65
66 #if defined(DEBUG) && !defined(NDEBUG)
67 #include "cs_dbg.h"
68 #endif
69
70 /*----------------------------------------------------------------------------
71 * Header for the current file
72 *----------------------------------------------------------------------------*/
73
74 #include "cs_gwf.h"
75
76 /*----------------------------------------------------------------------------*/
77
78 BEGIN_C_DECLS
79
80 /*!
81 \file cs_gwf.c
82
83 \brief Main functions dedicated to groundwater flows when using CDO schemes
84
85 */
86
87 /*============================================================================
88 * Local macro definitions
89 *============================================================================*/
90
91 /* Redefined names of function from cs_math to get shorter names */
92 #define _dp3 cs_math_3_dot_product
93
94 #define CS_GWF_DBG 0
95
96 /*============================================================================
97 * Local definitions
98 *============================================================================*/
99
100 static const char
101 cs_gwf_model_name[CS_GWF_N_MODEL_TYPES][CS_BASE_STRING_LEN] =
102 { N_("Saturated single-phase model"),
103 N_("Unsaturated single-phase model"),
104 N_("Two-phase model (capillary/gas pressure)")
105 };
106
107 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
108
109 /*============================================================================
110 * Static global variables
111 *============================================================================*/
112
113 static const char _err_empty_gw[] =
114 " Stop execution. The structure related to the groundwater module is empty.\n"
115 " Please check your settings.\n";
116
117 static cs_gwf_t *cs_gwf_main_structure = NULL;
118
119 /*============================================================================
120 * Private function prototypes
121 *============================================================================*/
122
123 /*----------------------------------------------------------------------------*/
124 /*!
125 * \brief Estimate the time at which the evaluation of properties has to be
126 * done
127 *
128 * \param[in] ts pointer to a cs_time_step_t structure
129 * \param[in] eq pointer to an equation structure
130 *
131 * \return the time value at which one has to perform evaluation
132 */
133 /*----------------------------------------------------------------------------*/
134
135 static cs_real_t
_get_time_eval(const cs_time_step_t * ts,cs_equation_t * eq)136 _get_time_eval(const cs_time_step_t *ts,
137 cs_equation_t *eq)
138
139 {
140 cs_gwf_t *gw = cs_gwf_main_structure;
141 cs_real_t time_eval = ts->t_cur;
142
143 const cs_real_t dt_cur = ts->dt[0];
144
145 /* Define the time at which one evaluates the properties */
146
147 cs_param_time_scheme_t time_scheme = cs_equation_get_time_scheme(eq);
148 cs_real_t theta = -1;
149
150 switch (time_scheme) {
151
152 case CS_TIME_SCHEME_STEADY:
153 case CS_TIME_N_SCHEMES:
154
155 /* Look for tracer equations */
156 for (int ieq = 0; ieq < gw->n_tracers; ieq++) {
157
158 cs_equation_t *_eq = gw->tracers[ieq]->eq;
159
160 theta = fmax(theta, cs_equation_get_theta_time_val(_eq));
161
162 }
163
164 if (theta > 0)
165 time_eval = ts->t_cur + theta*dt_cur;
166 else
167 time_eval = ts->t_cur;
168 break;
169
170 default:
171 theta = cs_equation_get_theta_time_val(eq);
172 time_eval = ts->t_cur + theta*dt_cur;
173 break;
174
175 } /* End of switch on the time scheme for the main equation */
176
177 return time_eval;
178 }
179
180 /*----------------------------------------------------------------------------*/
181 /*!
182 * \brief Retrieve the advection field related to the Darcy flux in the liquid
183 * phase
184 *
185 * \param[in] gw pointer to the main (high-level) GWF structure
186 *
187 * \return a pointer to a cs_adv_field_t structure or NULL
188 */
189 /*----------------------------------------------------------------------------*/
190
191 static cs_adv_field_t *
_get_l_adv_field(const cs_gwf_t * gw)192 _get_l_adv_field(const cs_gwf_t *gw)
193 {
194 switch (gw->model) {
195
196 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
197 {
198 cs_gwf_saturated_single_phase_t *mc = gw->model_context;
199
200 return mc->darcy->adv_field;
201 }
202 break;
203
204 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
205 {
206 cs_gwf_unsaturated_single_phase_t *mc = gw->model_context;
207
208 return mc->darcy->adv_field;
209 }
210 break;
211
212 case CS_GWF_MODEL_TWO_PHASE:
213 {
214 cs_gwf_miscible_two_phase_t *mc = gw->model_context;
215
216 if (mc->l_darcy != NULL)
217 return mc->l_darcy->adv_field;
218 }
219 break;
220
221 default:
222 break;
223
224 } /* Switch on the model */
225
226 return NULL;
227 }
228
229 /*----------------------------------------------------------------------------*/
230 /*!
231 * \brief Update head values (pressure head or head values for laws)
232 * Case of single-phase flows in porous media (saturated or not).
233 *
234 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
235 * \param[in] connect pointer to a cs_cdo_connect_t structure
236 * \param[in] richards pointer to the Richards equation
237 * \param[in, out] pressure_head pressure head field
238 * \param[in, out] head_in_law values of the head used in law
239 * \param[in] cur2prev true or false
240 */
241 /*----------------------------------------------------------------------------*/
242
243 static void
_spf_update_head(const cs_cdo_quantities_t * cdoq,const cs_cdo_connect_t * connect,const cs_equation_t * richards,cs_field_t * pressure_head,cs_real_t head_in_law[],bool cur2prev)244 _spf_update_head(const cs_cdo_quantities_t *cdoq,
245 const cs_cdo_connect_t *connect,
246 const cs_equation_t *richards,
247 cs_field_t *pressure_head,
248 cs_real_t head_in_law[],
249 bool cur2prev)
250 {
251 cs_gwf_t *gw = cs_gwf_main_structure;
252 assert(gw != NULL);
253
254 cs_param_space_scheme_t r_scheme = cs_equation_get_space_scheme(richards);
255 cs_field_t *hydraulic_head = cs_equation_get_field(richards);
256
257 if (gw->flag & CS_GWF_RESCALE_HEAD_TO_ZERO_MEAN_VALUE) {
258
259 switch (r_scheme) {
260
261 case CS_SPACE_SCHEME_CDOVB:
262 {
263 assert(hydraulic_head->location_id ==
264 cs_mesh_location_get_id_by_name("vertices"));
265
266 cs_real_t domain_integral =
267 cs_evaluate_scal_domain_integral_by_array(cs_flag_primal_vtx,
268 hydraulic_head->val);
269
270 const cs_real_t mean_value = domain_integral / cdoq->vol_tot;
271
272 # pragma omp parallel for if (cdoq->n_vertices > CS_THR_MIN)
273 for (cs_lnum_t i = 0; i < cdoq->n_vertices; i++)
274 hydraulic_head->val[i] -= mean_value;
275 }
276 break;
277
278 case CS_SPACE_SCHEME_CDOFB:
279 {
280 assert(hydraulic_head->location_id ==
281 cs_mesh_location_get_id_by_name("cells"));
282
283 cs_real_t domain_integral =
284 cs_evaluate_scal_domain_integral_by_array(cs_flag_primal_cell,
285 hydraulic_head->val);
286
287 const cs_real_t mean_value = domain_integral / cdoq->vol_tot;
288 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
289 for (cs_lnum_t i = 0; i < cdoq->n_cells; i++)
290 hydraulic_head->val[i] -= mean_value;
291 }
292 break;
293
294 default:
295 break; /* Nothing to do */
296
297 }
298
299 } /* Rescale hydraulic_head */
300
301 if (gw->flag & CS_GWF_GRAVITATION) { /* Update the pressure head (and if
302 needed head_in_law) */
303
304 cs_physical_constants_t *phys = cs_get_glob_physical_constants();
305
306 /* Sanity checks */
307 if (pressure_head == NULL)
308 bft_error(__FILE__, __LINE__, 0,
309 " The field related to the pressure head is not allocated.");
310
311 /* Copy current field values to previous values */
312 if (cur2prev)
313 cs_field_current_to_previous(pressure_head);
314
315 switch (r_scheme) {
316
317 case CS_SPACE_SCHEME_CDOVB:
318
319 # pragma omp parallel for if (cdoq->n_vertices > CS_THR_MIN)
320 for (cs_lnum_t i = 0; i < cdoq->n_vertices; i++) {
321 const cs_real_t gpot = _dp3(cdoq->vtx_coord + 3*i, phys->gravity);
322 pressure_head->val[i] = hydraulic_head->val[i] - gpot;
323 }
324
325 /* Update head_in_law */
326 if (head_in_law != NULL)
327 cs_reco_pv_at_cell_centers(connect->c2v, cdoq, pressure_head->val,
328 head_in_law);
329 break;
330
331 case CS_SPACE_SCHEME_CDOVCB:
332 {
333 assert(hydraulic_head->location_id ==
334 cs_mesh_location_get_id_by_name("vertices"));
335
336 # pragma omp parallel for if (cdoq->n_vertices > CS_THR_MIN)
337 for (cs_lnum_t i = 0; i < cdoq->n_vertices; i++) {
338 const cs_real_t gpot = _dp3(cdoq->vtx_coord + 3*i, phys->gravity);
339 pressure_head->val[i] = hydraulic_head->val[i] - gpot;
340 }
341
342 /* Update head_in_law */
343 if (head_in_law != NULL) {
344
345 const cs_real_t *hydraulic_head_cells =
346 cs_equation_get_cell_values(richards, false); /* current values */
347
348 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
349 for (cs_lnum_t i = 0; i < cdoq->n_cells; i++) {
350 const cs_real_t gpot = _dp3(cdoq->cell_centers + 3*i,
351 phys->gravity);
352 head_in_law[i] = hydraulic_head_cells[i] - gpot;
353 }
354
355 } /* head_in_law */
356
357 }
358 break;
359
360 case CS_SPACE_SCHEME_CDOFB:
361 case CS_SPACE_SCHEME_HHO_P0:
362
363 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
364 for (cs_lnum_t i = 0; i < cdoq->n_cells; i++) {
365 const cs_real_t gpot = _dp3(cdoq->cell_centers + 3*i, phys->gravity);
366 pressure_head->val[i] = hydraulic_head->val[i] - gpot;
367 }
368 break;
369
370 /* Head in law points either to hydraulic_head->val or pressure_head->val */
371
372 default:
373 bft_error(__FILE__, __LINE__, 0, " Invalid space scheme.");
374
375 } /* Switch on space scheme */
376
377 }
378 else { /* No gravity effect is taken into account */
379
380 if (head_in_law == NULL)
381 return;
382
383 /* Update head_in_law */
384 switch(r_scheme) {
385
386 case CS_SPACE_SCHEME_CDOVB:
387 cs_reco_pv_at_cell_centers(connect->c2v,
388 cdoq,
389 hydraulic_head->val,
390 head_in_law);
391 break;
392
393 case CS_SPACE_SCHEME_CDOVCB:
394 {
395 const cs_real_t *hydraulic_head_cells =
396 cs_equation_get_cell_values(richards, false); /* current values */
397
398 memcpy(head_in_law, hydraulic_head_cells,
399 sizeof(cs_real_t)*cdoq->n_cells);
400 }
401 break;
402
403 default:
404 break; /* Nothing to do for CDO-Fb schemes and HHO schemes */
405
406 } /* Switch on the space scheme related to the Richards equation */
407
408 } /* Gravity is activated or not */
409
410 }
411
412 /*----------------------------------------------------------------------------*/
413 /*!
414 * \brief Compute the steady-state of the groundwater flows module.
415 * Nothing is done if all equations are unsteady.
416 * Case of unstaturated/saturated single-phase flows in porous media.
417 *
418 * \param[in] mesh pointer to a cs_mesh_t structure
419 * \param[in] time_step pointer to a cs_time_step_t structure
420 * \param[in] connect pointer to a cs_cdo_connect_t structure
421 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
422 * \param[in, out] richards pointer to the Richards equation
423 */
424 /*----------------------------------------------------------------------------*/
425
426 static void
_spf_compute_steady_state(const cs_mesh_t * mesh,const cs_time_step_t * time_step,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_equation_t * richards)427 _spf_compute_steady_state(const cs_mesh_t *mesh,
428 const cs_time_step_t *time_step,
429 const cs_cdo_connect_t *connect,
430 const cs_cdo_quantities_t *cdoq,
431 cs_equation_t *richards)
432 {
433 cs_gwf_t *gw = cs_gwf_main_structure;
434 assert(gw != NULL && richards != NULL);
435
436 /* Sanity check */
437
438 assert(gw->model == CS_GWF_MODEL_SATURATED_SINGLE_PHASE);
439 assert(cs_equation_get_type(richards) == CS_EQUATION_TYPE_GROUNDWATER);
440
441 /* Build and solve the linear system related to the Richards equations */
442
443 if (cs_equation_is_steady(richards) ||
444 gw->flag & CS_GWF_FORCE_RICHARDS_ITERATIONS) {
445
446 /* Solve the algebraic system */
447
448 cs_equation_solve_steady_state(mesh, richards);
449
450 /* Update the variables related to the groundwater flow system */
451
452 cs_gwf_update(mesh, connect, cdoq, time_step, true);
453
454 }
455
456 }
457
458 /*----------------------------------------------------------------------------*/
459 /*!
460 * \brief Compute the new state for the groundwater flows module.
461 * Case of unstaturated/saturated single-phase flows in porous media.
462 *
463 * \param[in] mesh pointer to a cs_mesh_t structure
464 * \param[in] time_step pointer to a cs_time_step_t structure
465 * \param[in] connect pointer to a cs_cdo_connect_t structure
466 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
467 * \param[in, out] richards pointer to the Richards equation
468 */
469 /*----------------------------------------------------------------------------*/
470
471 static void
_spf_compute(const cs_mesh_t * mesh,const cs_time_step_t * time_step,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_equation_t * richards)472 _spf_compute(const cs_mesh_t *mesh,
473 const cs_time_step_t *time_step,
474 const cs_cdo_connect_t *connect,
475 const cs_cdo_quantities_t *cdoq,
476 cs_equation_t *richards)
477 {
478 cs_gwf_t *gw = cs_gwf_main_structure;
479
480 assert(gw != NULL && richards != NULL);
481
482 /* Sanity checks */
483
484 assert(gw->model == CS_GWF_MODEL_SATURATED_SINGLE_PHASE ||
485 gw->model == CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE);
486 assert(cs_equation_get_type(richards) == CS_EQUATION_TYPE_GROUNDWATER);
487
488 bool cur2prev = true;
489
490 /* Build and solve the linear system related to the Richards equations */
491
492 if (!cs_equation_is_steady(richards)) {
493
494 /* Solve the algebraic system. By default, a current to previous operation
495 is performed */
496
497 cs_equation_solve(cur2prev, mesh, richards);
498
499 /* Update the variables related to the groundwater flow system */
500
501 cs_gwf_update(mesh, connect, cdoq, time_step, cur2prev);
502
503 }
504
505 }
506
507 /* ==========================================================================
508 * Functions related to the model of saturated single phase flows (SSPF)
509 * ========================================================================== */
510
511 /*----------------------------------------------------------------------------*/
512 /*!
513 * \brief Allocate and initialize the modelling context for the model of
514 * saturated single-phase flows
515 *
516 * \return a pointer to a new allocated cs_gwf_saturated_single_phase_t struct
517 */
518 /*----------------------------------------------------------------------------*/
519
520 static cs_gwf_saturated_single_phase_t *
_sspf_init_context(void)521 _sspf_init_context(void)
522 {
523 cs_gwf_saturated_single_phase_t *mc = NULL;
524
525 BFT_MALLOC(mc, 1, cs_gwf_saturated_single_phase_t);
526
527 mc->pressure_head = NULL;
528
529 /* Create a new equation structure for Richards' equation */
530
531 mc->richards = cs_equation_add("Richards", /* equation name */
532 "hydraulic_head", /* variable name */
533 CS_EQUATION_TYPE_GROUNDWATER,
534 1,
535 CS_PARAM_BC_HMG_NEUMANN);
536
537 /* Define the Darcy flux structure
538 Add an advection field related to the darcian flux stemming from the
539 Richards equation. This advection field is steady since the head is
540 steady in this model.
541 */
542
543 mc->darcy = cs_gwf_darcy_flux_create(cs_flag_dual_face_byc);
544
545 cs_advection_field_status_t adv_status = CS_ADVECTION_FIELD_GWF |
546 CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX | CS_ADVECTION_FIELD_STEADY;
547
548 mc->darcy->adv_field = cs_advection_field_add("darcy_field",
549 adv_status);
550
551 /* Add a property related to the moisture content (should be a constant
552 in each soil) */
553 mc->moisture_content = cs_property_add("moisture_content", CS_PROPERTY_ISO);
554
555 /* Add the diffusion term to the Richards equation by associating the
556 absolute permeability to the diffusion property of the Richards eq. */
557
558 cs_equation_param_t *eqp = cs_equation_get_param(mc->richards);
559 cs_gwf_t *gw = cs_gwf_main_structure;
560
561 cs_equation_add_diffusion(eqp, gw->abs_permeability);
562
563 return mc;
564 }
565
566 /*----------------------------------------------------------------------------*/
567 /*!
568 * \brief Free the modelling context for the model of saturated single-phase
569 * flows
570 *
571 * \param[in, out] p_mc pointer of pointer to the model context structure
572 */
573 /*----------------------------------------------------------------------------*/
574
575 static void
_sspf_free_context(cs_gwf_saturated_single_phase_t ** p_mc)576 _sspf_free_context(cs_gwf_saturated_single_phase_t **p_mc)
577 {
578 if (p_mc == NULL)
579 return;
580 if ((*p_mc) == NULL)
581 return;
582
583 cs_gwf_saturated_single_phase_t *mc = *p_mc;
584
585 cs_gwf_darcy_flux_free(&(mc->darcy));
586
587 BFT_FREE(mc);
588 *p_mc = NULL;
589 }
590
591 /*----------------------------------------------------------------------------*/
592 /*!
593 * \brief Log the setup related to the modelling context of saturated
594 * single-phase flows
595 *
596 * \param[in] mc pointer to the model context structure
597 */
598 /*----------------------------------------------------------------------------*/
599
600 static void
_sspf_log_context(cs_gwf_saturated_single_phase_t * mc)601 _sspf_log_context(cs_gwf_saturated_single_phase_t *mc)
602 {
603 if (mc == NULL)
604 return;
605
606 cs_gwf_darcy_flux_log(mc->darcy);
607 }
608
609 /*----------------------------------------------------------------------------*/
610 /*!
611 * \brief Setup related to the modelling context of saturated single-phase
612 * flows. At this stage, all soils have been defined and equation
613 * parameters are set.
614 *
615 * \param[in, out] mc pointer to the model context structure
616 */
617 /*----------------------------------------------------------------------------*/
618
619 static void
_sspf_init_setup(cs_gwf_saturated_single_phase_t * mc)620 _sspf_init_setup(cs_gwf_saturated_single_phase_t *mc)
621 {
622 cs_gwf_t *gw = cs_gwf_main_structure;
623
624 if (mc == NULL)
625 return;
626 if (mc->richards == NULL)
627 bft_error(__FILE__, __LINE__, 0,
628 "%s: The Richards equation is not defined. Stop execution.\n",
629 __func__);
630 assert(cs_equation_is_steady(mc->richards) == true);
631
632 const int field_mask = CS_FIELD_INTENSIVE | CS_FIELD_VARIABLE | CS_FIELD_CDO;
633 const int c_loc_id = cs_mesh_location_get_id_by_name("cells");
634 const int v_loc_id = cs_mesh_location_get_id_by_name("vertices");
635 const int log_key = cs_field_key_id("log");
636 const int post_key = cs_field_key_id("post_vis");
637 const cs_param_space_scheme_t space_scheme =
638 cs_equation_get_space_scheme(mc->richards);
639
640 /* Handle gravity effects */
641 if (gw->flag & CS_GWF_GRAVITATION) {
642
643 switch (space_scheme) {
644 case CS_SPACE_SCHEME_CDOVB:
645 case CS_SPACE_SCHEME_CDOVCB:
646 mc->pressure_head = cs_field_create("pressure_head",
647 field_mask,
648 v_loc_id,
649 1,
650 false); /* has_previous */
651 break;
652
653 case CS_SPACE_SCHEME_CDOFB:
654 case CS_SPACE_SCHEME_HHO_P0:
655 mc->pressure_head = cs_field_create("pressure_head",
656 field_mask,
657 c_loc_id,
658 1,
659 false); /* has_previous */
660 break;
661
662 default:
663 bft_error(__FILE__, __LINE__, 0, " %s: Invalid space scheme.", __func__);
664 }
665
666 cs_field_set_key_int(mc->pressure_head, log_key, 1);
667 cs_field_set_key_int(mc->pressure_head, post_key, 1);
668
669 } /* Gravitation effect is activated */
670
671 /* Add default post-processing related to groundwater flow module */
672 cs_post_add_time_mesh_dep_output(cs_gwf_extra_post_sspf, gw);
673 }
674
675 /*----------------------------------------------------------------------------*/
676 /*!
677 * \brief Perform the last setup step in the case of a saturated single-phase
678 * flow model in porous media
679 *
680 * \param[in] connect pointer to a cs_cdo_connect_t structure
681 * \param[in] quant pointer to a cs_cdo_quantities_t structure
682 * \param[in, out] mc pointer to the casted model context
683 */
684 /*----------------------------------------------------------------------------*/
685
686 static void
_sspf_finalize_setup(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_gwf_saturated_single_phase_t * mc)687 _sspf_finalize_setup(const cs_cdo_connect_t *connect,
688 const cs_cdo_quantities_t *quant,
689 cs_gwf_saturated_single_phase_t *mc)
690 {
691 cs_gwf_t *gw = cs_gwf_main_structure;
692 assert(gw != NULL && mc != NULL);
693
694 const cs_param_space_scheme_t richards_scheme =
695 cs_equation_get_space_scheme(mc->richards);
696
697 /* Set the Darcian flux (in the volume and at the boundary) */
698
699 cs_gwf_darcy_flux_define(connect, quant, richards_scheme, mc->darcy);
700
701 /* Set permeability, moisture content and soil capacity according to the soil
702 settings */
703
704 cs_gwf_soil_saturated_set_properties(gw->abs_permeability,
705 mc->moisture_content);
706 }
707
708 /*----------------------------------------------------------------------------*/
709 /*!
710 * \brief Perform the update step in the case of a saturated single-phase
711 * flow model in porous media
712 *
713 * \param[in] connect pointer to a cs_cdo_connect_t structure
714 * \param[in] quant pointer to a cs_cdo_quantities_t structure
715 * \param[in] ts pointer to a cs_time_step_t structure
716 * \param[in] cur2prev true or false
717 * \param[in, out] mc pointer to the casted model context
718 *
719 *\return the value at which one has to evaluate the properties
720 */
721 /*----------------------------------------------------------------------------*/
722
723 static cs_real_t
_sspf_updates(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts,bool cur2prev,cs_gwf_saturated_single_phase_t * mc)724 _sspf_updates(const cs_cdo_connect_t *connect,
725 const cs_cdo_quantities_t *quant,
726 const cs_time_step_t *ts,
727 bool cur2prev,
728 cs_gwf_saturated_single_phase_t *mc)
729 {
730 cs_real_t time_eval = ts->t_cur;
731 if (cur2prev)
732 time_eval = _get_time_eval(ts, mc->richards);
733
734 /* Update head */
735
736 _spf_update_head(quant, connect, mc->richards,
737 mc->pressure_head,
738 NULL, /* there is no head_in_law for saturated soils */
739 cur2prev);
740
741 /* Update the advection field related to the groundwater flow module */
742
743 cs_gwf_darcy_flux_update(time_eval, mc->richards, cur2prev, mc->darcy);
744
745 /* Properties are constant in saturated soils. Therefore, there is no need fo
746 an update of the associated properties (permeability and moisture
747 content) */
748
749 return time_eval;
750 }
751
752 /*----------------------------------------------------------------------------*/
753 /*!
754 * \brief Predefined extra-operations for the groundwater flow module in case
755 * of saturated single phase flows in porous media
756 *
757 * \param[in] connect pointer to a cs_cdo_connect_t structure
758 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
759 * \param[in, out] mc pointer to the casted model context
760 */
761 /*----------------------------------------------------------------------------*/
762
763 static void
_sspf_extra_op(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_gwf_saturated_single_phase_t * mc)764 _sspf_extra_op(const cs_cdo_connect_t *connect,
765 const cs_cdo_quantities_t *cdoq,
766 cs_gwf_saturated_single_phase_t *mc)
767 {
768 cs_gwf_t *gw = cs_gwf_main_structure;
769 assert(gw != NULL && mc != NULL);
770
771 if (cs_flag_test(gw->post_flag, CS_GWF_POST_DARCY_FLUX_BALANCE) == false)
772 return; /* Nothing to do */
773
774 cs_gwf_darcy_flux_balance(connect, cdoq,
775 cs_equation_get_param(mc->richards),
776 mc->darcy);
777 }
778
779 /* ==========================================================================
780 * Functions related to the model of unsaturated single phase flows (USPF)
781 * ========================================================================== */
782
783 /*----------------------------------------------------------------------------*/
784 /*!
785 * \brief Allocate and initialize the modelling context for the model of
786 * unsaturated single-phase flows
787 *
788 * \param[in] pty_type type of permeability (iso, ortho...)
789 *
790 * \return a pointer to a new allocated cs_gwf_unsaturated_single_phase_t struct
791 */
792 /*----------------------------------------------------------------------------*/
793
794 static cs_gwf_unsaturated_single_phase_t *
_uspf_init_context(cs_property_type_t pty_type)795 _uspf_init_context(cs_property_type_t pty_type)
796 {
797 cs_gwf_unsaturated_single_phase_t *mc = NULL;
798
799 BFT_MALLOC(mc, 1, cs_gwf_unsaturated_single_phase_t);
800
801 mc->permeability_field = NULL;
802 mc->moisture_field = NULL;
803 mc->capacity_field = NULL;
804 mc->pressure_head = NULL;
805 mc->head_in_law = NULL;
806
807 /* Create a new equation structure for Richards' equation */
808
809 mc->richards = cs_equation_add("Richards", /* equation name */
810 "hydraulic_head", /* variable name */
811 CS_EQUATION_TYPE_GROUNDWATER,
812 1,
813 CS_PARAM_BC_HMG_NEUMANN);
814
815 /* Define the Darcy flux structure
816 Add an advection field related to the darcian flux stemming from the
817 Richards equation. This advection field is steady since the head is
818 steady in this model.
819 */
820
821 mc->darcy = cs_gwf_darcy_flux_create(cs_flag_dual_face_byc);
822
823 cs_advection_field_status_t adv_status = CS_ADVECTION_FIELD_GWF |
824 CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX;
825
826 mc->darcy->adv_field = cs_advection_field_add("darcy_field",
827 adv_status);
828
829 /* Add several properties:
830 * - the moisture content
831 * - the soil capacity
832 * - the (full) permeability = rel_permeability * abs_permeability
833 */
834
835 mc->moisture_content = cs_property_add("moisture_content", CS_PROPERTY_ISO);
836
837 mc->soil_capacity = cs_property_add("soil_capacity", CS_PROPERTY_ISO);
838
839 mc->permeability = cs_property_add("permeability", pty_type);
840
841 /* Add the time + diffusion terms to the Richards equation by associating the
842 full permeability to the diffusion property of the Richards eq. and the
843 soil capacity to the unsteady term */
844
845 cs_equation_param_t *eqp = cs_equation_get_param(mc->richards);
846
847 cs_equation_add_diffusion(eqp, mc->permeability);
848 cs_equation_add_time(eqp, mc->soil_capacity);
849
850 return mc;
851 }
852
853 /*----------------------------------------------------------------------------*/
854 /*!
855 * \brief Free the modelling context for the model of unsaturated single-phase
856 * flows
857 *
858 * \param[in, out] p_mc pointer of pointer to the model context structure
859 */
860 /*----------------------------------------------------------------------------*/
861
862 static void
_uspf_free_context(cs_gwf_unsaturated_single_phase_t ** p_mc)863 _uspf_free_context(cs_gwf_unsaturated_single_phase_t **p_mc)
864 {
865 if (p_mc == NULL)
866 return;
867 if (*p_mc == NULL)
868 return;
869
870 cs_gwf_unsaturated_single_phase_t *mc = *p_mc;
871
872 cs_gwf_darcy_flux_free(&(mc->darcy));
873
874 BFT_FREE(mc);
875 *p_mc = NULL;
876 }
877
878 /*----------------------------------------------------------------------------*/
879 /*!
880 * \brief Log the setup related to the modelling context of unsaturated
881 * single-phase flows
882 *
883 * \param[in] mc pointer to the model context structure
884 */
885 /*----------------------------------------------------------------------------*/
886
887 static void
_uspf_log_context(cs_gwf_unsaturated_single_phase_t * mc)888 _uspf_log_context(cs_gwf_unsaturated_single_phase_t *mc)
889 {
890 if (mc == NULL)
891 return;
892
893 cs_gwf_darcy_flux_log(mc->darcy);
894 }
895
896 /*----------------------------------------------------------------------------*/
897 /*!
898 * \brief Setup related to the modelling context of unsaturated single-phase
899 * flows. At this stage, all soils have been defined and equation
900 * parameters are set.
901 *
902 * \param[in, out] mc pointer to the model context structure
903 */
904 /*----------------------------------------------------------------------------*/
905
906 static void
_uspf_init_setup(cs_gwf_unsaturated_single_phase_t * mc)907 _uspf_init_setup(cs_gwf_unsaturated_single_phase_t *mc)
908 {
909 cs_gwf_t *gw = cs_gwf_main_structure;
910
911 if (mc == NULL)
912 return;
913 if (mc->richards == NULL)
914 bft_error(__FILE__, __LINE__, 0,
915 "%s: The Richards equation is not defined. Stop execution.\n",
916 __func__);
917 assert(cs_equation_is_steady(mc->richards) == false);
918
919 const int field_mask = CS_FIELD_INTENSIVE | CS_FIELD_VARIABLE | CS_FIELD_CDO;
920 const int c_loc_id = cs_mesh_location_get_id_by_name("cells");
921 const int v_loc_id = cs_mesh_location_get_id_by_name("vertices");
922 const int log_key = cs_field_key_id("log");
923 const int post_key = cs_field_key_id("post_vis");
924 const cs_param_space_scheme_t space_scheme =
925 cs_equation_get_space_scheme(mc->richards);
926
927 /* Handle gravity effects */
928 if (gw->flag & CS_GWF_GRAVITATION) {
929
930 switch (space_scheme) {
931 case CS_SPACE_SCHEME_CDOVB:
932 case CS_SPACE_SCHEME_CDOVCB:
933 mc->pressure_head = cs_field_create("pressure_head",
934 field_mask,
935 v_loc_id,
936 1,
937 true); /* has_previous */
938 break;
939
940 case CS_SPACE_SCHEME_CDOFB:
941 case CS_SPACE_SCHEME_HHO_P0:
942 mc->pressure_head = cs_field_create("pressure_head",
943 field_mask,
944 c_loc_id,
945 1,
946 true); /* has_previous */
947 break;
948
949 default:
950 bft_error(__FILE__, __LINE__, 0, " %s: Invalid space scheme.", __func__);
951 }
952
953 cs_field_set_key_int(mc->pressure_head, log_key, 1);
954 cs_field_set_key_int(mc->pressure_head, post_key, 1);
955
956 } /* Gravitation effect is activated */
957
958 /* Create fields at cells for properties */
959
960 int pty_mask = CS_FIELD_INTENSIVE | CS_FIELD_PROPERTY | CS_FIELD_CDO;
961 mc->moisture_field = cs_field_create("liquid_saturation",
962 pty_mask,
963 c_loc_id,
964 1, /* dimension */
965 true); /* has_previous */
966
967 cs_field_set_key_int(mc->moisture_field, log_key, 1);
968 if (gw->post_flag & CS_GWF_POST_LIQUID_SATURATION)
969 cs_field_set_key_int(mc->moisture_field, post_key, 1);
970
971 int permeability_dim = cs_property_get_dim(gw->abs_permeability);
972
973 mc->permeability_field = cs_field_create("permeability",
974 pty_mask,
975 c_loc_id,
976 permeability_dim,
977 true); /* has_previous */
978
979 if (gw->post_flag & CS_GWF_POST_PERMEABILITY) {
980 cs_field_set_key_int(mc->permeability_field, log_key, 1);
981 cs_field_set_key_int(mc->permeability_field, post_key, 1);
982 }
983
984 mc->capacity_field = cs_field_create("soil_capacity",
985 pty_mask,
986 c_loc_id,
987 1, /* dimension */
988 true);
989
990 cs_field_set_key_int(mc->capacity_field, log_key, 1);
991 if (gw->post_flag & CS_GWF_POST_SOIL_CAPACITY)
992 cs_field_set_key_int(mc->capacity_field, post_key, 1);
993
994 /* Add default post-processing related to groundwater flow module */
995 cs_post_add_time_mesh_dep_output(cs_gwf_extra_post_uspf, gw);
996 }
997
998 /*----------------------------------------------------------------------------*/
999 /*!
1000 * \brief Perform the last setup step in the case of a unsaturated single-phase
1001 * flow model in porous media
1002 *
1003 * \param[in] connect pointer to a cs_cdo_connect_t structure
1004 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1005 * \param[in, out] mc pointer to the casted model context
1006 */
1007 /*----------------------------------------------------------------------------*/
1008
1009 static void
_uspf_finalize_setup(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_gwf_unsaturated_single_phase_t * mc)1010 _uspf_finalize_setup(const cs_cdo_connect_t *connect,
1011 const cs_cdo_quantities_t *quant,
1012 cs_gwf_unsaturated_single_phase_t *mc)
1013 {
1014 cs_gwf_t *gw = cs_gwf_main_structure;
1015 assert(gw != NULL && mc != NULL);
1016
1017 const cs_field_t *hydraulic_head = cs_equation_get_field(mc->richards);
1018 const cs_param_space_scheme_t richards_scheme =
1019 cs_equation_get_space_scheme(mc->richards);
1020 const cs_lnum_t n_cells = connect->n_cells;
1021
1022 /* Set the Darcian flux (in the volume and at the boundary) */
1023
1024 cs_gwf_darcy_flux_define(connect, quant, richards_scheme, mc->darcy);
1025
1026 /* Allocate a head array defined at cells and used to update the soil
1027 properties */
1028
1029 switch (richards_scheme) {
1030
1031 case CS_SPACE_SCHEME_CDOVB:
1032 case CS_SPACE_SCHEME_CDOVCB:
1033 BFT_MALLOC(mc->head_in_law, n_cells, cs_real_t);
1034 #if defined(DEBUG) && !defined(NDEBUG)
1035 memset(mc->head_in_law, 0, sizeof(cs_real_t)*n_cells);
1036 #endif
1037 break;
1038
1039 case CS_SPACE_SCHEME_CDOFB:
1040
1041 if (gw->flag & CS_GWF_GRAVITATION)
1042 mc->head_in_law = mc->pressure_head->val;
1043 else
1044 mc->head_in_law = hydraulic_head->val;
1045
1046 bft_error(__FILE__, __LINE__, 0,
1047 "%s: Fb space scheme not fully implemented.", __func__);
1048 break;
1049
1050 default:
1051 bft_error(__FILE__, __LINE__, 0, "%s: Invalid space scheme.", __func__);
1052 break;
1053
1054 } /* Switch on Richards scheme */
1055
1056 /* Define the permeability property using a field */
1057
1058 cs_property_def_by_field(mc->permeability, mc->permeability_field);
1059
1060 /* Define the moisture content (liquid saturation) property using a field */
1061
1062 cs_property_def_by_field(mc->moisture_content, mc->moisture_field);
1063
1064 /* Define the soil capacity property using a field */
1065
1066 cs_property_def_by_field(mc->soil_capacity, mc->capacity_field);
1067
1068 /* Set soil context with array */
1069
1070 cs_gwf_soil_uspf_set_arrays(mc->head_in_law,
1071 mc->permeability_field->val,
1072 mc->moisture_field->val,
1073 mc->capacity_field->val);
1074
1075 }
1076
1077 /*----------------------------------------------------------------------------*/
1078 /*!
1079 * \brief Perform the update step in the case of a unsaturated single-phase
1080 * flow model in porous media
1081 *
1082 * \param[in] mesh pointer to a cs_mesh_t structure
1083 * \param[in] connect pointer to a cs_cdo_connect_t structure
1084 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1085 * \param[in] ts pointer to a cs_time_step_t structure
1086 * \param[in] cur2prev true or false
1087 * \param[in, out] mc pointer to the casted model context
1088 *
1089 *\return the value at which one has to evaluate the properties
1090 */
1091 /*----------------------------------------------------------------------------*/
1092
1093 static cs_real_t
_uspf_updates(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts,bool cur2prev,cs_gwf_unsaturated_single_phase_t * mc)1094 _uspf_updates(const cs_mesh_t *mesh,
1095 const cs_cdo_connect_t *connect,
1096 const cs_cdo_quantities_t *quant,
1097 const cs_time_step_t *ts,
1098 bool cur2prev,
1099 cs_gwf_unsaturated_single_phase_t *mc)
1100 {
1101 cs_real_t time_eval = ts->t_cur;
1102 if (cur2prev)
1103 time_eval = _get_time_eval(ts, mc->richards);
1104
1105 /* Update head */
1106
1107 _spf_update_head(quant, connect, mc->richards,
1108 mc->pressure_head,
1109 mc->head_in_law,
1110 cur2prev);
1111
1112 /* Update the advection field related to the groundwater flow module */
1113
1114 cs_gwf_darcy_flux_update(time_eval, mc->richards, cur2prev, mc->darcy);
1115
1116 /* Update properties related to soils.
1117 * Handle the permeability, the moisture content and the soil capacity
1118 */
1119
1120 if (cur2prev) {
1121
1122 cs_field_current_to_previous(mc->permeability_field);
1123 cs_field_current_to_previous(mc->moisture_field);
1124 cs_field_current_to_previous(mc->capacity_field);
1125
1126 }
1127
1128 /* Update soil properties with the new head values */
1129
1130 cs_gwf_soil_update(time_eval, mesh, connect, quant);
1131
1132 return time_eval;
1133 }
1134
1135 /*----------------------------------------------------------------------------*/
1136 /*!
1137 * \brief Predefined extra-operations for the groundwater flow module in case
1138 * of unsaturated single phase flows in porous media
1139 *
1140 * \param[in] connect pointer to a cs_cdo_connect_t structure
1141 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1142 * \param[in, out] mc pointer to the casted model context
1143 */
1144 /*----------------------------------------------------------------------------*/
1145
1146 static void
_uspf_extra_op(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_gwf_unsaturated_single_phase_t * mc)1147 _uspf_extra_op(const cs_cdo_connect_t *connect,
1148 const cs_cdo_quantities_t *cdoq,
1149 cs_gwf_unsaturated_single_phase_t *mc)
1150 {
1151 cs_gwf_t *gw = cs_gwf_main_structure;
1152 assert(gw != NULL && mc != NULL);
1153
1154 if (cs_flag_test(gw->post_flag, CS_GWF_POST_DARCY_FLUX_BALANCE) == false)
1155 return; /* Nothing to do */
1156
1157 cs_gwf_darcy_flux_balance(connect, cdoq,
1158 cs_equation_get_param(mc->richards),
1159 mc->darcy);
1160 }
1161
1162 /* ==========================================================================
1163 * Functions related to the model of two phase flows (TPF)
1164 * ========================================================================== */
1165
1166 /*----------------------------------------------------------------------------*/
1167 /*!
1168 * \brief Allocate and initialize the modelling context for the model of
1169 * (unsaturated) miscible two-phase flows
1170 *
1171 * \param[in] pty_type type of permeability (iso, ortho...)
1172 *
1173 * \return a pointer to a new allocated cs_gwf_miscible_two_phase_t structure
1174 */
1175 /*----------------------------------------------------------------------------*/
1176
1177 static cs_gwf_miscible_two_phase_t *
_mtpf_init_context(cs_property_type_t pty_type)1178 _mtpf_init_context(cs_property_type_t pty_type)
1179 {
1180 cs_gwf_miscible_two_phase_t *mc = NULL;
1181
1182 BFT_MALLOC(mc, 1, cs_gwf_miscible_two_phase_t);
1183
1184 /* Arrays */
1185
1186 mc->time_w_eq_array = NULL;
1187 mc->diff_w_eq_array = NULL;
1188 mc->time_h_eq_array = NULL;
1189 mc->diff_hl_eq_array = NULL;
1190 mc->diff_hg_eq_array = NULL;
1191 mc->l_rel_permeability = NULL;
1192
1193 /* Darcy flux (not used up to now) */
1194
1195 mc->l_darcy = NULL;
1196 mc->g_darcy = NULL;
1197
1198 /* Parameters (default values) */
1199
1200 mc->l_mass_density = 1000;
1201 mc->l_viscosity = 1e-3;
1202 mc->g_viscosity = 2e-5;
1203 mc->l_diffusivity_h = 1e-10;
1204 mc->w_molar_mass = 18e-3;
1205 mc->h_molar_mass = 3e-3;
1206 mc->ref_temperature = 280;
1207 mc->henry_constant = 1e-20; /* immiscible case */
1208
1209 /* Create a new equation for the water conservation */
1210
1211 mc->w_eq = cs_equation_add("w_conservation", /* equation name */
1212 "capillarity_pressure", /* variable name */
1213 CS_EQUATION_TYPE_GROUNDWATER,
1214 1,
1215 CS_PARAM_BC_HMG_NEUMANN);
1216
1217 /* Create a new equation for the hydrogen conservation */
1218
1219 mc->h_eq = cs_equation_add("h_conservation", /* equation name */
1220 "gas_pressure", /* variable name */
1221 CS_EQUATION_TYPE_GROUNDWATER,
1222 1,
1223 CS_PARAM_BC_HMG_NEUMANN);
1224
1225 /* Add properties:
1226 * - unsteady term for w_eq
1227 * - diffusion term for w_eq
1228 * - unsteady term for h_eq
1229 * - diffusion term for h_eq in the gas phase
1230 * - diffusion term for h_eq in the liquid phase (cross-term)
1231 */
1232
1233 mc->time_w_eq_pty = cs_property_add("time_w_eq_pty", CS_PROPERTY_ISO);
1234 mc->diff_w_eq_pty = cs_property_add("diff_w_eq_pty", pty_type);
1235 mc->time_h_eq_pty = cs_property_add("time_h_eq_pty", CS_PROPERTY_ISO);
1236 mc->diff_hg_eq_pty = cs_property_add("diff_hg_eq_pty", pty_type);
1237 mc->diff_hl_eq_pty = cs_property_add("diff_hl_eq_pty", pty_type);
1238
1239 /* Associate properties with equation to define new terms in these
1240 equations */
1241
1242 cs_equation_param_t *w_eqp = cs_equation_get_param(mc->w_eq);
1243
1244 cs_equation_add_time(w_eqp, mc->time_w_eq_pty);
1245 cs_equation_add_diffusion(w_eqp, mc->diff_w_eq_pty);
1246
1247 cs_equation_param_t *h_eqp = cs_equation_get_param(mc->h_eq);
1248
1249 cs_equation_add_time(h_eqp, mc->time_h_eq_pty);
1250 cs_equation_add_diffusion(h_eqp, mc->diff_hg_eq_pty);
1251
1252 return mc;
1253 }
1254
1255 /*----------------------------------------------------------------------------*/
1256 /*!
1257 * \brief Free the modelling context for the model of miscible two-phase
1258 * flows
1259 *
1260 * \param[in, out] p_mc pointer of pointer to the model context structure
1261 */
1262 /*----------------------------------------------------------------------------*/
1263
1264 static void
_mtpf_free_context(cs_gwf_miscible_two_phase_t ** p_mc)1265 _mtpf_free_context(cs_gwf_miscible_two_phase_t **p_mc)
1266 {
1267 if (p_mc == NULL)
1268 return;
1269 if (*p_mc == NULL)
1270 return;
1271
1272 cs_gwf_miscible_two_phase_t *mc = *p_mc;
1273
1274 cs_gwf_darcy_flux_free(&(mc->l_darcy));
1275 cs_gwf_darcy_flux_free(&(mc->g_darcy));
1276
1277 BFT_FREE(mc->time_w_eq_array);
1278 BFT_FREE(mc->diff_w_eq_array);
1279 BFT_FREE(mc->time_h_eq_array);
1280 BFT_FREE(mc->diff_hl_eq_array);
1281 BFT_FREE(mc->diff_hg_eq_array);
1282 BFT_FREE(mc->l_rel_permeability);
1283
1284 BFT_FREE(mc);
1285 *p_mc = NULL;
1286 }
1287
1288 /*----------------------------------------------------------------------------*/
1289 /*!
1290 * \brief Log the setup related to the modelling context of miscible
1291 * two-phase flows
1292 *
1293 * \param[in] mc pointer to the model context structure
1294 */
1295 /*----------------------------------------------------------------------------*/
1296
1297 static void
_mtpf_log_context(cs_gwf_miscible_two_phase_t * mc)1298 _mtpf_log_context(cs_gwf_miscible_two_phase_t *mc)
1299 {
1300 if (mc == NULL)
1301 return;
1302
1303 cs_gwf_darcy_flux_log(mc->l_darcy);
1304 cs_gwf_darcy_flux_log(mc->g_darcy);
1305
1306 cs_log_printf(CS_LOG_SETUP,
1307 " * GWF | Water mass density: %5.3e, viscosity: %5.3e,"
1308 " molar mass: %5.3e\n", mc->l_mass_density, mc->l_viscosity,
1309 mc->w_molar_mass);
1310 cs_log_printf(CS_LOG_SETUP,
1311 " * GWF | Main gas component: viscosity: %5.3e, diffusivity"
1312 " in the liquid phase: %5.3e, molar mass: %5.3e\n",
1313 mc->g_viscosity, mc->l_diffusivity_h, mc->h_molar_mass);
1314 cs_log_printf(CS_LOG_SETUP,
1315 " * GWF | Reference temperature: %5.3e K\n",
1316 mc->ref_temperature);
1317 cs_log_printf(CS_LOG_SETUP,
1318 " * GWF | Henry constant: %5.3e\n",
1319 mc->henry_constant);
1320
1321 }
1322
1323 /*----------------------------------------------------------------------------*/
1324 /*!
1325 * \brief Perform the initial setup step in the case of a miscible two-phase
1326 * flows model in porous media
1327 *
1328 * \param[in, out] mc pointer to the casted model context
1329 */
1330 /*----------------------------------------------------------------------------*/
1331
1332 static void
_mtpf_init_setup(cs_gwf_miscible_two_phase_t * mc)1333 _mtpf_init_setup(cs_gwf_miscible_two_phase_t *mc)
1334 {
1335 cs_gwf_t *gw = cs_gwf_main_structure;
1336 assert(gw != NULL && mc != NULL);
1337
1338 if (mc->w_eq == NULL || mc->h_eq == NULL)
1339 bft_error(__FILE__, __LINE__, 0,
1340 "%s: Equations are not defined for this model. Stop execution.\n",
1341 __func__);
1342
1343 assert(cs_equation_is_steady(mc->w_eq) == false);
1344 assert(cs_equation_is_steady(mc->h_eq) == false);
1345
1346 const int field_mask = CS_FIELD_INTENSIVE | CS_FIELD_VARIABLE | CS_FIELD_CDO;
1347 const int c_loc_id = cs_mesh_location_get_id_by_name("cells");
1348 const int v_loc_id = cs_mesh_location_get_id_by_name("vertices");
1349 const int log_key = cs_field_key_id("log");
1350 const int post_key = cs_field_key_id("post_vis");
1351 const cs_param_space_scheme_t space_scheme =
1352 cs_equation_get_space_scheme(mc->w_eq);
1353 assert(space_scheme == cs_equation_get_space_scheme(mc->h_eq));
1354
1355 /* One has to be consistent with the location of DoFs for the w_eq and h_eq
1356 * which are respectively related to the c_pressure and g_pressure */
1357
1358 int loc_id = c_loc_id;
1359 if (space_scheme == CS_SPACE_SCHEME_CDOVB ||
1360 space_scheme == CS_SPACE_SCHEME_CDOVCB)
1361 loc_id = v_loc_id;
1362
1363 mc->l_pressure = cs_field_create("liquid_pressure",
1364 field_mask,
1365 loc_id,
1366 1,
1367 true); /* has_previous */
1368
1369 cs_field_set_key_int(mc->l_pressure, log_key, 1);
1370 cs_field_set_key_int(mc->l_pressure, post_key, 1);
1371
1372 /* Create a liquid saturation field attached to cells: S_l */
1373
1374 int pty_mask = CS_FIELD_INTENSIVE | CS_FIELD_PROPERTY | CS_FIELD_CDO;
1375
1376 mc->l_saturation = cs_field_create("liquid_saturation",
1377 pty_mask,
1378 c_loc_id,
1379 1, /* dimension */
1380 true); /* has_previous */
1381
1382 cs_field_set_key_int(mc->l_saturation, log_key, 1);
1383 if (gw->post_flag & CS_GWF_POST_LIQUID_SATURATION)
1384 cs_field_set_key_int(mc->l_saturation, post_key, 1);
1385
1386 /* Add default post-processing related to groundwater flow module */
1387
1388 cs_post_add_time_mesh_dep_output(cs_gwf_extra_post_mtpf, gw);
1389 }
1390
1391 /*----------------------------------------------------------------------------*/
1392 /*!
1393 * \brief Perform the last setup step in the case of a (miscible) two-phase
1394 * flow model in porous media
1395 *
1396 * \param[in] connect pointer to a cs_cdo_connect_t structure
1397 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1398 * \param[in, out] mc pointer to the casted model context
1399 */
1400 /*----------------------------------------------------------------------------*/
1401
1402 static void
_mtpf_finalize_setup(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_gwf_miscible_two_phase_t * mc)1403 _mtpf_finalize_setup(const cs_cdo_connect_t *connect,
1404 const cs_cdo_quantities_t *quant,
1405 cs_gwf_miscible_two_phase_t *mc)
1406 {
1407 CS_UNUSED(quant);
1408
1409 cs_gwf_t *gw = cs_gwf_main_structure;
1410 assert(gw != NULL && mc != NULL);
1411
1412 const cs_lnum_t n_cells = connect->n_cells;
1413 const size_t csize = n_cells*sizeof(cs_real_t);
1414
1415 /* Set the Darcian flux (in the volume and at the boundary) -- Not done up to
1416 now */
1417
1418 /* Allocate and initialize the relative permeability in the liquid phase */
1419
1420 BFT_MALLOC(mc->l_rel_permeability, n_cells, cs_real_t);
1421 #pragma omp parallel for if (n_cells > CS_THR_MIN)
1422 for (cs_lnum_t i = 0; i < n_cells; i++)
1423 mc->l_rel_permeability[i] = 1; /* saturated by default */
1424
1425 /* Define the array storing the time property for the water eq. */
1426
1427 BFT_MALLOC(mc->time_w_eq_array, n_cells, cs_real_t);
1428 memset(mc->time_w_eq_array, 0, csize);
1429
1430 cs_property_def_by_array(mc->time_w_eq_pty,
1431 cs_flag_primal_cell, /* where data are located */
1432 mc->time_w_eq_array,
1433 false, /* not owner of the array */
1434 NULL); /* no index */
1435
1436 /* Define the array storing the diffusion property for the water eq. */
1437
1438 BFT_MALLOC(mc->diff_w_eq_array, n_cells, cs_real_t);
1439 memset(mc->diff_w_eq_array, 0, csize);
1440
1441 cs_property_def_by_array(mc->diff_w_eq_pty,
1442 cs_flag_primal_cell, /* where data are located */
1443 mc->diff_w_eq_array,
1444 false, /* not owner of the array */
1445 NULL); /* no index */
1446
1447 /* Define the array storing the time property for the hydrogen eq. */
1448
1449 BFT_MALLOC(mc->time_h_eq_array, n_cells, cs_real_t);
1450 memset(mc->time_h_eq_array, 0, csize);
1451
1452 cs_property_def_by_array(mc->time_h_eq_pty,
1453 cs_flag_primal_cell, /* where data are located */
1454 mc->time_h_eq_array,
1455 false, /* not owner of the array */
1456 NULL); /* no index */
1457
1458 /* Define the array storing the diffusion property for the liquid part in the
1459 hydrogen eq. */
1460
1461 BFT_MALLOC(mc->diff_hl_eq_array, n_cells, cs_real_t);
1462 memset(mc->diff_hl_eq_array, 0, csize);
1463
1464 cs_property_def_by_array(mc->diff_hl_eq_pty,
1465 cs_flag_primal_cell, /* where data are located */
1466 mc->diff_hl_eq_array,
1467 false, /* not owner of the array */
1468 NULL); /* no index */
1469
1470 /* Define the array storing the diffusion property for the gaz part in the
1471 hydrogen eq. */
1472
1473 BFT_MALLOC(mc->diff_hg_eq_array, n_cells, cs_real_t);
1474 memset(mc->diff_hg_eq_array, 0, csize);
1475
1476 cs_property_def_by_array(mc->diff_hg_eq_pty,
1477 cs_flag_primal_cell, /* where data are located */
1478 mc->diff_hg_eq_array,
1479 false, /* not owner of the array */
1480 NULL); /* no index */
1481
1482 /* Set soil context with array */
1483
1484 /* TODO */
1485 }
1486
1487 /*----------------------------------------------------------------------------*/
1488 /*!
1489 * \brief Update the pressure in the liquid phase
1490 * Case of two-phase flows in porous media.
1491 *
1492 * \param[in, out] mc pointer to a modelling context structure
1493 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1494 * \param[in] connect pointer to a cs_cdo_connect_t structure
1495 * \param[in] cur2prev true or false
1496 */
1497 /*----------------------------------------------------------------------------*/
1498
1499 /* static void */
1500 /* _tpf_water_pressure(cs_gwf_miscible_two_phase_t *mc, */
1501 /* const cs_cdo_quantities_t *cdoq, */
1502 /* const cs_cdo_connect_t *connect, */
1503 /* bool cur2prev) */
1504 /* { */
1505 /* cs_gwf_t *gw = cs_gwf_main_structure; */
1506 /* assert(mc != NULL && gw != NULL); */
1507 /* const cs_equation_t *w_eq = mc->w_eq; */
1508 /* const cs_equation_t *h_eq = mc->h_eq; */
1509
1510 /* cs_field_t *capillary_pressure = cs_equation_get_field(w_eq); */
1511 /* cs_field_t *gcomp_pressure = cs_equation_get_field(h_eq); */
1512 /* cs_field_t *water_pressure = mc->water_pressure; */
1513
1514 /* cs_param_space_scheme_t r_scheme = cs_equation_get_space_scheme(w_eq); */
1515
1516 /* switch (r_scheme) { */
1517
1518 /* case CS_SPACE_SCHEME_CDOVB: */
1519 /* { */
1520 /* assert(capillary_pressure->location_id == */
1521 /* cs_mesh_location_get_id_by_name("vertices")); */
1522 /* assert(gcomp_pressure->location_id == */
1523 /* cs_mesh_location_get_id_by_name("vertices")); */
1524 /* assert(water_pressure->location_id == */
1525 /* cs_mesh_location_get_id_by_name("vertices")); */
1526
1527 /* # pragma omp parallel for if (cdoq->n_vertices > CS_THR_MIN) */
1528 /* for (cs_lnum_t i = 0; i < cdoq->n_vertices; i++) */
1529 /* water_pressure->val[i] = gcomp_pressure->val[i] - capillary_pressure->val[i]; */
1530
1531 /* } */
1532 /* break; */
1533
1534 /* case CS_SPACE_SCHEME_CDOFB: */
1535 /* default: */
1536 /* break; /\* Nothing to do *\/ */
1537 /* } */
1538
1539 /* } */
1540
1541 /*----------------------------------------------------------------------------*/
1542 /*!
1543 * \brief Perform the update step in the case of a miscible two-phase flow
1544 * model in porous media
1545 *
1546 * \param[in] mesh pointer to a cs_mesh_t structure
1547 * \param[in] connect pointer to a cs_cdo_connect_t structure
1548 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1549 * \param[in] ts pointer to a cs_time_step_t structure
1550 * \param[in] cur2prev true or false
1551 * \param[in, out] mc pointer to the casted model context
1552 *
1553 *\return the value at which one has to evaluate the properties
1554 */
1555 /*----------------------------------------------------------------------------*/
1556
1557 static cs_real_t
_mtpf_updates(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts,bool cur2prev,cs_gwf_miscible_two_phase_t * mc)1558 _mtpf_updates(const cs_mesh_t *mesh,
1559 const cs_cdo_connect_t *connect,
1560 const cs_cdo_quantities_t *quant,
1561 const cs_time_step_t *ts,
1562 bool cur2prev,
1563 cs_gwf_miscible_two_phase_t *mc)
1564 {
1565 CS_UNUSED(mesh);
1566 CS_UNUSED(connect);
1567 CS_UNUSED(quant);
1568
1569 cs_gwf_t *gw = cs_gwf_main_structure;
1570 assert(gw != NULL && mc != NULL);
1571
1572 cs_real_t time_eval = ts->t_cur;
1573 if (cur2prev)
1574 time_eval = _get_time_eval(ts, mc->w_eq);
1575
1576 /* /\* Update head *\/ */
1577 /* _tpf_water_pressure(mc, quant, connect, cur2prev); */
1578
1579 /* /\* Update properties related to soils (water_saturation, water/gaz permeabilities, */
1580 /* compressibility coefficients *\/ */
1581
1582 /* if (cur2prev) { */
1583 /* cs_field_current_to_previous(gw->permea_field); */
1584 /* cs_field_current_to_previous(mc->gcomp_permea_field); */
1585 /* cs_field_current_to_previous(mc->water_pressure); */
1586 /* cs_field_current_to_previous(mc->l_saturation); */
1587 /* } */
1588
1589 /* const cs_equation_t *w_eq = mc->w_eq; */
1590 /* const cs_field_t *capillary_pressure = cs_equation_get_field(w_eq); */
1591
1592 /* const int n_soils = cs_gwf_get_n_soils(); */
1593 /* for (int i = 0; i < n_soils; i++) { */
1594
1595 /* cs_gwf_soil_t *soil = cs_gwf_soil_by_id(i); */
1596 /* const cs_zone_t *zone = cs_volume_zone_by_id(soil->zone_id); */
1597
1598 /* // Update soil properties in soil struct TO BE EDITED */
1599 /* soil->update_properties(time_eval, mesh, connect, quant, */
1600 /* capillary_pressure->val, */
1601 /* zone, */
1602 /* soil->input); */
1603
1604 /* } /\* Loop on soils *\/ */
1605
1606 /* /\* #if defined(DEBUG) && !defined(NDEBUG) && CS_GWF_DBG > 1 *\/ */
1607 /* /\* cs_dbg_darray_to_listing("MOISTURE_CONTENT", *\/ */
1608 /* /\* quant->n_cells, *\/ */
1609 /* /\* mc->moisture_field->val, 8); *\/ */
1610 /* /\* #endif *\/ */
1611
1612 return time_eval;
1613 }
1614
1615 /*----------------------------------------------------------------------------*/
1616 /*!
1617 * \brief Predefined extra-operations for the groundwater flow module in case
1618 * of (unsaturated) miscible two-phase flows in porous media
1619 *
1620 * \param[in] connect pointer to a cs_cdo_connect_t structure
1621 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1622 * \param[in, out] mc pointer to the casted model context
1623 */
1624 /*----------------------------------------------------------------------------*/
1625
1626 static void
_mtpf_extra_op(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_gwf_unsaturated_single_phase_t * mc)1627 _mtpf_extra_op(const cs_cdo_connect_t *connect,
1628 const cs_cdo_quantities_t *cdoq,
1629 cs_gwf_unsaturated_single_phase_t *mc)
1630 {
1631 CS_UNUSED(connect);
1632 CS_UNUSED(cdoq);
1633
1634 cs_gwf_t *gw = cs_gwf_main_structure;
1635 assert(gw != NULL && mc != NULL);
1636
1637 if (cs_flag_test(gw->post_flag, CS_GWF_POST_DARCY_FLUX_BALANCE) == false)
1638 return; /* Nothing to do */
1639
1640 /* TO BE DONE (not useful up to now) */
1641 }
1642
1643 /*----------------------------------------------------------------------------*/
1644 /*!
1645 * \brief Compute the new (unsteady) state for the groundwater flows module.
1646 * Case of (miscible) two-phase flows in porous media.
1647 *
1648 * \param[in] mesh pointer to a cs_mesh_t structure
1649 * \param[in] time_step pointer to a cs_time_step_t structure
1650 * \param[in] connect pointer to a cs_cdo_connect_t structure
1651 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1652 * \param[in, out] mc pointer to the casted model context
1653 */
1654 /*----------------------------------------------------------------------------*/
1655
1656 static void
_mtpf_compute(const cs_mesh_t * mesh,const cs_time_step_t * time_step,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_gwf_miscible_two_phase_t * mc)1657 _mtpf_compute(const cs_mesh_t *mesh,
1658 const cs_time_step_t *time_step,
1659 const cs_cdo_connect_t *connect,
1660 const cs_cdo_quantities_t *cdoq,
1661 cs_gwf_miscible_two_phase_t *mc)
1662 {
1663 cs_gwf_t *gw = cs_gwf_main_structure;
1664
1665 assert(gw != NULL && mc != NULL);
1666 assert(gw->model == CS_GWF_MODEL_TWO_PHASE);
1667
1668 cs_equation_t *w_eq = mc->w_eq;
1669 cs_equation_t *h_eq = mc->h_eq;
1670
1671 /* Sanity checks */
1672
1673 assert(w_eq != NULL && h_eq != NULL);
1674 assert(cs_equation_get_type(w_eq) == CS_EQUATION_TYPE_GROUNDWATER);
1675 assert(cs_equation_get_type(h_eq) == CS_EQUATION_TYPE_GROUNDWATER);
1676
1677 /* TODO */
1678
1679 bool cur2prev = true;
1680
1681 /* Build and solve the linear system related to the Richards equations */
1682
1683 if (!cs_equation_is_steady(w_eq) || !cs_equation_is_steady(h_eq)) {
1684
1685 /* Solve the algebraic system. By default, a current to previous operation
1686 is performed */
1687
1688 /* (@_@) Not clear how to do this: to check with J. Bonelle for this part */
1689
1690 cs_equation_solve(cur2prev, mesh, w_eq);
1691 cs_equation_solve(cur2prev, mesh, h_eq);
1692
1693 /* Update the variables related to the groundwater flow system */
1694
1695 cs_gwf_update(mesh, connect, cdoq, time_step, cur2prev);
1696
1697 }
1698 }
1699
1700 /*----------------------------------------------------------------------------*/
1701 /*!
1702 * \brief Create a structure dedicated to manage groundwater flows
1703 *
1704 * \return a pointer to a new allocated cs_gwf_t structure
1705 */
1706 /*----------------------------------------------------------------------------*/
1707
1708 static cs_gwf_t *
_gwf_create(void)1709 _gwf_create(void)
1710 {
1711 cs_gwf_t *gw = NULL;
1712
1713 BFT_MALLOC(gw, 1, cs_gwf_t);
1714
1715 /* Default initialization */
1716
1717 gw->model = CS_GWF_N_MODEL_TYPES;
1718 gw->flag = 0;
1719 gw->post_flag = CS_GWF_POST_DARCY_FLUX_BALANCE;
1720
1721 /* Property common to all model */
1722
1723 gw->abs_permeability = NULL;
1724
1725 /* Modelling context */
1726
1727 gw->model_context = NULL;
1728
1729 /* Tracer part */
1730
1731 gw->n_tracers = 0;
1732 gw->tracers = NULL;
1733 gw->finalize_tracer_setup = NULL;
1734 gw->add_tracer_terms = NULL;
1735
1736 return gw;
1737 }
1738
1739 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1740
1741 /*============================================================================
1742 * Public function prototypes
1743 *============================================================================*/
1744
1745 /*----------------------------------------------------------------------------*/
1746 /*!
1747 * \brief Check if the groundwater flow module has been activated
1748 *
1749 * \return true or false
1750 */
1751 /*----------------------------------------------------------------------------*/
1752
1753 bool
cs_gwf_is_activated(void)1754 cs_gwf_is_activated(void)
1755 {
1756 if (cs_gwf_main_structure == NULL)
1757 return false;
1758 else
1759 return true;
1760 }
1761
1762 /*----------------------------------------------------------------------------*/
1763 /*!
1764 * \brief Initialize the module dedicated to groundwater flows
1765 *
1766 * \param[in] pty_type type of permeability (iso, ortho...)
1767 * \param[in] model type of physical modelling
1768 * \param[in] option_flag optional flag to specify this module
1769 *
1770 * \return a pointer to a new allocated groundwater flow structure
1771 */
1772 /*----------------------------------------------------------------------------*/
1773
1774 cs_gwf_t *
cs_gwf_activate(cs_property_type_t pty_type,cs_gwf_model_type_t model,cs_gwf_option_flag_t option_flag)1775 cs_gwf_activate(cs_property_type_t pty_type,
1776 cs_gwf_model_type_t model,
1777 cs_gwf_option_flag_t option_flag)
1778 {
1779 cs_gwf_t *gw = _gwf_create();
1780
1781 /* Store the pointer to the groundawater flow structure */
1782
1783 cs_gwf_main_structure = gw;
1784
1785 gw->model = model;
1786 gw->flag = option_flag;
1787
1788 /* Add the absolute permeability property (used in the definition of the
1789 * diffusion term in the conservation equations and in the definition of the
1790 * Darcy flux. According to the type of soil model, the full permeability is
1791 * weigthed by a relative permeability (e.g. in a Van Genuchten-Mualen
1792 * model). */
1793
1794 gw->abs_permeability = cs_property_add("absolute_permeability", pty_type);
1795
1796 /* Allocate and initialize each model context (mc) */
1797
1798 switch (model) {
1799
1800 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
1801 gw->model_context = _sspf_init_context();
1802 break;
1803
1804 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
1805 gw->model_context = _uspf_init_context(pty_type);
1806 break;
1807
1808 case CS_GWF_MODEL_TWO_PHASE:
1809 gw->post_flag |= CS_GWF_POST_LIQUID_SATURATION;
1810 gw->model_context = _mtpf_init_context(pty_type);
1811 break;
1812
1813 default:
1814 bft_error(__FILE__, __LINE__, 0,
1815 " %s: Invalid model type for the GroundWater Flow module.\n",
1816 __func__);
1817
1818 }
1819
1820 return gw;
1821 }
1822
1823 /*----------------------------------------------------------------------------*/
1824 /*!
1825 * \brief Free all structures related to groundwater flows
1826 *
1827 * \return a NULL pointer
1828 */
1829 /*----------------------------------------------------------------------------*/
1830
1831 cs_gwf_t *
cs_gwf_destroy_all(void)1832 cs_gwf_destroy_all(void)
1833 {
1834 if (cs_gwf_main_structure == NULL)
1835 return NULL;
1836
1837 cs_gwf_t *gw = cs_gwf_main_structure;
1838
1839 switch (gw->model) {
1840
1841 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
1842 {
1843 cs_gwf_saturated_single_phase_t *mc = gw->model_context;
1844
1845 _sspf_free_context(&(mc));
1846 }
1847 break;
1848
1849 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
1850 {
1851 cs_gwf_unsaturated_single_phase_t *mc = gw->model_context;
1852
1853 _uspf_free_context(&(mc));
1854 }
1855 break;
1856
1857 case CS_GWF_MODEL_TWO_PHASE:
1858 {
1859 cs_gwf_miscible_two_phase_t *mc = gw->model_context;
1860
1861 _mtpf_free_context(&(mc));
1862 }
1863 break;
1864
1865 default:
1866 /* Nothing to do */
1867 break;
1868 }
1869
1870 /* Free all soils */
1871
1872 cs_gwf_soil_free_all();
1873
1874 /* Manage the tracer-related members */
1875
1876 for (int i = 0; i < gw->n_tracers; i++)
1877 gw->tracers[i] = cs_gwf_tracer_free(gw->tracers[i]);
1878 BFT_FREE(gw->tracers);
1879 BFT_FREE(gw->finalize_tracer_setup);
1880 BFT_FREE(gw->add_tracer_terms);
1881
1882 BFT_FREE(gw);
1883
1884 /* Fields, equations, advection fields and properties are freed elsewhere */
1885
1886 return NULL;
1887 }
1888
1889 /*----------------------------------------------------------------------------*/
1890 /*!
1891 * \brief Summary of the main cs_gwf_t structure
1892 */
1893 /*----------------------------------------------------------------------------*/
1894
1895 void
cs_gwf_log_setup(void)1896 cs_gwf_log_setup(void)
1897 {
1898 cs_gwf_t *gw = cs_gwf_main_structure;
1899
1900 if (gw == NULL)
1901 return;
1902 assert(gw->model < CS_GWF_N_MODEL_TYPES);
1903
1904 cs_log_printf(CS_LOG_SETUP, "\nSummary of the groundwater module\n");
1905 cs_log_printf(CS_LOG_SETUP, "%s", cs_sep_h1);
1906
1907 /* Display information on the general options */
1908
1909 if (gw->flag & CS_GWF_GRAVITATION)
1910 cs_log_printf(CS_LOG_SETUP, " * GWF | Gravitation: **True**\n");
1911 else
1912 cs_log_printf(CS_LOG_SETUP, " * GWF | Gravitation: **False**\n");
1913
1914 if (gw->flag & CS_GWF_ENFORCE_DIVERGENCE_FREE)
1915 cs_log_printf(CS_LOG_SETUP,
1916 " * GWF | Enforce the divergence-free constraint"
1917 " for the Darcy flux\n");
1918 if (gw->flag & CS_GWF_FORCE_RICHARDS_ITERATIONS)
1919 cs_log_printf(CS_LOG_SETUP,
1920 " * GWF | Force to solve Richards equation"
1921 " at each time step\n");
1922 if (gw->flag & CS_GWF_RESCALE_HEAD_TO_ZERO_MEAN_VALUE)
1923 cs_log_printf(CS_LOG_SETUP,
1924 " * GWF | Rescale head w.r.t zero mean value\n");
1925
1926 /* Display information on the post-processing options */
1927
1928 bool post_capacity =
1929 (gw->post_flag & CS_GWF_POST_SOIL_CAPACITY) ? true : false;
1930 bool post_liquid_saturation =
1931 (gw->post_flag & CS_GWF_POST_LIQUID_SATURATION) ? true : false;
1932 bool post_permeability =
1933 (gw->post_flag & CS_GWF_POST_PERMEABILITY) ? true : false;
1934
1935 cs_log_printf(CS_LOG_SETUP, " * GWF | Post:"
1936 " Soil capacity %s Liquid saturation %s Permeability %s\n",
1937 cs_base_strtf(post_capacity),
1938 cs_base_strtf(post_liquid_saturation),
1939 cs_base_strtf(post_permeability));
1940
1941 bool do_balance =
1942 (gw->post_flag & CS_GWF_POST_DARCY_FLUX_BALANCE) ? true : false;
1943 bool do_divergence =
1944 (gw->post_flag & CS_GWF_POST_DARCY_FLUX_DIVERGENCE) ? true : false;
1945 bool post_boundary =
1946 (gw->post_flag & CS_GWF_POST_DARCY_FLUX_AT_BOUNDARY) ? true : false;
1947
1948 cs_log_printf(CS_LOG_SETUP,
1949 " * GWF | Darcy Flux:"
1950 " Balance %s Divergence %s At boundary faces: %s\n",
1951 cs_base_strtf(do_balance),
1952 cs_base_strtf(do_divergence),
1953 cs_base_strtf(post_boundary));
1954
1955 /* Main options */
1956
1957 switch(gw->model) {
1958
1959 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
1960 cs_log_printf(CS_LOG_SETUP,
1961 " * GWF | Model: %s\n", cs_gwf_model_name[gw->model]);
1962 _sspf_log_context(gw->model_context);
1963 break;
1964
1965 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
1966 cs_log_printf(CS_LOG_SETUP,
1967 " * GWF | Model: %s\n", cs_gwf_model_name[gw->model]);
1968 _uspf_log_context(gw->model_context);
1969 break;
1970
1971 case CS_GWF_MODEL_TWO_PHASE:
1972 cs_log_printf(CS_LOG_SETUP,
1973 " * GWF | Model: %s\n", cs_gwf_model_name[gw->model]);
1974 _mtpf_log_context(gw->model_context);
1975 break;
1976
1977 default:
1978 break;
1979
1980 }
1981
1982 /* Soils */
1983
1984 cs_gwf_soil_log_setup();
1985
1986 /* Tracers */
1987
1988 cs_log_printf(CS_LOG_SETUP,
1989 " * GWF | Number of tracer equations: %d\n", gw->n_tracers);
1990
1991 for (int i = 0; i < gw->n_tracers; i++)
1992 cs_gwf_tracer_log_setup(gw->tracers[i]);
1993 }
1994
1995 /*----------------------------------------------------------------------------*/
1996 /*!
1997 * \brief Set the parameters defining the two-phase flow model.
1998 * Use SI unit if not prescribed otherwise.
1999 *
2000 * \param[in] l_mass_density mass density of the main liquid component
2001 * \param[in] l_viscosity viscosity in the liquid phase (Pa.s)
2002 * \param[in] g_viscosity viscosity in the gas phase (Pa.s)
2003 * \param[in] l_diffusivity_h diffusivity of the main gas component in the
2004 * liquid phase
2005 * \param[in] w_molar_mass molar mass of the main liquid component
2006 * \param[in] h_molar_mass molar mass of the main gas component
2007 * \param[in] ref_temperature reference temperature
2008 * \param[in] henry_constant constant in the Henry law
2009 */
2010 /*----------------------------------------------------------------------------*/
2011
2012 void
cs_gwf_set_two_phase_model(cs_real_t l_mass_density,cs_real_t l_viscosity,cs_real_t g_viscosity,cs_real_t l_diffusivity_h,cs_real_t w_molar_mass,cs_real_t h_molar_mass,cs_real_t ref_temperature,cs_real_t henry_constant)2013 cs_gwf_set_two_phase_model(cs_real_t l_mass_density,
2014 cs_real_t l_viscosity,
2015 cs_real_t g_viscosity,
2016 cs_real_t l_diffusivity_h,
2017 cs_real_t w_molar_mass,
2018 cs_real_t h_molar_mass,
2019 cs_real_t ref_temperature,
2020 cs_real_t henry_constant)
2021 {
2022 cs_gwf_t *gw = cs_gwf_main_structure;
2023
2024 /* Sanity checks */
2025 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2026 if (gw->model != CS_GWF_MODEL_TWO_PHASE)
2027 bft_error(__FILE__, __LINE__, 0,
2028 "%s: Invalid model. One expects a two-phase flow model.\n",
2029 __func__);
2030
2031 cs_gwf_miscible_two_phase_t *mc = gw->model_context;
2032
2033 /* Sanity checks */
2034
2035 assert(mc != NULL);
2036 assert(l_mass_density > 0);
2037 assert(ref_temperature > 0); /* In Kelvin */
2038 assert(w_molar_mass > 0 && h_molar_mass > 0);
2039 assert(l_viscosity > 0 && g_viscosity > 0);
2040
2041 /* Set the parameters */
2042
2043 mc->l_mass_density = l_mass_density;
2044 mc->l_viscosity = l_viscosity;
2045 mc->g_viscosity = g_viscosity;
2046 mc->l_diffusivity_h = l_diffusivity_h;
2047 mc->w_molar_mass = w_molar_mass;
2048 mc->h_molar_mass = h_molar_mass;
2049 mc->ref_temperature = ref_temperature;
2050 mc->henry_constant = henry_constant;
2051 }
2052
2053 /*----------------------------------------------------------------------------*/
2054 /*!
2055 * \brief Set the flag dedicated to the post-processing of the GWF module
2056 *
2057 * \param[in] post_flag flag to set
2058 * \param[in] reset reset post flag before
2059 */
2060 /*----------------------------------------------------------------------------*/
2061
2062 void
cs_gwf_set_post_options(cs_flag_t post_flag,bool reset)2063 cs_gwf_set_post_options(cs_flag_t post_flag,
2064 bool reset)
2065 {
2066 cs_gwf_t *gw = cs_gwf_main_structure;
2067 if (gw == NULL)
2068 return;
2069
2070 if (reset)
2071 gw->post_flag = post_flag;
2072 else
2073 gw->post_flag |= post_flag;
2074
2075 if (gw->post_flag & CS_GWF_POST_DARCY_FLUX_AT_BOUNDARY) {
2076
2077 cs_adv_field_t *adv = _get_l_adv_field(gw);
2078 if (adv != NULL)
2079 adv->status |= CS_ADVECTION_FIELD_DEFINE_AT_BOUNDARY_FACES;
2080
2081 if (gw->model == CS_GWF_MODEL_TWO_PHASE) {
2082
2083 cs_gwf_miscible_two_phase_t *mc =
2084 (cs_gwf_miscible_two_phase_t *)gw->model_context;
2085
2086 if (mc->g_darcy != NULL) {
2087 adv = mc->g_darcy->adv_field;
2088 if (adv != NULL)
2089 adv->status |= CS_ADVECTION_FIELD_DEFINE_AT_BOUNDARY_FACES;
2090 }
2091
2092 }
2093
2094 } /* Darcy flux at boundary */
2095 }
2096
2097 /*----------------------------------------------------------------------------*/
2098 /*!
2099 * \brief Retrieve the advection field related to the Darcy flux in the liquid
2100 * phase
2101 *
2102 * \return a pointer to a cs_adv_field_t structure or NULL
2103 */
2104 /*----------------------------------------------------------------------------*/
2105
2106 cs_adv_field_t *
cs_gwf_get_adv_field(void)2107 cs_gwf_get_adv_field(void)
2108 {
2109 cs_gwf_t *gw = cs_gwf_main_structure;
2110
2111 if (gw == NULL)
2112 return NULL;
2113
2114 return _get_l_adv_field(gw);
2115 }
2116
2117 /*----------------------------------------------------------------------------*/
2118 /*!
2119 * \brief Retrieve the head used in soil updates when an unsaturated
2120 * single-phase flow is considered. These values are located at cells.
2121 *
2122 * \return a pointer to the requested array of values or NULL if not defined
2123 */
2124 /*----------------------------------------------------------------------------*/
2125
2126 cs_real_t *
cs_gwf_get_uspf_head_in_law(void)2127 cs_gwf_get_uspf_head_in_law(void)
2128 {
2129 cs_gwf_t *gw = cs_gwf_main_structure;
2130
2131 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2132
2133 if (gw->model == CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE) {
2134
2135 cs_gwf_unsaturated_single_phase_t *mc = gw->model_context;
2136
2137 return mc->head_in_law;
2138 }
2139
2140 /* In all other cases */
2141
2142 return NULL;
2143 }
2144
2145 /*----------------------------------------------------------------------------*/
2146 /*!
2147 * \brief Create and add a new cs_gwf_soil_t structure. An initialization by
2148 * default of all members is performed.
2149 *
2150 * \param[in] z_name name of the volume zone corresponding to the soil
2151 * \param[in] bulk_density value of the mass density
2152 * \param[in] sat_moisture value of the saturated moisture content
2153 * \param[in] soil_model type of modelling for the hydraulic behavior
2154 *
2155 * \return a pointer to the new allocated soil structure
2156 */
2157 /*----------------------------------------------------------------------------*/
2158
2159 cs_gwf_soil_t *
cs_gwf_add_soil(const char * z_name,double bulk_density,double sat_moisture,cs_gwf_soil_hydraulic_model_t soil_model)2160 cs_gwf_add_soil(const char *z_name,
2161 double bulk_density,
2162 double sat_moisture,
2163 cs_gwf_soil_hydraulic_model_t soil_model)
2164 {
2165 cs_gwf_t *gw = cs_gwf_main_structure;
2166
2167 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2168
2169 const cs_zone_t *zone = cs_volume_zone_by_name_try(z_name);
2170
2171 /* Sanity checks */
2172 if (zone == NULL)
2173 bft_error(__FILE__, __LINE__, 0,
2174 " Zone %s related to the same soil is not defined.\n"
2175 " Stop adding a new soil.", z_name);
2176
2177 if (gw->model == CS_GWF_MODEL_SATURATED_SINGLE_PHASE)
2178 if (soil_model != CS_GWF_SOIL_SATURATED)
2179 bft_error(__FILE__, __LINE__, 0,
2180 "%s: Invalid type of soil with the general model used.\n"
2181 " In a saturated single-phase model, all soils have to be"
2182 " of type CS_GWF_SOIL_SATURATED.\n", __func__);
2183
2184 assert(bulk_density > 0);
2185 assert(sat_moisture > 0);
2186
2187 cs_property_type_t perm_type = cs_property_get_type(gw->abs_permeability);
2188
2189 cs_gwf_soil_t *soil = cs_gwf_soil_create(zone,
2190 soil_model,
2191 perm_type,
2192 sat_moisture,
2193 bulk_density);
2194
2195 return soil;
2196 }
2197
2198 /*----------------------------------------------------------------------------*/
2199 /*!
2200 * \brief Add a new equation related to the groundwater flow module
2201 * This equation is a particular type of unsteady advection-diffusion
2202 * reaction eq.
2203 * Tracer is advected thanks to the darcian velocity and
2204 * diffusion/reaction parameters result from a physical modelling.
2205 * Terms solved in the equation are activated according to the settings.
2206 * The advection field corresponds to that of the liquid phase.
2207 *
2208 * \param[in] tr_model physical modelling to consider (0 = default settings)
2209 * \param[in] eq_name name of the tracer equation
2210 * \param[in] var_name name of the related variable
2211 */
2212 /*----------------------------------------------------------------------------*/
2213
2214 cs_gwf_tracer_t *
cs_gwf_add_tracer(cs_gwf_tracer_model_t tr_model,const char * eq_name,const char * var_name)2215 cs_gwf_add_tracer(cs_gwf_tracer_model_t tr_model,
2216 const char *eq_name,
2217 const char *var_name)
2218 {
2219 cs_gwf_t *gw = cs_gwf_main_structure;
2220
2221 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2222
2223 if (tr_model & CS_GWF_TRACER_USER)
2224 bft_error(__FILE__, __LINE__, 0,
2225 "%s: User-defined is not allowed in this context.\n"
2226 " Please consider cs_gwf_add_user_tracer() instead.", __func__);
2227
2228 /* Set the advection field structure */
2229
2230 cs_adv_field_t *adv = _get_l_adv_field(gw);
2231 int tr_id = gw->n_tracers;
2232 cs_gwf_tracer_t *tracer = cs_gwf_tracer_init(tr_id,
2233 eq_name,
2234 var_name,
2235 adv,
2236 tr_model);
2237
2238 gw->n_tracers += 1;
2239 BFT_REALLOC(gw->tracers, gw->n_tracers, cs_gwf_tracer_t *);
2240 BFT_REALLOC(gw->finalize_tracer_setup,
2241 gw->n_tracers, cs_gwf_tracer_setup_t *);
2242 BFT_REALLOC(gw->add_tracer_terms,
2243 gw->n_tracers, cs_gwf_tracer_add_terms_t *);
2244
2245 gw->tracers[tr_id] = tracer;
2246
2247 if (gw->model == CS_GWF_MODEL_SATURATED_SINGLE_PHASE)
2248 gw->finalize_tracer_setup[tr_id] = cs_gwf_tracer_saturated_setup;
2249 else
2250 gw->finalize_tracer_setup[tr_id] = cs_gwf_tracer_unsaturated_setup;
2251
2252 gw->add_tracer_terms[tr_id] = cs_gwf_tracer_add_terms;
2253
2254 return tracer;
2255 }
2256
2257 /*----------------------------------------------------------------------------*/
2258 /*!
2259 * \brief Add a new equation related to the groundwater flow module
2260 * This equation is a particular type of unsteady advection-diffusion
2261 * reaction eq.
2262 * Tracer is advected thanks to the darcian velocity and
2263 * diffusion/reaction parameters result from a physical modelling.
2264 * Terms are activated according to the settings.
2265 * Modelling of the tracer parameters are left to the user
2266 *
2267 * \param[in] eq_name name of the tracer equation
2268 * \param[in] var_name name of the related variable
2269 * \param[in] setup function pointer (predefined prototype)
2270 * \param[in] add_terms function pointer (predefined prototype)
2271 */
2272 /*----------------------------------------------------------------------------*/
2273
2274 cs_gwf_tracer_t *
cs_gwf_add_user_tracer(const char * eq_name,const char * var_name,cs_gwf_tracer_setup_t * setup,cs_gwf_tracer_add_terms_t * add_terms)2275 cs_gwf_add_user_tracer(const char *eq_name,
2276 const char *var_name,
2277 cs_gwf_tracer_setup_t *setup,
2278 cs_gwf_tracer_add_terms_t *add_terms)
2279 {
2280 cs_gwf_t *gw = cs_gwf_main_structure;
2281
2282 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2283
2284 /* Set the advection field structure */
2285 cs_adv_field_t *adv = _get_l_adv_field(gw);
2286 int tr_id = gw->n_tracers;
2287 cs_gwf_tracer_t *tracer = cs_gwf_tracer_init(tr_id,
2288 eq_name,
2289 var_name,
2290 adv,
2291 CS_GWF_TRACER_USER);
2292
2293 gw->n_tracers += 1;
2294 BFT_REALLOC(gw->tracers, gw->n_tracers, cs_gwf_tracer_t *);
2295 BFT_REALLOC(gw->finalize_tracer_setup,
2296 gw->n_tracers, cs_gwf_tracer_setup_t *);
2297 BFT_REALLOC(gw->add_tracer_terms,
2298 gw->n_tracers, cs_gwf_tracer_add_terms_t *);
2299
2300 gw->tracers[tr_id] = tracer;
2301 gw->finalize_tracer_setup[tr_id] = setup;
2302 gw->add_tracer_terms[tr_id] = add_terms;
2303
2304 return tracer;
2305 }
2306
2307 /*----------------------------------------------------------------------------*/
2308 /*!
2309 * \brief Retrieve the pointer to the cs_gwf_tracer_t structure associated to
2310 * the name given as parameter
2311 *
2312 * \param[in] eq_name name of the tracer equation
2313 *
2314 * \return the pointer to a cs_gwf_tracer_t structure
2315 */
2316 /*----------------------------------------------------------------------------*/
2317
2318 cs_gwf_tracer_t *
cs_gwf_tracer_by_name(const char * eq_name)2319 cs_gwf_tracer_by_name(const char *eq_name)
2320 {
2321 cs_gwf_t *gw = cs_gwf_main_structure;
2322
2323 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2324
2325 if (eq_name == NULL)
2326 return NULL;
2327
2328 for (int i = 0; i < gw->n_tracers; i++) {
2329 cs_gwf_tracer_t *tracer = gw->tracers[i];
2330 const char *name_to_cmp = cs_equation_get_name(tracer->eq);
2331 if (strcmp(eq_name, name_to_cmp) == 0)
2332 return tracer;
2333 }
2334
2335 return NULL;
2336 }
2337
2338 /*----------------------------------------------------------------------------*/
2339 /*!
2340 * \brief Add new terms if needed (such as diffusion or reaction) to tracer
2341 * equations according to the settings
2342 */
2343 /*----------------------------------------------------------------------------*/
2344
2345 void
cs_gwf_add_tracer_terms(void)2346 cs_gwf_add_tracer_terms(void)
2347 {
2348 cs_gwf_t *gw = cs_gwf_main_structure;
2349
2350 /* Sanity checks */
2351 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2352
2353 /* Loop on tracer equations */
2354 for (int i = 0; i < gw->n_tracers; i++)
2355 gw->add_tracer_terms[i](gw->tracers[i]);
2356 }
2357
2358 /*----------------------------------------------------------------------------*/
2359 /*!
2360 * \brief Predefined settings for the Richards equation and the related
2361 * equations defining the groundwater flow module.
2362 * At this stage, all soils have been defined and equatino parameters
2363 * are set. Create new cs_field_t structures according to the setting.
2364 */
2365 /*----------------------------------------------------------------------------*/
2366
2367 void
cs_gwf_init_setup(void)2368 cs_gwf_init_setup(void)
2369 {
2370 cs_gwf_t *gw = cs_gwf_main_structure;
2371
2372 /* Sanity checks */
2373 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2374 cs_gwf_soil_check();
2375
2376 switch (gw->model) {
2377
2378 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
2379 _sspf_init_setup(gw->model_context);
2380 break;
2381
2382 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
2383 _uspf_init_setup(gw->model_context);
2384 break;
2385
2386 case CS_GWF_MODEL_TWO_PHASE:
2387 _mtpf_init_setup(gw->model_context);
2388 break;
2389
2390 default:
2391 bft_error(__FILE__, __LINE__, 0,
2392 " %s: Invalid model type for the GroundWater Flow module.\n",
2393 __func__);
2394 }
2395 }
2396
2397 /*----------------------------------------------------------------------------*/
2398 /*!
2399 * \brief Last initialization step of the groundwater flow module. At this
2400 * stage, the mesh quantities are defined.
2401 *
2402 * \param[in] connect pointer to a cs_cdo_connect_t structure
2403 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2404 */
2405 /*----------------------------------------------------------------------------*/
2406
2407 void
cs_gwf_finalize_setup(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)2408 cs_gwf_finalize_setup(const cs_cdo_connect_t *connect,
2409 const cs_cdo_quantities_t *quant)
2410 {
2411 cs_gwf_t *gw = cs_gwf_main_structure;
2412 if (gw == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_gw));
2413
2414 switch (gw->model) {
2415
2416 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
2417 _sspf_finalize_setup(connect, quant, gw->model_context);
2418 break;
2419
2420 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
2421 _uspf_finalize_setup(connect, quant, gw->model_context);
2422 break;
2423
2424 case CS_GWF_MODEL_TWO_PHASE:
2425 _mtpf_finalize_setup(connect, quant, gw->model_context);
2426 break;
2427
2428 default:
2429 bft_error(__FILE__, __LINE__, 0,
2430 " %s: Invalid model type for the GroundWater Flow module.\n",
2431 __func__);
2432 }
2433
2434 /* Store the soil id for each cell */
2435
2436 cs_gwf_soil_check();
2437 cs_gwf_build_cell2soil(quant->n_cells);
2438
2439 /* Loop on tracer equations. Link the advection field to each tracer
2440 equation */
2441
2442 cs_adv_field_t *adv = _get_l_adv_field(gw);
2443
2444 for (int i = 0; i < gw->n_tracers; i++) {
2445
2446 if (gw->model == CS_GWF_MODEL_SATURATED_SINGLE_PHASE)
2447 gw->finalize_tracer_setup[i](connect, quant, adv, NULL, gw->tracers[i]);
2448
2449 else {
2450
2451 cs_field_t *f = cs_field_by_name("liquid_saturation");
2452 assert(f != NULL);
2453
2454 gw->finalize_tracer_setup[i](connect, quant, adv, f->val, gw->tracers[i]);
2455
2456 }
2457
2458 } /* Loop on tracers */
2459 }
2460
2461 /*----------------------------------------------------------------------------*/
2462 /*!
2463 * \brief Update the groundwater system (pressure head, head in law, moisture
2464 * content, darcian velocity, soil capacity or permeability if needed)
2465 *
2466 * \param[in] mesh pointer to a cs_mesh_t structure
2467 * \param[in] connect pointer to a cs_cdo_connect_t structure
2468 * \param[in] quant pointer to a cs_cdo_quantities_t structure
2469 * \param[in] ts pointer to a cs_time_step_t structure
2470 * \param[in] cur2prev true or false
2471 */
2472 /*----------------------------------------------------------------------------*/
2473
2474 void
cs_gwf_update(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts,bool cur2prev)2475 cs_gwf_update(const cs_mesh_t *mesh,
2476 const cs_cdo_connect_t *connect,
2477 const cs_cdo_quantities_t *quant,
2478 const cs_time_step_t *ts,
2479 bool cur2prev)
2480 {
2481 cs_gwf_t *gw = cs_gwf_main_structure;
2482
2483 /* Sanity checks */
2484
2485 if (gw == NULL)
2486 bft_error(__FILE__, __LINE__, 0,
2487 "%s: Groundwater module is not allocated.", __func__);
2488
2489 cs_real_t time_eval = 0.;
2490 switch (gw->model) {
2491
2492 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
2493 time_eval = _sspf_updates(connect, quant, ts, cur2prev, gw->model_context);
2494 break;
2495
2496 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
2497 time_eval = _uspf_updates(mesh, connect, quant, ts, cur2prev,
2498 gw->model_context);
2499 break;
2500
2501 case CS_GWF_MODEL_TWO_PHASE:
2502 time_eval = _mtpf_updates(mesh, connect, quant, ts, cur2prev,
2503 gw->model_context);
2504 break;
2505
2506 default:
2507 bft_error(__FILE__, __LINE__, 0,
2508 " %s: Invalid model type for the GroundWater Flow module.\n",
2509 __func__);
2510 }
2511
2512 /* Update the diffusivity tensor associated to each tracer equation since the
2513 Darcy velocity may have changed */
2514 for (int i = 0; i < gw->n_tracers; i++) {
2515
2516 cs_gwf_tracer_t *tracer = gw->tracers[i];
2517 if (tracer->update_diff_tensor != NULL)
2518 tracer->update_diff_tensor(tracer, time_eval, mesh, connect, quant);
2519
2520 }
2521
2522 }
2523
2524 /*----------------------------------------------------------------------------*/
2525 /*!
2526 * \brief Compute the steady-state of the groundwater flows module.
2527 * Nothing is done if all equations are unsteady.
2528 *
2529 * \param[in] mesh pointer to a cs_mesh_t structure
2530 * \param[in] time_step pointer to a cs_time_step_t structure
2531 * \param[in] connect pointer to a cs_cdo_connect_t structure
2532 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
2533 */
2534 /*----------------------------------------------------------------------------*/
2535
2536 void
cs_gwf_compute_steady_state(const cs_mesh_t * mesh,const cs_time_step_t * time_step,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq)2537 cs_gwf_compute_steady_state(const cs_mesh_t *mesh,
2538 const cs_time_step_t *time_step,
2539 const cs_cdo_connect_t *connect,
2540 const cs_cdo_quantities_t *cdoq)
2541 {
2542 cs_gwf_t *gw = cs_gwf_main_structure;
2543
2544 switch (gw->model) {
2545
2546 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
2547 {
2548 cs_gwf_saturated_single_phase_t *mc = gw->model_context;
2549
2550 _spf_compute_steady_state(mesh, time_step, connect, cdoq, mc->richards);
2551 }
2552 break;
2553
2554 default:
2555 bft_error(__FILE__, __LINE__, 0,
2556 " %s: Invalid model type for steady-state computations with the"
2557 " GroundWater Flow module.\n",
2558 __func__);
2559 }
2560
2561 /* Compute tracers */
2562 /* --------------- */
2563
2564 for (int i = 0; i < gw->n_tracers; i++) {
2565
2566 cs_gwf_tracer_t *tracer = gw->tracers[i];
2567
2568 if (cs_equation_is_steady(tracer->eq)) {
2569
2570 /* Solve the algebraic system */
2571 cs_equation_solve_steady_state(mesh, tracer->eq);
2572
2573 if (tracer->update_precipitation != NULL)
2574 tracer->update_precipitation(tracer,
2575 time_step->t_cur,
2576 mesh, connect, cdoq);
2577
2578 } /* Solve this equation which is steady */
2579
2580 } /* Loop on tracer equations */
2581
2582 }
2583
2584 /*----------------------------------------------------------------------------*/
2585 /*!
2586 * \brief Compute the system related to groundwater flows module
2587 *
2588 * \param[in] mesh pointer to a cs_mesh_t structure
2589 * \param[in] time_step pointer to a cs_time_step_t structure
2590 * \param[in] connect pointer to a cs_cdo_connect_t structure
2591 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
2592 */
2593 /*----------------------------------------------------------------------------*/
2594
2595 void
cs_gwf_compute(const cs_mesh_t * mesh,const cs_time_step_t * time_step,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq)2596 cs_gwf_compute(const cs_mesh_t *mesh,
2597 const cs_time_step_t *time_step,
2598 const cs_cdo_connect_t *connect,
2599 const cs_cdo_quantities_t *cdoq)
2600 {
2601 cs_gwf_t *gw = cs_gwf_main_structure;
2602
2603 switch (gw->model) {
2604
2605 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
2606 {
2607 cs_gwf_saturated_single_phase_t *mc = gw->model_context;
2608
2609 _spf_compute(mesh, time_step, connect, cdoq, mc->richards);
2610 }
2611 break;
2612
2613 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
2614 {
2615 cs_gwf_unsaturated_single_phase_t *mc = gw->model_context;
2616
2617 _spf_compute(mesh, time_step, connect, cdoq, mc->richards);
2618 }
2619 break;
2620
2621 case CS_GWF_MODEL_TWO_PHASE:
2622 _mtpf_compute(mesh, time_step, connect, cdoq, gw->model_context);
2623 break;
2624
2625 default:
2626 bft_error(__FILE__, __LINE__, 0,
2627 " %s: Invalid model type for the GroundWater Flow module.\n",
2628 __func__);
2629 }
2630
2631 /* Compute tracers */
2632 /* --------------- */
2633 bool cur2prev = true;
2634
2635 for (int i = 0; i < gw->n_tracers; i++) {
2636
2637 cs_gwf_tracer_t *tracer = gw->tracers[i];
2638
2639 if (!cs_equation_is_steady(tracer->eq)) { /* unsteady ? */
2640
2641 /* Solve the algebraic system. By default, a current to previous operation
2642 is performed */
2643 cs_equation_solve(cur2prev, mesh, tracer->eq);
2644
2645 if (tracer->update_precipitation != NULL)
2646 tracer->update_precipitation(tracer,
2647 time_step->t_cur,
2648 mesh, connect, cdoq);
2649
2650 } /* Solve this equation which is unsteady */
2651
2652 } /* Loop on tracer equations */
2653
2654 }
2655
2656 /*----------------------------------------------------------------------------*/
2657 /*!
2658 * \brief Compute the integral over a given set of cells of the field related
2659 * to a tracer equation. This integral turns out to be exact for linear
2660 * functions.
2661 *
2662 * \param[in] connect pointer to a \ref cs_cdo_connect_t structure
2663 * \param[in] cdoq pointer to a \ref cs_cdo_quantities_t structure
2664 * \param[in] tracer pointer to a \ref cs_gwf_tracer_t structure
2665 * \param[in] z_name name of the volumic zone where the integral is done
2666 * (if NULL or "" all cells are considered)
2667 *
2668 * \return the value of the integral
2669 */
2670 /*----------------------------------------------------------------------------*/
2671
2672 cs_real_t
cs_gwf_integrate_tracer(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,const cs_gwf_tracer_t * tracer,const char * z_name)2673 cs_gwf_integrate_tracer(const cs_cdo_connect_t *connect,
2674 const cs_cdo_quantities_t *cdoq,
2675 const cs_gwf_tracer_t *tracer,
2676 const char *z_name)
2677 {
2678 cs_gwf_t *gw = cs_gwf_main_structure;
2679
2680 /* Sanity checks */
2681
2682 if (gw == NULL)
2683 bft_error(__FILE__, __LINE__, 0,
2684 "%s: Groundwater module is not allocated.", __func__);
2685
2686 const int z_id = cs_get_vol_zone_id(z_name);
2687 const cs_zone_t *zone = cs_volume_zone_by_id(z_id);
2688
2689 if (gw->model == CS_GWF_MODEL_SATURATED_SINGLE_PHASE)
2690 return cs_gwf_tracer_integrate_sat(connect, cdoq, tracer, zone);
2691 else
2692 return cs_gwf_tracer_integrate(connect, cdoq, tracer, zone);
2693
2694 }
2695
2696 /*----------------------------------------------------------------------------*/
2697 /*!
2698 * \brief Predefined extra-operations for the groundwater flow module
2699 *
2700 * \param[in] connect pointer to a cs_cdo_connect_t structure
2701 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
2702 */
2703 /*----------------------------------------------------------------------------*/
2704
2705 void
cs_gwf_extra_op(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq)2706 cs_gwf_extra_op(const cs_cdo_connect_t *connect,
2707 const cs_cdo_quantities_t *cdoq)
2708 {
2709 cs_gwf_t *gw = cs_gwf_main_structure;
2710
2711 if (gw == NULL)
2712 return;
2713
2714 switch (gw->model) {
2715
2716 case CS_GWF_MODEL_SATURATED_SINGLE_PHASE:
2717 _sspf_extra_op(connect, cdoq, gw->model_context);
2718 break;
2719
2720 case CS_GWF_MODEL_UNSATURATED_SINGLE_PHASE:
2721 _uspf_extra_op(connect, cdoq, gw->model_context);
2722 break;
2723
2724 case CS_GWF_MODEL_TWO_PHASE:
2725 _mtpf_extra_op(connect, cdoq, gw->model_context);
2726 break;
2727
2728 default:
2729 bft_error(__FILE__, __LINE__, 0,
2730 " %s: Invalid model type for the GroundWater Flow module.\n",
2731 __func__);
2732 }
2733
2734 }
2735
2736 /*----------------------------------------------------------------------------*/
2737 /*!
2738 * \brief Predefined post-processing output for the groundwater flow module
2739 * in case of saturated single-phase flows (sspf) in porous media.
2740 * Prototype of this function is given since it is a function pointer
2741 * defined in cs_post.h (\ref cs_post_time_mesh_dep_output_t)
2742 *
2743 * \param[in, out] input pointer to a optional structure (here a
2744 * cs_gwf_t structure)
2745 * \param[in] mesh_id id of the output mesh for the current call
2746 * \param[in] cat_id category id of the output mesh for this call
2747 * \param[in] ent_flag indicate global presence of cells (ent_flag[0]),
2748 * interior faces (ent_flag[1]), boundary faces
2749 * (ent_flag[2]), particles (ent_flag[3]) or probes
2750 * (ent_flag[4])
2751 * \param[in] n_cells local number of cells of post_mesh
2752 * \param[in] n_i_faces local number of interior faces of post_mesh
2753 * \param[in] n_b_faces local number of boundary faces of post_mesh
2754 * \param[in] cell_ids list of cells (0 to n-1)
2755 * \param[in] i_face_ids list of interior faces (0 to n-1)
2756 * \param[in] b_face_ids list of boundary faces (0 to n-1)
2757 * \param[in] time_step pointer to a cs_time_step_t struct.
2758 */
2759 /*----------------------------------------------------------------------------*/
2760
2761 void
cs_gwf_extra_post_sspf(void * input,int mesh_id,int cat_id,int ent_flag[5],cs_lnum_t n_cells,cs_lnum_t n_i_faces,cs_lnum_t n_b_faces,const cs_lnum_t cell_ids[],const cs_lnum_t i_face_ids[],const cs_lnum_t b_face_ids[],const cs_time_step_t * time_step)2762 cs_gwf_extra_post_sspf(void *input,
2763 int mesh_id,
2764 int cat_id,
2765 int ent_flag[5],
2766 cs_lnum_t n_cells,
2767 cs_lnum_t n_i_faces,
2768 cs_lnum_t n_b_faces,
2769 const cs_lnum_t cell_ids[],
2770 const cs_lnum_t i_face_ids[],
2771 const cs_lnum_t b_face_ids[],
2772 const cs_time_step_t *time_step)
2773 {
2774 CS_UNUSED(cat_id);
2775 CS_UNUSED(ent_flag);
2776 CS_UNUSED(n_cells);
2777 CS_UNUSED(n_i_faces);
2778 CS_UNUSED(n_b_faces);
2779 CS_UNUSED(cell_ids);
2780 CS_UNUSED(i_face_ids);
2781 CS_UNUSED(b_face_ids);
2782 CS_UNUSED(time_step);
2783
2784 if (input == NULL)
2785 return;
2786
2787 const cs_gwf_t *gw = (const cs_gwf_t *)input;
2788
2789 if (mesh_id == CS_POST_MESH_VOLUME) {
2790
2791 const cs_gwf_saturated_single_phase_t *mc = gw->model_context;
2792
2793 if (gw->post_flag & CS_GWF_POST_DARCY_FLUX_DIVERGENCE) {
2794
2795 cs_adv_field_t *adv = _get_l_adv_field(gw);
2796
2797 /* Only case avalaible up to now */
2798
2799 if (cs_advection_field_get_deftype(adv) == CS_XDEF_BY_ARRAY) {
2800
2801 cs_real_t *divergence =
2802 cs_advection_field_divergence_at_vertices(adv, time_step->t_cur);
2803
2804 cs_post_write_vertex_var(mesh_id,
2805 CS_POST_WRITER_DEFAULT,
2806 "darcy_flux_divergence",
2807 1,
2808 false,
2809 false,
2810 CS_POST_TYPE_cs_real_t,
2811 divergence,
2812 time_step);
2813
2814 BFT_FREE(divergence);
2815
2816 }
2817
2818 } /* Post-processing of the divergence is requested */
2819
2820 if (gw->post_flag & CS_GWF_POST_LIQUID_SATURATION) {
2821
2822 cs_real_t *liquid_saturation = NULL;
2823 BFT_MALLOC(liquid_saturation, n_cells, cs_real_t);
2824
2825 cs_property_eval_at_cells(time_step->t_cur,
2826 mc->moisture_content,
2827 liquid_saturation);
2828
2829 cs_post_write_var(mesh_id,
2830 CS_POST_WRITER_DEFAULT,
2831 "liquid_saturation",
2832 1,
2833 false,
2834 false,
2835 CS_POST_TYPE_cs_real_t,
2836 liquid_saturation,
2837 NULL,
2838 NULL,
2839 time_step);
2840
2841 BFT_FREE(liquid_saturation);
2842
2843 } /* Post-processing of the liquid saturation */
2844
2845 if (gw->post_flag & CS_GWF_POST_PERMEABILITY) {
2846
2847 int dim = cs_property_get_dim(gw->abs_permeability);
2848 cs_real_t *permeability = NULL;
2849 BFT_MALLOC(permeability, dim*n_cells, cs_real_t);
2850
2851 cs_property_eval_at_cells(time_step->t_cur,
2852 gw->abs_permeability,
2853 permeability);
2854
2855 cs_post_write_var(mesh_id,
2856 CS_POST_WRITER_DEFAULT,
2857 "permeability",
2858 dim,
2859 true,
2860 false,
2861 CS_POST_TYPE_cs_real_t,
2862 permeability,
2863 NULL,
2864 NULL,
2865 time_step);
2866
2867 BFT_FREE(permeability);
2868
2869 } /* Post-processing of the permeability field */
2870
2871 } /* volume mesh id */
2872
2873 }
2874
2875 /*----------------------------------------------------------------------------*/
2876 /*!
2877 * \brief Predefined post-processing output for the groundwater flow module
2878 * in case of unsaturated single-phase flows (uspf) in porous media.
2879 * Prototype of this function is given since it is a function pointer
2880 * defined in cs_post.h (\ref cs_post_time_mesh_dep_output_t)
2881 *
2882 * \param[in, out] input pointer to a optional structure (here a
2883 * cs_gwf_t structure)
2884 * \param[in] mesh_id id of the output mesh for the current call
2885 * \param[in] cat_id category id of the output mesh for this call
2886 * \param[in] ent_flag indicate global presence of cells (ent_flag[0]),
2887 * interior faces (ent_flag[1]), boundary faces
2888 * (ent_flag[2]), particles (ent_flag[3]) or probes
2889 * (ent_flag[4])
2890 * \param[in] n_cells local number of cells of post_mesh
2891 * \param[in] n_i_faces local number of interior faces of post_mesh
2892 * \param[in] n_b_faces local number of boundary faces of post_mesh
2893 * \param[in] cell_ids list of cells (0 to n-1)
2894 * \param[in] i_face_ids list of interior faces (0 to n-1)
2895 * \param[in] b_face_ids list of boundary faces (0 to n-1)
2896 * \param[in] time_step pointer to a cs_time_step_t struct.
2897 */
2898 /*----------------------------------------------------------------------------*/
2899
2900 void
cs_gwf_extra_post_uspf(void * input,int mesh_id,int cat_id,int ent_flag[5],cs_lnum_t n_cells,cs_lnum_t n_i_faces,cs_lnum_t n_b_faces,const cs_lnum_t cell_ids[],const cs_lnum_t i_face_ids[],const cs_lnum_t b_face_ids[],const cs_time_step_t * time_step)2901 cs_gwf_extra_post_uspf(void *input,
2902 int mesh_id,
2903 int cat_id,
2904 int ent_flag[5],
2905 cs_lnum_t n_cells,
2906 cs_lnum_t n_i_faces,
2907 cs_lnum_t n_b_faces,
2908 const cs_lnum_t cell_ids[],
2909 const cs_lnum_t i_face_ids[],
2910 const cs_lnum_t b_face_ids[],
2911 const cs_time_step_t *time_step)
2912 {
2913 CS_UNUSED(cat_id);
2914 CS_UNUSED(ent_flag);
2915 CS_UNUSED(n_cells);
2916 CS_UNUSED(n_i_faces);
2917 CS_UNUSED(n_b_faces);
2918 CS_UNUSED(cell_ids);
2919 CS_UNUSED(i_face_ids);
2920 CS_UNUSED(b_face_ids);
2921 CS_UNUSED(time_step);
2922
2923 if (input == NULL)
2924 return;
2925
2926 const cs_gwf_t *gw = (const cs_gwf_t *)input;
2927
2928 if (mesh_id == CS_POST_MESH_VOLUME) {
2929 if (gw->post_flag & CS_GWF_POST_DARCY_FLUX_DIVERGENCE) {
2930
2931 cs_adv_field_t *adv = _get_l_adv_field(gw);
2932
2933 /* Only case avalaible up to now */
2934 if (cs_advection_field_get_deftype(adv) == CS_XDEF_BY_ARRAY) {
2935
2936 cs_real_t *divergence =
2937 cs_advection_field_divergence_at_vertices(adv, time_step->t_cur);
2938
2939 cs_post_write_vertex_var(mesh_id,
2940 CS_POST_WRITER_DEFAULT,
2941 "darcy_flux_divergence",
2942 1,
2943 false,
2944 false,
2945 CS_POST_TYPE_cs_real_t,
2946 divergence,
2947 time_step);
2948
2949 BFT_FREE(divergence);
2950
2951 }
2952
2953 } /* Post-processing of the divergence is requested */
2954 } /* volume mesh id */
2955
2956 }
2957
2958 /*----------------------------------------------------------------------------*/
2959 /*!
2960 * \brief Predefined post-processing output for the groundwater flow module
2961 * in case of miscible two-phase flows (mtpf) in porous media.
2962 * Prototype of this function is given since it is a function pointer
2963 * defined in cs_post.h (\ref cs_post_time_mesh_dep_output_t)
2964 *
2965 * \param[in, out] input pointer to a optional structure (here a
2966 * cs_gwf_t structure)
2967 * \param[in] mesh_id id of the output mesh for the current call
2968 * \param[in] cat_id category id of the output mesh for this call
2969 * \param[in] ent_flag indicate global presence of cells (ent_flag[0]),
2970 * interior faces (ent_flag[1]), boundary faces
2971 * (ent_flag[2]), particles (ent_flag[3]) or probes
2972 * (ent_flag[4])
2973 * \param[in] n_cells local number of cells of post_mesh
2974 * \param[in] n_i_faces local number of interior faces of post_mesh
2975 * \param[in] n_b_faces local number of boundary faces of post_mesh
2976 * \param[in] cell_ids list of cells (0 to n-1)
2977 * \param[in] i_face_ids list of interior faces (0 to n-1)
2978 * \param[in] b_face_ids list of boundary faces (0 to n-1)
2979 * \param[in] time_step pointer to a cs_time_step_t struct.
2980 */
2981 /*----------------------------------------------------------------------------*/
2982
2983 void
cs_gwf_extra_post_mtpf(void * input,int mesh_id,int cat_id,int ent_flag[5],cs_lnum_t n_cells,cs_lnum_t n_i_faces,cs_lnum_t n_b_faces,const cs_lnum_t cell_ids[],const cs_lnum_t i_face_ids[],const cs_lnum_t b_face_ids[],const cs_time_step_t * time_step)2984 cs_gwf_extra_post_mtpf(void *input,
2985 int mesh_id,
2986 int cat_id,
2987 int ent_flag[5],
2988 cs_lnum_t n_cells,
2989 cs_lnum_t n_i_faces,
2990 cs_lnum_t n_b_faces,
2991 const cs_lnum_t cell_ids[],
2992 const cs_lnum_t i_face_ids[],
2993 const cs_lnum_t b_face_ids[],
2994 const cs_time_step_t *time_step)
2995 {
2996 CS_UNUSED(cat_id);
2997 CS_UNUSED(ent_flag);
2998 CS_UNUSED(n_cells);
2999 CS_UNUSED(n_i_faces);
3000 CS_UNUSED(n_b_faces);
3001 CS_UNUSED(cell_ids);
3002 CS_UNUSED(i_face_ids);
3003 CS_UNUSED(b_face_ids);
3004 CS_UNUSED(time_step);
3005
3006 if (input == NULL)
3007 return;
3008
3009 const cs_gwf_t *gw = (const cs_gwf_t *)input;
3010 const cs_gwf_miscible_two_phase_t *mc =
3011 (const cs_gwf_miscible_two_phase_t *)gw->model_context;
3012
3013 if (mesh_id == CS_POST_MESH_VOLUME) {
3014
3015 if (gw->post_flag & CS_GWF_POST_DARCY_FLUX_DIVERGENCE) {
3016
3017 cs_adv_field_t *adv = _get_l_adv_field(gw);
3018
3019 /* Only case avalaible up to now */
3020 if (cs_advection_field_get_deftype(adv) == CS_XDEF_BY_ARRAY) {
3021
3022 cs_real_t *divergence =
3023 cs_advection_field_divergence_at_vertices(adv, time_step->t_cur);
3024
3025 cs_post_write_vertex_var(mesh_id,
3026 CS_POST_WRITER_DEFAULT,
3027 "darcy_flux_divergence",
3028 1,
3029 false,
3030 false,
3031 CS_POST_TYPE_cs_real_t,
3032 divergence,
3033 time_step);
3034
3035 BFT_FREE(divergence);
3036
3037 }
3038
3039 } /* Post-processing of the divergence is requested */
3040
3041 if (gw->post_flag & CS_GWF_POST_PERMEABILITY) {
3042
3043 int dim = cs_property_get_dim(gw->abs_permeability);
3044 cs_real_t *permeability = NULL;
3045 BFT_MALLOC(permeability, dim*n_cells, cs_real_t);
3046
3047 cs_property_eval_at_cells(time_step->t_cur,
3048 gw->abs_permeability,
3049 permeability);
3050
3051 /* permeability = abs_permeability * l_rel_permeability */
3052
3053 assert(mc->l_rel_permeability != NULL);
3054 for (cs_lnum_t c = 0; c < n_cells; c++)
3055 for (int k = 0; k < dim; k++)
3056 permeability[dim*c+k] *= mc->l_rel_permeability[c];
3057
3058 cs_post_write_var(mesh_id,
3059 CS_POST_WRITER_DEFAULT,
3060 "permeability",
3061 dim,
3062 true,
3063 false,
3064 CS_POST_TYPE_cs_real_t,
3065 permeability,
3066 NULL,
3067 NULL,
3068 time_step);
3069
3070 BFT_FREE(permeability);
3071
3072 } /* Post-processing of the permeability field */
3073
3074 } /* volume mesh id */
3075
3076 }
3077
3078 /*----------------------------------------------------------------------------*/
3079
3080 #undef _dp3
3081
3082 END_C_DECLS
3083