1 /*============================================================================
2 * Build an algebraic CDO vertex+cell-based system for unsteady convection
3 * diffusion reaction of scalar-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 <stdio.h>
35 #include <stdlib.h>
36 #include <math.h>
37 #include <assert.h>
38 #include <string.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include <bft_mem.h>
45
46 #include "cs_cdo_advection.h"
47 #include "cs_cdo_bc.h"
48 #include "cs_cdo_diffusion.h"
49 #include "cs_cdo_local.h"
50 #include "cs_defs.h"
51 #include "cs_equation_assemble.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_post.h"
60 #include "cs_quadrature.h"
61 #include "cs_reco.h"
62 #include "cs_scheme_geometry.h"
63 #include "cs_search.h"
64 #include "cs_sles.h"
65 #include "cs_source_term.h"
66 #include "cs_static_condensation.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_cdovcb_scaleq.h"
78
79 /*----------------------------------------------------------------------------*/
80
81 BEGIN_C_DECLS
82
83 /*!
84 \file cs_cdovcb_scaleq.c
85
86 \brief Build an algebraic CDO vertex+cell-based system for unsteady
87 convection diffusion reaction of scalar-valued equations with source
88 terms
89 */
90
91 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
92
93 /*=============================================================================
94 * Local Macro definitions and structure definitions
95 *============================================================================*/
96
97 #define CS_CDOVCB_SCALEQ_DBG 0
98
99 /* Redefined the name of functions from cs_math to get shorter names */
100
101 #define _dp3 cs_math_3_dot_product
102
103 /* Algebraic system for CDO vertex+cell-based discretization */
104
105 struct _cs_cdovcb_scaleq_t {
106
107 /* Ids related to the variable field and to the boundary flux field */
108
109 int var_field_id;
110 int bflux_field_id;
111
112 /* System size */
113 cs_lnum_t n_dofs; /* n_vertices + n_cells */
114
115 /* Store the values of the field at cell centers and the data needed to
116 compute the cell values from the vertex values. No need to synchronize
117 all these quantities since they are only cellwise quantities. */
118
119 cs_real_t *cell_values;
120 cs_real_t *cell_values_pre;
121
122 cs_real_t *cell_rhs; /* right-hand side related to cell dofs */
123
124 /* Assembly process */
125
126 cs_equation_assembly_t *assemble;
127
128 /* Members related to the static condensation */
129
130 cs_real_t *rc_tilda; /* Acc^-1 * RHS_cell */
131 cs_real_t *acv_tilda; /* Acc^-1 * Acv
132 Cell-vertices lower-Left block of the full matrix
133 Access to the values thanks to the c2v
134 connectivity */
135
136 /* Array storing the value of the contribution of all source terms */
137
138 cs_real_t *source_terms;
139
140 /* Boundary conditions */
141
142 cs_cdo_enforce_bc_t *enforce_dirichlet;
143 cs_cdo_enforce_bc_t *enforce_robin_bc;
144 cs_flag_t *vtx_bc_flag;
145
146 /* Pointers of function to build the diffusion term */
147
148 cs_hodge_t **diffusion_hodge;
149 cs_hodge_compute_t *get_stiffness_matrix;
150
151 /* Pointers of function to build the advection term */
152
153 cs_cdovb_advection_t *get_advection_matrix;
154 cs_cdovb_advection_bc_t *add_advection_bc;
155
156 /* If one needs to build a local hodge op. for time and reaction */
157
158 cs_hodge_param_t mass_hodgep;
159 cs_hodge_t **mass_hodge;
160 cs_hodge_compute_t *get_mass_matrix;
161
162 };
163
164 /*============================================================================
165 * Type definitions
166 *============================================================================*/
167
168 /* Algebraic system for CDO vertex-based discretization */
169
170 typedef struct _cs_cdovcb_scaleq_t cs_cdovcb_scaleq_t;
171
172 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
173
174 /*============================================================================
175 * Private variables
176 *============================================================================*/
177
178 /* Size = 1 if openMP is not used */
179
180 static cs_cell_sys_t **_vcbs_cell_system = NULL;
181 static cs_cell_builder_t **_vcbs_cell_builder = NULL;
182
183 /* Pointer to shared structures (owned by a cs_domain_t structure) */
184
185 static const cs_cdo_quantities_t *cs_shared_quant;
186 static const cs_cdo_connect_t *cs_shared_connect;
187 static const cs_time_step_t *cs_shared_time_step;
188 static const cs_matrix_structure_t *cs_shared_ms;
189
190 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
191
192 /*============================================================================
193 * Private function prototypes
194 *============================================================================*/
195
196 /*----------------------------------------------------------------------------*/
197 /*!
198 * \brief Initialize the local builder structure used for building the system
199 * cellwise
200 *
201 * \param[in] connect pointer to a cs_cdo_connect_t structure
202 *
203 * \return a pointer to a new allocated cs_cell_builder_t structure
204 */
205 /*----------------------------------------------------------------------------*/
206
207 static cs_cell_builder_t *
_cell_builder_create(const cs_cdo_connect_t * connect)208 _cell_builder_create(const cs_cdo_connect_t *connect)
209 {
210 const int n_vc = connect->n_max_vbyc;
211 const int n_ec = connect->n_max_ebyc;
212 const int n_fc = connect->n_max_fbyc;
213
214 cs_cell_builder_t *cb = cs_cell_builder_create();
215
216 int size = n_vc + 1;
217 BFT_MALLOC(cb->ids, size, int);
218 memset(cb->ids, 0, size*sizeof(int));
219
220 size = 2*n_vc + 3*n_ec + n_fc;
221 BFT_MALLOC(cb->values, size, double);
222 memset(cb->values, 0, size*sizeof(cs_real_t));
223
224 size = 2*n_ec + n_vc;
225 BFT_MALLOC(cb->vectors, size, cs_real_3_t);
226 memset(cb->vectors, 0, size*sizeof(cs_real_3_t));
227
228 /* Local square dense matrices used during the construction of
229 operators */
230
231 cb->loc = cs_sdm_square_create(n_vc + 1);
232 cb->aux = cs_sdm_square_create(n_vc + 1);
233
234 return cb;
235 }
236
237 /*----------------------------------------------------------------------------*/
238 /*!
239 * \brief Set the boundary conditions known from the settings
240 *
241 * \param[in] t_eval time at which one evaluates BCs
242 * \param[in] mesh pointer to a cs_mesh_t structure
243 * \param[in] connect pointer to a cs_cdo_connect_t structure
244 * \param[in] eqp pointer to a cs_equation_param_t structure
245 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
246 * \param[in, out] vtx_bc_flag pointer to an array of BC flag for each vertex
247 */
248 /*----------------------------------------------------------------------------*/
249
250 static void
_setup_vcb(cs_real_t t_eval,const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,cs_flag_t * vtx_bc_flag)251 _setup_vcb(cs_real_t t_eval,
252 const cs_mesh_t *mesh,
253 const cs_cdo_connect_t *connect,
254 const cs_equation_param_t *eqp,
255 cs_equation_builder_t *eqb,
256 cs_flag_t *vtx_bc_flag)
257 {
258 assert(vtx_bc_flag != NULL); /* Sanity check */
259 const cs_cdo_quantities_t *quant = cs_shared_quant;
260
261 /* Compute the values of the Dirichlet BC */
262
263 BFT_MALLOC(eqb->dir_values, quant->n_vertices, cs_real_t);
264
265 cs_equation_compute_dirichlet_vb(t_eval,
266 mesh,
267 quant,
268 connect,
269 eqp,
270 eqb->face_bc,
271 _vcbs_cell_builder[0], /* static variable */
272 vtx_bc_flag,
273 eqb->dir_values);
274
275 /* Internal enforcement of DoFs */
276
277 if (cs_equation_param_has_internal_enforcement(eqp))
278 eqb->enforced_values =
279 cs_enforcement_define_at_vertices(connect,
280 eqp->n_enforcements,
281 eqp->enforcement_params);
282 }
283
284 /*----------------------------------------------------------------------------*/
285 /*!
286 * \brief Initialize the local structure for the current cell
287 *
288 * \param[in] cm pointer to a cellwise view of the mesh
289 * \param[in] eqp pointer to a cs_equation_param_t structure
290 * \param[in] eqb pointer to a cs_equation_builder_t structure
291 * \param[in] eqc pointer to a cs_cdovcb_scaleq_t structure
292 * \param[in] vtx_bc_flag Flag related to BC associated to each vertex
293 * \param[in] field_tn values of the field at the last computed time
294 * \param[in, out] csys pointer to a cellwise view of the system
295 * \param[in, out] cb pointer to a cellwise builder
296 */
297 /*----------------------------------------------------------------------------*/
298
299 static void
_svcb_init_cell_system(const cs_cell_mesh_t * cm,const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_cdovcb_scaleq_t * eqc,const cs_flag_t vtx_bc_flag[],const cs_real_t field_tn[],cs_cell_sys_t * csys,cs_cell_builder_t * cb)300 _svcb_init_cell_system(const cs_cell_mesh_t *cm,
301 const cs_equation_param_t *eqp,
302 const cs_equation_builder_t *eqb,
303 const cs_cdovcb_scaleq_t *eqc,
304 const cs_flag_t vtx_bc_flag[],
305 const cs_real_t field_tn[],
306 cs_cell_sys_t *csys,
307 cs_cell_builder_t *cb)
308 {
309 /* Cell-wise view of the linear system to build */
310
311 const short int n_dofs = cm->n_vc + 1;
312
313 csys->c_id = cm->c_id;
314 csys->n_dofs = n_dofs;
315
316 /* Initialize the local system */
317
318 cs_cell_sys_reset(cm->n_fc, csys);
319
320 cs_sdm_square_init(n_dofs, csys->mat);
321
322 for (short int v = 0; v < cm->n_vc; v++) {
323 csys->dof_ids[v] = cm->v_ids[v];
324 csys->val_n[v] = field_tn[cm->v_ids[v]];
325 }
326 csys->dof_ids[cm->n_vc] = cm->c_id;
327 csys->val_n[cm->n_vc] = eqc->cell_values[cm->c_id];
328
329 /* Store the local values attached to Dirichlet values if the current cell
330 has at least one border face */
331
332 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
333
334 cs_equation_vb_set_cell_bc(cm,
335 eqp,
336 eqb->face_bc,
337 vtx_bc_flag,
338 eqb->dir_values,
339 cb->t_bc_eval,
340 csys,
341 cb);
342
343 #if defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
344 cs_dbg_check_hmg_dirichlet_cw(__func__, csys);
345 #endif
346 } /* Border cell */
347
348 /* Special case to handle if enforcement by penalization or algebraic
349 * This situation may happen with a tetrahedron with one vertex or an edge
350 * lying on the boundary (but no face)
351 */
352
353 if (cb->cell_flag == CS_FLAG_BOUNDARY_CELL_BY_VERTEX) {
354
355 assert(vtx_bc_flag != NULL);
356
357 for (short int v = 0; v < cm->n_vc; v++) {
358 csys->dof_flag[v] = vtx_bc_flag[cm->v_ids[v]];
359 if (cs_cdo_bc_is_dirichlet(csys->dof_flag[v])) {
360 csys->has_dirichlet = true;
361 csys->dir_values[v] = eqb->dir_values[cm->v_ids[v]];
362 }
363 }
364
365 }
366
367 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 2
368 if (cs_dbg_cw_test(eqp, cm, csys)) cs_cell_mesh_dump(cm);
369 #endif
370 }
371
372 /*----------------------------------------------------------------------------*/
373 /*!
374 * \brief Build the local matrices arising from the diffusion, advection,
375 * reaction terms in CDO-VCb schemes. If asked, a mass matrix is also
376 * computed and stored in mass_hodge->matrix
377 *
378 * \param[in] eqp pointer to a cs_equation_param_t structure
379 * \param[in] eqb pointer to a cs_equation_builder_t structure
380 * \param[in] eqc context for this kind of discretization
381 * \param[in] cm pointer to a cellwise view of the mesh
382 * \param[in, out] fm pointer to a facewise view of the mesh
383 * \param[in, out] mass_hodge pointer to a discrete Hodge op. (mass matrix)
384 * \param[in, out] diff_hodge pointer to a discrete Hodge op. for diffusion
385 * \param[in, out] csys pointer to a cellwise view of the system
386 * \param[in, out] cb pointer to a cellwise builder
387 */
388 /*----------------------------------------------------------------------------*/
389
390 static void
_svcb_conv_diff_reac(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_cdovcb_scaleq_t * eqc,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * mass_hodge,cs_hodge_t * diff_hodge,cs_cell_sys_t * csys,cs_cell_builder_t * cb)391 _svcb_conv_diff_reac(const cs_equation_param_t *eqp,
392 const cs_equation_builder_t *eqb,
393 const cs_cdovcb_scaleq_t *eqc,
394 const cs_cell_mesh_t *cm,
395 cs_face_mesh_t *fm,
396 cs_hodge_t *mass_hodge,
397 cs_hodge_t *diff_hodge,
398 cs_cell_sys_t *csys,
399 cs_cell_builder_t *cb)
400 {
401 if (cs_equation_param_has_diffusion(eqp)) { /* DIFFUSION TERM
402 * ============== */
403 assert(diff_hodge != NULL);
404
405 /* Set the diffusion property */
406
407 if (!(eqb->diff_pty_uniform))
408 cs_hodge_set_property_value_cw(cm, cb->t_pty_eval, cb->cell_flag,
409 diff_hodge);
410
411 /* Define the local stiffness matrix: local matrix owned by the cellwise
412 builder (store in cb->loc) */
413
414 eqc->get_stiffness_matrix(cm, diff_hodge, cb);
415
416 /* Add the local diffusion operator to the local system */
417
418 cs_sdm_add(csys->mat, cb->loc);
419
420 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
421 if (cs_dbg_cw_test(eqp, cm, csys))
422 cs_cell_sys_dump("\n>> Cell system after diffusion", csys);
423 #endif
424 }
425
426 if (cs_equation_param_has_convection(eqp) &&
427 ((cb->cell_flag & CS_FLAG_SOLID_CELL) == 0)) { /* ADVECTION TERM
428 * ============== */
429
430 cs_property_data_t *diff_pty =
431 (diff_hodge == NULL) ? NULL : diff_hodge->pty_data;
432
433 /* Define the local advection matrix */
434
435 eqc->get_advection_matrix(eqp, cm, diff_pty, fm, cb);
436
437 /* Add it to the local system */
438
439 if (eqp->adv_scaling_property == NULL)
440 cs_sdm_add(csys->mat, cb->loc);
441
442 else {
443
444 if (cs_property_is_uniform(eqp->adv_scaling_property))
445 cs_sdm_add_mult(csys->mat,
446 eqp->adv_scaling_property->ref_value, cb->loc);
447 else {
448 cs_real_t scaling = cs_property_value_in_cell(cm,
449 eqp->adv_scaling_property,
450 cb->t_pty_eval);
451 cs_sdm_add_mult(csys->mat, scaling, cb->loc);
452 }
453
454 }
455
456 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
457 if (cs_dbg_cw_test(eqp, cm, csys))
458 cs_cell_sys_dump("\n>> Cell system after advection", csys);
459 #endif
460 }
461
462 if (eqb->sys_flag & CS_FLAG_SYS_MASS_MATRIX) { /* MASS MATRIX
463 * =========== */
464 assert(mass_hodge != NULL);
465
466 /* Build the mass matrix and store it in mass_hodge->matrix */
467
468 eqc->get_mass_matrix(cm, mass_hodge, cb);
469
470 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
471 if (cs_dbg_cw_test(eqp, cm, csys)) {
472 cs_log_printf(CS_LOG_DEFAULT, ">> Cell mass matrix");
473 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, mass_hodge->matrix);
474 }
475 #endif
476 }
477
478 if (cs_equation_param_has_reaction(eqp)) { /* REACTION TERM
479 * ============= */
480
481 /* Update the value of the reaction property(ies) if needed */
482
483 cs_equation_set_reaction_properties_cw(eqp, eqb, cm, cb);
484
485 if (eqb->sys_flag & CS_FLAG_SYS_REAC_DIAG) {
486
487 /* |c|*wvc = |dual_cell(v) cap c| */
488
489 assert(cs_eflag_test(eqb->msh_flag, CS_FLAG_COMP_PVQ));
490 const double ptyc = cb->rpty_val * cm->vol_c;
491 for (short int i = 0; i < cm->n_vc; i++)
492 csys->mat->val[i*(cm->n_vc + 1)] += 0.75 * cm->wvc[i] * ptyc;
493
494 /* Cell DoF */
495
496 csys->mat->val[cm->n_vc*(cm->n_vc + 1)] += 0.25 * ptyc;
497
498 }
499 else {
500
501 assert(cs_flag_test(eqb->sys_flag, CS_FLAG_SYS_MASS_MATRIX));
502
503 /* Update local system matrix with the reaction term */
504
505 cs_sdm_add_mult(csys->mat, cb->rpty_val, mass_hodge->matrix);
506
507 }
508
509 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
510 if (cs_dbg_cw_test(eqp, cm, csys))
511 cs_cell_sys_dump("\n>> Cell system after reaction", csys);
512 #endif
513 }
514 }
515
516 /*----------------------------------------------------------------------------*/
517 /*!
518 * \brief First pass to apply boundary conditions enforced weakly in CDO-VCb
519 * schemes. Update the local system before applying the time scheme.
520 *
521 * \param[in] eqp pointer to a cs_equation_param_t structure
522 * \param[in] eqc context for this kind of discretization
523 * \param[in] cm pointer to a cellwise view of the mesh
524 * \param[in, out] fm pointer to a facewise view of the mesh
525 * \param[in, out] diff_hodge pointer to a discrete Hodge op. for diffusion
526 * \param[in, out] csys pointer to a cellwise view of the system
527 * \param[in, out] cb pointer to a cellwise builder
528 */
529 /*----------------------------------------------------------------------------*/
530
531 static void
_svcb_apply_weak_bc(const cs_equation_param_t * eqp,const cs_cdovcb_scaleq_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)532 _svcb_apply_weak_bc(const cs_equation_param_t *eqp,
533 const cs_cdovcb_scaleq_t *eqc,
534 const cs_cell_mesh_t *cm,
535 cs_face_mesh_t *fm,
536 cs_hodge_t *diff_hodge,
537 cs_cell_sys_t *csys,
538 cs_cell_builder_t *cb)
539 {
540 /* BOUNDARY CONDITION CONTRIBUTION TO THE ALGEBRAIC SYSTEM
541 * Operations that have to be performed BEFORE the static condensation */
542
543 if (cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
544
545 if (cs_equation_param_has_convection(eqp) &&
546 ((cb->cell_flag & CS_FLAG_SOLID_CELL) == 0))
547
548 /* Apply boundary conditions related to the advection term
549 csys is updated inside (matrix and rhs) */
550
551 eqc->add_advection_bc(cm, eqp, cb->t_bc_eval, fm, cb, csys);
552
553 /* Weakly enforced Dirichlet BCs for cells attached to the boundary
554 csys is updated inside (matrix and rhs) */
555
556 if (cs_equation_param_has_diffusion(eqp)) {
557
558 if (csys->has_robin) {
559 assert(eqc->enforce_robin_bc != NULL);
560 eqc->enforce_robin_bc(eqp, cm, fm, diff_hodge, cb, csys);
561 }
562
563 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
564 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM)
565 eqc->enforce_dirichlet(eqp, cm, fm, diff_hodge, cb, csys);
566
567 }
568
569 /* Neumann boundary conditions:
570 * The common practice is to define Phi_neu = - lambda * grad(u) . n_fc
571 * An outward flux is a positive flux whereas an inward flux is negative
572 * The minus just above implies the minus just below */
573
574 if (csys->has_nhmg_neumann) {
575 for (short int v = 0; v < cm->n_vc; v++)
576 csys->rhs[v] -= csys->neu_values[v];
577 }
578
579 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
580 if (cs_dbg_cw_test(eqp, cm, csys))
581 cs_cell_sys_dump(">> Cell system matrix after weak BC treatment", csys);
582 #endif
583 } /* Cell with at least one boundary face */
584 }
585
586 /*----------------------------------------------------------------------------*/
587 /*!
588 * \brief Second pass to apply boundary conditions. Only Dirichlet BCs which
589 * are enforced strongly. Apply also the enforcement of internal DoFs.
590 * Update the local system after applying the time scheme and the
591 * static condensation.
592 *
593 * \param[in] eqp pointer to a cs_equation_param_t structure
594 * \param[in] eqb pointer to a cs_equation_builder_t structure
595 * \param[in] eqc context for this kind of discretization
596 * \param[in] cm pointer to a cellwise view of the mesh
597 * \param[in, out] fm pointer to a facewise view of the mesh
598 * \param[in, out] diff_hodge pointer to a discrete Hodge op. for diffusion
599 * \param[in, out] csys pointer to a cellwise view of the system
600 * \param[in, out] cb pointer to a cellwise builder
601 */
602 /*----------------------------------------------------------------------------*/
603
604 static void
_svcb_enforce_values(const cs_equation_param_t * eqp,const cs_equation_builder_t * eqb,const cs_cdovcb_scaleq_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)605 _svcb_enforce_values(const cs_equation_param_t *eqp,
606 const cs_equation_builder_t *eqb,
607 const cs_cdovcb_scaleq_t *eqc,
608 const cs_cell_mesh_t *cm,
609 cs_face_mesh_t *fm,
610 cs_hodge_t *diff_hodge,
611 cs_cell_sys_t *csys,
612 cs_cell_builder_t *cb)
613 {
614 /* Internal enforcement of DoFs: Update csys (matrix and rhs) */
615
616 if (cs_equation_param_has_internal_enforcement(eqp)) {
617
618 cs_equation_enforced_internal_dofs(eqb, cb, csys);
619
620 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_SCALEQ_DBG > 2
621 if (cs_dbg_cw_test(eqp, cm, csys))
622 cs_cell_sys_dump("\n>> Cell system after the internal enforcement",
623 csys);
624 #endif
625 }
626
627 /* BOUNDARY CONDITION CONTRIBUTION TO THE ALGEBRAIC SYSTEM
628 * Operations that have to be performed AFTER the static condensation */
629
630 if (cs_cell_has_boundary_elements(cb) && csys->has_dirichlet) {
631
632 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED ||
633 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_ALGEBRAIC) {
634
635 /* Strongly enforced Dirichlet BCs for cells attached to the boundary
636 csys is updated inside (matrix and rhs) */
637
638 eqc->enforce_dirichlet(eqp, cm, fm, diff_hodge, cb, csys);
639
640 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 2
641 if (cs_dbg_cw_test(eqp, cm, csys))
642 cs_cell_sys_dump("\n>> Cell system after strong BC treatment", csys);
643 #endif
644 }
645
646 } /* Boundary cell */
647
648 }
649
650 /*----------------------------------------------------------------------------*/
651 /*!
652 * \brief Compute the residual normalization at the cellwise level according
653 * to the requested type of renormalization
654 * This function should be called after the static condensation.
655 * Case of CDO scalar-valued vertex+cell-based scheme.
656 *
657 * \param[in] type type of renormalization
658 * \param[in] cm pointer to a cs_cell_mesh_t structure
659 * \param[in] csys pointer to a cs_cell_sys_t structure
660 *
661 * \return the value of the cellwise contribution to the normalization of
662 * the residual
663 */
664 /*----------------------------------------------------------------------------*/
665
666 static double
_svcb_cw_rhs_normalization(cs_param_resnorm_type_t type,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys)667 _svcb_cw_rhs_normalization(cs_param_resnorm_type_t type,
668 const cs_cell_mesh_t *cm,
669 const cs_cell_sys_t *csys)
670 {
671 double _rhs_norm = 0;
672
673 switch (type) {
674
675 case CS_PARAM_RESNORM_WEIGHTED_RHS:
676 for (short int i = 0; i < cm->n_vc; i++)
677 _rhs_norm += cm->wvc[i] * csys->rhs[i]*csys->rhs[i];
678 _rhs_norm = _rhs_norm * cm->vol_c;
679 break;
680
681 case CS_PARAM_RESNORM_FILTERED_RHS:
682 for (short int i = 0; i < cm->n_vc; i++) {
683 if (csys->dof_flag[i] & CS_CDO_BC_DIRICHLET)
684 continue;
685 else if (csys->dof_is_forced[i])
686 continue;
687 else
688 _rhs_norm += csys->rhs[i]*csys->rhs[i];
689 }
690 break;
691
692 case CS_PARAM_RESNORM_NORM2_RHS:
693 for (short int i = 0; i < cm->n_vc; i++)
694 _rhs_norm += csys->rhs[i]*csys->rhs[i];
695 break;
696
697 default:
698 break; /* Nothing to do */
699
700 } /* Type of residual normalization */
701
702 return _rhs_norm;
703 }
704
705 /*----------------------------------------------------------------------------*/
706 /*!
707 * \brief Perform the assembly step
708 *
709 * \param[in] eqc context for this kind of discretization
710 * \param[in] cm pointer to a cellwise view of the mesh
711 * \param[in] csys pointer to a cellwise view of the system
712 * \param[in] rs pointer to a cs_range_set_t structure
713 * \param[in, out] eqa pointer to a cs_equation_assemble_t structure
714 * \param[in, out] mav pointer to a cs_matrix_assembler_values_t structure
715 * \param[in, out] rhs right-hand side array
716 */
717 /*----------------------------------------------------------------------------*/
718
719 inline static void
_assemble(const cs_cdovcb_scaleq_t * eqc,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,const cs_range_set_t * rs,cs_equation_assemble_t * eqa,cs_matrix_assembler_values_t * mav,cs_real_t * rhs)720 _assemble(const cs_cdovcb_scaleq_t *eqc,
721 const cs_cell_mesh_t *cm,
722 const cs_cell_sys_t *csys,
723 const cs_range_set_t *rs,
724 cs_equation_assemble_t *eqa,
725 cs_matrix_assembler_values_t *mav,
726 cs_real_t *rhs)
727 {
728 /* Matrix assembly */
729
730 eqc->assemble(csys->mat, csys->dof_ids, rs, eqa, mav);
731
732 /* RHS assembly */
733
734 #if CS_CDO_OMP_SYNC_SECTIONS > 0
735 # pragma omp critical
736 {
737 for (short int v = 0; v < cm->n_vc; v++)
738 rhs[cm->v_ids[v]] += csys->rhs[v];
739 }
740
741 if (eqc->source_terms != NULL) {
742
743 /* Source term assembly */
744
745 # pragma omp critical
746 {
747 for (short int v = 0; v < cm->n_vc; v++)
748 eqc->source_terms[cm->v_ids[v]] += csys->source[v];
749 }
750
751 }
752
753 #else /* Use atomic barrier */
754
755 for (short int v = 0; v < cm->n_vc; v++)
756 # pragma omp atomic
757 rhs[cm->v_ids[v]] += csys->rhs[v];
758
759 if (eqc->source_terms != NULL) {
760
761 /* Source term assembly */
762
763 for (int v = 0; v < cm->n_vc; v++)
764 # pragma omp atomic
765 eqc->source_terms[cm->v_ids[v]] += csys->source[v];
766
767 }
768
769 #endif
770
771 if (eqc->source_terms != NULL) {
772 cs_real_t *cell_sources = eqc->source_terms + cs_shared_quant->n_vertices;
773
774 cell_sources[cm->c_id] = csys->source[cm->n_vc];
775 }
776 }
777
778 /*----------------------------------------------------------------------------*/
779 /*!
780 * \brief Update the variables related to CDO-VCb system after a resolution
781 *
782 * \param[in, out] tce pointer to a timer counter
783 * \param[in, out] fld pointer to a cs_field_t structure
784 * \param[in, out] eqc pointer to a context structure
785 * \param[in] cur2prev true if one performs "current to previous" op.
786 */
787 /*----------------------------------------------------------------------------*/
788
789 static inline void
_update_cell_fields(cs_timer_counter_t * tce,cs_field_t * fld,cs_cdovcb_scaleq_t * eqc,bool cur2prev)790 _update_cell_fields(cs_timer_counter_t *tce,
791 cs_field_t *fld,
792 cs_cdovcb_scaleq_t *eqc,
793 bool cur2prev)
794 {
795 cs_timer_t t0 = cs_timer_time();
796
797 /* Copy current field values to previous values */
798
799 if (cur2prev && eqc->cell_values_pre != NULL) {
800 size_t bsize = cs_shared_connect->n_cells*sizeof(cs_real_t);
801 memcpy(eqc->cell_values_pre, eqc->cell_values, bsize);
802 }
803
804 /* Compute values at cells pc = acc^-1*(RHS - Acv*pv) */
805
806 cs_static_condensation_recover_scalar(cs_shared_connect->c2v,
807 eqc->rc_tilda,
808 eqc->acv_tilda,
809 fld->val,
810 eqc->cell_values);
811
812 cs_timer_t t1 = cs_timer_time();
813 cs_timer_counter_add_diff(tce, &t0, &t1);
814 }
815
816 /*----------------------------------------------------------------------------*/
817 /*!
818 * \brief Update the value of stabilization coefficient in given situation
819 *
820 * \param[in] eqp pointer to a cs_equation_param_t structure
821 */
822 /*----------------------------------------------------------------------------*/
823
824 static void
_set_cip_coef(const cs_equation_param_t * eqp)825 _set_cip_coef(const cs_equation_param_t *eqp)
826 {
827 const double gseed = 1e-2; /* Default value to multiply according to the
828 problem and the ratio of diameters */
829
830 const cs_cdo_quantities_t *quant = cs_shared_quant;
831 const double hc_max = quant->cell_info.h_max;
832 const double hc_min = quant->cell_info.h_min;
833 const double hf_max = quant->face_info.h_max;
834 const double hf_min = quant->face_info.h_min;
835 const double hcMm = hc_max * hc_min;
836 const double hfMm = hf_min * hf_max;
837 const double rho_fc = hcMm / hfMm;
838
839 double gamma = 10 * gseed * hc_max * hc_max * rho_fc;
840
841 /* If not pure convection */
842
843 if (cs_equation_param_has_diffusion(eqp) ||
844 cs_equation_param_has_reaction(eqp) ||
845 cs_equation_param_has_time(eqp))
846 gamma *= 0.1;
847
848 cs_cdo_advection_set_cip_coef(gamma);
849 }
850
851 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
852
853 /*============================================================================
854 * Public function prototypes
855 *============================================================================*/
856
857 /*----------------------------------------------------------------------------*/
858 /*!
859 * \brief Check if the generic structures for building a CDO-vertex+cell
860 * based scheme are allocated
861 *
862 * \return true or false
863 */
864 /*----------------------------------------------------------------------------*/
865
866 bool
cs_cdovcb_scaleq_is_initialized(void)867 cs_cdovcb_scaleq_is_initialized(void)
868 {
869 if (_vcbs_cell_system == NULL || _vcbs_cell_builder == NULL)
870 return false;
871 else
872 return true;
873 }
874
875 /*----------------------------------------------------------------------------*/
876 /*!
877 * \brief Allocate work buffer and general structures related to CDO
878 * vertex+cell-based schemes
879 * Set shared pointers.
880 *
881 * \param[in] quant additional mesh quantities struct.
882 * \param[in] connect pointer to a cs_cdo_connect_t struct.
883 * \param[in] time_step pointer to a time step structure
884 * \param[in] ms pointer to a cs_matrix_structure_t structure
885 */
886 /*----------------------------------------------------------------------------*/
887
888 void
cs_cdovcb_scaleq_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)889 cs_cdovcb_scaleq_init_common(const cs_cdo_quantities_t *quant,
890 const cs_cdo_connect_t *connect,
891 const cs_time_step_t *time_step,
892 const cs_matrix_structure_t *ms)
893 {
894 /* Assign static const pointers */
895
896 cs_shared_quant = quant;
897 cs_shared_connect = connect;
898 cs_shared_time_step = time_step;
899 cs_shared_ms = ms;
900
901 /* Specific treatment for handling openMP */
902
903 BFT_MALLOC(_vcbs_cell_system, cs_glob_n_threads, cs_cell_sys_t *);
904 BFT_MALLOC(_vcbs_cell_builder, cs_glob_n_threads, cs_cell_builder_t *);
905
906 for (int i = 0; i < cs_glob_n_threads; i++) {
907 _vcbs_cell_system[i] = NULL;
908 _vcbs_cell_builder[i] = NULL;
909 }
910
911 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
912 #pragma omp parallel
913 {
914 int t_id = omp_get_thread_num();
915 assert(t_id < cs_glob_n_threads);
916
917 _vcbs_cell_system[t_id] = cs_cell_sys_create(connect->n_max_vbyc + 1,
918 connect->n_max_fbyc,
919 1, NULL);
920 _vcbs_cell_builder[t_id] = _cell_builder_create(connect);
921 }
922 #else
923 assert(cs_glob_n_threads == 1);
924 _vcbs_cell_system[0] = cs_cell_sys_create(connect->n_max_vbyc + 1,
925 connect->n_max_fbyc,
926 1, NULL);
927 _vcbs_cell_builder[0] = _cell_builder_create(connect);
928 #endif /* openMP */
929 }
930
931 /*----------------------------------------------------------------------------*/
932 /*!
933 * \brief Retrieve work buffers used for building a CDO system cellwise
934 *
935 * \param[out] csys pointer to a pointer on a cs_cell_sys_t structure
936 * \param[out] cb pointer to a pointer on a cs_cell_builder_t structure
937 */
938 /*----------------------------------------------------------------------------*/
939
940 void
cs_cdovcb_scaleq_get(cs_cell_sys_t ** csys,cs_cell_builder_t ** cb)941 cs_cdovcb_scaleq_get(cs_cell_sys_t **csys,
942 cs_cell_builder_t **cb)
943 {
944 int t_id = 0;
945
946 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
947 t_id = omp_get_thread_num();
948 assert(t_id < cs_glob_n_threads);
949 #endif /* openMP */
950
951 *csys = _vcbs_cell_system[t_id];
952 *cb = _vcbs_cell_builder[t_id];
953 }
954
955 /*----------------------------------------------------------------------------*/
956 /*!
957 * \brief Free work buffer and general structure related to CDO vertex-based
958 * schemes
959 */
960 /*----------------------------------------------------------------------------*/
961
962 void
cs_cdovcb_scaleq_finalize_common(void)963 cs_cdovcb_scaleq_finalize_common(void)
964 {
965 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
966 #pragma omp parallel
967 {
968 int t_id = omp_get_thread_num();
969 cs_cell_sys_free(&(_vcbs_cell_system[t_id]));
970 cs_cell_builder_free(&(_vcbs_cell_builder[t_id]));
971 }
972 #else
973 assert(cs_glob_n_threads == 1);
974 cs_cell_sys_free(&(_vcbs_cell_system[0]));
975 cs_cell_builder_free(&(_vcbs_cell_builder[0]));
976 #endif /* openMP */
977
978 BFT_FREE(_vcbs_cell_system);
979 BFT_FREE(_vcbs_cell_builder);
980 _vcbs_cell_builder = NULL;
981 _vcbs_cell_system = NULL;
982 }
983
984 /*----------------------------------------------------------------------------*/
985 /*!
986 * \brief Initialize a cs_cdovcb_scaleq_t structure storing data useful
987 * for building and managing such a scheme
988 *
989 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
990 * \param[in] var_id id of the variable field
991 * \param[in] bflux_id id of the boundary flux field
992 * \param[in, out] eqb pointer to a \ref cs_equation_builder_t structure
993 *
994 * \return a pointer to a new allocated cs_cdovcb_scaleq_t structure
995 */
996 /*----------------------------------------------------------------------------*/
997
998 void *
cs_cdovcb_scaleq_init_context(const cs_equation_param_t * eqp,int var_id,int bflux_id,cs_equation_builder_t * eqb)999 cs_cdovcb_scaleq_init_context(const cs_equation_param_t *eqp,
1000 int var_id,
1001 int bflux_id,
1002 cs_equation_builder_t *eqb)
1003 {
1004 assert(eqp != NULL && eqb != NULL);
1005
1006 if (eqp->space_scheme != CS_SPACE_SCHEME_CDOVCB && eqp->dim != 1)
1007 bft_error(__FILE__, __LINE__, 0, " Invalid type of equation.\n"
1008 " Expected: scalar-valued CDO vertex+cell-based equation.");
1009
1010 const cs_cdo_connect_t *connect = cs_shared_connect;
1011 const cs_lnum_t n_vertices = connect->n_vertices;
1012 const cs_lnum_t n_cells = connect->n_cells;
1013
1014 cs_cdovcb_scaleq_t *eqc = NULL;
1015
1016 BFT_MALLOC(eqc, 1, cs_cdovcb_scaleq_t);
1017
1018 eqc->var_field_id = var_id;
1019 eqc->bflux_field_id = bflux_id;
1020
1021 /* System dimension */
1022
1023 eqc->n_dofs = n_vertices + n_cells;
1024
1025 /* Store the last computed values of the field at cell centers and the data
1026 needed to compute the cell values from the vertex values.
1027 No need to synchronize all these quantities since they are only cellwise
1028 quantities. */
1029
1030 BFT_MALLOC(eqc->cell_values, n_cells, cs_real_t);
1031 BFT_MALLOC(eqc->rc_tilda, n_cells, cs_real_t);
1032 BFT_MALLOC(eqc->acv_tilda, connect->c2v->idx[n_cells], cs_real_t);
1033
1034 eqc->cell_values_pre = NULL; /* For a future usage */
1035 memset(eqc->cell_values, 0, sizeof(cs_real_t)*n_cells);
1036 memset(eqc->rc_tilda, 0, sizeof(cs_real_t)*n_cells);
1037 memset(eqc->acv_tilda, 0, sizeof(cs_real_t)*connect->c2v->idx[n_cells]);
1038
1039 /* Flag to indicate what to build in a cell mesh */
1040
1041 eqb->msh_flag = CS_FLAG_COMP_PV | CS_FLAG_COMP_PVQ | CS_FLAG_COMP_DEQ |
1042 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_EV |
1043 CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_PFC |
1044 CS_FLAG_COMP_HFQ;
1045 eqb->bd_msh_flag = 0;
1046
1047 bool need_eigen =
1048 (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_NITSCHE ||
1049 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_WEAK_SYM) ? true : false;
1050
1051 /* Diffusion term */
1052
1053 eqc->diffusion_hodge = NULL;
1054 eqc->get_stiffness_matrix = NULL;
1055
1056 if (cs_equation_param_has_diffusion(eqp)) {
1057
1058 eqc->diffusion_hodge = cs_hodge_init_context(connect,
1059 eqp->diffusion_property,
1060 &(eqp->diffusion_hodgep),
1061 true, /* tensor ? */
1062 need_eigen); /* eigen ? */
1063
1064 eqc->get_stiffness_matrix = cs_hodge_vcb_get_stiffness;
1065
1066 }
1067
1068 /* Dirichlet boundary condition enforcement */
1069
1070 BFT_MALLOC(eqc->vtx_bc_flag, n_vertices, cs_flag_t);
1071 cs_equation_set_vertex_bc_flag(connect, eqb->face_bc, eqc->vtx_bc_flag);
1072
1073 /* Boundary conditions (no other choice for Robin boundary conditions */
1074
1075 eqc->enforce_robin_bc = cs_cdo_diffusion_svb_wbs_robin;
1076 eqc->enforce_dirichlet = NULL;
1077 switch (eqp->default_enforcement) {
1078
1079 case CS_PARAM_BC_ENFORCE_ALGEBRAIC:
1080 eqc->enforce_dirichlet = cs_cdo_diffusion_alge_dirichlet;
1081 break;
1082 case CS_PARAM_BC_ENFORCE_PENALIZED:
1083 eqc->enforce_dirichlet = cs_cdo_diffusion_pena_dirichlet;
1084 break;
1085 case CS_PARAM_BC_ENFORCE_WEAK_NITSCHE:
1086 eqc->enforce_dirichlet = cs_cdo_diffusion_vcb_weak_dirichlet;
1087 break;
1088 case CS_PARAM_BC_ENFORCE_WEAK_SYM:
1089 eqc->enforce_dirichlet = cs_cdo_diffusion_vcb_wsym_dirichlet;
1090 break;
1091
1092 default:
1093 bft_error(__FILE__, __LINE__, 0,
1094 " %s: Invalid type of algorithm to enforce Dirichlet BC.",
1095 __func__);
1096 break;
1097 }
1098
1099 /* Non-homogeneous Neumann BCs */
1100
1101 if (eqb->face_bc->n_nhmg_neu_faces > 0)
1102 eqb->bd_msh_flag = CS_FLAG_COMP_FV;
1103
1104 /* Advection term */
1105
1106 eqc->get_advection_matrix = NULL;
1107 eqc->add_advection_bc = NULL;
1108
1109 if (cs_equation_param_has_convection(eqp)) {
1110
1111 switch (eqp->adv_scheme) {
1112 case CS_PARAM_ADVECTION_SCHEME_CIP:
1113
1114 eqb->msh_flag |= CS_FLAG_COMP_EF;
1115 _set_cip_coef(eqp);
1116 eqc->add_advection_bc = cs_cdo_advection_vcb_bc;
1117
1118 if (cs_advection_field_is_cellwise(eqp->adv_field))
1119 eqc->get_advection_matrix = cs_cdo_advection_vcb_cw_cst;
1120 else
1121 eqc->get_advection_matrix = cs_cdo_advection_vcb;
1122 break;
1123
1124 case CS_PARAM_ADVECTION_SCHEME_CIP_CW:
1125 eqb->msh_flag |= CS_FLAG_COMP_EF;
1126 _set_cip_coef(eqp);
1127 eqc->add_advection_bc = cs_cdo_advection_vcb_bc;
1128 eqc->get_advection_matrix = cs_cdo_advection_vcb_cw_cst;
1129 break;
1130
1131 default:
1132 bft_error(__FILE__, __LINE__, 0,
1133 " Invalid advection scheme for vertex-based discretization");
1134
1135 } /* Scheme */
1136
1137 }
1138 else {
1139
1140 if (eqp->default_enforcement != CS_PARAM_BC_ENFORCE_WEAK_NITSCHE)
1141 eqb->sys_flag |= CS_FLAG_SYS_SYM; /* Algebraic system is symmetric */
1142
1143 }
1144
1145 /* A mass matrix can be requested either for the reaction term, the unsteady
1146 term or for the source term */
1147
1148 cs_hodge_algo_t reac_hodge_algo = CS_HODGE_N_ALGOS;
1149 cs_hodge_algo_t time_hodge_algo = CS_HODGE_N_ALGOS;
1150 cs_hodge_algo_t srct_hodge_algo = CS_HODGE_N_ALGOS;
1151
1152 /* Reaction term */
1153
1154 if (cs_equation_param_has_reaction(eqp)) {
1155
1156 if (eqp->do_lumping) {
1157
1158 eqb->sys_flag |= CS_FLAG_SYS_REAC_DIAG;
1159 reac_hodge_algo = CS_HODGE_ALGO_VORONOI;
1160
1161 }
1162 else {
1163
1164 switch (eqp->reaction_hodgep.algo) {
1165
1166 case CS_HODGE_ALGO_VORONOI:
1167 eqb->sys_flag |= CS_FLAG_SYS_REAC_DIAG;
1168 reac_hodge_algo = CS_HODGE_ALGO_VORONOI;
1169 break;
1170 case CS_HODGE_ALGO_WBS:
1171 eqb->sys_flag |= CS_FLAG_SYS_MASS_MATRIX;
1172 reac_hodge_algo = CS_HODGE_ALGO_WBS;
1173 break;
1174 default:
1175 bft_error(__FILE__, __LINE__, 0,
1176 "%s: Invalid choice of algorithm for the reaction term.",
1177 __func__);
1178 break;
1179 }
1180
1181 } /* Lumping or not lumping */
1182
1183 } /* Reaction */
1184
1185 /* Unsteady term */
1186
1187 if (cs_equation_param_has_time(eqp)) {
1188
1189 if (eqp->do_lumping) {
1190
1191 eqb->sys_flag |= CS_FLAG_SYS_TIME_DIAG;
1192 time_hodge_algo = CS_HODGE_ALGO_VORONOI;
1193
1194 }
1195 else {
1196
1197 switch (eqp->time_hodgep.algo) {
1198
1199 case CS_HODGE_ALGO_VORONOI:
1200 eqb->sys_flag |= CS_FLAG_SYS_TIME_DIAG;
1201 time_hodge_algo = CS_HODGE_ALGO_VORONOI;
1202 break;
1203
1204 case CS_HODGE_ALGO_WBS:
1205 eqb->sys_flag |= CS_FLAG_SYS_MASS_MATRIX;
1206 time_hodge_algo = CS_HODGE_ALGO_WBS;
1207 break;
1208 default:
1209 bft_error(__FILE__, __LINE__, 0,
1210 "%s: Invalid choice of algorithm for the time term.",
1211 __func__);
1212 break;
1213 }
1214
1215 } /* Lumping or not lumping */
1216
1217 } /* Unsteady term requested */
1218
1219 /* Source term */
1220
1221 eqc->source_terms = NULL;
1222
1223 if (cs_equation_param_has_sourceterm(eqp)) {
1224
1225 /* This is a mandatory choice for this kind of scheme */
1226
1227 srct_hodge_algo = CS_HODGE_ALGO_WBS;
1228
1229 if (cs_equation_param_has_time(eqp)) {
1230
1231 if (eqp->time_scheme == CS_TIME_SCHEME_THETA ||
1232 eqp->time_scheme == CS_TIME_SCHEME_CRANKNICO) {
1233
1234 BFT_MALLOC(eqc->source_terms, eqc->n_dofs, cs_real_t);
1235 # pragma omp parallel for if (eqc->n_dofs > CS_THR_MIN)
1236 for (cs_lnum_t i = 0; i < eqc->n_dofs; i++)
1237 eqc->source_terms[i] = 0;
1238
1239 } /* Theta scheme */
1240 } /* Time-dependent equation */
1241
1242 } /* There is at least one source term */
1243
1244 /* Pre-defined a cs_hodge_param_t structure */
1245
1246 eqc->mass_hodgep.inv_pty = false;
1247 eqc->mass_hodgep.coef = 1.0; /* not useful in this case */
1248 eqc->mass_hodgep.type = CS_HODGE_TYPE_VC;
1249 eqc->mass_hodgep.algo = cs_hodge_set_mass_algo(eqp->name,
1250 reac_hodge_algo,
1251 time_hodge_algo,
1252 srct_hodge_algo);
1253
1254 if (eqc->mass_hodgep.algo == CS_HODGE_ALGO_WBS)
1255 eqb->msh_flag |= CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PEQ
1256 | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_PFC;
1257
1258 /* Initialize the hodge structure for the mass matrix */
1259
1260 eqc->mass_hodge = cs_hodge_init_context(connect,
1261 NULL,
1262 &(eqc->mass_hodgep),
1263 false, /* tensor ? */
1264 false); /* eigen ? */
1265
1266 if (eqp->verbosity > 1 && eqb->sys_flag & CS_FLAG_SYS_MASS_MATRIX) {
1267 cs_log_printf(CS_LOG_SETUP,
1268 "#### Parameters of the mass matrix of the equation %s\n",
1269 eqp->name);
1270 cs_hodge_param_log("Mass matrix", NULL, eqc->mass_hodgep);
1271 }
1272
1273 /* Set the function pointer */
1274
1275 eqc->get_mass_matrix = cs_hodge_get_func(__func__, eqc->mass_hodgep);
1276
1277 /* Assembly process */
1278
1279 eqc->assemble = cs_equation_assemble_set(CS_SPACE_SCHEME_CDOVCB,
1280 CS_CDO_CONNECT_VTX_SCAL);
1281
1282 return eqc;
1283 }
1284
1285 /*----------------------------------------------------------------------------*/
1286 /*!
1287 * \brief Destroy a cs_cdovcb_scaleq_t structure
1288 *
1289 * \param[in, out] data pointer to a cs_cdovcb_scaleq_t structure
1290 *
1291 * \return a NULL pointer
1292 */
1293 /*----------------------------------------------------------------------------*/
1294
1295 void *
cs_cdovcb_scaleq_free_context(void * data)1296 cs_cdovcb_scaleq_free_context(void *data)
1297 {
1298 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)data;
1299
1300 if (eqc == NULL)
1301 return eqc;
1302
1303 cs_hodge_free_context(&(eqc->diffusion_hodge));
1304 cs_hodge_free_context(&(eqc->mass_hodge));
1305
1306 BFT_FREE(eqc->cell_values);
1307 BFT_FREE(eqc->rc_tilda);
1308 BFT_FREE(eqc->acv_tilda);
1309
1310 BFT_FREE(eqc->vtx_bc_flag);
1311 BFT_FREE(eqc->source_terms);
1312
1313 /* Last free */
1314
1315 BFT_FREE(eqc);
1316
1317 return NULL;
1318 }
1319
1320 /*----------------------------------------------------------------------------*/
1321 /*!
1322 * \brief Set the initial values of the variable field taking into account
1323 * the boundary conditions.
1324 * Case of scalar-valued CDO-VCb schemes.
1325 *
1326 * \param[in] t_eval time at which one evaluates BCs
1327 * \param[in] field_id id related to the variable field of this equation
1328 * \param[in] mesh pointer to a cs_mesh_t structure
1329 * \param[in] eqp pointer to a cs_equation_param_t structure
1330 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1331 * \param[in, out] context pointer to the scheme context (cast on-the-fly)
1332 */
1333 /*----------------------------------------------------------------------------*/
1334
1335 void
cs_cdovcb_scaleq_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)1336 cs_cdovcb_scaleq_init_values(cs_real_t t_eval,
1337 const int field_id,
1338 const cs_mesh_t *mesh,
1339 const cs_equation_param_t *eqp,
1340 cs_equation_builder_t *eqb,
1341 void *context)
1342 {
1343 const cs_cdo_quantities_t *quant = cs_shared_quant;
1344 const cs_cdo_connect_t *connect = cs_shared_connect;
1345
1346 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
1347 cs_field_t *fld = cs_field_by_id(field_id);
1348 cs_real_t *v_vals = fld->val;
1349 cs_real_t *c_vals = eqc->cell_values;
1350
1351 /* By default, 0 is set as initial condition for the computational domain.
1352
1353 Warning: This operation has to be done after the settings of the
1354 Dirichlet boundary conditions where an interface sum is performed
1355 for vertex-based schemes
1356 */
1357
1358 memset(v_vals, 0, quant->n_vertices*sizeof(cs_real_t));
1359 memset(c_vals, 0, quant->n_cells*sizeof(cs_real_t));
1360
1361 if (eqp->n_ic_defs > 0) {
1362
1363 cs_lnum_t *def2v_ids = (cs_lnum_t *)cs_equation_get_tmpbuf();
1364 cs_lnum_t *def2v_idx = NULL;
1365 BFT_MALLOC(def2v_idx, eqp->n_ic_defs + 1, cs_lnum_t);
1366
1367 cs_equation_sync_vol_def_at_vertices(connect,
1368 eqp->n_ic_defs,
1369 eqp->ic_defs,
1370 def2v_idx,
1371 def2v_ids);
1372
1373 /* Initialize values at mesh vertices */
1374
1375 cs_flag_t v_dof_flag = CS_FLAG_SCALAR | cs_flag_primal_vtx;
1376 cs_flag_t c_dof_flag = CS_FLAG_SCALAR | cs_flag_primal_cell;
1377
1378 for (int def_id = 0; def_id < eqp->n_ic_defs; def_id++) {
1379
1380 /* Get and then set the definition of the initial condition */
1381
1382 const cs_xdef_t *def = eqp->ic_defs[def_id];
1383 const cs_lnum_t n_v_selected = def2v_idx[def_id+1] - def2v_idx[def_id];
1384 const cs_lnum_t *selected_lst = def2v_ids + def2v_idx[def_id];
1385
1386 switch(def->type) {
1387
1388 case CS_XDEF_BY_VALUE:
1389 cs_evaluate_potential_at_vertices_by_value(def,
1390 n_v_selected,
1391 selected_lst,
1392 v_vals);
1393 cs_evaluate_potential_at_cells_by_value(def, c_vals);
1394 break;
1395
1396 case CS_XDEF_BY_QOV:
1397 cs_evaluate_potential_by_qov(v_dof_flag | c_dof_flag, def,
1398 v_vals, c_vals);
1399 break;
1400
1401 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1402 if (eqp->dof_reduction != CS_PARAM_REDUCTION_DERHAM)
1403 bft_error(__FILE__, __LINE__, 0,
1404 " %s: Incompatible reduction for equation %s.\n",
1405 __func__, eqp->name);
1406 cs_evaluate_potential_at_vertices_by_analytic(def,
1407 t_eval,
1408 n_v_selected,
1409 selected_lst,
1410 v_vals);
1411 cs_evaluate_potential_at_cells_by_analytic(def, t_eval, c_vals);
1412 break;
1413
1414 default:
1415 bft_error(__FILE__, __LINE__, 0,
1416 " %s: Invalid way to initialize field values for eq. %s.\n",
1417 __func__, eqp->name);
1418
1419 } /* Switch on possible type of definition */
1420
1421 } /* Loop on definitions */
1422
1423 BFT_FREE(def2v_idx);
1424
1425 } /* Initial values to set */
1426
1427 /* Set the boundary values as initial values: Compute the values of the
1428 Dirichlet BC */
1429
1430 cs_equation_compute_dirichlet_vb(t_eval,
1431 mesh,
1432 quant,
1433 connect,
1434 eqp,
1435 eqb->face_bc,
1436 _vcbs_cell_builder[0], /* static variable */
1437 eqc->vtx_bc_flag,
1438 v_vals);
1439 }
1440
1441 /*----------------------------------------------------------------------------*/
1442 /*!
1443 * \brief Build and solve the linear system arising from a scalar steady-state
1444 * equation with a CDO-VCb scheme. Use for interpolation purpose from
1445 * cell values to vertex values.
1446 * One works cellwise and then process to the assembly.
1447 *
1448 * \param[in] mesh pointer to a cs_mesh_t structure
1449 * \param[in] cell_values array of cell values
1450 * \param[in] field_id id of the variable field
1451 * \param[in] eqp pointer to a cs_equation_param_t structure
1452 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1453 * \param[in, out] context pointer to cs_cdovcb_scaleq_t structure
1454 */
1455 /*----------------------------------------------------------------------------*/
1456
1457 void
cs_cdovcb_scaleq_interpolate(const cs_mesh_t * mesh,const cs_real_t * cell_values,const int field_id,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)1458 cs_cdovcb_scaleq_interpolate(const cs_mesh_t *mesh,
1459 const cs_real_t *cell_values,
1460 const int field_id,
1461 const cs_equation_param_t *eqp,
1462 cs_equation_builder_t *eqb,
1463 void *context)
1464 {
1465 const cs_cdo_connect_t *connect = cs_shared_connect;
1466 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_VTX_SCAL];
1467 const cs_cdo_quantities_t *quant = cs_shared_quant;
1468 const cs_lnum_t n_vertices = quant->n_vertices;
1469 const cs_time_step_t *ts = cs_shared_time_step;
1470 const cs_real_t time_eval = ts->t_cur; /* Interpolation case */
1471
1472 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
1473 cs_field_t *fld = cs_field_by_id(field_id);
1474
1475 cs_timer_t t0 = cs_timer_time();
1476
1477 /* Build an array storing the Dirichlet values at vertices
1478 * First argument is set to t_cur even if this is a steady computation since
1479 * one can call this function to compute a steady-state solution at each time
1480 * step of an unsteady computation. Dirichlet boundary conditions are always
1481 * evaluated at t_cur in case of interpolation
1482 */
1483
1484 _setup_vcb(time_eval, mesh, connect, eqp, eqb, eqc->vtx_bc_flag);
1485
1486 if (eqb->init_step)
1487 eqb->init_step = false;
1488
1489 /* Initialize the local system: matrix and rhs */
1490
1491 cs_matrix_t *matrix = cs_matrix_create(cs_shared_ms);
1492 cs_real_t *rhs = NULL;
1493
1494 BFT_MALLOC(rhs, n_vertices, cs_real_t);
1495 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
1496 for (cs_lnum_t i = 0; i < n_vertices; i++) rhs[i] = 0.0;
1497
1498 /* Initialize the structure to assemble values */
1499
1500 cs_matrix_assembler_values_t *mav
1501 = cs_matrix_assembler_values_init(matrix, NULL, NULL);
1502
1503 /* -------------------------
1504 * Main OpenMP block on cell
1505 * ------------------------- */
1506
1507 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
1508 shared(quant, connect, eqp, eqb, eqc, rhs, matrix, mav, fld, rs, \
1509 _vcbs_cell_system, _vcbs_cell_builder, cell_values) \
1510 firstprivate(time_eval)
1511 {
1512 /* Set variables and structures inside the OMP section so that each thread
1513 has its own value */
1514
1515 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
1516 int t_id = omp_get_thread_num();
1517 #else
1518 int t_id = 0;
1519 #endif
1520
1521 /* Each thread get back its related structures:
1522 Get the cell-wise view of the mesh and the algebraic system */
1523
1524 cs_face_mesh_t *fm = cs_cdo_local_get_face_mesh(t_id);
1525 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
1526 cs_cell_sys_t *csys = _vcbs_cell_system[t_id];
1527 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
1528 cs_equation_assemble_t *eqa = cs_equation_assemble_get(t_id);
1529 cs_hodge_t *diff_hodge =
1530 (eqc->diffusion_hodge == NULL) ? NULL : eqc->diffusion_hodge[t_id];
1531 cs_hodge_t *mass_hodge =
1532 (eqc->mass_hodge == NULL) ? NULL : eqc->mass_hodge[t_id];
1533
1534 /* Set times at which one evaluates quantities if needed */
1535
1536 cb->t_pty_eval = time_eval;
1537 cb->t_bc_eval = time_eval;
1538 cb->t_st_eval = time_eval;
1539
1540 /* Initialization of the values of properties */
1541
1542 cs_equation_init_properties(eqp, eqb, diff_hodge, cb);
1543
1544 /* ---------------------------------------------
1545 * Main loop on cells to build the linear system
1546 * --------------------------------------------- */
1547
1548 # pragma omp for CS_CDO_OMP_SCHEDULE
1549 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1550
1551 /* Set the current cell flag */
1552
1553 cb->cell_flag = connect->cell_flag[c_id];
1554
1555 /* Set the local mesh structure for the current cell */
1556
1557 cs_cell_mesh_build(c_id,
1558 cs_equation_cell_mesh_flag(cb->cell_flag, eqb),
1559 connect, quant, cm);
1560
1561 /* Set the local (i.e. cellwise) structures for the current cell */
1562
1563 _svcb_init_cell_system(cm, eqp, eqb, eqc, eqc->vtx_bc_flag, fld->val,
1564 csys, cb);
1565
1566 /* Build and add the diffusion/advection/reaction term to the local
1567 system. A mass matrix is also built if needed. */
1568
1569 _svcb_conv_diff_reac(eqp, eqb, eqc, cm, fm,
1570 mass_hodge, diff_hodge, csys, cb);
1571
1572 if (cs_equation_param_has_sourceterm(eqp)) { /* SOURCE TERM
1573 * =========== */
1574
1575 /* Reset the local contribution */
1576
1577 memset(csys->source, 0, csys->n_dofs*sizeof(cs_real_t));
1578
1579 /* Source term contribution to the algebraic system
1580 If the equation is steady, the source term has already been computed
1581 and is added to the right-hand side during its initialization. */
1582
1583 cs_source_term_compute_cellwise(eqp->n_source_terms,
1584 (cs_xdef_t *const *)eqp->source_terms,
1585 cm,
1586 eqb->source_mask,
1587 eqb->compute_source,
1588 cb->t_st_eval,
1589 mass_hodge,
1590 cb,
1591 csys->source);
1592
1593 for (short int v = 0; v < cm->n_vc; v++)
1594 csys->rhs[v] += csys->source[v];
1595 csys->rhs[cm->n_vc] += csys->source[cm->n_vc];
1596
1597 } /* End of term source */
1598
1599 /* Apply boundary conditions (those which are weakly enforced) */
1600
1601 _svcb_apply_weak_bc(eqp, eqc, cm, fm, diff_hodge, csys, cb);
1602
1603 { /* Reduce the system size since one has the knowledge of the cell
1604 value */
1605
1606 /* Reshape the local system */
1607
1608 for (short int i = 0; i < cm->n_vc; i++) {
1609
1610 double *old_i = csys->mat->val + csys->n_dofs*i; /* Old "i" row */
1611 double *new_i = csys->mat->val + cm->n_vc*i; /* New "i" row */
1612
1613 for (short int j = 0; j < cm->n_vc; j++)
1614 new_i[j] = old_i[j];
1615
1616 /* Update RHS: RHS = RHS - Avc*pc */
1617
1618 csys->rhs[i] -= cell_values[csys->c_id] * old_i[cm->n_vc];
1619
1620 }
1621 csys->n_dofs = cm->n_vc;
1622 csys->mat->n_rows = csys->mat->n_cols = cm->n_vc;
1623
1624 }
1625
1626 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
1627 if (cs_dbg_cw_test(eqp, cm, csys))
1628 cs_cell_sys_dump(">> Cell system matrix after condensation", csys);
1629 #endif
1630
1631 /* Enforce values if needed (internal or Dirichlet) */
1632
1633 _svcb_enforce_values(eqp, eqb, eqc, cm, fm, diff_hodge, csys, cb);
1634
1635 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 0
1636 if (cs_dbg_cw_test(eqp, cm, csys))
1637 cs_cell_sys_dump(">> (FINAL) Cell system matrix", csys);
1638 #endif
1639
1640 /* ASSEMBLY PROCESS
1641 * ================ */
1642
1643 _assemble(eqc, cm, csys, rs, eqa, mav, rhs);
1644
1645 } /* Main loop on cells */
1646
1647 } /* OPENMP Block */
1648
1649 cs_matrix_assembler_values_done(mav); /* optional */
1650
1651 /* Free temporary buffers and structures */
1652
1653 cs_equation_builder_reset(eqb);
1654 cs_matrix_assembler_values_finalize(&mav);
1655
1656 /* End of the system building */
1657
1658 cs_timer_t t1 = cs_timer_time();
1659 cs_timer_counter_add_diff(&(eqb->tcb), &t0, &t1);
1660
1661 /* Copy current field values to previous values */
1662
1663 cs_field_current_to_previous(fld);
1664
1665 /* Solve the linear system */
1666 cs_real_t normalization = 1.0;
1667 cs_sles_t *sles = cs_sles_find_or_add(eqp->sles_param->field_id, NULL);
1668
1669 cs_equation_solve_scalar_system(n_vertices,
1670 eqp->sles_param,
1671 matrix,
1672 rs,
1673 normalization,
1674 true, /* rhs_redux */
1675 sles,
1676 fld->val,
1677 rhs);
1678
1679 /* Compute values at cells pc = acc^-1*(RHS - Acv*pv) */
1680
1681 memcpy(eqc->cell_values, cell_values, quant->n_cells*sizeof(cs_real_t));
1682
1683 /* Free remaining buffers */
1684
1685 BFT_FREE(rhs);
1686 cs_sles_free(sles);
1687 cs_matrix_destroy(&matrix);
1688 }
1689
1690 /*----------------------------------------------------------------------------*/
1691 /*!
1692 * \brief Build and solve the linear system arising from a scalar steady-state
1693 * convection/diffusion/reaction equation with a CDO-VCb scheme
1694 * One works cellwise and then process to the assembly.
1695 *
1696 * \param[in] cur2prev true="current to previous" operation is performed
1697 * \param[in] mesh pointer to a cs_mesh_t structure
1698 * \param[in] field_id id of the variable field related to this equation
1699 * \param[in] eqp pointer to a cs_equation_param_t structure
1700 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1701 * \param[in, out] context pointer to cs_cdovcb_scaleq_t structure
1702 */
1703 /*----------------------------------------------------------------------------*/
1704
1705 void
cs_cdovcb_scaleq_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)1706 cs_cdovcb_scaleq_solve_steady_state(bool cur2prev,
1707 const cs_mesh_t *mesh,
1708 const int field_id,
1709 const cs_equation_param_t *eqp,
1710 cs_equation_builder_t *eqb,
1711 void *context)
1712 {
1713 cs_timer_t t0 = cs_timer_time();
1714
1715 const cs_cdo_connect_t *connect = cs_shared_connect;
1716 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_VTX_SCAL];
1717 const cs_cdo_quantities_t *quant = cs_shared_quant;
1718 const cs_lnum_t n_vertices = quant->n_vertices;
1719 const cs_time_step_t *ts = cs_shared_time_step;
1720
1721 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
1722 cs_field_t *fld = cs_field_by_id(field_id);
1723
1724 /* Build an array storing the Dirichlet values at vertices
1725 * First argument is set to t_cur even if this is a steady computation since
1726 * one can call this function to compute a steady-state solution at each time
1727 * step of an unsteady computation. Dirichlet boundary conditions are always
1728 * evaluated at t + dt */
1729
1730 _setup_vcb(ts->t_cur + ts->dt[0], mesh, connect, eqp, eqb, eqc->vtx_bc_flag);
1731
1732 if (eqb->init_step)
1733 eqb->init_step = false;
1734
1735 /* Initialize the local system: matrix and rhs */
1736
1737 cs_matrix_t *matrix = cs_matrix_create(cs_shared_ms);
1738 cs_real_t *rhs = NULL;
1739 double rhs_norm = 0.;
1740
1741 BFT_MALLOC(rhs, n_vertices, cs_real_t);
1742 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
1743 for (cs_lnum_t i = 0; i < n_vertices; i++) rhs[i] = 0.0;
1744
1745 /* Initialize the structure to assemble values */
1746
1747 cs_matrix_assembler_values_t *mav
1748 = cs_matrix_assembler_values_init(matrix, NULL, NULL);
1749
1750 /* -------------------------
1751 * Main OpenMP block on cell
1752 * ------------------------- */
1753
1754 # pragma omp parallel if (quant->n_cells > CS_THR_MIN)
1755 {
1756 /* Set variables and structures inside the OMP section so that each thread
1757 has its own value */
1758
1759 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
1760 int t_id = omp_get_thread_num();
1761 #else
1762 int t_id = 0;
1763 #endif
1764
1765 /* Each thread get back its related structures:
1766 Get the cell-wise view of the mesh and the algebraic system */
1767
1768 cs_face_mesh_t *fm = cs_cdo_local_get_face_mesh(t_id);
1769 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
1770 cs_cell_sys_t *csys = _vcbs_cell_system[t_id];
1771 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
1772 cs_equation_assemble_t *eqa = cs_equation_assemble_get(t_id);
1773 cs_hodge_t *diff_hodge =
1774 (eqc->diffusion_hodge == NULL) ? NULL : eqc->diffusion_hodge[t_id];
1775 cs_hodge_t *mass_hodge =
1776 (eqc->mass_hodge == NULL) ? NULL : eqc->mass_hodge[t_id];
1777
1778 /* Set times at which one evaluates quantities if needed */
1779
1780 const cs_real_t time_eval = ts->t_cur + ts->dt[0]; /* dummy parameter */
1781 cb->t_pty_eval = time_eval;
1782 cb->t_bc_eval = time_eval;
1783 cb->t_st_eval = time_eval;
1784
1785 /* Initialization of the values of properties */
1786
1787 cs_equation_init_properties(eqp, eqb, diff_hodge, cb);
1788
1789 /* ---------------------------------------------
1790 * Main loop on cells to build the linear system
1791 * --------------------------------------------- */
1792
1793 # pragma omp for CS_CDO_OMP_SCHEDULE reduction(+:rhs_norm)
1794 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
1795
1796 /* Set the current cell flag */
1797
1798 cb->cell_flag = connect->cell_flag[c_id];
1799
1800 /* Set the local mesh structure for the current cell */
1801
1802 cs_cell_mesh_build(c_id,
1803 cs_equation_cell_mesh_flag(cb->cell_flag, eqb),
1804 connect, quant, cm);
1805
1806 /* Set the local (i.e. cellwise) structures for the current cell */
1807
1808 _svcb_init_cell_system(cm, eqp, eqb, eqc, eqc->vtx_bc_flag, fld->val,
1809 csys, cb);
1810
1811 /* Build and add the diffusion/advection/reaction term to the local
1812 system. A mass matrix is also built if needed. */
1813
1814 _svcb_conv_diff_reac(eqp, eqb, eqc, cm, fm,
1815 mass_hodge, diff_hodge, csys, cb);
1816
1817 if (cs_equation_param_has_sourceterm(eqp)) { /* SOURCE TERM
1818 * =========== */
1819 /* Reset the local contribution */
1820
1821 memset(csys->source, 0, csys->n_dofs*sizeof(cs_real_t));
1822
1823 /* Source term contribution to the algebraic system
1824 If the equation is steady, the source term has already been computed
1825 and is added to the right-hand side during its initialization. */
1826
1827 cs_source_term_compute_cellwise(eqp->n_source_terms,
1828 (cs_xdef_t *const *)eqp->source_terms,
1829 cm,
1830 eqb->source_mask,
1831 eqb->compute_source,
1832 cb->t_st_eval,
1833 mass_hodge,
1834 cb,
1835 csys->source);
1836
1837 for (short int v = 0; v < cm->n_vc; v++)
1838 csys->rhs[v] += csys->source[v];
1839 csys->rhs[cm->n_vc] += csys->source[cm->n_vc];
1840
1841 } /* End of term source */
1842
1843 /* Apply boundary conditions (those which are weakly enforced) */
1844
1845 _svcb_apply_weak_bc(eqp, eqc, cm, fm, diff_hodge, csys, cb);
1846
1847 /* STATIC CONDENSATION
1848 * ===================
1849 * of the local system matrix of size n_vc + 1 into a matrix of size n_vc.
1850 * Store data in rc_tilda and acv_tilda to compute the values at cell
1851 * centers after solving the system */
1852
1853 cs_static_condensation_scalar_eq(connect->c2v,
1854 eqc->rc_tilda, eqc->acv_tilda,
1855 cb, csys);
1856
1857 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
1858 if (cs_dbg_cw_test(eqp, cm, csys))
1859 cs_cell_sys_dump(">> Cell system matrix after condensation", csys);
1860 #endif
1861
1862 /* Compute a cellwise norm of the RHS for the normalization of the
1863 residual during the resolution of the linear system */
1864
1865 rhs_norm += _svcb_cw_rhs_normalization(eqp->sles_param->resnorm_type,
1866 cm, csys);
1867
1868 /* Enforce values if needed (internal or Dirichlet) */
1869
1870 _svcb_enforce_values(eqp, eqb, eqc, cm, fm, diff_hodge, csys, cb);
1871
1872 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 0
1873 if (cs_dbg_cw_test(eqp, cm, csys))
1874 cs_cell_sys_dump(">> (FINAL) Cell system matrix", csys);
1875 #endif
1876
1877 /* ASSEMBLY PROCESS
1878 * ================ */
1879
1880 _assemble(eqc, cm, csys, rs, eqa, mav, rhs);
1881
1882 } /* Main loop on cells */
1883
1884 } /* OPENMP Block */
1885
1886 cs_matrix_assembler_values_done(mav); /* optional */
1887
1888 /* Free temporary buffers and structures */
1889
1890 cs_equation_builder_reset(eqb);
1891 cs_matrix_assembler_values_finalize(&mav);
1892
1893 /* Copy current field values to previous values */
1894
1895 if (cur2prev)
1896 cs_field_current_to_previous(fld);
1897
1898 /* End of the system building */
1899
1900 cs_timer_t t1 = cs_timer_time();
1901 cs_timer_counter_add_diff(&(eqb->tcb), &t0, &t1);
1902
1903 /* Solve the linear system */
1904 /* ======================= */
1905
1906 /* Last step in the computation of the renormalization coefficient */
1907
1908 cs_equation_sync_rhs_normalization(eqp->sles_param->resnorm_type,
1909 n_vertices,
1910 rhs,
1911 &rhs_norm);
1912
1913 cs_sles_t *sles = cs_sles_find_or_add(eqp->sles_param->field_id, NULL);
1914
1915 cs_equation_solve_scalar_system(n_vertices,
1916 eqp->sles_param,
1917 matrix,
1918 rs,
1919 rhs_norm,
1920 true, /* rhs_redux */
1921 sles,
1922 fld->val,
1923 rhs);
1924
1925 cs_timer_t t2 = cs_timer_time();
1926 cs_timer_counter_add_diff(&(eqb->tcs), &t1, &t2);
1927
1928 /* Update fields */
1929
1930 _update_cell_fields(&(eqb->tce), fld, eqc, cur2prev);
1931
1932 /* Free remaining buffers */
1933
1934 BFT_FREE(rhs);
1935 cs_sles_free(sles);
1936 cs_matrix_destroy(&matrix);
1937 }
1938
1939 /*----------------------------------------------------------------------------*/
1940 /*!
1941 * \brief Build and solve the linear system arising from a scalar unsteady
1942 * convection/diffusion/reaction equation with a CDO-VCb scheme
1943 * Time scheme is an implicit Euler
1944 * One works cellwise and then process to the assembly.
1945 *
1946 * \param[in] cur2prev true="current to previous" operation is performed
1947 * \param[in] mesh pointer to a cs_mesh_t structure
1948 * \param[in] field_id id of the variable field related to this equation
1949 * \param[in] eqp pointer to a cs_equation_param_t structure
1950 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
1951 * \param[in, out] context pointer to cs_cdovcb_scaleq_t structure
1952 */
1953 /*----------------------------------------------------------------------------*/
1954
1955 void
cs_cdovcb_scaleq_solve_implicit(bool cur2prev,const cs_mesh_t * mesh,const int field_id,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)1956 cs_cdovcb_scaleq_solve_implicit(bool cur2prev,
1957 const cs_mesh_t *mesh,
1958 const int field_id,
1959 const cs_equation_param_t *eqp,
1960 cs_equation_builder_t *eqb,
1961 void *context)
1962 {
1963 cs_timer_t t0 = cs_timer_time();
1964
1965 const cs_cdo_connect_t *connect = cs_shared_connect;
1966 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_VTX_SCAL];
1967 const cs_cdo_quantities_t *quant = cs_shared_quant;
1968 const cs_lnum_t n_vertices = quant->n_vertices;
1969 const cs_time_step_t *ts = cs_shared_time_step;
1970
1971 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
1972 cs_field_t *fld = cs_field_by_id(field_id);
1973
1974 assert(cs_equation_param_has_time(eqp) == true);
1975 assert(eqp->time_scheme == CS_TIME_SCHEME_EULER_IMPLICIT);
1976
1977 /* Build an array storing the Dirichlet values at vertices.
1978 * Dirichlet boundary conditions are always evaluated at t + dt */
1979
1980 _setup_vcb(ts->t_cur + ts->dt[0], mesh, connect, eqp, eqb, eqc->vtx_bc_flag);
1981
1982 if (eqb->init_step)
1983 eqb->init_step = false;
1984
1985 /* Initialize the local system: matrix and rhs */
1986
1987 cs_matrix_t *matrix = cs_matrix_create(cs_shared_ms);
1988 cs_real_t *rhs = NULL;
1989 double rhs_norm = 0.;
1990
1991 BFT_MALLOC(rhs, n_vertices, cs_real_t);
1992 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
1993 for (cs_lnum_t i = 0; i < n_vertices; i++) rhs[i] = 0.0;
1994
1995 /* Initialize the structure to assemble values */
1996
1997 cs_matrix_assembler_values_t *mav
1998 = cs_matrix_assembler_values_init(matrix, NULL, NULL);
1999
2000 /* -------------------------
2001 * Main OpenMP block on cell
2002 * ------------------------- */
2003
2004 # pragma omp parallel if (quant->n_cells > CS_THR_MIN)
2005 {
2006 /* Set variables and structures inside the OMP section so that each thread
2007 has its own value */
2008
2009 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
2010 int t_id = omp_get_thread_num();
2011 #else
2012 int t_id = 0;
2013 #endif
2014
2015 /* Each thread get back its related structures:
2016 Get the cell-wise view of the mesh and the algebraic system */
2017
2018 cs_face_mesh_t *fm = cs_cdo_local_get_face_mesh(t_id);
2019 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
2020 cs_cell_sys_t *csys = _vcbs_cell_system[t_id];
2021 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
2022 cs_equation_assemble_t *eqa = cs_equation_assemble_get(t_id);
2023 cs_hodge_t *diff_hodge =
2024 (eqc->diffusion_hodge == NULL) ? NULL : eqc->diffusion_hodge[t_id];
2025 cs_hodge_t *mass_hodge =
2026 (eqc->mass_hodge == NULL) ? NULL : eqc->mass_hodge[t_id];
2027
2028 /* Set times at which one evaluates quantities if needed */
2029
2030 const cs_real_t time_eval = ts->t_cur + ts->dt[0];
2031 const cs_real_t inv_dtcur = 1./ts->dt[0];
2032
2033 cb->t_pty_eval = time_eval;
2034 cb->t_bc_eval = time_eval;
2035 cb->t_st_eval = time_eval;
2036
2037 /* Initialization of the values of properties */
2038
2039 cs_equation_init_properties(eqp, eqb, diff_hodge, cb);
2040
2041 /* ---------------------------------------------
2042 * Main loop on cells to build the linear system
2043 * --------------------------------------------- */
2044
2045 # pragma omp for CS_CDO_OMP_SCHEDULE reduction(+:rhs_norm)
2046 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
2047
2048 /* Set the current cell flag */
2049
2050 cb->cell_flag = connect->cell_flag[c_id];
2051
2052 /* Set the local mesh structure for the current cell */
2053
2054 cs_cell_mesh_build(c_id,
2055 cs_equation_cell_mesh_flag(cb->cell_flag, eqb),
2056 connect, quant, cm);
2057
2058 /* Set the local (i.e. cellwise) structures for the current cell */
2059
2060 _svcb_init_cell_system(cm, eqp, eqb, eqc, eqc->vtx_bc_flag, fld->val,
2061 csys, cb);
2062
2063 /* Build and add the diffusion/advection/reaction term to the local
2064 system. A mass matrix is also built if needed. */
2065
2066 _svcb_conv_diff_reac(eqp, eqb, eqc, cm, fm,
2067 mass_hodge, diff_hodge, csys, cb);
2068
2069 if (cs_equation_param_has_sourceterm(eqp)) { /* SOURCE TERM
2070 * =========== */
2071 /* Reset the local contribution */
2072
2073 memset(csys->source, 0, csys->n_dofs*sizeof(cs_real_t));
2074
2075 /* Source term contribution to the algebraic system
2076 If the equation is steady, the source term has already been computed
2077 and is added to the right-hand side during its initialization. */
2078
2079 cs_source_term_compute_cellwise(eqp->n_source_terms,
2080 (cs_xdef_t *const *)eqp->source_terms,
2081 cm,
2082 eqb->source_mask,
2083 eqb->compute_source,
2084 cb->t_st_eval,
2085 mass_hodge,
2086 cb,
2087 csys->source);
2088
2089 for (short int v = 0; v < cm->n_vc; v++)
2090 csys->rhs[v] += csys->source[v];
2091 csys->rhs[cm->n_vc] += csys->source[cm->n_vc];
2092
2093 } /* End of term source */
2094
2095 /* Apply boundary conditions (those which are weakly enforced) */
2096
2097 _svcb_apply_weak_bc(eqp, eqc, cm, fm, diff_hodge, csys, cb);
2098
2099 /* UNSTEADY TERM + TIME SCHEME
2100 * =========================== */
2101
2102 if (!(eqb->time_pty_uniform))
2103 cb->tpty_val = cs_property_value_in_cell(cm, eqp->time_property,
2104 cb->t_pty_eval);
2105
2106 if (eqb->sys_flag & CS_FLAG_SYS_TIME_DIAG) { /* Mass lumping */
2107
2108 /* |c|*wvc = |dual_cell(v) cap c| */
2109
2110 const double ptyc = cb->tpty_val * cm->vol_c * inv_dtcur;
2111
2112 /* STEPS >> Compute the time contribution to the RHS: Mtime*pn
2113 >> Update the cellwise system with the time matrix */
2114
2115 for (short int i = 0; i < cm->n_vc; i++) {
2116
2117 const double dval = 0.75 * ptyc * cm->wvc[i];
2118
2119 /* Update the RHS with values at time t_n */
2120
2121 csys->rhs[i] += dval * csys->val_n[i];
2122
2123 /* Add the diagonal contribution from time matrix */
2124
2125 csys->mat->val[i*(csys->n_dofs + 1)] += dval;
2126
2127 }
2128
2129 /* Cell row */
2130
2131 const double dvalc = 0.25 * ptyc;
2132
2133 /* Update the RHS with values at time t_n */
2134
2135 csys->rhs[cm->n_vc] += dvalc * csys->val_n[cm->n_vc];
2136
2137 /* Add the diagonal contribution from time matrix */
2138
2139 csys->mat->val[cm->n_vc*(csys->n_dofs + 1)] += dvalc;
2140
2141 }
2142 else { /* Use the mass matrix */
2143
2144 const double tpty_coef = cb->tpty_val * inv_dtcur;
2145 const cs_sdm_t *mass_mat = mass_hodge->matrix;
2146
2147 /* STEPS >> Compute the time contribution to the RHS: Mtime*pn
2148 >> Update the cellwise system with the time matrix */
2149
2150 /* Update rhs with csys->mat*p^n */
2151
2152 double *time_pn = cb->values;
2153 cs_sdm_square_matvec(mass_mat, csys->val_n, time_pn);
2154 for (short int i = 0; i < csys->n_dofs; i++)
2155 csys->rhs[i] += tpty_coef*time_pn[i];
2156
2157 /* Update the cellwise system with the time matrix */
2158
2159 cs_sdm_add_mult(csys->mat, tpty_coef, mass_mat);
2160
2161 }
2162
2163 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
2164 if (cs_dbg_cw_test(eqp, cm, csys))
2165 cs_cell_sys_dump("\n>> Cell system after adding time", csys);
2166 #endif
2167
2168 /* STATIC CONDENSATION
2169 * ===================
2170 * of the local system matrix of size n_vc + 1 into a matrix of size n_vc.
2171 * Store data in rc_tilda and acv_tilda to compute the values at cell
2172 * centers after solving the system */
2173
2174 cs_static_condensation_scalar_eq(connect->c2v,
2175 eqc->rc_tilda, eqc->acv_tilda,
2176 cb, csys);
2177
2178 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
2179 if (cs_dbg_cw_test(eqp, cm, csys))
2180 cs_cell_sys_dump(">> Cell system matrix after condensation", csys);
2181 #endif
2182
2183 /* Enforce values if needed (internal or Dirichlet) */
2184
2185 _svcb_enforce_values(eqp, eqb, eqc, cm, fm, diff_hodge, csys, cb);
2186
2187 /* Compute a cellwise norm of the RHS for the normalization of the
2188 residual during the resolution of the linear system */
2189
2190 rhs_norm += _svcb_cw_rhs_normalization(eqp->sles_param->resnorm_type,
2191 cm, csys);
2192
2193 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 0
2194 if (cs_dbg_cw_test(eqp, cm, csys))
2195 cs_cell_sys_dump(">> (FINAL) Cell system matrix", csys);
2196 #endif
2197
2198 /* ASSEMBLY PROCESS
2199 * ================ */
2200
2201 _assemble(eqc, cm, csys, rs, eqa, mav, rhs);
2202
2203 } /* Main loop on cells */
2204
2205 } /* OpenMP Block */
2206
2207 cs_matrix_assembler_values_done(mav); /* optional */
2208
2209 /* Free temporary buffers and structures */
2210
2211 cs_equation_builder_reset(eqb);
2212 cs_matrix_assembler_values_finalize(&mav);
2213
2214 /* Copy current field values to previous values */
2215
2216 if (cur2prev)
2217 cs_field_current_to_previous(fld);
2218
2219 /* End of the system building */
2220
2221 cs_timer_t t1 = cs_timer_time();
2222 cs_timer_counter_add_diff(&(eqb->tcb), &t0, &t1);
2223
2224 /* Solve the linear system */
2225 /* ======================= */
2226
2227 /* Last step in the computation of the renormalization coefficient */
2228
2229 cs_equation_sync_rhs_normalization(eqp->sles_param->resnorm_type,
2230 n_vertices,
2231 rhs,
2232 &rhs_norm);
2233
2234 cs_sles_t *sles = cs_sles_find_or_add(eqp->sles_param->field_id, NULL);
2235
2236 cs_equation_solve_scalar_system(n_vertices,
2237 eqp->sles_param,
2238 matrix,
2239 rs,
2240 rhs_norm,
2241 true, /* rhs_redux */
2242 sles,
2243 fld->val,
2244 rhs);
2245
2246 cs_timer_t t2 = cs_timer_time();
2247 cs_timer_counter_add_diff(&(eqb->tcs), &t1, &t2);
2248
2249 /* Update fields */
2250
2251 _update_cell_fields(&(eqb->tce), fld, eqc, cur2prev);
2252
2253 /* Free remaining buffers */
2254
2255 BFT_FREE(rhs);
2256 cs_sles_free(sles);
2257 cs_matrix_destroy(&matrix);
2258 }
2259
2260 /*----------------------------------------------------------------------------*/
2261 /*!
2262 * \brief Build and solve the linear system arising from a scalar unsteady
2263 * convection/diffusion/reaction equation with a CDO-VCb scheme
2264 * Time scheme is a theta scheme.
2265 * One works cellwise and then process to the assembly.
2266 *
2267 * \param[in] cur2prev true="current to previous" operation is performed
2268 * \param[in] mesh pointer to a cs_mesh_t structure
2269 * \param[in] field_id id of the variable field related to this equation
2270 * \param[in] eqp pointer to a cs_equation_param_t structure
2271 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
2272 * \param[in, out] context pointer to cs_cdovcb_scaleq_t structure
2273 */
2274 /*----------------------------------------------------------------------------*/
2275
2276 void
cs_cdovcb_scaleq_solve_theta(bool cur2prev,const cs_mesh_t * mesh,const int field_id,const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)2277 cs_cdovcb_scaleq_solve_theta(bool cur2prev,
2278 const cs_mesh_t *mesh,
2279 const int field_id,
2280 const cs_equation_param_t *eqp,
2281 cs_equation_builder_t *eqb,
2282 void *context)
2283 {
2284 cs_timer_t t0 = cs_timer_time();
2285
2286 const cs_cdo_connect_t *connect = cs_shared_connect;
2287 const cs_range_set_t *rs = connect->range_sets[CS_CDO_CONNECT_VTX_SCAL];
2288 const cs_cdo_quantities_t *quant = cs_shared_quant;
2289 const cs_lnum_t n_vertices = quant->n_vertices;
2290 const cs_time_step_t *ts = cs_shared_time_step;
2291 const cs_real_t tcoef = 1 - eqp->theta;
2292
2293 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
2294 cs_field_t *fld = cs_field_by_id(field_id);
2295
2296 assert(cs_equation_param_has_time(eqp) == true);
2297 assert(eqp->time_scheme == CS_TIME_SCHEME_THETA ||
2298 eqp->time_scheme == CS_TIME_SCHEME_CRANKNICO);
2299
2300 /* Build an array storing the Dirichlet values at vertices.
2301 * Always evaluated at t + dt */
2302
2303 _setup_vcb(ts->t_cur + ts->dt[0], mesh, connect, eqp, eqb, eqc->vtx_bc_flag);
2304
2305 /* Initialize the local system: matrix and rhs */
2306
2307 cs_matrix_t *matrix = cs_matrix_create(cs_shared_ms);
2308 cs_real_t *rhs = NULL;
2309 double rhs_norm = 0.;
2310
2311 BFT_MALLOC(rhs, n_vertices, cs_real_t);
2312 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
2313 for (cs_lnum_t i = 0; i < n_vertices; i++) rhs[i] = 0.0;
2314
2315 /* Initialize the structure to assemble values */
2316
2317 cs_matrix_assembler_values_t *mav
2318 = cs_matrix_assembler_values_init(matrix, NULL, NULL);
2319
2320 /* Detect the first call (in this case, we compute the initial source term) */
2321
2322 bool compute_initial_source = false;
2323 if (eqb->init_step) {
2324
2325 eqb->init_step = false;
2326 if (cs_equation_param_has_sourceterm(eqp))
2327 compute_initial_source = true; /* First iteration */
2328
2329 }
2330 else {
2331
2332 if (cs_equation_param_has_sourceterm(eqp)) {
2333
2334 /* Add contribution of the previous computed source term
2335 * Only vertices (and not cells) since there is an assembly process
2336 * on vertices (and not on cells) */
2337
2338 for (cs_lnum_t v = 0; v < n_vertices; v++)
2339 rhs[v] += tcoef * eqc->source_terms[v];
2340 memset(eqc->source_terms, 0, n_vertices * sizeof(cs_real_t));
2341
2342 if (eqp->default_enforcement == CS_PARAM_BC_ENFORCE_ALGEBRAIC ||
2343 eqp->default_enforcement == CS_PARAM_BC_ENFORCE_PENALIZED) {
2344
2345 assert(eqc->vtx_bc_flag != NULL);
2346 for (cs_lnum_t v = 0; v < n_vertices; v++) {
2347 if (cs_cdo_bc_is_dirichlet(eqc->vtx_bc_flag[v]))
2348 rhs[v] = 0.;
2349 }
2350
2351 } /* Algebraic or penalized enforcement is set */
2352 } /* At least one source term is defined */
2353
2354 } /* Not the first time step */
2355
2356 /* -------------------------
2357 * Main OpenMP block on cell
2358 * ------------------------- */
2359
2360 # pragma omp parallel if (quant->n_cells > CS_THR_MIN)
2361 {
2362 /* Set variables and structures inside the OMP section so that each thread
2363 has its own value */
2364
2365 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
2366 int t_id = omp_get_thread_num();
2367 #else
2368 int t_id = 0;
2369 #endif
2370
2371 /* Each thread get back its related structures:
2372 Get the cell-wise view of the mesh and the algebraic system */
2373
2374 cs_equation_assemble_t *eqa = cs_equation_assemble_get(t_id);
2375 cs_face_mesh_t *fm = cs_cdo_local_get_face_mesh(t_id);
2376 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
2377 cs_cell_sys_t *csys = _vcbs_cell_system[t_id];
2378 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
2379 cs_real_t *cell_sources = eqc->source_terms + quant->n_vertices;
2380 cs_hodge_t *diff_hodge =
2381 (eqc->diffusion_hodge == NULL) ? NULL : eqc->diffusion_hodge[t_id];
2382 cs_hodge_t *mass_hodge =
2383 (eqc->mass_hodge == NULL) ? NULL : eqc->mass_hodge[t_id];
2384
2385 const cs_real_t t_cur = ts->t_cur;
2386 const cs_real_t dt_cur = ts->dt[0];
2387 const cs_real_t inv_dtcur = 1/dt_cur;
2388
2389 /* Set times at which one evaluates quantities if needed.
2390 * Time_eval = (1-theta).t^n + theta.t^(n+1) = t^n + theta.dt
2391 * since t^(n+1) = t^n + dt */
2392
2393 cb->t_pty_eval = t_cur + eqp->theta*dt_cur;
2394 cb->t_bc_eval = t_cur + dt_cur;
2395 cb->t_st_eval = t_cur + dt_cur;
2396
2397 /* Initialization of the values of properties */
2398
2399 cs_equation_init_properties(eqp, eqb, diff_hodge, cb);
2400
2401 /* ---------------------------------------------
2402 * Main loop on cells to build the linear system
2403 * --------------------------------------------- */
2404
2405 # pragma omp for CS_CDO_OMP_SCHEDULE reduction(+:rhs_norm)
2406 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
2407
2408 /* Set the current cell flag */
2409
2410 cb->cell_flag = connect->cell_flag[c_id];
2411
2412 /* Set the local mesh structure for the current cell */
2413
2414 cs_cell_mesh_build(c_id,
2415 cs_equation_cell_mesh_flag(cb->cell_flag, eqb),
2416 connect, quant, cm);
2417
2418 /* Set the local (i.e. cellwise) structures for the current cell */
2419
2420 _svcb_init_cell_system(cm, eqp, eqb, eqc, eqc->vtx_bc_flag, fld->val,
2421 csys, cb);
2422
2423 /* Build and add the diffusion/advection/reaction term to the local
2424 system. A mass matrix is also built if needed. */
2425
2426 _svcb_conv_diff_reac(eqp, eqb, eqc, cm, fm,
2427 mass_hodge, diff_hodge, csys, cb);
2428
2429 if (cs_equation_param_has_sourceterm(eqp)) { /* SOURCE TERM
2430 * =========== */
2431
2432 if (compute_initial_source) { /* First time step */
2433
2434 /* Reset the local contribution */
2435
2436 memset(csys->source, 0, csys->n_dofs*sizeof(cs_real_t));
2437
2438 cs_source_term_compute_cellwise(eqp->n_source_terms,
2439 (cs_xdef_t *const *)eqp->source_terms,
2440 cm,
2441 eqb->source_mask,
2442 eqb->compute_source,
2443 t_cur,
2444 mass_hodge,
2445 cb,
2446 csys->source);
2447
2448 for (short int v = 0; v < cm->n_vc; v++)
2449 csys->rhs[v] += tcoef * csys->source[v];
2450 csys->rhs[cm->n_vc] += tcoef * csys->source[cm->n_vc];
2451
2452 }
2453 else { /* Add the contribution of the previous time step */
2454
2455 /* Contribution at vertices has been done before the loop on cells */
2456
2457 csys->rhs[cm->n_vc] += tcoef*cell_sources[cm->c_id];
2458
2459 }
2460
2461 /* Reset the local contribution */
2462
2463 memset(csys->source, 0, csys->n_dofs*sizeof(cs_real_t));
2464
2465 /* Source term contribution to the algebraic system
2466 If the equation is steady, the source term has already been computed
2467 and is added to the right-hand side during its initialization. */
2468
2469 cs_source_term_compute_cellwise(eqp->n_source_terms,
2470 (cs_xdef_t *const *)eqp->source_terms,
2471 cm,
2472 eqb->source_mask,
2473 eqb->compute_source,
2474 cb->t_st_eval,
2475 mass_hodge,
2476 cb,
2477 csys->source);
2478
2479 for (short int v = 0; v < cm->n_vc; v++)
2480 csys->rhs[v] += eqp->theta * csys->source[v];
2481 csys->rhs[cm->n_vc] += eqp->theta * csys->source[cm->n_vc];
2482
2483 } /* End of term source */
2484
2485 /* Apply boundary conditions (those which are weakly enforced) */
2486
2487 _svcb_apply_weak_bc(eqp, eqc, cm, fm, diff_hodge, csys, cb);
2488
2489 /* UNSTEADY TERM + TIME SCHEME
2490 * ===========================
2491 *
2492 * STEP.1 >> Compute the contribution of the "adr" to the RHS:
2493 * tcoef*adr_pn where adr_pn = csys->mat * p_n */
2494
2495 double *adr_pn = cb->values;
2496 cs_sdm_square_matvec(csys->mat, csys->val_n, adr_pn);
2497 for (short int i = 0; i < csys->n_dofs; i++)
2498 csys->rhs[i] -= tcoef*adr_pn[i];
2499
2500 /* STEP.2 >> Multiply csys->mat by theta */
2501
2502 for (int i = 0; i < csys->n_dofs*csys->n_dofs; i++)
2503 csys->mat->val[i] *= eqp->theta;
2504
2505 /* STEP.3 >> Handle the mass matrix
2506 * Two contributions for the mass matrix
2507 * a) add to csys->mat
2508 * b) add to rhs: mass_mat * p_n */
2509
2510 if (!(eqb->time_pty_uniform))
2511 cb->tpty_val = cs_property_value_in_cell(cm, eqp->time_property,
2512 cb->t_pty_eval);
2513
2514 if (eqb->sys_flag & CS_FLAG_SYS_TIME_DIAG) { /* Mass lumping */
2515
2516 /* |c|*wvc = |dual_cell(v) cap c| */
2517
2518 const double ptyc = cb->tpty_val * cm->vol_c * inv_dtcur;
2519
2520 /* STEPS >> Compute the time contribution to the RHS: Mtime*pn
2521 * >> Update the cellwise system with the time matrix */
2522
2523 for (short int i = 0; i < cm->n_vc; i++) {
2524
2525 const double dval = 0.75 * ptyc * cm->wvc[i];
2526
2527 /* Update the RHS with values at time t_n */
2528
2529 csys->rhs[i] += dval * csys->val_n[i];
2530
2531 /* Add the diagonal contribution from time matrix */
2532
2533 csys->mat->val[i*(csys->n_dofs + 1)] += dval;
2534
2535 }
2536
2537 /* Cell row */
2538
2539 const double dvalc = 0.25 * ptyc;
2540
2541 /* Update the RHS with values at time t_n */
2542
2543 csys->rhs[cm->n_vc] += dvalc * csys->val_n[cm->n_vc];
2544
2545 /* Add the diagonal contribution from time matrix */
2546
2547 csys->mat->val[cm->n_vc*(csys->n_dofs + 1)] += dvalc;
2548
2549 }
2550 else { /* Use the mass matrix */
2551
2552 const double tpty_coef = cb->tpty_val * inv_dtcur;
2553 const cs_sdm_t *mass_mat = mass_hodge->matrix;
2554
2555 /* Update rhs with mass_mat*p^n */
2556
2557 double *time_pn = cb->values;
2558 cs_sdm_square_matvec(mass_mat, csys->val_n, time_pn);
2559 for (short int i = 0; i < csys->n_dofs; i++)
2560 csys->rhs[i] += tpty_coef*time_pn[i];
2561
2562 /* Update the cellwise system with the time matrix */
2563
2564 cs_sdm_add_mult(csys->mat, tpty_coef, mass_mat);
2565
2566 }
2567
2568 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
2569 if (cs_dbg_cw_test(eqp, cm, csys))
2570 cs_cell_sys_dump("\n>> Cell system after adding time", csys);
2571 #endif
2572
2573 /* STATIC CONDENSATION
2574 * ===================
2575 * of the local system matrix of size n_vc + 1 into a matrix of size n_vc.
2576 * Store data in rc_tilda and acv_tilda to compute the values at cell
2577 * centers after solving the system */
2578
2579 cs_static_condensation_scalar_eq(connect->c2v,
2580 eqc->rc_tilda, eqc->acv_tilda,
2581 cb, csys);
2582
2583 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 1
2584 if (cs_dbg_cw_test(eqp, cm, csys))
2585 cs_cell_sys_dump(">> Cell system matrix after condensation", csys);
2586 #endif
2587
2588 /* Compute a cellwise norm of the RHS for the normalization of the
2589 residual during the resolution of the linear system */
2590
2591 rhs_norm += _svcb_cw_rhs_normalization(eqp->sles_param->resnorm_type,
2592 cm, csys);
2593
2594 /* Enforce values if needed (internal or Dirichlet) */
2595
2596 _svcb_enforce_values(eqp, eqb, eqc, cm, fm, diff_hodge, csys, cb);
2597
2598 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVCB_SCALEQ_DBG > 0
2599 if (cs_dbg_cw_test(eqp, cm, csys))
2600 cs_cell_sys_dump(">> (FINAL) Cell system matrix", csys);
2601 #endif
2602
2603 /* ASSEMBLY PROCESS
2604 * ================ */
2605
2606 _assemble(eqc, cm, csys, rs, eqa, mav, rhs);
2607
2608 } /* Main loop on cells */
2609
2610 } /* OpenMP Block */
2611
2612 cs_matrix_assembler_values_done(mav); /* optional */
2613
2614 /* Free temporary buffers and structures */
2615
2616 cs_equation_builder_reset(eqb);
2617 cs_matrix_assembler_values_finalize(&mav);
2618
2619 /* Copy current field values to previous values */
2620
2621 if (cur2prev)
2622 cs_field_current_to_previous(fld);
2623
2624 /* End of the system building */
2625
2626 cs_timer_t t1 = cs_timer_time();
2627 cs_timer_counter_add_diff(&(eqb->tcb), &t0, &t1);
2628
2629 /* Solve the linear system */
2630 /* ======================= */
2631
2632 /* Last step in the computation of the renormalization coefficient */
2633
2634 cs_equation_sync_rhs_normalization(eqp->sles_param->resnorm_type,
2635 n_vertices,
2636 rhs,
2637 &rhs_norm);
2638
2639 cs_sles_t *sles = cs_sles_find_or_add(eqp->sles_param->field_id, NULL);
2640
2641 cs_equation_solve_scalar_system(n_vertices,
2642 eqp->sles_param,
2643 matrix,
2644 rs,
2645 rhs_norm,
2646 true, /* rhs_redux */
2647 sles,
2648 fld->val,
2649 rhs);
2650
2651 cs_timer_t t2 = cs_timer_time();
2652 cs_timer_counter_add_diff(&(eqb->tcs), &t1, &t2);
2653
2654 /* Update fields */
2655
2656 _update_cell_fields(&(eqb->tce), fld, eqc, cur2prev);
2657
2658 /* Free remaining buffers */
2659
2660 BFT_FREE(rhs);
2661 cs_sles_free(sles);
2662 cs_matrix_destroy(&matrix);
2663 }
2664
2665 /*----------------------------------------------------------------------------*/
2666 /*!
2667 * \brief Retrieve an array of values at mesh vertices for the variable field
2668 * associated to the given context
2669 * The lifecycle of this array is managed by the code. So one does not
2670 * have to free the return pointer.
2671 *
2672 * \param[in, out] context pointer to a data structure cast on-the-fly
2673 * \param[in] previous retrieve the previous state (true/false)
2674 *
2675 * \return a pointer to an array of cs_real_t (size: n_vertices)
2676 */
2677 /*----------------------------------------------------------------------------*/
2678
2679 cs_real_t *
cs_cdovcb_scaleq_get_vertex_values(void * context,bool previous)2680 cs_cdovcb_scaleq_get_vertex_values(void *context,
2681 bool previous)
2682 {
2683 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
2684
2685 if (eqc == NULL)
2686 return NULL;
2687
2688 cs_field_t *pot = cs_field_by_id(eqc->var_field_id);
2689
2690 if (previous)
2691 return pot->val_pre;
2692 else
2693 return pot->val;
2694 }
2695
2696 /*----------------------------------------------------------------------------*/
2697 /*!
2698 * \brief Get the computed values at mesh cells from the inverse operation
2699 * w.r.t. the static condensation (DoF used in the linear system are
2700 * located at primal vertices and field related to the structure
2701 * equation is also attached to primal vertices)
2702 * The lifecycle of this array is managed by the code. So one does not
2703 * have to free the return pointer.
2704 *
2705 * \param[in, out] context pointer to a data structure cast on-the-fly
2706 * \param[in] previous retrieve the previous state (true/false)
2707 *
2708 * \return a pointer to an array of cs_real_t (size: n_cells)
2709 */
2710 /*----------------------------------------------------------------------------*/
2711
2712 cs_real_t *
cs_cdovcb_scaleq_get_cell_values(void * context,bool previous)2713 cs_cdovcb_scaleq_get_cell_values(void *context,
2714 bool previous)
2715 {
2716 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
2717
2718 if (eqc == NULL)
2719 return NULL;
2720
2721 if (previous)
2722 return eqc->cell_values_pre;
2723 else
2724 return eqc->cell_values;
2725 }
2726
2727 /*----------------------------------------------------------------------------*/
2728 /*!
2729 * \brief Compute for each vertex of a boundary face, the portion of diffusive
2730 * flux across the boundary face. The surface attached to each vertex
2731 * corresponds to the intersection of its dual cell (associated to
2732 * a vertex of the face) with the face.
2733 * Case of scalar-valued CDO-VCb schemes
2734 *
2735 * \param[in] t_eval time at which one performs the evaluation
2736 * \param[in] eqp pointer to a cs_equation_param_t structure
2737 * \param[in] pot_v pointer to an array of field values at vertices
2738 * \param[in] pot_c pointer to an array of field values at cells
2739 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
2740 * \param[in, out] context pointer to a scheme builder structure
2741 * \param[in, out] vf_flux pointer to the values of the diffusive flux
2742 */
2743 /*----------------------------------------------------------------------------*/
2744
2745 void
cs_cdovcb_scaleq_boundary_diff_flux(const cs_real_t t_eval,const cs_equation_param_t * eqp,const cs_real_t * pot_v,const cs_real_t * pot_c,cs_equation_builder_t * eqb,void * context,cs_real_t * vf_flux)2746 cs_cdovcb_scaleq_boundary_diff_flux(const cs_real_t t_eval,
2747 const cs_equation_param_t *eqp,
2748 const cs_real_t *pot_v,
2749 const cs_real_t *pot_c,
2750 cs_equation_builder_t *eqb,
2751 void *context,
2752 cs_real_t *vf_flux)
2753 {
2754 if (vf_flux == NULL)
2755 return;
2756
2757 cs_timer_t t0 = cs_timer_time();
2758
2759 const cs_cdo_quantities_t *quant = cs_shared_quant;
2760 const cs_cdo_connect_t *connect = cs_shared_connect;
2761
2762 if (cs_equation_param_has_diffusion(eqp) == false) {
2763 memset(vf_flux, 0, connect->bf2v->idx[quant->n_b_faces]*sizeof(cs_real_t));
2764
2765 cs_timer_t t1 = cs_timer_time();
2766 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
2767 return;
2768 }
2769
2770 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
2771
2772 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
2773 shared(quant, connect, eqp, eqb, eqc, vf_flux, pot_v, pot_c, \
2774 _vcbs_cell_builder) \
2775 firstprivate(t_eval)
2776 {
2777 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
2778 int t_id = omp_get_thread_num();
2779 #else
2780 int t_id = 0;
2781 #endif
2782
2783 const cs_cdo_bc_face_t *face_bc = eqb->face_bc;
2784 const cs_adjacency_t *bf2v = connect->bf2v;
2785 const cs_adjacency_t *f2c = connect->f2c;
2786 const cs_lnum_t fidx_shift = f2c->idx[quant->n_i_faces];
2787
2788 /* Set inside the OMP section so that each thread has its own value */
2789
2790 double *tmp = cs_cdo_local_get_d_buffer(t_id);
2791 cs_real_t *pot = tmp,
2792 *flux = tmp + connect->n_max_vbyc + 1;
2793
2794 /* Each thread get back its related structures:
2795 Get the cellwise view of the mesh and the algebraic system */
2796
2797 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
2798 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
2799 cs_hodge_t *hodge = eqc->diffusion_hodge[t_id];
2800 cs_property_data_t *diff_pty = hodge->pty_data;
2801
2802 /* msh_flag for Neumann and Robin BCs. Add add_flag for the other cases
2803 when one has to reconstruct a flux */
2804
2805 cs_eflag_t msh_flag = CS_FLAG_COMP_PV | CS_FLAG_COMP_FV;
2806 cs_eflag_t add_flag = CS_FLAG_COMP_EV | CS_FLAG_COMP_FE |
2807 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PVQ |
2808 CS_FLAG_COMP_DEQ | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_HFQ;
2809
2810 if (eqb->diff_pty_uniform) /* c_id = 0 */
2811 cs_property_get_cell_tensor(0, t_eval,
2812 eqp->diffusion_property,
2813 hodge->param->inv_pty,
2814 diff_pty->tensor);
2815
2816 # pragma omp for CS_CDO_OMP_SCHEDULE
2817 for (cs_lnum_t bf_id = 0; bf_id < quant->n_b_faces; bf_id++) {
2818
2819 const cs_lnum_t f_id = bf_id + quant->n_i_faces;
2820 const cs_lnum_t c_id = f2c->ids[bf_id + fidx_shift];
2821 const cs_lnum_t *idx = bf2v->idx + bf_id;
2822
2823 cs_real_t *_flx = vf_flux + idx[0];
2824
2825 switch (face_bc->flag[bf_id]) {
2826
2827 case CS_CDO_BC_HMG_NEUMANN:
2828 memset(_flx, 0, (idx[1]-idx[0])*sizeof(cs_real_t));
2829 break;
2830
2831 case CS_CDO_BC_NEUMANN:
2832 {
2833 cs_real_t *neu_values = cb->values;
2834
2835 /* Set the local mesh structure for the current cell */
2836
2837 cs_cell_mesh_build(c_id, msh_flag, connect, quant, cm);
2838
2839 const short int f = cs_cell_mesh_get_f(f_id, cm);
2840
2841 cs_equation_compute_neumann_sv(t_eval,
2842 face_bc->def_ids[bf_id],
2843 f,
2844 eqp,
2845 cm,
2846 neu_values);
2847
2848 short int n_vf = 0;
2849 for (int i = cm->f2v_idx[f]; i < cm->f2v_idx[f+1]; i++)
2850 _flx[n_vf++] = neu_values[cm->f2v_ids[i]];
2851
2852 }
2853 break;
2854
2855 case CS_CDO_BC_ROBIN:
2856 {
2857 cs_real_t *robin_values = cb->values;
2858 cs_real_t *wvf = cb->values + 3;
2859
2860 /* Set the local mesh structure for the current cell */
2861
2862 cs_cell_mesh_build(c_id, msh_flag, connect, quant, cm);
2863
2864 const short int f = cs_cell_mesh_get_f(f_id, cm);
2865 const cs_real_t f_area = quant->b_face_surf[bf_id];
2866
2867 /* Robin BC expression: K du/dn + alpha*(p - p0) = g */
2868
2869 cs_equation_compute_robin(t_eval,
2870 face_bc->def_ids[bf_id],
2871 f,
2872 eqp,
2873 cm,
2874 robin_values);
2875
2876 const cs_real_t alpha = robin_values[0];
2877 const cs_real_t p0 = robin_values[1];
2878 const cs_real_t g = robin_values[2];
2879
2880 cs_cdo_quantities_compute_b_wvf(connect, quant, bf_id, wvf);
2881
2882 short int n_vf = 0;
2883 for (int i = cm->f2v_idx[f]; i < cm->f2v_idx[f+1]; i++) {
2884 const cs_real_t pv = pot_v[cm->v_ids[cm->f2v_ids[i]]];
2885 _flx[n_vf] = f_area*wvf[n_vf]*(alpha*(p0 - pv) + g);
2886 n_vf++;
2887 }
2888 }
2889 break;
2890
2891 default:
2892 { /* Reconstruct a normal flux at the boundary face */
2893
2894 /* Set the local mesh structure for the current cell */
2895
2896 cs_cell_mesh_build(c_id, msh_flag | add_flag, connect, quant, cm);
2897
2898 const short int f = cs_cell_mesh_get_f(f_id, cm);
2899
2900 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDOVB_SCALEQ_DBG > 2
2901 if (cs_dbg_cw_test(eqp, cm, NULL)) cs_cell_mesh_dump(cm);
2902 #endif
2903 if (!eqb->diff_pty_uniform)
2904 cs_property_tensor_in_cell(cm,
2905 eqp->diffusion_property,
2906 t_eval,
2907 hodge->param->inv_pty,
2908 diff_pty->tensor);
2909
2910 /* Define a local buffer keeping the value of the discrete potential
2911 for the current cell */
2912
2913 for (short int v = 0; v < cm->n_vc; v++)
2914 pot[v] = pot_v[cm->v_ids[v]];
2915 pot[cm->n_vc] = pot_c[c_id];
2916
2917 cs_cdo_diffusion_wbs_vbyf_flux(f, cm, pot, hodge, cb, flux);
2918
2919 /* Fill the global flux array */
2920
2921 short int n_vf = 0;
2922 for (int i = cm->f2v_idx[f]; i < cm->f2v_idx[f+1]; i++)
2923 _flx[n_vf++] = flux[cm->f2v_ids[i]];
2924
2925 }
2926 break;
2927
2928 } /* End of switch */
2929
2930 } /* End of loop on boundary faces */
2931
2932 } /* End of Open block */
2933
2934 cs_timer_t t1 = cs_timer_time();
2935 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
2936 }
2937
2938 /*----------------------------------------------------------------------------*/
2939 /*!
2940 * \brief Compute the diffusive and convective flux across a list of faces
2941 * Case of scalar-valued CDO-VCb schemes
2942 *
2943 * \param[in] normal indicate in which direction flux is > 0
2944 * \param[in] pdi pointer to an array of field values
2945 * \param[in] eqp pointer to a cs_equation_param_t structure
2946 * \param[in] ml_id id related to a cs_mesh_location_t struct.
2947 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
2948 * \param[in, out] context pointer to data specific for this scheme
2949 * \param[in, out] d_flux pointer to the value of the diffusive flux
2950 * \param[in, out] c_flux pointer to the value of the convective flux
2951 */
2952 /*----------------------------------------------------------------------------*/
2953
2954 void
cs_cdovcb_scaleq_flux_across_plane(const cs_real_t normal[],const cs_real_t * pdi,const cs_equation_param_t * eqp,int ml_id,cs_equation_builder_t * eqb,void * context,double * d_flux,double * c_flux)2955 cs_cdovcb_scaleq_flux_across_plane(const cs_real_t normal[],
2956 const cs_real_t *pdi,
2957 const cs_equation_param_t *eqp,
2958 int ml_id,
2959 cs_equation_builder_t *eqb,
2960 void *context,
2961 double *d_flux,
2962 double *c_flux)
2963 {
2964 *d_flux = 0.;
2965 *c_flux = 0.;
2966
2967 if (pdi == NULL)
2968 return;
2969
2970 cs_mesh_location_type_t ml_t = cs_mesh_location_get_type(ml_id);
2971
2972 if (ml_t != CS_MESH_LOCATION_INTERIOR_FACES &&
2973 ml_t != CS_MESH_LOCATION_BOUNDARY_FACES) {
2974 cs_base_warn(__FILE__, __LINE__);
2975 cs_log_printf(CS_LOG_DEFAULT,
2976 _(" Mesh location type is incompatible with the computation\n"
2977 " of the flux across faces.\n"));
2978 return;
2979 }
2980
2981 const cs_timer_t t0 = cs_timer_time();
2982 const cs_lnum_t *n_elts = cs_mesh_location_get_n_elts(ml_id);
2983 const cs_lnum_t *elt_ids = cs_mesh_location_get_elt_list(ml_id);
2984
2985 if (n_elts[0] > 0 && elt_ids == NULL)
2986 bft_error(__FILE__, __LINE__, 0,
2987 _(" Computing the flux across all interior or border faces is not"
2988 " managed yet."));
2989
2990 const cs_cdo_connect_t *connect = cs_shared_connect;
2991 const cs_adjacency_t *f2c = connect->f2c;
2992 const cs_cdo_quantities_t *quant = cs_shared_quant;
2993 const cs_real_t t_cur = cs_shared_time_step->t_cur;
2994
2995 double flx, p_f;
2996 cs_real_33_t pty_tens;
2997 cs_nvec3_t adv_c;
2998
2999 /* To be modified for an fully integration of openMP */
3000
3001 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
3002 cs_face_mesh_t *fm = cs_cdo_local_get_face_mesh(0);
3003 cs_cell_builder_t *cb = _vcbs_cell_builder[0];
3004
3005 /* Set inside the OMP section so that each thread has its own value */
3006
3007 double *p_v = cs_cdo_local_get_d_buffer(0);
3008
3009 if (ml_t == CS_MESH_LOCATION_BOUNDARY_FACES) {
3010
3011 const cs_lnum_t n_i_faces = connect->n_faces[CS_INT_FACES];
3012 const cs_lnum_t *cell_ids = f2c->ids + f2c->idx[n_i_faces];
3013
3014 for (cs_lnum_t id = 0; id < n_elts[0]; id++) { /* Loop on selected faces */
3015
3016 const cs_lnum_t bf_id = elt_ids[id];
3017 const cs_lnum_t f_id = n_i_faces + bf_id;
3018 const cs_lnum_t c_id = cell_ids[bf_id];
3019
3020 /* Build a face-wise view of the mesh */
3021
3022 cs_face_mesh_build(c_id, f_id, connect, quant, fm);
3023
3024 const short int sgn = (_dp3(fm->face.unitv, normal) < 0) ? -1 : 1;
3025
3026 /* Store values of the potential related to this face */
3027
3028 for (short int v = 0; v < fm->n_vf; v++)
3029 p_v[v] = pdi[fm->v_ids[v]];
3030
3031 /* Interpolate a value at the face center */
3032
3033 p_f = cs_reco_fw_scalar_pv_at_face_center(fm, p_v);
3034
3035 /* Compute the local diffusive flux */
3036
3037 if (cs_equation_param_has_diffusion(eqp)) {
3038
3039 cs_property_get_cell_tensor(c_id, t_cur,
3040 eqp->diffusion_property,
3041 eqp->diffusion_hodgep.inv_pty,
3042 pty_tens);
3043
3044 flx = cs_cdo_diffusion_wbs_face_flux(fm,
3045 (const cs_real_3_t (*))pty_tens,
3046 p_v, p_f, eqc->cell_values[c_id],
3047 cb);
3048 *d_flux += sgn * flx;
3049
3050 } /* Diffusive flux */
3051
3052 /* Compute the local advective flux */
3053
3054 if (cs_equation_param_has_convection(eqp)) {
3055
3056 const double coef = sgn * fm->face.meas * p_f;
3057
3058 cs_advection_field_get_cell_vector(c_id, eqp->adv_field, &adv_c);
3059 *c_flux += coef * adv_c.meas * _dp3(adv_c.unitv, fm->face.unitv);
3060
3061 }
3062
3063 } /* Loop on selected border faces */
3064
3065 }
3066 else if (ml_t == CS_MESH_LOCATION_INTERIOR_FACES) {
3067
3068 for (cs_lnum_t i = 0; i < n_elts[0]; i++) {
3069
3070 const cs_lnum_t f_id = elt_ids[i];
3071
3072 for (cs_lnum_t j = f2c->idx[f_id]; j < f2c->idx[f_id+1]; j++) {
3073
3074 const cs_lnum_t c_id = f2c->ids[j];
3075
3076 /* Build a face-wise view of the mesh */
3077
3078 cs_face_mesh_build(c_id, f_id, connect, quant, fm);
3079
3080 const short int sgn = (_dp3(fm->face.unitv, normal) < 0) ? -1 : 1;
3081
3082 /* Store values related to this face */
3083
3084 for (short int v = 0; v < fm->n_vf; v++)
3085 p_v[v] = pdi[fm->v_ids[v]];
3086
3087 p_f = cs_reco_fw_scalar_pv_at_face_center(fm, p_v);
3088
3089 /* Compute the diffusive flux seen from cell c1 */
3090
3091 if (cs_equation_param_has_diffusion(eqp)) {
3092
3093 cs_property_get_cell_tensor(c_id, t_cur,
3094 eqp->diffusion_property,
3095 eqp->diffusion_hodgep.inv_pty,
3096 pty_tens);
3097
3098 flx = cs_cdo_diffusion_wbs_face_flux(fm,
3099 (const cs_real_3_t (*))pty_tens,
3100 p_v, p_f, eqc->cell_values[c_id],
3101 cb);
3102
3103 *d_flux += sgn * 0.5 * flx; /* Centered approximation */
3104
3105 } /* Diffusive flux */
3106
3107 if (cs_equation_param_has_convection(eqp)) {
3108
3109 /* Compute the local advective flux seen from cell */
3110
3111 cs_advection_field_get_cell_vector(c_id,
3112 eqp->adv_field,
3113 &adv_c);
3114 flx = adv_c.meas * _dp3(adv_c.unitv, fm->face.unitv);
3115
3116 /* Centered approximation of the advective flux */
3117
3118 *c_flux += sgn * 0.5 * flx * p_f * fm->face.meas;
3119
3120 } /* Advective flux */
3121
3122 } /* Loop on cells attached to this interior face */
3123
3124 } /* Loop on selected interior faces */
3125
3126 } /* Set of interior or border faces */
3127
3128 cs_timer_t t1 = cs_timer_time();
3129 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
3130 }
3131
3132 /*----------------------------------------------------------------------------*/
3133 /*!
3134 * \brief Cellwise computation of the diffusive flux in each cells.
3135 * Case of scalar-valued CDO-VCb schemes
3136 *
3137 * \param[in] values discrete values for the potential
3138 * \param[in] eqp pointer to a cs_equation_param_t structure
3139 * \param[in] t_eval time at which one performs the evaluation
3140 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
3141 * \param[in, out] context pointer to data structure cast on-the-fly
3142 * \param[in, out] diff_flux value of the diffusive flux
3143 */
3144 /*----------------------------------------------------------------------------*/
3145
3146 void
cs_cdovcb_scaleq_diff_flux_in_cells(const cs_real_t * values,const cs_equation_param_t * eqp,cs_real_t t_eval,cs_equation_builder_t * eqb,void * context,cs_real_t * diff_flux)3147 cs_cdovcb_scaleq_diff_flux_in_cells(const cs_real_t *values,
3148 const cs_equation_param_t *eqp,
3149 cs_real_t t_eval,
3150 cs_equation_builder_t *eqb,
3151 void *context,
3152 cs_real_t *diff_flux)
3153 {
3154 if (diff_flux == NULL)
3155 return ;
3156
3157 const cs_cdo_quantities_t *quant = cs_shared_quant;
3158 const cs_cdo_connect_t *connect = cs_shared_connect;
3159
3160 if (cs_equation_param_has_diffusion(eqp) == false) {
3161 memset(diff_flux, 0, 3*quant->n_cells);
3162 return;
3163 }
3164
3165 assert(eqp->diffusion_hodgep.algo == CS_HODGE_ALGO_WBS);
3166
3167 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
3168
3169 cs_timer_t t0 = cs_timer_time();
3170
3171 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
3172 shared(quant, connect, eqp, eqb, eqc, diff_flux, values, \
3173 t_eval, _vcbs_cell_builder)
3174 {
3175 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
3176 int t_id = omp_get_thread_num();
3177 #else
3178 int t_id = 0;
3179 #endif
3180
3181 cs_eflag_t msh_flag = CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ |
3182 CS_FLAG_COMP_HFQ | CS_FLAG_COMP_DEQ | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV;
3183
3184 /* Each thread get back its related structures:
3185 Get the cellwise view of the mesh and the algebraic system */
3186
3187 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
3188 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
3189 cs_hodge_t *hodge = eqc->diffusion_hodge[t_id];
3190
3191 /* Set inside the OMP section so that each thread has its own value */
3192
3193 double *pot = cs_cdo_local_get_d_buffer(t_id);
3194
3195 cb->t_pty_eval = t_eval;
3196 if (eqb->diff_pty_uniform) /* cell_id =0, cell_flag = 0 */
3197 cs_hodge_set_property_value(0, cb->t_pty_eval, 0, hodge);
3198
3199 /* Define the flux by cellwise contributions */
3200
3201 # pragma omp for CS_CDO_OMP_SCHEDULE
3202 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
3203
3204 cb->cell_flag = 0; /* No need to make a distincition between cells */
3205
3206 /* Set the local mesh structure for the current cell */
3207
3208 cs_cell_mesh_build(c_id, msh_flag, connect, quant, cm);
3209
3210 if (!eqb->diff_pty_uniform) /* cell_flag is always set to 0 */
3211 cs_hodge_set_property_value_cw(cm, cb->t_pty_eval, cb->cell_flag,
3212 hodge);
3213
3214 /* Define a local buffer keeping the value of the discrete potential
3215 for the current cell */
3216
3217 for (short int v = 0; v < cm->n_vc; v++)
3218 pot[v] = values[cm->v_ids[v]];
3219 pot[cm->n_vc] = eqc->cell_values[c_id];
3220
3221 cs_cdo_diffusion_wbs_get_cell_flux(cm, pot, hodge,
3222 cb, diff_flux + 3*c_id);
3223
3224 } /* Loop on cells */
3225
3226 } /* OMP Section */
3227
3228 cs_timer_t t1 = cs_timer_time();
3229 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
3230 }
3231
3232 /*----------------------------------------------------------------------------*/
3233 /*!
3234 * \brief Cellwise computation of the diffusive flux across dual faces
3235 * Case of scalar-valued CDO-VCb schemes
3236 *
3237 * \param[in] values discrete values for the potential
3238 * \param[in] eqp pointer to a cs_equation_param_t structure
3239 * \param[in] t_eval time at which one performs the evaluation
3240 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
3241 * \param[in, out] context pointer to data structure cast on-the-fly
3242 * \param[in, out] diff_flux value of the diffusive flux
3243 */
3244 /*----------------------------------------------------------------------------*/
3245
3246 void
cs_cdovcb_scaleq_diff_flux_dfaces(const cs_real_t * values,const cs_equation_param_t * eqp,cs_real_t t_eval,cs_equation_builder_t * eqb,void * context,cs_real_t * diff_flux)3247 cs_cdovcb_scaleq_diff_flux_dfaces(const cs_real_t *values,
3248 const cs_equation_param_t *eqp,
3249 cs_real_t t_eval,
3250 cs_equation_builder_t *eqb,
3251 void *context,
3252 cs_real_t *diff_flux)
3253 {
3254 if (diff_flux == NULL)
3255 return ;
3256
3257 const cs_cdo_quantities_t *quant = cs_shared_quant;
3258 const cs_cdo_connect_t *connect = cs_shared_connect;
3259
3260 if (cs_equation_param_has_diffusion(eqp) == false) {
3261 memset(diff_flux, 0, connect->c2e->idx[quant->n_cells]*sizeof(cs_real_t));
3262 return;
3263 }
3264
3265 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
3266
3267 assert(eqp->diffusion_hodgep.algo == CS_HODGE_ALGO_WBS);
3268 assert(eqc->diffusion_hodge != NULL);
3269
3270 cs_timer_t t0 = cs_timer_time();
3271
3272 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
3273 shared(quant, connect, eqp, eqb, eqc, diff_flux, values, \
3274 t_eval, _vcbs_cell_builder)
3275 {
3276 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
3277 int t_id = omp_get_thread_num();
3278 #else
3279 int t_id = 0;
3280 #endif
3281
3282 cs_eflag_t msh_flag = CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ |
3283 CS_FLAG_COMP_SEF | CS_FLAG_COMP_DEQ | CS_FLAG_COMP_FEQ |
3284 CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ;
3285
3286 /* Each thread get back its related structures:
3287 Get the cellwise view of the mesh and the algebraic system as well
3288 as the Hodge operator for the diffusion term */
3289
3290 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
3291 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
3292 cs_hodge_t *hodge = eqc->diffusion_hodge[t_id];
3293
3294 /* Set inside the OMP section so that each thread has its own value */
3295
3296 double *pot = cs_cdo_local_get_d_buffer(t_id);
3297
3298 cb->t_pty_eval = t_eval;
3299 if (eqb->diff_pty_uniform) /* cell_id =0, cell_flag = 0 */
3300 cs_hodge_set_property_value(0, cb->t_pty_eval, 0, hodge);
3301
3302 /* Define the flux by cellwise contributions */
3303
3304 # pragma omp for CS_CDO_OMP_SCHEDULE
3305 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
3306
3307 cb->cell_flag = 0; /* No need to make a distincition between cells */
3308
3309 /* Set the local mesh structure for the current cell */
3310
3311 cs_cell_mesh_build(c_id, msh_flag, connect, quant, cm);
3312
3313 if (!eqb->diff_pty_uniform) /* cell_flag is always set to 0 */
3314 cs_hodge_set_property_value_cw(cm, cb->t_pty_eval, cb->cell_flag,
3315 hodge);
3316
3317 /* Define a local buffer keeping the value of the discrete potential
3318 for the current cell */
3319
3320 for (short int v = 0; v < cm->n_vc; v++)
3321 pot[v] = values[cm->v_ids[v]];
3322 pot[cm->n_vc] = eqc->cell_values[c_id];
3323
3324 /* Only hodge->pty_data is used */
3325
3326 cs_cdo_diffusion_wbs_get_dfbyc_flux(cm, pot, hodge,
3327 cb,
3328 diff_flux + connect->c2e->idx[c_id]);
3329
3330 } /* Loop on cells */
3331
3332 } /* OMP Section */
3333
3334 cs_timer_t t1 = cs_timer_time();
3335 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
3336 }
3337
3338 /*----------------------------------------------------------------------------*/
3339 /*!
3340 * \brief Cellwise computation of the discrete gradient at vertices
3341 *
3342 * \param[in] v_values discrete values for the potential at vertices
3343 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
3344 * \param[in, out] context pointer to data structure cast on-the-fly
3345 * \param[in, out] v_gradient gradient at vertices
3346 */
3347 /*----------------------------------------------------------------------------*/
3348
3349 void
cs_cdovcb_scaleq_vtx_gradient(const cs_real_t * v_values,cs_equation_builder_t * eqb,void * context,cs_real_t * v_gradient)3350 cs_cdovcb_scaleq_vtx_gradient(const cs_real_t *v_values,
3351 cs_equation_builder_t *eqb,
3352 void *context,
3353 cs_real_t *v_gradient)
3354 {
3355 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
3356
3357 const cs_cdo_quantities_t *quant = cs_shared_quant;
3358 const cs_cdo_connect_t *connect = cs_shared_connect;
3359
3360 if (v_gradient == NULL)
3361 bft_error(__FILE__, __LINE__, 0,
3362 " Result array has to be allocated prior to the call.");
3363
3364 cs_real_t *dualcell_vol = NULL;
3365 BFT_MALLOC(dualcell_vol, quant->n_vertices, cs_real_t);
3366
3367 # pragma omp parallel for if (3*quant->n_vertices > CS_THR_MIN)
3368 for (cs_lnum_t i = 0; i < 3*quant->n_vertices; i++)
3369 v_gradient[i] = 0;
3370 # pragma omp parallel for if (quant->n_vertices > CS_THR_MIN)
3371 for (cs_lnum_t i = 0; i < quant->n_vertices; i++)
3372 dualcell_vol[i] = 0;
3373
3374 cs_timer_t t0 = cs_timer_time();
3375
3376 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
3377 shared(quant, connect, eqc, v_gradient, v_values, dualcell_vol, \
3378 _vcbs_cell_builder, cs_glob_n_ranks)
3379 {
3380 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
3381 int t_id = omp_get_thread_num();
3382 #else
3383 int t_id = 0;
3384 #endif
3385 cs_real_3_t cgrd;
3386
3387 /* Set inside the OMP section so that each thread has its own value */
3388
3389 double *pot = cs_cdo_local_get_d_buffer(t_id);
3390
3391 cs_eflag_t msh_flag = CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ |
3392 CS_FLAG_COMP_DEQ | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ;
3393
3394 /* Each thread get back its related structures:
3395 Get the cellwise view of the mesh and the algebraic system */
3396
3397 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
3398 cs_cell_builder_t *cb = _vcbs_cell_builder[t_id];
3399
3400 /* Define the flux by cellwise contributions */
3401
3402 # pragma omp for CS_CDO_OMP_SCHEDULE
3403 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
3404
3405 /* Set the local mesh structure for the current cell */
3406
3407 cs_cell_mesh_build(c_id, msh_flag, connect, quant, cm);
3408
3409 /* Define a local buffer keeping the value of the discrete potential
3410 for the current cell */
3411
3412 for (short int v = 0; v < cm->n_vc; v++)
3413 pot[v] = v_values[cm->v_ids[v]];
3414 pot[cm->n_vc] = eqc->cell_values[c_id];
3415
3416 cs_reco_cw_cgrd_wbs_from_pvc(cm, pot, cb, cgrd);
3417
3418 for (short int v = 0; v < cm->n_vc; v++) {
3419 const double dvol = cm->wvc[v] * cm->vol_c;
3420 # pragma omp atomic
3421 dualcell_vol[cm->v_ids[v]] += dvol;
3422 for (int k = 0; k < 3; k++)
3423 # pragma omp atomic
3424 v_gradient[3*cm->v_ids[v] + k] += dvol*cgrd[k];
3425 }
3426
3427 } /* Loop on cells */
3428
3429 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL) {
3430
3431 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
3432 connect->n_vertices,
3433 1,
3434 true, /* interlace */
3435 CS_REAL_TYPE,
3436 dualcell_vol);
3437
3438 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
3439 connect->n_vertices,
3440 3,
3441 true, /* interlace */
3442 CS_REAL_TYPE,
3443 v_gradient);
3444 }
3445
3446 # pragma omp for CS_CDO_OMP_SCHEDULE
3447 for (cs_lnum_t i = 0; i < quant->n_vertices; i++) {
3448 cs_real_t inv_dualcell_vol = 1/dualcell_vol[i];
3449 for (int k = 0; k < 3; k++)
3450 v_gradient[3*i + k] *= inv_dualcell_vol;
3451 }
3452
3453 } /* OMP Section */
3454
3455 BFT_FREE(dualcell_vol);
3456
3457 cs_timer_t t1 = cs_timer_time();
3458 cs_timer_counter_add_diff(&(eqb->tce), &t0, &t1);
3459 }
3460
3461 /*----------------------------------------------------------------------------*/
3462 /*!
3463 * \brief Read additional arrays (not defined as fields) but useful for the
3464 * checkpoint/restart process
3465 *
3466 * \param[in, out] restart pointer to \ref cs_restart_t structure
3467 * \param[in] eqname name of the related equation
3468 * \param[in] scheme_context pointer to a data structure cast on-the-fly
3469 */
3470 /*----------------------------------------------------------------------------*/
3471
3472 void
cs_cdovcb_scaleq_read_restart(cs_restart_t * restart,const char * eqname,void * scheme_context)3473 cs_cdovcb_scaleq_read_restart(cs_restart_t *restart,
3474 const char *eqname,
3475 void *scheme_context)
3476 {
3477 /* Only the cell values are handled. Vertex values are stored in a cs_field_t
3478 structure and thus are handled automatically. */
3479
3480 if (restart == NULL)
3481 return;
3482 if (eqname == NULL)
3483 bft_error(__FILE__, __LINE__, 0, " %s: Name is NULL", __func__);
3484 if (scheme_context == NULL)
3485 bft_error(__FILE__, __LINE__, 0, " %s: Scheme context is NULL", __func__);
3486
3487 int retcode = CS_RESTART_SUCCESS;
3488 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)scheme_context;
3489
3490 const int cell_loc_id = cs_mesh_location_get_id_by_name("cells");
3491
3492 /* Define the section name */
3493
3494 char sec_name[128];
3495 snprintf(sec_name, 127, "%s::cell_vals", eqname);
3496
3497 /* Check section */
3498
3499 retcode = cs_restart_check_section(restart,
3500 sec_name,
3501 cell_loc_id,
3502 1, /* scalar-valued */
3503 CS_TYPE_cs_real_t);
3504
3505 if (retcode == CS_RESTART_SUCCESS)
3506 retcode = cs_restart_read_section(restart,
3507 sec_name,
3508 cell_loc_id,
3509 1, /* scalar-valued */
3510 CS_TYPE_cs_real_t,
3511 eqc->cell_values);
3512 }
3513
3514 /*----------------------------------------------------------------------------*/
3515 /*!
3516 * \brief Write additional arrays (not defined as fields) but useful for the
3517 * checkpoint/restart process
3518 *
3519 * \param[in, out] restart pointer to \ref cs_restart_t structure
3520 * \param[in] eqname name of the related equation
3521 * \param[in] scheme_context pointer to a data structure cast on-the-fly
3522 */
3523 /*----------------------------------------------------------------------------*/
3524
3525 void
cs_cdovcb_scaleq_write_restart(cs_restart_t * restart,const char * eqname,void * scheme_context)3526 cs_cdovcb_scaleq_write_restart(cs_restart_t *restart,
3527 const char *eqname,
3528 void *scheme_context)
3529 {
3530 /* Only the cell values are handled. Vertex values are stored in a cs_field_t
3531 structure and thus are handled automatically. */
3532
3533 if (restart == NULL)
3534 return;
3535 if (eqname == NULL)
3536 bft_error(__FILE__, __LINE__, 0, " %s: Name is NULL", __func__);
3537
3538 const cs_cdovcb_scaleq_t *eqc = (const cs_cdovcb_scaleq_t *)scheme_context;
3539 const int cell_loc_id = cs_mesh_location_get_id_by_name("cells");
3540
3541 /* Define the section name */
3542
3543 char sec_name[128];
3544 snprintf(sec_name, 127, "%s::cell_vals", eqname);
3545
3546 /* Write section */
3547
3548 cs_restart_write_section(restart,
3549 sec_name,
3550 cell_loc_id,
3551 1, /* scalar-valued */
3552 CS_TYPE_cs_real_t,
3553 eqc->cell_values);
3554 }
3555
3556 /*----------------------------------------------------------------------------*/
3557 /*!
3558 * \brief Operate a current to previous operation for the field associated to
3559 * this equation and potentially for related fields/arrays.
3560 *
3561 * \param[in] eqp pointer to a cs_equation_param_t structure
3562 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
3563 * \param[in, out] context pointer to cs_cdovcb_scaleq_t structure
3564 */
3565 /*----------------------------------------------------------------------------*/
3566
3567 void
cs_cdovcb_scaleq_current_to_previous(const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)3568 cs_cdovcb_scaleq_current_to_previous(const cs_equation_param_t *eqp,
3569 cs_equation_builder_t *eqb,
3570 void *context)
3571 {
3572 CS_UNUSED(eqp);
3573 CS_UNUSED(eqb);
3574
3575 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
3576 cs_field_t *fld = cs_field_by_id(eqc->var_field_id);
3577
3578 /* Cell values */
3579
3580 if (eqc->cell_values_pre != NULL)
3581 memcpy(eqc->cell_values_pre, eqc->cell_values,
3582 sizeof(cs_real_t)*cs_shared_quant->n_cells);
3583
3584 /* Vertex values */
3585
3586 cs_field_current_to_previous(fld);
3587 }
3588
3589 /*----------------------------------------------------------------------------*/
3590 /*!
3591 * \brief Predefined extra-operations related to this equation
3592 *
3593 * \param[in] eqp pointer to a cs_equation_param_t structure
3594 * \param[in, out] eqb pointer to a cs_equation_builder_t structure
3595 * \param[in, out] context pointer to cs_cdovcb_scaleq_t structure
3596 */
3597 /*----------------------------------------------------------------------------*/
3598
3599 void
cs_cdovcb_scaleq_extra_post(const cs_equation_param_t * eqp,cs_equation_builder_t * eqb,void * context)3600 cs_cdovcb_scaleq_extra_post(const cs_equation_param_t *eqp,
3601 cs_equation_builder_t *eqb,
3602 void *context)
3603 {
3604 cs_cdovcb_scaleq_t *eqc = (cs_cdovcb_scaleq_t *)context;
3605
3606 /* TODO */
3607 CS_UNUSED(eqp);
3608 CS_UNUSED(eqb);
3609 CS_UNUSED(eqc);
3610 }
3611
3612 /*----------------------------------------------------------------------------*/
3613
3614 #undef _dp3
3615
3616 END_C_DECLS
3617