1 /*============================================================================
2 * Build an algebraic CDO vertex-based system for unsteady convection diffusion
3 * reaction of vector-valued equations with source terms
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <assert.h>
35 #include <float.h>
36 #include <math.h>
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <string.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include <bft_mem.h>
46
47 #include "cs_cdo_advection.h"
48 #include "cs_cdo_bc.h"
49 #include "cs_cdo_diffusion.h"
50 #include "cs_cdo_local.h"
51 #include "cs_cdovb_priv.h"
52 #include "cs_equation_bc.h"
53 #include "cs_equation_common.h"
54 #include "cs_evaluate.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_post.h"
61 #include "cs_quadrature.h"
62 #include "cs_reco.h"
63 #include "cs_scheme_geometry.h"
64 #include "cs_search.h"
65 #include "cs_sles.h"
66 #include "cs_source_term.h"
67 #include "cs_timer.h"
68
69 #if defined(DEBUG) && !defined(NDEBUG)
70 #include "cs_dbg.h"
71 #endif
72
73 /*----------------------------------------------------------------------------
74 * Header for the current file
75 *----------------------------------------------------------------------------*/
76
77 #include "cs_cdovb_vecteq.h"
78
79 /*----------------------------------------------------------------------------*/
80
81 BEGIN_C_DECLS
82
83 /*!
84 * \file cs_cdovb_vecteq.c
85 *
86 * \brief Build an algebraic CDO vertex-based system for unsteady
87 * convection-diffusion-reaction of vector-valued equations with
88 * source terms
89 *
90 */
91
92 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
93
94 /*=============================================================================
95 * Local Macro definitions and structure definitions
96 *============================================================================*/
97
98 #define CS_CDOVB_VECTEQ_DBG 0
99
100 /* Redefined the name of functions from cs_math to get shorter names */
101 #define _dp3 cs_math_3_dot_product
102
103 /*============================================================================
104 * Type definitions
105 *============================================================================*/
106
107 /* Algebraic system for CDO vertex-based discretization */
108
109 typedef struct _cs_cdovb_t cs_cdovb_vecteq_t;
110
111 /*============================================================================
112 * Private variables
113 *============================================================================*/
114
115 /* Structure to enable a full cellwise strategy during the system building */
116
117 static cs_cell_sys_t **_vvb_cell_system = NULL;
118 static cs_cell_builder_t **_vvb_cell_builder = NULL;
119
120 /* Pointer to shared structures */
121
122 static const cs_cdo_quantities_t *cs_shared_quant;
123 static const cs_cdo_connect_t *cs_shared_connect;
124 static const cs_time_step_t *cs_shared_time_step;
125 static const cs_matrix_structure_t *cs_shared_ms;
126
127 /*============================================================================
128 * Private function prototypes
129 *============================================================================*/
130
131 /*----------------------------------------------------------------------------*/
132 /*!
133 * \brief Initialize the local builder structure used for building the system
134 * cellwise
135 *
136 * \param[in] connect pointer to a cs_cdo_connect_t structure
137 *
138 * \return a pointer to a new allocated cs_cell_builder_t structure
139 */
140 /*----------------------------------------------------------------------------*/
141
142 static cs_cell_builder_t *
_vvb_create_cell_builder(const cs_cdo_connect_t * connect)143 _vvb_create_cell_builder(const cs_cdo_connect_t *connect)
144 {
145 const int n_vc = connect->n_max_vbyc;
146 const int n_ec = connect->n_max_ebyc;
147
148 cs_cell_builder_t *cb = cs_cell_builder_create();
149
150 BFT_MALLOC(cb->ids, n_vc, int);
151 memset(cb->ids, 0, n_vc*sizeof(int));
152
153 int size = n_ec*(n_ec+1);
154 size = CS_MAX(4*n_ec + 3*n_vc, size);
155 BFT_MALLOC(cb->values, size, double);
156 memset(cb->values, 0, size*sizeof(cs_real_t));
157
158 size = 2*n_ec;
159 BFT_MALLOC(cb->vectors, size, cs_real_3_t);
160 memset(cb->vectors, 0, size*sizeof(cs_real_3_t));
161
162 /* Local square dense matrices used during the construction of
163 operators */
164
165 cb->aux = cs_sdm_square_create(n_ec);
166 cb->loc = cs_sdm_block33_create(n_vc, n_vc);
167
168 return cb;
169 }
170
171 /*----------------------------------------------------------------------------*/
172 /*!
173 * \brief Set the boundary conditions known from the settings
174 * Define an indirection array for the enforcement of internal DoFs
175 * only if needed
176 *
177 * \param[in] t_eval time at which one evaluates BCs
178 * \param[in] mesh pointer to a cs_mesh_t structure
179 * \param[in] eqp pointer to a cs_equation_param_t structure
180 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
181 * \param[in, out] vtx_bc_flag pointer to an array of BC flag for each vertex
182 */
183 /*----------------------------------------------------------------------------*/
184
185 static void
_vvb_setup(cs_real_t t_eval,const cs_mesh_t * mesh,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,cs_flag_t vtx_bc_flag[])186 _vvb_setup(cs_real_t t_eval,
187 const cs_mesh_t *mesh,
188 const cs_equation_param_t *eqp,
189 cs_equation_builder_t *eqb,
190 cs_flag_t vtx_bc_flag[])
191 {
192 assert(vtx_bc_flag != NULL); /* Sanity check */
193 const cs_cdo_quantities_t *quant = cs_shared_quant;
194 const cs_cdo_connect_t *connect = cs_shared_connect;
195
196 /* Compute the values of the Dirichlet BC */
197
198 BFT_MALLOC(eqb->dir_values, 3*quant->n_vertices, cs_real_t);
199
200 cs_equation_compute_dirichlet_vb(t_eval,
201 mesh,
202 quant,
203 connect,
204 eqp,
205 eqb->face_bc,
206 _vvb_cell_builder[0], /* static variable */
207 vtx_bc_flag,
208 eqb->dir_values);
209
210 /* Internal enforcement of DoFs */
211
212 if (cs_equation_param_has_internal_enforcement(eqp))
213 eqb->enforced_values =
214 cs_enforcement_define_at_vertices(connect,
215 eqp->n_enforcements,
216 eqp->enforcement_params);
217 }
218
219 /*----------------------------------------------------------------------------*/
220 /*!
221 * \brief Initialize the local structure for the current cell.
222 * Case of vector-valued CDO-Vb schemes
223 *
224 * \param[in] cm pointer to a cellwise view of the mesh
225 * \param[in] eqp pointer to a cs_equation_param_t structure
226 * \param[in] eqb pointer to a cs_equation_builder_t structure
227 * \param[in] vtx_bc_flag Flag related to BC associated to each vertex
228 * \param[in] field_tn values of the field at the last computed time
229 * \param[in, out] csys pointer to a cellwise view of the system
230 * \param[in, out] cb pointer to a cellwise builder
231 */
232 /*----------------------------------------------------------------------------*/
233
234 static void
_vvb_init_cell_system(const cs_cell_mesh_t * cm,const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_flag_t vtx_bc_flag[],const cs_real_t field_tn[],cs_cell_sys_t * csys,cs_cell_builder_t * cb)235 _vvb_init_cell_system(const cs_cell_mesh_t *cm,
236 const cs_equation_param_t *eqp,
237 const cs_equation_builder_t *eqb,
238 const cs_flag_t vtx_bc_flag[],
239 const cs_real_t field_tn[],
240 cs_cell_sys_t *csys,
241 cs_cell_builder_t *cb)
242 {
243 csys->c_id = cm->c_id;
244 csys->n_dofs = 3*cm->n_vc;
245
246 /* Cell-wise view of the linear system to build:
247 Initialize the local system */
248
249 cs_cell_sys_reset(cm->n_fc, csys); /* Generic part */
250
251 cs_sdm_block33_init(csys->mat, cm->n_vc, cm->n_vc);
252
253 for (int v = 0; v < cm->n_vc; v++) {
254 const cs_lnum_t v_id = cm->v_ids[v];
255 const cs_real_t *val_p = field_tn + 3*v_id;
256 for (int k = 0; k < 3; k++) {
257 csys->dof_ids[3*v + k] = 3*v_id + k;
258 csys->val_n[3*v + k] = val_p[k];
259 }
260 }
261
262 /* Store the local values attached to Dirichlet values if the current cell
263 has at least one border face */
264
265 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
266
267 /* Set the bc (specific part) */
268
269 cs_equation_vb_set_cell_bc(cm,
270 eqp,
271 eqb->face_bc,
272 vtx_bc_flag,
273 eqb->dir_values,
274 cb->t_bc_eval,
275 csys,
276 cb);
277
278 #if defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
279 cs_dbg_check_hmg_dirichlet_cw(__func__, csys);
280 #endif
281 } /* Border cell */
282
283 /* Special case to handle if there is an enforcement by penalization or
284 * algebraic. This situation may happen with a tetrahedron with one vertex
285 * or an edge lying on the boundary (but no face)
286 */
287
288 if (cb->cell_flag == CS_FLAG_BOUNDARY_CELL_BY_VERTEX) {
289
290 assert(vtx_bc_flag != NULL);
291
292 for (int v = 0; v < cm->n_vc; v++) {
293
294 for (int k = 0; k < 3; k++)
295 csys->dof_flag[3*v+k] = vtx_bc_flag[cm->v_ids[v]];
296 if (cs_cdo_bc_is_dirichlet(vtx_bc_flag[cm->v_ids[v]])) {
297 csys->has_dirichlet = true;
298 const cs_real_t *bc_val = eqb->dir_values + 3*cm->v_ids[v];
299 for (int k = 0; k < 3; k++)
300 csys->dir_values[3*v+k] = bc_val[k];
301 }
302
303 }
304
305 }
306
307 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 2
308 if (cs_dbg_cw_test(eqp, cm, csys)) cs_cell_mesh_dump(cm);
309 #endif
310 }
311
312 /*----------------------------------------------------------------------------*/
313 /*!
314 * \brief Build the local matrices arising from the diffusion, advection,
315 * reaction terms. If asked, a mass matrix is also computed and
316 * stored in mass_hodge->matrix
317 * Case of vector-valued CDO-Vb schemes
318 *
319 * \param[in] eqp pointer to a cs_equation_param_t structure
320 * \param[in] eqb pointer to a cs_equation_builder_t structure
321 * \param[in] eqc context for this kind of discretization
322 * \param[in] cm pointer to a cellwise view of the mesh
323 * \param[in, out] mass_hodge pointer to a cs_hodge_t structure (mass matrix)
324 * \param[in, out] diff_hodge pointer to a cs_hodge_t structure (diffusion)
325 * \param[in, out] csys pointer to a cellwise view of the system
326 * \param[in, out] cb pointer to a cellwise builder
327 */
328 /*----------------------------------------------------------------------------*/
329
330 static void
_vvb_conv_diff_reac(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_cdovb_vecteq_t * eqc,const cs_cell_mesh_t * cm,cs_hodge_t * mass_hodge,cs_hodge_t * diff_hodge,cs_cell_sys_t * csys,cs_cell_builder_t * cb)331 _vvb_conv_diff_reac(const cs_equation_param_t *eqp,
332 const cs_equation_builder_t *eqb,
333 const cs_cdovb_vecteq_t *eqc,
334 const cs_cell_mesh_t *cm,
335 cs_hodge_t *mass_hodge,
336 cs_hodge_t *diff_hodge,
337 cs_cell_sys_t *csys,
338 cs_cell_builder_t *cb)
339 {
340 if (cs_equation_param_has_diffusion(eqp)) { /* DIFFUSION TERM
341 * ============== */
342 assert(diff_hodge != NULL);
343
344 /* Set the diffusion property */
345
346 if (!(eqb->diff_pty_uniform))
347 cs_hodge_set_property_value_cw(cm, cb->t_pty_eval, cb->cell_flag,
348 diff_hodge);
349
350 /* Define the local stiffness matrix: local matrix owned by the cellwise
351 builder (store in cb->loc) */
352
353 eqc->get_stiffness_matrix(cm, diff_hodge, cb);
354
355 /* Add the local diffusion operator to the local system */
356
357 const cs_real_t *sval = cb->loc->val;
358 for (int bi = 0; bi < cm->n_vc; bi++) {
359 for (int bj = 0; bj < cm->n_vc; bj++) {
360
361 /* Retrieve the 3x3 matrix */
362
363 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bi, bj);
364 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
365
366 const cs_real_t _val = sval[cm->n_vc*bi+bj];
367 bij->val[0] += _val;
368 bij->val[4] += _val;
369 bij->val[8] += _val;
370
371 }
372 }
373
374 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 1
375 if (cs_dbg_cw_test(eqp, cm, csys))
376 cs_cell_sys_dump("\n>> Cell system after adding diffusion", csys);
377 #endif
378 }
379
380 if (eqb->sys_flag & CS_FLAG_SYS_MASS_MATRIX) { /* MASS MATRIX
381 * =========== */
382 assert(mass_hodge != NULL);
383
384 /* Build the mass matrix and store it in mass_hodge->matrix */
385
386 eqc->get_mass_matrix(cm, mass_hodge, cb);
387
388 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 1
389 if (cs_dbg_cw_test(eqp, cm, csys)) {
390 cs_log_printf(CS_LOG_DEFAULT, ">> Cell mass matrix");
391 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, mass_hodge->matrix);
392 }
393 #endif
394 }
395
396 if (cs_equation_param_has_reaction(eqp)) { /* REACTION TERM
397 * ============= */
398
399 /* Update the value of the reaction property(ies) if needed */
400
401 cs_equation_set_reaction_properties_cw(eqp, eqb, cm, cb);
402
403 if (eqb->sys_flag & CS_FLAG_SYS_REAC_DIAG) {
404
405 /* |c|*wvc = |dual_cell(v) cap c| */
406
407 assert(cs_eflag_test(eqb->msh_flag, CS_FLAG_COMP_PVQ));
408 const double ptyc = cb->rpty_val * cm->vol_c;
409
410 /* Only the diagonal block and its diagonal entries are modified */
411
412 for (int bi = 0; bi < cm->n_vc; bi++) {
413
414 const double vpty = cm->wvc[bi] * ptyc;
415
416 /* Retrieve the 3x3 matrix */
417
418 cs_sdm_t *bii = cs_sdm_get_block(csys->mat, bi, bi);
419 assert(bii->n_rows == 3 && bii->n_cols == 3);
420
421 bii->val[0] += vpty;
422 bii->val[4] += vpty;
423 bii->val[8] += vpty;
424
425 }
426
427 }
428 else {
429
430 assert(cs_flag_test(eqb->sys_flag, CS_FLAG_SYS_MASS_MATRIX));
431
432 /* Add the local mass matrix to the local system */
433
434 const cs_real_t *mval = mass_hodge->matrix->val;
435 for (int bi = 0; bi < cm->n_vc; bi++) {
436 for (int bj = 0; bj < cm->n_vc; bj++) {
437
438 /* Retrieve the 3x3 matrix */
439
440 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bi, bj);
441 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
442
443 const cs_real_t _val = cb->rpty_val * mval[cm->n_vc*bi+bj];
444 bij->val[0] += _val;
445 bij->val[4] += _val;
446 bij->val[8] += _val;
447
448 }
449 }
450
451 } /* Lumping or not lumping */
452
453 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 1
454 if (cs_dbg_cw_test(eqp, cm, csys))
455 cs_cell_sys_dump("\n>> Cell system after adding reaction", csys);
456 #endif
457 } /* Reaction term */
458 }
459
460 /*----------------------------------------------------------------------------*/
461 /*!
462 * \brief First pass to apply boundary conditions enforced weakly. Update
463 * the local system before applying the time scheme.
464 * Case of vector-valued CDO-Vb schemes
465 *
466 * \param[in] eqp pointer to a cs_equation_param_t structure
467 * \param[in] eqc context for this kind of discretization
468 * \param[in] cm pointer to a cellwise view of the mesh
469 * \param[in, out] fm pointer to a facewise view of the mesh
470 * \param[in, out] diff_hodge pointer to a discrete Hodge op. for diffusion
471 * \param[in, out] csys pointer to a cellwise view of the system
472 * \param[in, out] cb pointer to a cellwise builder
473 */
474 /*----------------------------------------------------------------------------*/
475
476 static void
_vvb_apply_weak_bc(const cs_equation_param_t * eqp,const cs_cdovb_vecteq_t * eqc,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * diff_hodge,cs_cell_sys_t * csys,cs_cell_builder_t * cb)477 _vvb_apply_weak_bc(const cs_equation_param_t *eqp,
478 const cs_cdovb_vecteq_t *eqc,
479 const cs_cell_mesh_t *cm,
480 cs_face_mesh_t *fm,
481 cs_hodge_t *diff_hodge,
482 cs_cell_sys_t *csys,
483 cs_cell_builder_t *cb)
484 {
485 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
486
487 /* Neumann boundary conditions:
488 * The common practice is to define Phi_neu = - lambda * \grad(u) . n_fc
489 * An outward flux is a positive flux whereas an inward flux is negative
490 * The minus just above implies the minus just below */
491
492 if (csys->has_nhmg_neumann) {
493 for (short int i = 0; i < 3*cm->n_vc; i++)
494 csys->rhs[i] -= csys->neu_values[i];
495 }
496
497 /* The enforcement of the Dirichlet has to be done after all
498 other contributions */
499
500 if (cs_equation_param_has_diffusion(eqp)) {
501
502 if (csys->has_dirichlet) /* csys is updated inside (matrix and rhs) */
503 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
504 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM)
505 eqc->enforce_dirichlet(eqp, cm, fm, diff_hodge, cb, csys);
506
507 }
508
509 if (csys->has_sliding)
510 eqc->enforce_sliding(eqp, cm, fm, diff_hodge, cb, csys);
511
512 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 1
513 if (cs_dbg_cw_test(eqp, cm, csys))
514 cs_cell_sys_dump("\n>> Cell system after BC weak treatment", csys);
515 #endif
516 } /* Cell with at least one boundary face */
517 }
518
519 /*----------------------------------------------------------------------------*/
520 /*!
521 * \brief Second pass to apply boundary conditions. Only Dirichlet BCs which
522 * are enforced strongly. Apply also the enforcement of internal DoFs.
523 * Update the local system after applying the time scheme.
524 * Case of vector-valued CDO-Vb schemes
525 *
526 * \param[in] eqp pointer to a cs_equation_param_t structure
527 * \param[in] eqb pointer to a cs_equation_builder_t structure
528 * \param[in] eqc context for this kind of discretization
529 * \param[in] cm pointer to a cellwise view of the mesh
530 * \param[in, out] fm pointer to a facewise view of the mesh
531 * \param[in, out] diff_hodge pointer to a discrete Hodge op. for diffusion
532 * \param[in, out] csys pointer to a cellwise view of the system
533 * \param[in, out] cb pointer to a cellwise builder
534 */
535 /*----------------------------------------------------------------------------*/
536
537 static void
_vvb_enforce_values(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_cdovb_vecteq_t * eqc,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * diff_hodge,cs_cell_sys_t * csys,cs_cell_builder_t * cb)538 _vvb_enforce_values(const cs_equation_param_t *eqp,
539 const cs_equation_builder_t *eqb,
540 const cs_cdovb_vecteq_t *eqc,
541 const cs_cell_mesh_t *cm,
542 cs_face_mesh_t *fm,
543 cs_hodge_t *diff_hodge,
544 cs_cell_sys_t *csys,
545 cs_cell_builder_t *cb)
546 {
547 if (cs_cell_has_boundary_elements(cb) && csys->has_dirichlet) {
548
549 /* Boundary element (through either vertices or faces) */
550
551 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_ALGEBRAIC ||
552 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED) {
553
554 /* csys is updated inside (matrix and rhs) */
555
556 eqc->enforce_dirichlet(eqp, cm, fm, diff_hodge, cb, csys);
557
558 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 2
559 if (cs_dbg_cw_test(eqp, cm, csys))
560 cs_cell_sys_dump("\n>> Cell system after strong BC treatment", csys);
561 #endif
562 }
563 }
564
565 if (cs_equation_param_has_internal_enforcement(eqp)) {
566
567 /* Internal enforcement of DoFs: Update csys (matrix and rhs) */
568
569 cs_equation_enforced_internal_block_dofs(eqb, cb, csys);
570
571 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 2
572 if (cs_dbg_cw_test(eqp, cm, csys))
573 cs_cell_sys_dump("\n>> Cell system after the internal enforcement",
574 csys);
575 #endif
576 }
577 }
578
579 /*----------------------------------------------------------------------------*/
580 /*!
581 * \brief Compute the residual normalization at the cellwise level according
582 * to the requested type of renormalization
583 * Case of vector-valued system.
584 *
585 * \param[in] type type of renormalization
586 * \param[in] cm pointer to a cs_cell_mesh_t structure
587 * \param[in] csys pointer to a cs_cell_sys_t structure
588 *
589 * \return the value of the cellwise contribution to the normalization of
590 * the residual
591 */
592 /*----------------------------------------------------------------------------*/
593
594 static double
_vvb_cw_rhs_normalization(cs_param_resnorm_type_t type,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys)595 _vvb_cw_rhs_normalization(cs_param_resnorm_type_t type,
596 const cs_cell_mesh_t *cm,
597 const cs_cell_sys_t *csys)
598 {
599 double _rhs_norm = 0;
600
601 switch (type) {
602
603 case CS_PARAM_RESNORM_WEIGHTED_RHS:
604 for (short int i = 0; i < cm->n_vc; i++) {
605 const cs_real_t w = cm->wvc[i];
606 const cs_real_t *_rhs = csys->rhs + 3*i;
607 _rhs_norm += w * (_rhs[0]*_rhs[0] + _rhs[1]*_rhs[1] + _rhs[2]*_rhs[2]);
608 }
609 _rhs_norm = cm->vol_c * _rhs_norm;
610 break;
611
612 case CS_PARAM_RESNORM_FILTERED_RHS:
613 for (short int i = 0; i < csys->n_dofs; i++) {
614 if (csys->dof_flag[i] & CS_CDO_BC_DIRICHLET)
615 continue;
616 else if (csys->dof_is_forced[i])
617 continue;
618 else
619 _rhs_norm += csys->rhs[i]*csys->rhs[i];
620 }
621 break;
622
623 case CS_PARAM_RESNORM_NORM2_RHS:
624 for (short int i = 0; i < csys->n_dofs; i++)
625 _rhs_norm += csys->rhs[i]*csys->rhs[i];
626 break;
627
628 default:
629 break; /* Nothing to do */
630
631 } /* Type of residual normalization */
632
633 return _rhs_norm;
634 }
635
636 /*! \endcond DOXYGEN_SHOULD_SKIP_THIS */
637
638 /*============================================================================
639 * Public function prototypes
640 *============================================================================*/
641
642 /*----------------------------------------------------------------------------*/
643 /*!
644 * \brief Check if the generic structures for building a CDO-Vb scheme are
645 * allocated
646 *
647 * \return true or false
648 */
649 /*----------------------------------------------------------------------------*/
650
651 bool
cs_cdovb_vecteq_is_initialized(void)652 cs_cdovb_vecteq_is_initialized(void)
653 {
654 if (_vvb_cell_system == NULL || _vvb_cell_builder == NULL)
655 return false;
656 else
657 return true;
658 }
659
660 /*----------------------------------------------------------------------------*/
661 /*!
662 * \brief Allocate work buffer and general structures related to CDO
663 * vector-valued vertex-based schemes
664 * Set shared pointers.
665 *
666 * \param[in] quant additional mesh quantities struct.
667 * \param[in] connect pointer to a cs_cdo_connect_t struct.
668 * \param[in] time_step pointer to a time step structure
669 * \param[in] ms pointer to a cs_matrix_structure_t structure
670 */
671 /*----------------------------------------------------------------------------*/
672
673 void
cs_cdovb_vecteq_init_common(const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,const cs_time_step_t * time_step,const cs_matrix_structure_t * ms)674 cs_cdovb_vecteq_init_common(const cs_cdo_quantities_t *quant,
675 const cs_cdo_connect_t *connect,
676 const cs_time_step_t *time_step,
677 const cs_matrix_structure_t *ms)
678 {
679 /* Assign static const pointers */
680
681 cs_shared_quant = quant;
682 cs_shared_connect = connect;
683 cs_shared_time_step = time_step;
684 cs_shared_ms = ms;
685
686 /* Structure used to build the final system by a cell-wise process */
687
688 assert(cs_glob_n_threads > 0);
689 BFT_MALLOC(_vvb_cell_system, cs_glob_n_threads, cs_cell_sys_t *);
690 BFT_MALLOC(_vvb_cell_builder, cs_glob_n_threads, cs_cell_builder_t *);
691
692 for (int i = 0; i < cs_glob_n_threads; i++) {
693 _vvb_cell_system[i] = NULL;
694 _vvb_cell_builder[i] = NULL;
695 }
696
697 const int n_max_dofs = 3*connect->n_max_vbyc;
698
699 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
700 #pragma omp parallel
701 {
702 int t_id = omp_get_thread_num();
703 assert(t_id < cs_glob_n_threads);
704
705 cs_cell_builder_t *cb = _vvb_create_cell_builder(connect);
706 _vvb_cell_builder[t_id] = cb;
707
708 int block_size = 3;
709 _vvb_cell_system[t_id] = cs_cell_sys_create(n_max_dofs,
710 connect->n_max_fbyc,
711 1,
712 &block_size);
713 }
714 #else
715 assert(cs_glob_n_threads == 1);
716
717 cs_cell_builder_t *cb = _vvb_create_cell_builder(connect);
718 _vvb_cell_builder[0] = cb;
719
720 int block_size = 3;
721 _vvb_cell_system[0] = cs_cell_sys_create(n_max_dofs,
722 connect->n_max_fbyc,
723 1,
724 &block_size);
725 #endif /* openMP */
726 }
727
728 /*----------------------------------------------------------------------------*/
729 /*!
730 * \brief Retrieve work buffers used for building a CDO system cellwise
731 *
732 * \param[out] csys pointer to a pointer on a cs_cell_sys_t structure
733 * \param[out] cb pointer to a pointer on a cs_cell_builder_t structure
734 */
735 /*----------------------------------------------------------------------------*/
736
737 void
cs_cdovb_vecteq_get(cs_cell_sys_t ** csys,cs_cell_builder_t ** cb)738 cs_cdovb_vecteq_get(cs_cell_sys_t **csys,
739 cs_cell_builder_t **cb)
740 {
741 int t_id = 0;
742
743 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
744 t_id = omp_get_thread_num();
745 assert(t_id < cs_glob_n_threads);
746 #endif /* openMP */
747
748 *csys = _vvb_cell_system[t_id];
749 *cb = _vvb_cell_builder[t_id];
750 }
751
752 /*----------------------------------------------------------------------------*/
753 /*!
754 * \brief Free work buffer and general structure related to CDO vertex-based
755 * schemes
756 */
757 /*----------------------------------------------------------------------------*/
758
759 void
cs_cdovb_vecteq_finalize_common(void)760 cs_cdovb_vecteq_finalize_common(void)
761 {
762 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
763 #pragma omp parallel
764 {
765 int t_id = omp_get_thread_num();
766 cs_cell_sys_free(&(_vvb_cell_system[t_id]));
767 cs_cell_builder_free(&(_vvb_cell_builder[t_id]));
768 }
769 #else
770 assert(cs_glob_n_threads == 1);
771 cs_cell_sys_free(&(_vvb_cell_system[0]));
772 cs_cell_builder_free(&(_vvb_cell_builder[0]));
773 #endif /* openMP */
774
775 BFT_FREE(_vvb_cell_system);
776 BFT_FREE(_vvb_cell_builder);
777 _vvb_cell_builder = NULL;
778 _vvb_cell_system = NULL;
779 }
780
781 /*----------------------------------------------------------------------------*/
782 /*!
783 * \brief Initialize a cs_cdovb_vecteq_t structure storing data useful for
784 * building and managing such a scheme
785 *
786 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
787 * \param[in] var_id id of the variable field
788 * \param[in] bflux_id id of the boundary flux field
789 * \param[in, out] eqb pointer to a \ref cs_equation_builder_t structure
790 *
791 * \return a pointer to a new allocated cs_cdovb_vecteq_t structure
792 */
793 /*----------------------------------------------------------------------------*/
794
795 void *
cs_cdovb_vecteq_init_context(const cs_equation_param_t * eqp,int var_id,int bflux_id,cs_equation_builder_t * eqb)796 cs_cdovb_vecteq_init_context(const cs_equation_param_t *eqp,
797 int var_id,
798 int bflux_id,
799 cs_equation_builder_t *eqb)
800 {
801 assert(eqp != NULL && eqb != NULL);
802
803 if (eqp->space_scheme != CS_SPACE_SCHEME_CDOVB || eqp->dim != 3)
804 bft_error(__FILE__, __LINE__, 0,
805 " %s: Invalid type of equation.\n"
806 " Expected: vector-valued CDO vertex-based equation.", __func__);
807
808 eqb->sys_flag = CS_FLAG_SYS_VECTOR;
809
810 const cs_cdo_connect_t *connect = cs_shared_connect;
811 const cs_lnum_t n_vertices = connect->n_vertices;
812
813 cs_cdovb_vecteq_t *eqc = NULL;
814
815 BFT_MALLOC(eqc, 1, cs_cdovb_vecteq_t);
816
817 eqc->var_field_id = var_id;
818 eqc->bflux_field_id = bflux_id;
819
820 eqc->n_dofs = 3*n_vertices;
821
822 /* Flag to indicate the minimal set of quantities to build in a cell mesh
823 According to the situation, additional flags have to be set */
824
825 eqb->msh_flag = CS_FLAG_COMP_PV | CS_FLAG_COMP_PVQ | CS_FLAG_COMP_PE |
826 CS_FLAG_COMP_EV;
827
828 /* Store additional flags useful for building boundary operator.
829 Only activated on boundary cells */
830
831 eqb->bd_msh_flag = CS_FLAG_COMP_PF | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
832 CS_FLAG_COMP_FEQ;
833
834 bool need_eigen =
835 (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
836 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM) ? true : false;
837
838 /* Diffusion term */
839
840 eqc->get_stiffness_matrix = NULL;
841 eqc->get_stiffness_matrix = NULL;
842
843 if (cs_equation_param_has_diffusion(eqp)) {
844
845 eqc->diffusion_hodge = cs_hodge_init_context(connect,
846 eqp->diffusion_property,
847 &(eqp->diffusion_hodgep),
848 true, /* tensor ? */
849 need_eigen); /* eigen ? */
850
851 const cs_property_data_t *diff_pty = eqc->diffusion_hodge[0]->pty_data;
852
853 if (diff_pty->is_iso == false)
854 bft_error(__FILE__, __LINE__, 0, " %s: Case not handle yet\n", __func__);
855
856 switch (eqp->diffusion_hodgep.algo) {
857
858 case CS_HODGE_ALGO_COST:
859 eqb->msh_flag |= CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ;
860 eqb->bd_msh_flag |= CS_FLAG_COMP_DEQ;
861 if (diff_pty->is_iso || diff_pty->is_unity)
862 eqc->get_stiffness_matrix = cs_hodge_vb_cost_get_iso_stiffness;
863 else
864 eqc->get_stiffness_matrix = cs_hodge_vb_cost_get_aniso_stiffness;
865 break;
866
867 case CS_HODGE_ALGO_BUBBLE:
868 eqb->msh_flag |= CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ;
869 eqb->bd_msh_flag |= CS_FLAG_COMP_DEQ;
870 if (diff_pty->is_iso || diff_pty->is_unity)
871 eqc->get_stiffness_matrix = cs_hodge_vb_bubble_get_iso_stiffness;
872 else
873 eqc->get_stiffness_matrix = cs_hodge_vb_bubble_get_aniso_stiffness;
874 break;
875
876 case CS_HODGE_ALGO_OCS2:
877 eqb->msh_flag |= CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ | CS_FLAG_COMP_SEF;
878 eqb->bd_msh_flag |= CS_FLAG_COMP_DEQ;
879 eqc->get_stiffness_matrix = cs_hodge_vb_ocs2_get_aniso_stiffness;
880 break;
881
882 case CS_HODGE_ALGO_VORONOI:
883 eqb->msh_flag |= CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ;
884 eqb->bd_msh_flag |= CS_FLAG_COMP_DEQ;
885 eqc->get_stiffness_matrix = cs_hodge_vb_voro_get_stiffness;
886 break;
887
888 case CS_HODGE_ALGO_WBS:
889 eqb->msh_flag |= CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PEQ |
890 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_HFQ | CS_FLAG_COMP_PFC;
891 eqc->get_stiffness_matrix = cs_hodge_vb_wbs_get_stiffness;
892 break;
893
894 default:
895 bft_error(__FILE__, __LINE__, 0,
896 " %s: Invalid type of algorithm to build the diffusion term.",
897 __func__);
898
899 } /* Switch on Hodge algo. */
900
901 } /* Has diffusion */
902
903 /* Boundary conditions */
904
905 BFT_MALLOC(eqc->vtx_bc_flag, n_vertices, cs_flag_t);
906 cs_equation_set_vertex_bc_flag(connect, eqb->face_bc, eqc->vtx_bc_flag);
907
908 eqc->enforce_robin_bc = NULL;
909 if (cs_equation_param_has_robin_bc(eqp))
910 bft_error(__FILE__, __LINE__, 0,
911 (" %s: Robin boundary conditions are not handled yet."),
912 __func__);
913
914 eqc->enforce_dirichlet = NULL;
915 switch (eqp->default_enforcement) {
916
917 case CS_PARAM_BC_ENFORCE_ALGEBRAIC:
918 eqc->enforce_dirichlet = cs_cdo_diffusion_alge_block_dirichlet;
919 break;
920 case CS_PARAM_BC_ENFORCE_PENALIZED:
921 eqc->enforce_dirichlet = cs_cdo_diffusion_pena_block_dirichlet;
922 break;
923
924 case CS_PARAM_BC_ENFORCE_WEAK_NITSCHE:
925 eqb->bd_msh_flag |= CS_FLAG_COMP_DEQ | CS_FLAG_COMP_HFQ;
926 eqc->enforce_dirichlet = cs_cdo_diffusion_vvb_ocs_weak_dirichlet;
927 break;
928
929 case CS_PARAM_BC_ENFORCE_WEAK_SYM:
930 default:
931 bft_error(__FILE__, __LINE__, 0,
932 " %s: Invalid type of algorithm to enforce Dirichlet BC.",
933 __func__);
934
935 }
936
937 eqc->enforce_sliding = NULL;
938 if (eqb->face_bc->n_sliding_faces > 0) {
939
940 /* There is at least one face with a sliding condition to handle */
941
942 eqb->bd_msh_flag |= CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_HFQ;
943 eqc->enforce_sliding = cs_cdo_diffusion_vvb_ocs_sliding;
944
945 }
946
947 /* Advection term */
948
949 eqc->get_advection_matrix = NULL;
950 eqc->add_advection_bc = NULL;
951
952 /* A mass matrix can be requested either for the reaction term, the unsteady
953 term or for the source term */
954
955 cs_hodge_algo_t reac_hodge_algo = CS_HODGE_N_ALGOS;
956 cs_hodge_algo_t time_hodge_algo = CS_HODGE_N_ALGOS;
957 cs_hodge_algo_t srct_hodge_algo = CS_HODGE_N_ALGOS;
958
959 /* Reaction term */
960
961 if (cs_equation_param_has_reaction(eqp)) {
962
963 if (eqp->do_lumping) {
964
965 eqb->sys_flag |= CS_FLAG_SYS_REAC_DIAG;
966 reac_hodge_algo = CS_HODGE_ALGO_VORONOI;
967
968 }
969 else {
970
971 switch (eqp->reaction_hodgep.algo) {
972
973 case CS_HODGE_ALGO_VORONOI:
974 eqb->sys_flag |= CS_FLAG_SYS_REAC_DIAG;
975 reac_hodge_algo = CS_HODGE_ALGO_VORONOI;
976 break;
977 case CS_HODGE_ALGO_WBS:
978 eqb->sys_flag |= CS_FLAG_SYS_MASS_MATRIX;
979 reac_hodge_algo = CS_HODGE_ALGO_WBS;
980 break;
981 default:
982 bft_error(__FILE__, __LINE__, 0,
983 "%s: Invalid choice of algorithm for the reaction term.",
984 __func__);
985 break;
986 }
987
988 } /* Lumping or not lumping */
989
990 } /* Reaction term is requested */
991
992 /* Unsteady term */
993
994 if (cs_equation_param_has_time(eqp)) {
995
996 if (eqp->do_lumping) {
997
998 eqb->sys_flag |= CS_FLAG_SYS_TIME_DIAG;
999 time_hodge_algo = CS_HODGE_ALGO_VORONOI;
1000
1001 }
1002 else {
1003
1004 switch (eqp->time_hodgep.algo) {
1005
1006 case CS_HODGE_ALGO_VORONOI:
1007 eqb->sys_flag |= CS_FLAG_SYS_TIME_DIAG;
1008 time_hodge_algo = CS_HODGE_ALGO_VORONOI;
1009 break;
1010 case CS_HODGE_ALGO_WBS:
1011 eqb->sys_flag |= CS_FLAG_SYS_MASS_MATRIX;
1012 time_hodge_algo = CS_HODGE_ALGO_WBS;
1013 break;
1014 default:
1015 bft_error(__FILE__, __LINE__, 0,
1016 "%s: Invalid choice of algorithm for the unsteady term.",
1017 __func__);
1018 break;
1019 }
1020
1021 } /* Lumping or not lumping */
1022
1023 } /* Unsteady term is requested */
1024
1025 /* Source term */
1026
1027 eqc->source_terms = NULL;
1028
1029 if (cs_equation_param_has_sourceterm(eqp)) {
1030
1031 if (cs_equation_param_has_time(eqp)) {
1032 if (eqp->time_scheme == CS_TIME_SCHEME_THETA ||
1033 eqp->time_scheme == CS_TIME_SCHEME_CRANKNICO) {
1034
1035 BFT_MALLOC(eqc->source_terms, eqc->n_dofs, cs_real_t);
1036 memset(eqc->source_terms, 0, eqc->n_dofs*sizeof(cs_real_t));
1037
1038 } /* Theta scheme */
1039 } /* Time-dependent equation */
1040
1041 /* Need a mass matrix */
1042
1043 for (int st_id = 0; st_id < eqp->n_source_terms; st_id++) {
1044
1045 cs_xdef_t *st = eqp->source_terms[st_id];
1046 if (st->meta & CS_FLAG_PRIMAL) {
1047 eqb->sys_flag |= CS_FLAG_SYS_MASS_MATRIX;
1048 srct_hodge_algo = CS_HODGE_ALGO_WBS;
1049 }
1050
1051 }
1052
1053 } /* There is at least one source term */
1054
1055 /* Pre-defined a cs_hodge_builder_t structure */
1056
1057 eqc->mass_hodgep.inv_pty = false;
1058 eqc->mass_hodgep.coef = 1.0; /* not useful in this case */
1059 eqc->mass_hodgep.type = CS_HODGE_TYPE_VPCD;
1060
1061 eqc->mass_hodgep.algo = cs_hodge_set_mass_algo(eqp->name,
1062 reac_hodge_algo,
1063 time_hodge_algo,
1064 srct_hodge_algo);
1065
1066 if (eqc->mass_hodgep.algo == CS_HODGE_ALGO_WBS)
1067 eqb->msh_flag |= CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PEQ
1068 | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_PFC;
1069
1070 /* Initialize the hodge structure for the mass matrix */
1071
1072 eqc->mass_hodge = cs_hodge_init_context(connect,
1073 NULL,
1074 &(eqc->mass_hodgep),
1075 false, /* tensor ? */
1076 false); /* eigen ? */
1077
1078 if (eqp->verbosity > 1 && eqb->sys_flag & CS_FLAG_SYS_MASS_MATRIX) {
1079 cs_log_printf(CS_LOG_SETUP,
1080 "#### Parameters of the mass matrix of the equation %s\n",
1081 eqp->name);
1082 cs_hodge_param_log("Mass matrix", NULL, eqc->mass_hodgep);
1083 }
1084
1085 /* Set the function pointer */
1086
1087 eqc->get_mass_matrix = cs_hodge_get_func(__func__, eqc->mass_hodgep);
1088
1089 /* Assembly process */
1090
1091 eqc->assemble = cs_equation_assemble_set(CS_SPACE_SCHEME_CDOVB,
1092 CS_CDO_CONNECT_VTX_VECT);
1093
1094 /* Array used for extra-operations */
1095
1096 eqc->cell_values = NULL;
1097
1098 return eqc;
1099 }
1100
1101 /*----------------------------------------------------------------------------*/
1102 /*!
1103 * \brief Destroy a cs_cdovb_vecteq_t structure
1104 *
1105 * \param[in, out] builder pointer to a cs_cdovb_vecteq_t structure
1106 *
1107 * \return a NULL pointer
1108 */
1109 /*----------------------------------------------------------------------------*/
1110
1111 void *
cs_cdovb_vecteq_free_context(void * builder)1112 cs_cdovb_vecteq_free_context(void *builder)
1113 {
1114 cs_cdovb_vecteq_t *eqc = (cs_cdovb_vecteq_t *)builder;
1115
1116 if (eqc == NULL)
1117 return eqc;
1118
1119 BFT_FREE(eqc->source_terms);
1120 BFT_FREE(eqc->cell_values);
1121 BFT_FREE(eqc->vtx_bc_flag);
1122
1123 cs_hodge_free_context(&(eqc->diffusion_hodge));
1124 cs_hodge_free_context(&(eqc->mass_hodge));
1125
1126 /* Last free */
1127
1128 BFT_FREE(eqc);
1129
1130 return NULL;
1131 }
1132
1133 /*----------------------------------------------------------------------------*/
1134 /*!
1135 * \brief Set the initial values of the variable field taking into account
1136 * the boundary conditions.
1137 * Case of vector-valued CDO-Vb schemes.
1138 *
1139 * \param[in] t_eval time at which one evaluates BCs
1140 * \param[in] field_id id related to the variable field of this equation
1141 * \param[in] mesh pointer to a cs_mesh_t structure
1142 * \param[in] eqp pointer to a cs_equation_param_t structure
1143 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1144 * \param[in, out] context pointer to the scheme context (cast on-the-fly)
1145 */
1146 /*----------------------------------------------------------------------------*/
1147
1148 void
cs_cdovb_vecteq_init_values(cs_real_t t_eval,const int field_id,const cs_mesh_t * mesh,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)1149 cs_cdovb_vecteq_init_values(cs_real_t t_eval,
1150 const int field_id,
1151 const cs_mesh_t *mesh,
1152 const cs_equation_param_t *eqp,
1153 cs_equation_builder_t *eqb,
1154 void *context)
1155 {
1156 const cs_cdo_quantities_t *quant = cs_shared_quant;
1157 const cs_cdo_connect_t *connect = cs_shared_connect;
1158
1159 cs_cdovb_vecteq_t *eqc = (cs_cdovb_vecteq_t *)context;
1160 cs_field_t *fld = cs_field_by_id(field_id);
1161 cs_real_t *v_vals = fld->val;
1162
1163 /* By default, 0 is set as initial condition for the computational domain.
1164
1165 Warning: This operation has to be done after the settings of the
1166 Dirichlet boundary conditions where an interface sum is performed
1167 for vertex-based schemes
1168 */
1169
1170 memset(v_vals, 0, 3*quant->n_vertices*sizeof(cs_real_t));
1171
1172 if (eqp->n_ic_defs > 0) {
1173
1174 cs_lnum_t *def2v_ids = (cs_lnum_t *)cs_equation_get_tmpbuf();
1175 cs_lnum_t *def2v_idx = NULL;
1176 BFT_MALLOC(def2v_idx, eqp->n_ic_defs + 1, cs_lnum_t);
1177
1178 cs_equation_sync_vol_def_at_vertices(connect,
1179 eqp->n_ic_defs,
1180 eqp->ic_defs,
1181 def2v_idx,
1182 def2v_ids);
1183
1184 /* Initialize values at mesh vertices */
1185
1186 cs_flag_t dof_flag = CS_FLAG_VECTOR | cs_flag_primal_vtx;
1187
1188 for (int def_id = 0; def_id < eqp->n_ic_defs; def_id++) {
1189
1190 /* Get and then set the definition of the initial condition */
1191
1192 const cs_xdef_t *def = eqp->ic_defs[def_id];
1193 const cs_lnum_t n_v_selected = def2v_idx[def_id+1] - def2v_idx[def_id];
1194 const cs_lnum_t *selected_lst = def2v_ids + def2v_idx[def_id];
1195
1196 switch(def->type) {
1197
1198 case CS_XDEF_BY_VALUE:
1199 cs_evaluate_potential_at_vertices_by_value(def,
1200 n_v_selected,
1201 selected_lst,
1202 v_vals);
1203 break;
1204
1205 case CS_XDEF_BY_QOV:
1206 /* Synchronization is performed inside */
1207 cs_evaluate_potential_by_qov(dof_flag, def, v_vals, NULL);
1208 break;
1209
1210 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1211 assert(eqp->dof_reduction == CS_PARAM_REDUCTION_DERHAM);
1212 cs_evaluate_potential_at_vertices_by_analytic(def,
1213 t_eval,
1214 n_v_selected,
1215 selected_lst,
1216 v_vals);
1217 break;
1218
1219 default:
1220 bft_error(__FILE__, __LINE__, 0,
1221 " %s: Invalid way to initialize field values for eq. %s.\n",
1222 __func__, eqp->name);
1223
1224 } /* Switch on possible type of definition */
1225
1226 } /* Loop on definitions */
1227
1228 BFT_FREE(def2v_idx);
1229
1230 } /* Initial values to set */
1231
1232 /* Set the boundary values as initial values: Compute the values of the
1233 Dirichlet BC */
1234
1235 cs_equation_compute_dirichlet_vb(t_eval,
1236 mesh,
1237 quant,
1238 connect,
1239 eqp,
1240 eqb->face_bc,
1241 _vvb_cell_builder[0], /* static variable */
1242 eqc->vtx_bc_flag,
1243 v_vals);
1244 }
1245
1246 /*----------------------------------------------------------------------------*/
1247 /*!
1248 * \brief Build and solve the linear system arising from a vector steady-state
1249 * convection/diffusion/reaction equation with a CDO-Vb scheme.
1250 * One works cellwise and then process to the assembly.
1251 *
1252 * \param[in] cur2prev true="current to previous" operation is performed
1253 * \param[in] mesh pointer to a cs_mesh_t structure
1254 * \param[in] field_id id of the variable field related to this equation
1255 * \param[in] eqp pointer to a cs_equation_param_t structure
1256 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1257 * \param[in, out] context pointer to cs_cdovb_vecteq_t structure
1258 */
1259 /*----------------------------------------------------------------------------*/
1260
1261 void
cs_cdovb_vecteq_solve_steady_state(bool cur2prev,const cs_mesh_t * mesh,const int field_id,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)1262 cs_cdovb_vecteq_solve_steady_state(bool cur2prev,
1263 const cs_mesh_t *mesh,
1264 const int field_id,
1265 const cs_equation_param_t *eqp,
1266 cs_equation_builder_t *eqb,
1267 void *context)
1268 {
1269 cs_timer_t t0 = cs_timer_time();
1270
1271 const cs_cdo_connect_t *connect = cs_shared_connect;
1272 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_VTX_VECT];
1273 const cs_cdo_quantities_t *quant = cs_shared_quant;
1274 const cs_time_step_t *ts = cs_shared_time_step;
1275 const cs_real_t time_eval = ts->t_cur + ts->dt[0];
1276
1277 cs_cdovb_vecteq_t *eqc = (cs_cdovb_vecteq_t *)context;
1278 cs_field_t *fld = cs_field_by_id(field_id);
1279
1280 /* Build an array storing the Dirichlet values at vertices
1281 * First argument is set to t_cur + dt_cur even if this is a steady
1282 * computation since one can call this function to compute a steady-state
1283 * solution at each time step of an unsteady computation.
1284 */
1285
1286 _vvb_setup(time_eval, mesh, eqp, eqb, eqc->vtx_bc_flag);
1287
1288 /* Initialize the local system: matrix and rhs */
1289
1290 cs_matrix_t *matrix = cs_matrix_create(cs_shared_ms);
1291 cs_real_t *rhs = NULL;
1292 double rhs_norm = 0.0;
1293
1294 assert(3*quant->n_vertices == eqc->n_dofs);
1295 BFT_MALLOC(rhs, eqc->n_dofs, cs_real_t);
1296 # pragma omp parallel for if (quant->n_vertices > CS_THR_MIN)
1297 for (cs_lnum_t i = 0; i < eqc->n_dofs; i++) rhs[i] = 0.0;
1298
1299 /* Initialize the structure to assemble values */
1300
1301 cs_matrix_assembler_values_t *mav
1302 = cs_matrix_assembler_values_init(matrix, NULL, NULL);
1303
1304 /* ------------------------- */
1305 /* Main OpenMP block on cell */
1306 /* ------------------------- */
1307
1308 # pragma omp parallel if (quant->n_cells > CS_THR_MIN)
1309 {
1310 /* Set variables and structures inside the OMP section so that each thread
1311 has its own value */
1312
1313 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
1314 int t_id = omp_get_thread_num();
1315 #else
1316 int t_id = 0;
1317 #endif
1318
1319 /* Each thread get back its related structures:
1320 Get the cell-wise view of the mesh and the algebraic system */
1321
1322 cs_face_mesh_t *fm = cs_cdo_local_get_face_mesh(t_id);
1323 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
1324 cs_cell_sys_t *csys = _vvb_cell_system[t_id];
1325 cs_cell_builder_t *cb = _vvb_cell_builder[t_id];
1326 cs_equation_assemble_t *eqa = cs_equation_assemble_get(t_id);
1327 cs_hodge_t *diff_hodge =
1328 (eqc->diffusion_hodge == NULL) ? NULL : eqc->diffusion_hodge[t_id];
1329 cs_hodge_t *mass_hodge =
1330 (eqc->mass_hodge == NULL) ? NULL : eqc->mass_hodge[t_id];
1331
1332 /* Set times at which one evaluates quantities if needed */
1333
1334 cb->t_pty_eval = time_eval;
1335 cb->t_bc_eval = time_eval;
1336 cb->t_st_eval = time_eval;
1337
1338 /* Initialization of the values of properties */
1339
1340 cs_equation_init_properties(eqp, eqb, diff_hodge, cb);
1341
1342 /* --------------------------------------------- */
1343 /* Main loop on cells to build the linear system */
1344 /* --------------------------------------------- */
1345
1346 # pragma omp for CS_CDO_OMP_SCHEDULE reduction(+:rhs_norm)
1347 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1348
1349 /* Set the current cell flag */
1350
1351 cb->cell_flag = connect->cell_flag[c_id];
1352
1353 /* Set the local mesh structure for the current cell */
1354
1355 cs_cell_mesh_build(c_id,
1356 cs_equation_cell_mesh_flag(cb->cell_flag, eqb),
1357 connect, quant, cm);
1358
1359 /* Set the local (i.e. cellwise) structures for the current cell */
1360
1361 _vvb_init_cell_system(cm, eqp, eqb, eqc->vtx_bc_flag, fld->val,
1362 csys, cb);
1363
1364 /* Build and add the diffusion/advection/reaction term to the local
1365 system. A mass matrix is also built if needed (mass_hodge->matrix) */
1366
1367 _vvb_conv_diff_reac(eqp, eqb, eqc, cm, mass_hodge, diff_hodge, csys, cb);
1368
1369 if (cs_equation_param_has_sourceterm(eqp)) { /* SOURCE TERM
1370 * =========== */
1371
1372 /* Reset the local contribution */
1373
1374 memset(csys->source, 0, csys->n_dofs*sizeof(cs_real_t));
1375
1376 /* Source term contribution to the algebraic system */
1377
1378 cs_source_term_compute_cellwise(eqp->n_source_terms,
1379 (cs_xdef_t *const *)eqp->source_terms,
1380 cm,
1381 eqb->source_mask,
1382 eqb->compute_source,
1383 cb->t_st_eval,
1384 mass_hodge,
1385 cb,
1386 csys->source);
1387
1388 /* Update the RHS */
1389
1390 for (int v = 0; v < csys->n_dofs; v++)
1391 csys->rhs[v] += csys->source[v];
1392
1393 } /* End of source term */
1394
1395 /* Apply boundary conditions (those which are weakly enforced) */
1396
1397 _vvb_apply_weak_bc(eqp, eqc, cm, fm, diff_hodge, csys, cb);
1398
1399 /* Enforce values if needed (internal or Dirichlet) */
1400
1401 _vvb_enforce_values(eqp, eqb, eqc, cm, fm, diff_hodge, csys, cb);
1402
1403 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_VECTEQ_DBG > 0
1404 if (cs_dbg_cw_test(eqp, cm, csys))
1405 cs_cell_sys_dump(">> (FINAL) Cell system matrix", csys);
1406 #endif
1407
1408 /* Compute a norm of the RHS for the normalization of the residual
1409 of the linear system to solve */
1410
1411 rhs_norm += _vvb_cw_rhs_normalization(eqp->sles_param->resnorm_type,
1412 cm, csys);
1413
1414 /* ASSEMBLY PROCESS
1415 * ================ */
1416
1417 eqc->assemble(csys->mat, csys->dof_ids, rs, eqa, mav);
1418
1419 # pragma omp critical
1420 {
1421 for (int i = 0; i < csys->n_dofs; i++) /* RHS Assembly */
1422 rhs[csys->dof_ids[i]] += csys->rhs[i];
1423 }
1424
1425 /* ********************** END OF ASSEMBLY PROCESS ******************* */
1426
1427 } /* Main loop on cells */
1428
1429 } /* OPENMP Block */
1430
1431 cs_matrix_assembler_values_done(mav); /* optional */
1432
1433 /* Free temporary buffers and structures */
1434
1435 cs_equation_builder_reset(eqb);
1436 cs_matrix_assembler_values_finalize(&mav);
1437
1438 /* Last step in the computation of the renormalization coefficient */
1439
1440 cs_equation_sync_rhs_normalization(eqp->sles_param->resnorm_type,
1441 eqc->n_dofs, /* 3*n_vertices */
1442 rhs,
1443 &rhs_norm);
1444
1445 /* Copy current field values to previous values */
1446
1447 if (cur2prev)
1448 cs_field_current_to_previous(fld);
1449
1450 /* End of the system building */
1451
1452 cs_timer_t t1 = cs_timer_time();
1453 cs_timer_counter_add_diff(&(eqb->tcb), &t0, &t1);
1454
1455 /* Solve the linear system (treated as a scalar-valued system
1456 with 3 times more DoFs) */
1457
1458 cs_sles_t *sles = cs_sles_find_or_add(eqp->sles_param->field_id, NULL);
1459
1460 cs_equation_solve_scalar_system(eqc->n_dofs, /* 3*n_vertices */
1461 eqp->sles_param,
1462 matrix,
1463 rs,
1464 rhs_norm,
1465 true, /* rhs_redux */
1466 sles,
1467 fld->val,
1468 rhs);
1469
1470 cs_timer_t t2 = cs_timer_time();
1471 cs_timer_counter_add_diff(&(eqb->tcs), &t1, &t2);
1472
1473 BFT_FREE(rhs);
1474 cs_sles_free(sles);
1475 cs_matrix_destroy(&matrix);
1476 }
1477
1478 /*----------------------------------------------------------------------------*/
1479 /*!
1480 * \brief Store solution(s) of the linear system into a field structure
1481 * Update extra-field values if required (for hybrid discretization)
1482 *
1483 * \param[in] solu solution array
1484 * \param[in] rhs rhs associated to this solution array
1485 * \param[in] eqp pointer to a cs_equation_param_t structure
1486 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1487 * \param[in, out] data pointer to cs_cdovb_vecteq_t structure
1488 * \param[in, out] field_val pointer to the current value of the field
1489 */
1490 /*----------------------------------------------------------------------------*/
1491
1492 void
cs_cdovb_vecteq_update_field(const cs_real_t * solu,const cs_real_t * rhs,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * data,cs_real_t * field_val)1493 cs_cdovb_vecteq_update_field(const cs_real_t *solu,
1494 const cs_real_t *rhs,
1495 const cs_equation_param_t *eqp,
1496 cs_equation_builder_t *eqb,
1497 void *data,
1498 cs_real_t *field_val)
1499 {
1500 CS_UNUSED(rhs);
1501 CS_UNUSED(eqp);
1502
1503 cs_cdovb_vecteq_t *eqc = (cs_cdovb_vecteq_t *)data;
1504 cs_timer_t t0 = cs_timer_time();
1505
1506 /* Set the computed solution in field array */
1507
1508 # pragma omp parallel for if (eqc->n_dofs > CS_THR_MIN)
1509 for (cs_lnum_t i = 0; i < eqc->n_dofs; i++)
1510 field_val[i] = solu[i];
1511
1512 cs_timer_t t1 = cs_timer_time();
1513 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
1514 }
1515
1516 /*----------------------------------------------------------------------------*/
1517 /*!
1518 * \brief Retrieve an array of values at mesh vertices for the variable field
1519 * associated to the given context
1520 * The lifecycle of this array is managed by the code. So one does not
1521 * have to free the return pointer.
1522 *
1523 * \param[in, out] context pointer to a data structure cast on-the-fly
1524 * \param[in] previous retrieve the previous state (true/false)
1525 *
1526 * \return a pointer to an array of cs_real_t (size: 3*n_vertices)
1527 */
1528 /*----------------------------------------------------------------------------*/
1529
1530 cs_real_t *
cs_cdovb_vecteq_get_vertex_values(void * context,bool previous)1531 cs_cdovb_vecteq_get_vertex_values(void *context,
1532 bool previous)
1533 {
1534 cs_cdovb_vecteq_t *eqc = (cs_cdovb_vecteq_t *)context;
1535
1536 if (eqc == NULL)
1537 return NULL;
1538
1539 cs_field_t *pot = cs_field_by_id(eqc->var_field_id);
1540
1541 if (previous)
1542 return pot->val_pre;
1543 else
1544 return pot->val;
1545 }
1546
1547 /*----------------------------------------------------------------------------*/
1548 /*!
1549 * \brief Compute an array of values at mesh cells by interpolating the
1550 * variable field associated to the given context located at mesh
1551 * vertices
1552 * The lifecycle of this array is managed by the code. So one does not
1553 * have to free the return pointer.
1554 *
1555 * \param[in, out] context pointer to a data structure cast on-the-fly
1556 * \param[in] previous retrieve the previous state (true/false)
1557 *
1558 * \return a pointer to an array of cs_real_t (size: 3*n_cells)
1559 */
1560 /*----------------------------------------------------------------------------*/
1561
1562 cs_real_t *
cs_cdovb_vecteq_get_cell_values(void * context,bool previous)1563 cs_cdovb_vecteq_get_cell_values(void *context,
1564 bool previous)
1565 {
1566 cs_cdovb_vecteq_t *eqc = (cs_cdovb_vecteq_t *)context;
1567
1568 if (eqc == NULL)
1569 return NULL;
1570
1571 const cs_cdo_quantities_t *quant = cs_shared_quant;
1572 const cs_cdo_connect_t *connect = cs_shared_connect;
1573
1574 cs_field_t *pot = cs_field_by_id(eqc->var_field_id);
1575
1576 cs_real_t *vtx_values = NULL;
1577 if (previous)
1578 vtx_values = pot->val_pre;
1579 else
1580 vtx_values = pot->val;
1581 assert(vtx_values != NULL);
1582
1583 /* Reset buffer of values */
1584
1585 if (eqc->cell_values == NULL)
1586 BFT_MALLOC(eqc->cell_values, 3*quant->n_cells, cs_real_t);
1587 memset(eqc->cell_values, 0, 3*quant->n_cells*sizeof(cs_real_t));
1588
1589 /* Compute the values at cell centers from an interpolation of the field
1590 values defined at vertices */
1591
1592 cs_reco_vect_pv_at_cell_centers(connect->c2v, quant, vtx_values,
1593 eqc->cell_values);
1594
1595 return eqc->cell_values;
1596 }
1597
1598 /*----------------------------------------------------------------------------*/
1599 /*!
1600 * \brief Operate a current to previous operation for the field associated to
1601 * this equation and potentially for related fields/arrays.
1602 *
1603 * \param[in] eqp pointer to a cs_equation_param_t structure
1604 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1605 * \param[in, out] context pointer to cs_cdovb_vecteq_t structure
1606 */
1607 /*----------------------------------------------------------------------------*/
1608
1609 void
cs_cdovb_vecteq_current_to_previous(const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)1610 cs_cdovb_vecteq_current_to_previous(const cs_equation_param_t *eqp,
1611 cs_equation_builder_t *eqb,
1612 void *context)
1613 {
1614 CS_UNUSED(eqp);
1615 CS_UNUSED(eqb);
1616
1617 cs_cdovb_vecteq_t *eqc = (cs_cdovb_vecteq_t *)context;
1618 cs_field_t *fld = cs_field_by_id(eqc->var_field_id);
1619
1620 cs_field_current_to_previous(fld);
1621 }
1622
1623 /*----------------------------------------------------------------------------*/
1624 /*!
1625 * \brief Predefined extra-operations related to this equation
1626 *
1627 * \param[in] eqp pointer to a cs_equation_param_t structure
1628 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1629 * \param[in, out] context pointer to cs_cdovb_vecteq_t structure
1630 */
1631 /*----------------------------------------------------------------------------*/
1632
1633 void
cs_cdovb_vecteq_extra_post(const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)1634 cs_cdovb_vecteq_extra_post(const cs_equation_param_t *eqp,
1635 cs_equation_builder_t *eqb,
1636 void *context)
1637 {
1638 CS_UNUSED(context);
1639 CS_UNUSED(eqp);
1640
1641 const cs_timer_t t0 = cs_timer_time();
1642
1643 /* TODO */
1644
1645 cs_timer_t t1 = cs_timer_time();
1646 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
1647 }
1648
1649 /*----------------------------------------------------------------------------*/
1650
1651 #undef _dp3
1652
1653 END_C_DECLS
1654