1 /*============================================================================
2 * Helper 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 <stdlib.h>
37 #include <string.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include <bft_mem.h>
44
45 #include "cs_advection_field.h"
46 #include "cs_parall.h"
47
48 /*----------------------------------------------------------------------------
49 * Header for the current file
50 *----------------------------------------------------------------------------*/
51
52 #include "cs_gwf_priv.h"
53
54 /*----------------------------------------------------------------------------*/
55
56 BEGIN_C_DECLS
57
58 /*!
59 \file cs_gwf_priv.c
60
61 \brief Helper functions dedicated to groundwater flows when using CDO schemes
62
63 */
64
65 /*============================================================================
66 * Local macro definitions
67 *============================================================================*/
68
69 #define CS_GWF_PRIV_DBG 0
70
71 /*============================================================================
72 * Local definitions
73 *============================================================================*/
74
75 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
76
77 /*============================================================================
78 * Static global variables
79 *============================================================================*/
80
81 /*============================================================================
82 * Private function prototypes
83 *============================================================================*/
84
85 /*----------------------------------------------------------------------------*/
86 /*!
87 * \brief Compute the associated boundary Darcy flux for each vertex of
88 * boundary faces.
89 * Case of a vertex-based discretization and single-phase flows in
90 * porous media (saturated or not).
91 *
92 * \param[in] t_eval time at which one performs the evaluation
93 * \param[in] eq pointer to the equation related to this Darcy flux
94 * \param[in, out] adv pointer to the Darcy advection field
95 */
96 /*----------------------------------------------------------------------------*/
97
98 static void
_update_vb_darcy_flux_at_boundary(cs_real_t t_eval,const cs_equation_t * eq,cs_adv_field_t * adv)99 _update_vb_darcy_flux_at_boundary(cs_real_t t_eval,
100 const cs_equation_t *eq,
101 cs_adv_field_t *adv)
102 {
103 if (adv->n_bdy_flux_defs > 1 ||
104 adv->bdy_flux_defs[0]->type != CS_XDEF_BY_ARRAY)
105 bft_error(__FILE__, __LINE__, 0,
106 " %s: Invalid definition of the advection field at the boundary",
107 __func__);
108
109 cs_xdef_t *def = adv->bdy_flux_defs[0];
110 cs_xdef_array_context_t *actx = (cs_xdef_array_context_t *)def->context;
111 cs_real_t *nflx_val = actx->values;
112
113 if (cs_flag_test(actx->loc, cs_flag_dual_closure_byf) == false)
114 bft_error(__FILE__, __LINE__, 0,
115 " %s: Invalid definition of the advection field at the boundary",
116 __func__);
117
118 cs_equation_compute_boundary_diff_flux(t_eval, eq, nflx_val);
119 }
120
121 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
122
123 /*============================================================================
124 * Public function prototypes
125 *============================================================================*/
126
127 /*----------------------------------------------------------------------------*/
128 /*!
129 * \brief Allocate and initialize a \ref cs_gwf_darcy_flux_t structure
130 *
131 * \param[in] loc_flag flag to define where the flux is defined
132 *
133 * \return a pointer to the newly allocated structure
134 */
135 /*----------------------------------------------------------------------------*/
136
137 cs_gwf_darcy_flux_t *
cs_gwf_darcy_flux_create(cs_flag_t loc_flag)138 cs_gwf_darcy_flux_create(cs_flag_t loc_flag)
139 {
140 cs_gwf_darcy_flux_t *darcy = NULL;
141
142 BFT_MALLOC(darcy, 1, cs_gwf_darcy_flux_t);
143
144 darcy->flux_location = loc_flag;
145 darcy->adv_field = NULL;
146 darcy->flux_val = NULL;
147 darcy->boundary_flux_val = NULL;
148
149 return darcy;
150 }
151
152 /*----------------------------------------------------------------------------*/
153 /*!
154 * \brief Free a \ref cs_gwf_darcy_flux_t structure
155 *
156 * \param[in, out] p_darcy pointer of pointer to the darcy structure
157 */
158 /*----------------------------------------------------------------------------*/
159
160 void
cs_gwf_darcy_flux_free(cs_gwf_darcy_flux_t ** p_darcy)161 cs_gwf_darcy_flux_free(cs_gwf_darcy_flux_t **p_darcy)
162 {
163 if (p_darcy == NULL)
164 return;
165 if (*p_darcy == NULL)
166 return;
167
168 cs_gwf_darcy_flux_t *darcy = *p_darcy;
169
170 /* flux_val and boundary_flux_val are allocated only if the related advection
171 field is defined by array. At the definition step, the GWF module has kept
172 the ownership of the lifecycle of the flux_val and boundary_flux_val
173 arrays. In this case, the lifecycle is not managed by the definition and
174 thus one has to free the arrays now. */
175
176 BFT_FREE(darcy->boundary_flux_val);
177 BFT_FREE(darcy->flux_val);
178
179 BFT_FREE(darcy);
180 *p_darcy = NULL;
181 }
182
183 /*----------------------------------------------------------------------------*/
184 /*!
185 * \brief Log a \ref cs_gwf_darcy_flux_t structure
186 *
187 * \param[in, out] darcy pointer to the darcy structure
188 */
189 /*----------------------------------------------------------------------------*/
190
191 void
cs_gwf_darcy_flux_log(cs_gwf_darcy_flux_t * darcy)192 cs_gwf_darcy_flux_log(cs_gwf_darcy_flux_t *darcy)
193 {
194 if (darcy == NULL)
195 return;
196 assert(darcy->adv_field != NULL);
197
198 cs_log_printf(CS_LOG_SETUP,
199 " * GWF | Darcy name: %s with flux location: %s\n",
200 darcy->adv_field->name,
201 cs_flag_str_location(darcy->flux_location));
202 }
203
204 /*----------------------------------------------------------------------------*/
205 /*!
206 * \brief Set the definition of the advection field attached to a
207 * \ref cs_gwf_darcy_flux_t structure
208 *
209 * \param[in] connect pointer to a cs_cdo_connect_t structure
210 * \param[in] quant pointer to a cs_cdo_quantities_t structure
211 * \param[in] space_scheme space discretization using this structure
212 * \param[in, out] darcy pointer to the darcy structure
213 */
214 /*----------------------------------------------------------------------------*/
215
216 void
cs_gwf_darcy_flux_define(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_param_space_scheme_t space_scheme,cs_gwf_darcy_flux_t * darcy)217 cs_gwf_darcy_flux_define(const cs_cdo_connect_t *connect,
218 const cs_cdo_quantities_t *quant,
219 cs_param_space_scheme_t space_scheme,
220 cs_gwf_darcy_flux_t *darcy)
221 {
222 if (darcy == NULL)
223 return;
224
225 cs_adv_field_t *adv = darcy->adv_field;
226 assert(adv != NULL);
227
228 switch (space_scheme) {
229
230 case CS_SPACE_SCHEME_CDOVB:
231 case CS_SPACE_SCHEME_CDOVCB:
232 {
233 const cs_adjacency_t *bf2v = connect->bf2v;
234
235 /* Define the flux of the advection field at the boundary */
236
237 cs_flag_t array_location = CS_FLAG_SCALAR | cs_flag_dual_closure_byf;
238 size_t array_size = bf2v->idx[quant->n_b_faces];
239
240 BFT_MALLOC(darcy->boundary_flux_val, array_size, cs_real_t);
241 memset(darcy->boundary_flux_val, 0, array_size*sizeof(cs_real_t));
242
243 /* Do not transfer the ownership */
244
245 cs_advection_field_def_boundary_flux_by_array(adv,
246 NULL,
247 array_location,
248 darcy->boundary_flux_val,
249 false,
250 bf2v->idx);
251
252 /* Define the advection field in the volume */
253
254 if (cs_flag_test(darcy->flux_location, cs_flag_dual_face_byc)) {
255
256 const cs_adjacency_t *c2e = connect->c2e;
257
258 array_location = CS_FLAG_SCALAR | darcy->flux_location;
259 array_size = c2e->idx[quant->n_cells];
260 BFT_MALLOC(darcy->flux_val, array_size, cs_real_t);
261 memset(darcy->flux_val, 0, array_size*sizeof(cs_real_t));
262
263 /* Do not transfer the ownership */
264
265 cs_advection_field_def_by_array(adv,
266 array_location,
267 darcy->flux_val,
268 false, /* transfer ownership */
269 c2e->idx);
270
271 /* Reset the type of advection field */
272
273 if (adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR)
274 adv->status -= CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR;
275 adv->status |= CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX;
276
277 }
278 else if (cs_flag_test(darcy->flux_location, cs_flag_primal_cell)) {
279
280 cs_field_t *cell_adv_field =
281 cs_advection_field_get_field(adv, CS_MESH_LOCATION_CELLS);
282 assert(cell_adv_field != NULL);
283
284 cs_advection_field_def_by_field(adv, cell_adv_field);
285
286 /* Reset the type of advection field */
287
288 if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX)
289 adv->status -= CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX;
290 adv->status |= CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR;
291
292 }
293 else
294 bft_error(__FILE__, __LINE__, 0,
295 " %s: Invalid location for the definition of the Darcy flux.",
296 __func__);
297
298 }
299 break;
300
301 case CS_SPACE_SCHEME_CDOFB: /* TODO */
302 bft_error(__FILE__, __LINE__, 0,
303 " %s: CDO-Fb space scheme not fully implemented.", __func__);
304 break;
305
306 default:
307 bft_error(__FILE__, __LINE__, 0, " %s: Invalid space scheme.", __func__);
308 break;
309
310 } /* Switch on the space scheme */
311
312 }
313
314 /*----------------------------------------------------------------------------*/
315 /*!
316 * \brief Update the advection field/arrays related to the Darcy flux.
317 *
318 * \param[in] t_eval time at which one performs the evaluation
319 * \param[in] eq pointer to the equation related to this Darcy flux
320 * \param[in] cur2prev true or false
321 * \param[in, out] darcy pointer to the darcy structure
322 */
323 /*----------------------------------------------------------------------------*/
324
325 void
cs_gwf_darcy_flux_update(const cs_real_t t_eval,const cs_equation_t * eq,bool cur2prev,cs_gwf_darcy_flux_t * darcy)326 cs_gwf_darcy_flux_update(const cs_real_t t_eval,
327 const cs_equation_t *eq,
328 bool cur2prev,
329 cs_gwf_darcy_flux_t *darcy)
330 {
331 cs_adv_field_t *adv = darcy->adv_field;
332 assert(adv != NULL);
333
334 /* Update the velocity field at cell centers induced by the Darcy flux */
335
336 cs_field_t *vel = cs_advection_field_get_field(adv, CS_MESH_LOCATION_CELLS);
337
338 assert(vel != NULL); /* Sanity check */
339 if (cur2prev)
340 cs_field_current_to_previous(vel);
341
342 /* Update arrays related to the Darcy flux:
343 * Compute the new darcian flux and darcian velocity inside each cell */
344
345 switch (cs_equation_get_space_scheme(eq)) {
346
347 case CS_SPACE_SCHEME_CDOVB:
348 case CS_SPACE_SCHEME_CDOVCB:
349
350 /* Update the array of flux values associated to the advection field */
351
352 if (cs_flag_test(darcy->flux_location, cs_flag_dual_face_byc)) {
353
354 assert(darcy->flux_val != NULL);
355 if (adv->definition->type != CS_XDEF_BY_ARRAY)
356 bft_error(__FILE__, __LINE__, 0,
357 " %s: Invalid definition of the advection field", __func__);
358
359 cs_equation_compute_diff_flux_cellwise(eq,
360 darcy->flux_location,
361 t_eval,
362 darcy->flux_val);
363
364 #if defined(DEBUG) && !defined(NDEBUG) && CS_GWF_PRIV_DBG > 2
365 cs_dbg_darray_to_listing("DARCIAN_FLUX_DFbyC",
366 connect->c2e->idx[cdoq->n_cells],
367 darcy->flux_val, 8);
368 #endif
369
370 /* Set the new values of the vector field at cell centers */
371
372 cs_advection_field_in_cells(adv, t_eval, vel->val);
373
374 }
375 else if (cs_flag_test(darcy->flux_location, cs_flag_primal_cell))
376 cs_equation_compute_diff_flux_cellwise(eq,
377 darcy->flux_location,
378 t_eval,
379 vel->val);
380
381 /* Update the Darcy flux at the boundary */
382
383 _update_vb_darcy_flux_at_boundary(t_eval, eq, adv);
384 break;
385
386 case CS_SPACE_SCHEME_CDOFB:
387 case CS_SPACE_SCHEME_HHO_P0:
388 bft_error(__FILE__, __LINE__, 0, " TODO.");
389 break;
390
391 default:
392 bft_error(__FILE__, __LINE__, 0, " Invalid space scheme.");
393
394 } /* End of switch */
395
396 #if defined(DEBUG) && !defined(NDEBUG) && CS_GWF_PRIV_DBG > 1
397 cs_dbg_darray_to_listing("DARCIAN_FLUX_CELL", 3*cdoq->n_cells, vel->val, 3);
398 #endif
399
400 cs_field_t *bdy_nflx =
401 cs_advection_field_get_field(adv, CS_MESH_LOCATION_BOUNDARY_FACES);
402
403 if (bdy_nflx != NULL) { /* Values of the Darcy flux at boundary face exist */
404
405 if (cur2prev)
406 cs_field_current_to_previous(bdy_nflx);
407
408 /* Set the new values of the field related to the normal boundary flux */
409
410 cs_advection_field_across_boundary(adv, t_eval, bdy_nflx->val);
411
412 }
413
414 }
415
416 /*----------------------------------------------------------------------------*/
417 /*!
418 * \brief Operate the balance by zone (relying on the splitting arising from
419 * the boundary settings) for the advection field attached to a \ref
420 * cs_gwf_darcy_flux_t structure
421 *
422 * \param[in] connect pointer to a cs_cdo_connect_t structure
423 * \param[in] quant pointer to a cs_cdo_quantities_t structure
424 * \param[in] eqp pointer to the set of equation parameters
425 * \param[in, out] darcy pointer to the darcy structure
426 */
427 /*----------------------------------------------------------------------------*/
428
429 void
cs_gwf_darcy_flux_balance(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_equation_param_t * eqp,cs_gwf_darcy_flux_t * darcy)430 cs_gwf_darcy_flux_balance(const cs_cdo_connect_t *connect,
431 const cs_cdo_quantities_t *quant,
432 const cs_equation_param_t *eqp,
433 cs_gwf_darcy_flux_t *darcy)
434 {
435 if (darcy == NULL)
436 return;
437
438 const cs_lnum_t n_b_faces = quant->n_b_faces;
439 const cs_adv_field_t *adv = darcy->adv_field;
440 assert(adv != NULL);
441
442 /* Try to retrieve the boundary flux values */
443
444 const cs_field_t *nflx =
445 cs_advection_field_get_field(adv, CS_MESH_LOCATION_BOUNDARY_FACES);
446 cs_real_t *flux_val = (nflx == NULL) ? darcy->boundary_flux_val : nflx->val;
447
448 if (flux_val == NULL && n_b_faces > 0) /* No value on which operates */
449 return;
450
451 const cs_adjacency_t *bf2v = connect->bf2v;
452
453 /* Define the balance by zone */
454
455 bool *is_counted = NULL;
456 BFT_MALLOC(is_counted, n_b_faces, bool);
457 # pragma omp parallel for if (n_b_faces > CS_THR_MIN)
458 for (int i = 0; i < n_b_faces; i++) is_counted[i] = false;
459
460 /* n_bc_defs + 1 to take into account the default boundary condition */
461
462 cs_real_t *balances = NULL;
463 BFT_MALLOC(balances, eqp->n_bc_defs + 1, cs_real_t);
464
465 for (int ibc = 0; ibc < eqp->n_bc_defs; ibc++) {
466
467 const cs_xdef_t *def = eqp->bc_defs[ibc];
468 const cs_zone_t *z = cs_boundary_zone_by_id(def->z_id);
469
470 balances[ibc] = 0;
471
472 if (nflx == NULL) { /* The definition of the boundary flux relies on the
473 bf2v adjacency */
474
475 #if defined(DEBUG) && !defined(NDEBUG)
476 cs_xdef_t *_def = adv->bdy_flux_defs[0];
477 cs_xdef_array_context_t *actx = (cs_xdef_array_context_t *)_def->context;
478
479 assert(adv->n_bdy_flux_defs == 1 && _def->type == CS_XDEF_BY_ARRAY);
480 assert(cs_flag_test(actx->loc, cs_flag_dual_closure_byf) == true);
481 #endif
482
483 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
484 const cs_lnum_t bf_id = z->elt_ids[i];
485 is_counted[bf_id] = true;
486 for (cs_lnum_t j = bf2v->idx[bf_id]; j < bf2v->idx[bf_id+1]; j++)
487 balances[ibc] += flux_val[j];
488 }
489
490 }
491 else {
492
493 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
494 const cs_lnum_t bf_id = z->elt_ids[i];
495 is_counted[bf_id] = true;
496 balances[ibc] += flux_val[bf_id];
497 }
498
499 } /* nflux is NULL ? */
500
501 } /* Loop on BC definitions */
502
503 /* Manage the default boundary condition if there is at least one face not
504 considered in the previous loop */
505
506 bool display = false;
507 balances[eqp->n_bc_defs] = 0.;
508 for (cs_lnum_t bf_id = 0; bf_id < n_b_faces; bf_id++) {
509 if (is_counted[bf_id] == false) {
510
511 display = true;
512 if (nflx == NULL) {
513
514 for (cs_lnum_t j = bf2v->idx[bf_id]; j < bf2v->idx[bf_id+1]; j++)
515 balances[eqp->n_bc_defs] += flux_val[j];
516
517 }
518 else
519 balances[eqp->n_bc_defs] += flux_val[bf_id];
520
521 } /* Not already counted */
522 } /* Loop on boundary faces */
523
524 int display_flag = display ? 1 : 0;
525
526 /* Parallel synchronizations */
527
528 cs_parall_max(1, CS_INT_TYPE, &display_flag);
529 cs_parall_sum(eqp->n_bc_defs + 1, CS_REAL_TYPE, balances);
530
531 /* Output into the default log file */
532
533 cs_log_printf(CS_LOG_DEFAULT,
534 "-b- Balance of %s across the boundary zones:\n", adv->name);
535
536 for (int ibc = 0; ibc < eqp->n_bc_defs; ibc++) {
537 const cs_zone_t *z = cs_boundary_zone_by_id((eqp->bc_defs[ibc])->z_id);
538 cs_log_printf(CS_LOG_DEFAULT, "-b- %-32s: % -5.3e\n",
539 z->name, balances[ibc]);
540 }
541
542 if (display_flag > 0)
543 cs_log_printf(CS_LOG_DEFAULT, "-b- %-32s: % -5.3e\n",
544 "Remaining part of the boundary", balances[eqp->n_bc_defs]);
545
546 /* Free memory */
547
548 BFT_FREE(is_counted);
549 BFT_FREE(balances);
550 }
551
552 /*----------------------------------------------------------------------------*/
553
554 END_C_DECLS
555