1 /*============================================================================
2 * Convection diffusion equation solver.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------*/
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <stdarg.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39 #include <assert.h>
40 #include <math.h>
41
42 #if defined(HAVE_MPI)
43 #include <mpi.h>
44 #endif
45
46 /*----------------------------------------------------------------------------
47 * Local headers
48 *----------------------------------------------------------------------------*/
49
50 #include "bft_mem.h"
51 #include "bft_error.h"
52 #include "bft_printf.h"
53
54 #include "cs_balance.h"
55 #include "cs_blas.h"
56 #include "cs_convection_diffusion.h"
57 #include "cs_field.h"
58 #include "cs_field_pointer.h"
59 #include "cs_halo.h"
60 #include "cs_log.h"
61 #include "cs_math.h"
62 #include "cs_mesh.h"
63 #include "cs_field.h"
64 #include "cs_gradient.h"
65 #include "cs_mesh_quantities.h"
66 #include "cs_parameters.h"
67 #include "cs_porous_model.h"
68 #include "cs_prototypes.h"
69 #include "cs_timer.h"
70 #include "cs_parall.h"
71 #include "cs_matrix_building.h"
72 #include "cs_matrix_default.h"
73 #include "cs_sles.h"
74 #include "cs_sles_default.h"
75
76 /*----------------------------------------------------------------------------
77 * Header for the current file
78 *----------------------------------------------------------------------------*/
79
80 #include "cs_equation_iterative_solve.h"
81
82 /*----------------------------------------------------------------------------*/
83
84 BEGIN_C_DECLS
85
86 /*============================================================================
87 * Private function definitions
88 *============================================================================*/
89
90 /*============================================================================
91 * Public function definitions
92 *============================================================================*/
93
94 /*----------------------------------------------------------------------------*/
95 /*! \file cs_equation_iterative_solve.c
96 *
97 * \brief This file gathers functions that solve advection diffusion equations
98 * with source terms for one time step for a scalar, vector or tensor variable.
99 *
100 */
101 /*----------------------------------------------------------------------------*/
102
103 /*----------------------------------------------------------------------------*/
104 /*! \brief This function solves an advection diffusion equation with source
105 * terms for one time step for the variable \f$ a \f$.
106 *
107 * The equation reads:
108 *
109 * \f[
110 * f_s^{imp}(a^{n+1}-a^n)
111 * + \divs \left( a^{n+1} \rho \vect{u} - \mu \grad a^{n+1} \right)
112 * = Rhs
113 * \f]
114 *
115 * This equation is rewritten as:
116 *
117 * \f[
118 * f_s^{imp} \delta a
119 * + \divs \left( \delta a \rho \vect{u} - \mu \grad \delta a \right)
120 * = Rhs^1
121 * \f]
122 *
123 * where \f$ \delta a = a^{n+1} - a^n\f$ and
124 * \f$ Rhs^1 = Rhs - \divs( a^n \rho \vect{u} - \mu \grad a^n)\f$
125 *
126 *
127 * It is in fact solved with the following iterative process:
128 *
129 * \f[
130 * f_s^{imp} \delta a^k
131 * + \divs \left(\delta a^k \rho \vect{u}-\mu\grad\delta a^k \right)
132 * = Rhs^k
133 * \f]
134 *
135 * where \f$Rhs^k=Rhs-f_s^{imp}(a^k-a^n)
136 * - \divs \left( a^k\rho\vect{u}-\mu\grad a^k \right)\f$
137 *
138 * Be careful, it is forbidden to modify \f$ f_s^{imp} \f$ here!
139 *
140 * Please refer to the
141 * <a href="../../theory.pdf#codits"><b>codits</b></a> section of the
142 * theory guide for more informations.
143 *
144 * \param[in] idtvar indicator of the temporal scheme
145 * \param[in] iterns external sub-iteration number
146 * \param[in] f_id field id (or -1)
147 * \param[in] name associated name if f_id < 0, or NULL
148 * \param[in] iescap compute the predictor indicator if 1
149 * \param[in] imucpp indicator
150 * - 0 do not multiply the convectiv term by Cp
151 * - 1 do multiply the convectiv term by Cp
152 * \param[in] normp Reference norm to solve the system (optional)
153 * if negative: recomputed here
154 * \param[in] var_cal_opt pointer to a cs_var_cal_opt_t structure which
155 * contains variable calculation options
156 * \param[in] pvara variable at the previous time step
157 * \f$ a^n \f$
158 * \param[in] pvark variable at the previous sub-iteration
159 * \f$ a^k \f$.
160 * If you sub-iter on Navier-Stokes, then
161 * it allows to initialize by something else than
162 * pvara (usually pvar=pvara)
163 * \param[in] coefap boundary condition array for the variable
164 * (explicit part)
165 * \param[in] coefbp boundary condition array for the variable
166 * (implicit part)
167 * \param[in] cofafp boundary condition array for the diffusion
168 * of the variable (explicit part)
169 * \param[in] cofbfp boundary condition array for the diffusion
170 * of the variable (implicit part)
171 * \param[in] i_massflux mass flux at interior faces
172 * \param[in] b_massflux mass flux at boundary faces
173 * \param[in] i_viscm \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
174 * at interior faces for the matrix
175 * \param[in] b_viscm \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
176 * at boundary faces for the matrix
177 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
178 * at interior faces for the r.h.s.
179 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
180 * at boundary faces for the r.h.s.
181 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
182 * \param[in] weighf internal face weight between cells i j in case
183 * of tensor diffusion
184 * \param[in] weighb boundary face weight for cells i in case
185 * of tensor diffusion
186 * \param[in] icvflb global indicator of boundary convection flux
187 * - 0 upwind scheme at all boundary faces
188 * - 1 imposed flux at some boundary faces
189 * \param[in] icvfli boundary face indicator array of convection flux
190 * - 0 upwind scheme
191 * - 1 imposed flux
192 * \param[in] rovsdt \f$ f_s^{imp} \f$
193 * \param[in] smbrp Right hand side \f$ Rhs^k \f$
194 * \param[in,out] pvar current variable
195 * \param[out] dpvar last variable increment
196 * \param[in] xcpp array of specific heat (Cp)
197 * \param[out] eswork prediction-stage error estimator
198 * (if iescap > 0)
199 */
200 /*----------------------------------------------------------------------------*/
201
202 void
cs_equation_iterative_solve_scalar(int idtvar,int iterns,int f_id,const char * name,int iescap,int imucpp,cs_real_t normp,cs_var_cal_opt_t * var_cal_opt,const cs_real_t pvara[],const cs_real_t pvark[],const cs_real_t coefap[],const cs_real_t coefbp[],const cs_real_t cofafp[],const cs_real_t cofbfp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_viscm[],const cs_real_t b_viscm[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_6_t viscel[],const cs_real_2_t weighf[],const cs_real_t weighb[],int icvflb,const int icvfli[],const cs_real_t rovsdt[],cs_real_t smbrp[],cs_real_t pvar[],cs_real_t dpvar[],const cs_real_t xcpp[],cs_real_t eswork[])203 cs_equation_iterative_solve_scalar(int idtvar,
204 int iterns,
205 int f_id,
206 const char *name,
207 int iescap,
208 int imucpp,
209 cs_real_t normp,
210 cs_var_cal_opt_t *var_cal_opt,
211 const cs_real_t pvara[],
212 const cs_real_t pvark[],
213 const cs_real_t coefap[],
214 const cs_real_t coefbp[],
215 const cs_real_t cofafp[],
216 const cs_real_t cofbfp[],
217 const cs_real_t i_massflux[],
218 const cs_real_t b_massflux[],
219 const cs_real_t i_viscm[],
220 const cs_real_t b_viscm[],
221 const cs_real_t i_visc[],
222 const cs_real_t b_visc[],
223 cs_real_6_t viscel[],
224 const cs_real_2_t weighf[],
225 const cs_real_t weighb[],
226 int icvflb,
227 const int icvfli[],
228 const cs_real_t rovsdt[],
229 cs_real_t smbrp[],
230 cs_real_t pvar[],
231 cs_real_t dpvar[],
232 const cs_real_t xcpp[],
233 cs_real_t eswork[])
234 {
235 /* Local variables */
236
237 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
238
239 int iconvp = var_cal_opt->iconv;
240 int idiffp = var_cal_opt->idiff;
241 int iwarnp = var_cal_opt->verbosity;
242 int iswdyp = var_cal_opt->iswdyn;
243 int ndircp = var_cal_opt->ndircl;
244 double epsrsp = var_cal_opt->epsrsm;
245 double epsilp = var_cal_opt->epsilo;
246 double relaxp = var_cal_opt->relaxv;
247 double thetap = var_cal_opt->thetav;
248
249 const cs_real_t *cell_vol = cs_glob_mesh_quantities->cell_vol;
250 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
251 const cs_lnum_t n_i_faces = cs_glob_mesh->n_i_faces;
252 const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
253
254 int isym, inc, isweep, niterf, iccocg, nswmod;
255 int lvar, ibsize, iesize, imasac, key_sinfo_id;
256 double residu, rnorm, ressol, rnorm2;
257 double thetex, nadxkm1, nadxk, paxm1ax, paxm1rk, paxkrk, alph, beta;
258
259 cs_lnum_t eb_size[4], db_size[4];
260
261 cs_solving_info_t sinfo;
262
263 cs_field_t *f = NULL;
264 int coupling_id = -1;
265
266 cs_real_t *dam, *xam, *smbini;
267 cs_real_t *dam_conv = NULL, *xam_conv = NULL;
268 cs_real_t *dam_diff = NULL, *xam_diff = NULL;
269
270 cs_real_t *w1 = NULL;
271
272 bool conv_diff_mg = false;
273
274 /*============================================================================
275 * 0. Initialization
276 *==========================================================================*/
277
278 /* Name */
279 const char *var_name = cs_sles_name(f_id, name);
280
281 /* solving info */
282 key_sinfo_id = cs_field_key_id("solving_info");
283 if (f_id > -1) {
284 f = cs_field_by_id(f_id);
285 cs_field_get_key_struct(f, key_sinfo_id, &sinfo);
286 coupling_id = cs_field_get_key_int(f, cs_field_key_id("coupling_entity"));
287 }
288
289 /* Determine if we are in a case with special requirements */
290
291 if (coupling_id < 0 && iconvp > 0) {
292 cs_sles_t *sc = cs_sles_find_or_add(f_id, name);
293 const char *sles_type = cs_sles_get_type(sc);
294 if (strcmp(sles_type, "cs_multigrid_t") == 0)
295 conv_diff_mg = true;
296 }
297
298 /* Allocate temporary arrays */
299
300 BFT_MALLOC(dam, n_cells_ext, cs_real_t);
301 if (conv_diff_mg) {
302 BFT_MALLOC(dam_conv, n_cells_ext, cs_real_t);
303 BFT_MALLOC(dam_diff, n_cells_ext, cs_real_t);
304 }
305 BFT_MALLOC(smbini, n_cells_ext, cs_real_t);
306
307 cs_real_t *adxk = NULL, *adxkm1 = NULL, *dpvarm1 = NULL, *rhs0 = NULL;
308
309 if (iswdyp >= 1) {
310 BFT_MALLOC(adxk, n_cells_ext, cs_real_t);
311 BFT_MALLOC(adxkm1, n_cells_ext, cs_real_t);
312 BFT_MALLOC(dpvarm1, n_cells_ext, cs_real_t);
313 BFT_MALLOC(rhs0, n_cells_ext, cs_real_t);
314 }
315
316 /* Symmetric matrix, except if advection */
317 isym = 1;
318 if (iconvp > 0) isym = 2;
319
320 bool symmetric = (isym == 1) ? true : false;
321
322 BFT_MALLOC(xam, isym*n_i_faces, cs_real_t);
323 if (conv_diff_mg) {
324 BFT_MALLOC(xam_conv, 2*n_i_faces, cs_real_t);
325 BFT_MALLOC(xam_diff, n_i_faces, cs_real_t);
326 }
327
328 /* Matrix block size */
329 ibsize = 1;
330 iesize = 1;
331
332 db_size[0] = ibsize;
333 db_size[1] = ibsize;
334 db_size[2] = ibsize;
335 db_size[3] = ibsize*ibsize;
336
337 eb_size[0] = iesize;
338 eb_size[1] = iesize;
339 eb_size[2] = iesize;
340 eb_size[3] = iesize*iesize;
341
342 /* Periodicity has to be taken into account */
343
344 /* Initialization for test before matrix vector product
345 for computation of initial residual */
346
347 /*============================================================================
348 * 1. Building of the "simplified" matrix
349 *==========================================================================*/
350
351 if (conv_diff_mg) {
352 cs_matrix_wrapper_scalar_conv_diff(iconvp,
353 idiffp,
354 ndircp,
355 thetap,
356 imucpp,
357 coefbp,
358 cofbfp,
359 rovsdt,
360 i_massflux,
361 b_massflux,
362 i_viscm,
363 b_viscm,
364 xcpp,
365 dam,
366 xam,
367 dam_conv,
368 xam_conv,
369 dam_diff,
370 xam_diff);
371 }
372 else {
373 cs_matrix_wrapper_scalar(iconvp,
374 idiffp,
375 ndircp,
376 isym,
377 thetap,
378 imucpp,
379 coefbp,
380 cofbfp,
381 rovsdt,
382 i_massflux,
383 b_massflux,
384 i_viscm,
385 b_viscm,
386 xcpp,
387 dam,
388 xam);
389 }
390
391 /* For steady computations, the diagonal is relaxed */
392 if (idtvar < 0) {
393 # pragma omp parallel for
394 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
395 dam[iel] /= relaxp;
396 }
397
398 /*==========================================================================
399 * 2. Iterative process to handle non orthogonalities (starting from the
400 * second iteration).
401 *==========================================================================*/
402
403 /* Application of the theta-scheme */
404
405 /* On calcule le bilan explicite total */
406 thetex = 1. - thetap;
407
408 /* Compute the min/ max limiter */
409 inc = 1;
410
411 if (f_id > -1)
412 cs_beta_limiter_building(f_id, inc, rovsdt);
413
414 /* If thetex = 0, no need to do more */
415 if (fabs(thetex) > cs_math_epzero) {
416 inc = 1;
417 iccocg = 1;
418
419 /* The added convective scalar mass flux is:
420 (thetex*Y_\face-imasac*Y_\celli)*mf.
421 When building the explicit part of the rhs, one
422 has to impose 0 on mass accumulation. */
423 imasac = 0;
424
425 var_cal_opt->thetav = thetex;
426
427 cs_balance_scalar(idtvar,
428 f_id,
429 imucpp,
430 imasac,
431 inc,
432 iccocg,
433 var_cal_opt,
434 NULL, /* pvar == pvara */
435 pvara,
436 coefap,
437 coefbp,
438 cofafp,
439 cofbfp,
440 i_massflux,
441 b_massflux,
442 i_visc,
443 b_visc,
444 viscel,
445 xcpp,
446 weighf,
447 weighb,
448 icvflb,
449 icvfli,
450 smbrp);
451
452 var_cal_opt->thetav = thetap;
453
454 }
455
456 /* Before looping, the RHS without reconstruction is stored in smbini */
457
458 cs_lnum_t has_dc = mq->has_disable_flag;
459 # pragma omp parallel if(n_cells > CS_THR_MIN)
460 {
461 # pragma omp for nowait
462 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
463 smbini[cell_id] = smbrp[cell_id];
464 smbrp[cell_id] = 0.;
465 }
466
467 # pragma omp for nowait
468 for (cs_lnum_t iel = 0; iel < n_cells_ext; iel++)
469 pvar[iel] = pvark[iel];
470 }
471
472 /* In the following, cs_balance_scalar is called with inc=1,
473 except for Weight Matrix (nswrsp=-1) */
474 inc = 1;
475
476 if (var_cal_opt->nswrsm == -1) {
477 var_cal_opt->nswrsm = 1;
478 inc = 0;
479 }
480
481 /* Incrementation and rebuild of right hand side */
482
483 /* We entered with an explicit RHS based on PVARA.
484 If we initialize with PVAR with something else than PVARA
485 we must correc the RHS (this is the case when iterating on navsto) */
486
487 iccocg = 1;
488
489 /* The added convective scalar mass flux is:
490 (thetap*Y_\face-imasac*Y_\celli)*mf.
491 When building the implicit part of the rhs, one
492 has to impose 1 on mass accumulation. */
493 imasac = 1;
494
495 cs_balance_scalar(idtvar,
496 f_id,
497 imucpp,
498 imasac,
499 inc,
500 iccocg,
501 var_cal_opt,
502 pvar,
503 pvara,
504 coefap,
505 coefbp,
506 cofafp,
507 cofbfp,
508 i_massflux,
509 b_massflux,
510 i_visc,
511 b_visc,
512 viscel,
513 xcpp,
514 weighf,
515 weighb,
516 icvflb,
517 icvfli,
518 smbrp);
519
520 if (iswdyp >= 1) {
521 # pragma omp parallel for
522 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
523 rhs0[iel] = smbrp[iel];
524 smbini[iel] -= rovsdt[iel]*(pvar[iel] - pvara[iel]);
525 smbrp[iel] += smbini[iel];
526
527 adxkm1[iel] = 0.;
528 adxk[iel] = 0.;
529 dpvar[iel] = 0.;
530 }
531
532 /* ||A.dx^0||^2 = 0 */
533 nadxk = 0.;
534 }
535 else {
536 # pragma omp parallel for
537 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
538 smbini[iel] -= rovsdt[iel]*(pvar[iel] - pvara[iel]);
539 smbrp[iel] += smbini[iel];
540 }
541 }
542
543 /* --- Right hand side residual */
544 residu = sqrt(cs_gdot(n_cells, smbrp, smbrp));
545
546 if (normp > 0.)
547 rnorm = normp;
548
549 else {
550 /* Normalization residual
551 (L2-norm of B.C. + source terms + non-orthogonality terms)
552
553 Caution: when calling a matrix-vector product, here for a variable
554 which is not "by increments" and is assumed initialized, including
555 for ghost values. */
556
557 /* Allocate a temporary array */
558 BFT_MALLOC(w1, n_cells_ext, cs_real_t);
559
560 cs_real_t *w2;
561 BFT_MALLOC(w2, n_cells_ext, cs_real_t);
562
563 cs_real_t p_mean = cs_gmean(n_cells, mq->cell_vol, pvar);
564
565 if (iwarnp >= 2)
566 bft_printf("Spatial average of X^n = %f\n", p_mean);
567 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
568 w2[cell_id] = (pvar[cell_id]-p_mean);
569
570 cs_matrix_vector_native_multiply(symmetric,
571 db_size,
572 eb_size,
573 f_id,
574 dam,
575 xam,
576 w2,
577 w1);
578
579 BFT_FREE(w2);
580
581 if (iwarnp >= 2) {
582 bft_printf("L2 norm ||AX^n|| = %f\n", sqrt(cs_gdot(n_cells, w1, w1)));
583 bft_printf("L2 norm ||B^n|| = %f\n", sqrt(cs_gdot(n_cells, smbrp, smbrp)));
584 }
585 # pragma omp parallel for
586 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
587 w1[cell_id] += smbrp[cell_id];
588 /* Remove contributions from penalized cells */
589 if (has_dc * mq->c_disable_flag[has_dc * cell_id] != 0)
590 w1[cell_id] = 0.;
591 }
592
593 rnorm2 = cs_gdot(n_cells, w1, w1);
594 rnorm = sqrt(rnorm2);
595 }
596
597 sinfo.rhs_norm = rnorm;
598
599 /* Free memory */
600 BFT_FREE(w1);
601
602 /* Warning: for weight matrix, one and only one sweep is done. */
603 nswmod = CS_MAX(var_cal_opt->nswrsm, 1);
604
605 /* Reconstruction loop (beginning) */
606 if (iterns <= 1)
607 sinfo.n_it = 0;
608 isweep = 1;
609
610 while ((isweep <= nswmod && residu > epsrsp*rnorm) || isweep == 1) {
611
612 /* Solving on the increment: dpvar */
613
614 if (iswdyp >= 1) {
615 # pragma omp parallel for
616 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
617 dpvarm1[iel] = dpvar[iel];
618 dpvar[iel] = 0.;
619 }
620 }
621 else {
622 # pragma omp parallel for
623 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
624 dpvar[iel] = 0.;
625 }
626
627 /* Solver residual */
628 ressol = residu;
629
630 if (conv_diff_mg)
631 cs_sles_setup_native_conv_diff(f_id,
632 var_name,
633 db_size,
634 eb_size,
635 dam,
636 xam,
637 dam_conv,
638 xam_conv,
639 dam_diff,
640 xam_diff);
641
642 cs_sles_solve_native(f_id,
643 var_name,
644 symmetric,
645 db_size,
646 eb_size,
647 dam,
648 xam,
649 epsilp,
650 rnorm,
651 &niterf,
652 &ressol,
653 smbrp,
654 dpvar);
655
656 /* Dynamic relaxation of the system */
657 if (iswdyp >= 1) {
658
659 /* Computation of the variable relaxation coefficient */
660 lvar = -1;
661
662 # pragma omp parallel for
663 for (cs_lnum_t iel = 0; iel < n_cells_ext; iel++) {
664 adxkm1[iel] = adxk[iel];
665 adxk[iel] = - rhs0[iel];
666 }
667
668 cs_balance_scalar(idtvar,
669 lvar,
670 imucpp,
671 imasac,
672 inc,
673 iccocg,
674 var_cal_opt,
675 dpvar,
676 NULL, /* dpvara == dpvar */
677 coefap,
678 coefbp,
679 cofafp,
680 cofbfp,
681 i_massflux,
682 b_massflux,
683 i_visc,
684 b_visc,
685 viscel,
686 xcpp,
687 weighf,
688 weighb,
689 icvflb,
690 icvfli,
691 adxk);
692
693 /* ||E.dx^(k-1)-E.0||^2 */
694 nadxkm1 = nadxk;
695
696 /* ||E.dx^k-E.0||^2 */
697 nadxk = cs_gdot(n_cells, adxk, adxk);
698
699 /* < E.dx^k-E.0; r^k > */
700 paxkrk = cs_gdot(n_cells, smbrp, adxk);
701
702 /* Relaxation with respect to dx^k and dx^(k-1) */
703 if (iswdyp >= 2) {
704
705 /* < E.dx^(k-1)-E.0; r^k > */
706 paxm1rk = cs_gdot(n_cells, smbrp, adxkm1);
707
708 /* < E.dx^(k-1)-E.0; E.dx^k-E.0 > */
709 paxm1ax = cs_gdot(n_cells, adxk, adxkm1);
710
711 if (nadxkm1 > 1e-30*rnorm2 && (nadxk*nadxkm1-pow(paxm1ax,2)) > 1e-30*rnorm2)
712 beta = (paxkrk*paxm1ax - nadxk*paxm1rk)/(nadxk*nadxkm1-pow(paxm1ax,2));
713 else
714 beta = 0.;
715
716 } else {
717 beta = 0.;
718 paxm1ax = 1.;
719 paxm1rk = 0.;
720 paxm1ax = 0.;
721 }
722
723 /* The first sweep is not relaxed */
724 if (isweep == 1) {
725 alph = 1.;
726 beta = 0.;
727 } else if (isweep == 2) {
728 beta = 0.;
729 alph = -paxkrk/CS_MAX(nadxk, 1e-30*rnorm2);
730 } else {
731 alph = -(paxkrk + beta*paxm1ax)/CS_MAX(nadxk, 1e-30*rnorm2);
732 }
733
734 /* Writing */
735 if (iwarnp >= 3)
736 bft_printf("%s Sweep: %d Dynamic relaxation: alpha = %12.5e, "
737 "beta = %12.5e,\n< dI^k : R^k > = %12.5e, "
738 "||dI^k ||^2 = %12.5e,\n< dI^k-1 : R^k > = %12.5e, "
739 "||dI^k-1||^2 = %12.5e,\n< dI^k-1 : dI^k > = %12.5e\n",
740 var_name, isweep, alph, beta, paxkrk, nadxk, paxm1rk,
741 nadxkm1, paxm1ax);
742 }
743
744 /* --- Update the solution with the increment */
745
746 if (iswdyp <= 0) {
747 # pragma omp parallel for
748 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
749 pvar[iel] += dpvar[iel];
750 } else if (iswdyp == 1) {
751 if (alph < 0.) break;
752 # pragma omp parallel for
753 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
754 pvar[iel] += alph*dpvar[iel];
755 } else if (iswdyp >= 2) {
756 # pragma omp parallel for
757 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
758 pvar[iel] += alph*dpvar[iel] + beta*dpvarm1[iel];
759 }
760
761 /* ---> Handle parallelism and periodicity */
762 if (cs_glob_rank_id >= 0 || cs_glob_mesh->n_init_perio > 0) {
763 cs_mesh_sync_var_scal(pvar);
764 }
765
766 /* --- Update the right hand side And compute the new residual */
767
768 iccocg = 0;
769
770 if (iswdyp <= 0) {
771 # pragma omp parallel for
772 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
773 /* smbini already contains unsteady terms and mass source terms
774 of the RHS updated at each sweep */
775 smbini[iel] -= rovsdt[iel]*dpvar[iel];
776 smbrp[iel] = smbini[iel];
777 }
778 } else if (iswdyp == 1) {
779 # pragma omp parallel for
780 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
781 /* smbini already contains unsteady terms and mass source terms
782 of the RHS updated at each sweep */
783 smbini[iel] -= rovsdt[iel]*alph*dpvar[iel];
784 smbrp[iel] = smbini[iel];
785 }
786 } else if (iswdyp >= 2) {
787 # pragma omp parallel for
788 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
789 /* smbini already contains unsteady terms and mass source terms
790 of the RHS updated at each sweep */
791 smbini[iel] -= rovsdt[iel]*(alph*dpvar[iel]+beta*dpvarm1[iel]);
792 smbrp[iel] = smbini[iel];
793 }
794 }
795
796 /* Compute the beta (min/max) limiter */
797 if (f_id > -1)
798 cs_beta_limiter_building(f_id, inc, rovsdt);
799
800 /* The added convective scalar mass flux is:
801 (thetex*Y_\face-imasac*Y_\celli)*mf.
802 When building the implicit part of the rhs, one
803 has to impose 1 on mass accumulation. */
804 imasac = 1;
805
806 cs_balance_scalar(idtvar,
807 f_id,
808 imucpp,
809 imasac,
810 inc,
811 iccocg,
812 var_cal_opt,
813 pvar,
814 pvara,
815 coefap,
816 coefbp,
817 cofafp,
818 cofbfp,
819 i_massflux,
820 b_massflux,
821 i_visc,
822 b_visc,
823 viscel,
824 xcpp,
825 weighf,
826 weighb,
827 icvflb,
828 icvfli,
829 smbrp);
830
831 /* --- Convergence test */
832 residu = sqrt(cs_gdot(n_cells, smbrp, smbrp));
833
834 /* Writing */
835 sinfo.n_it = sinfo.n_it + niterf;
836
837 /* Writing */
838 if (iwarnp >= 2) {
839 bft_printf("%s: CV_DIF_TS, IT: %d, Res: %12.5e, Norm: %12.5e\n",
840 var_name, isweep, residu, rnorm);
841 bft_printf("%s: Current reconstruction sweep: %d, "
842 "Iterations for solver: %d\n", var_name, isweep, niterf);
843 }
844
845 isweep++;
846
847 }
848 /* --- Reconstruction loop (end) */
849
850 /* Writing: convergence */
851 if (fabs(rnorm) > cs_math_epzero)
852 sinfo.res_norm = residu/rnorm;
853 else
854 sinfo.res_norm = 0.;
855
856 /* Save convergence info */
857 if (f_id > -1) {
858 f = cs_field_by_id(f_id);
859 cs_field_set_key_struct(f, key_sinfo_id, &sinfo);
860 }
861
862 if (iwarnp >= 1) {
863 if (residu <= epsrsp*rnorm)
864 bft_printf("%s: CV_DIF_TS, IT : %d, Res : %12.5e, Norm : %12.5e\n",
865 var_name, isweep-1, residu, rnorm);
866 /* Writing: non-convergence */
867 else if (isweep > nswmod)
868 bft_printf("@\n@ @@ WARNING: %s CONVECTION-DIFFUSION-SOURCE-TERMS\n@"
869 "=========\n@ Maximum number of iterations %d reached\n@",
870 var_name,nswmod);
871 }
872
873
874 /* Rebuild final flux at faces */
875
876 if (f_id > -1) {
877 f = cs_field_by_id(f_id);
878
879 const int kiflux = cs_field_key_id("inner_flux_id");
880 const int kbflux = cs_field_key_id("boundary_flux_id");
881
882 int i_flux_id = cs_field_get_key_int(f, kiflux);
883 int b_flux_id = cs_field_get_key_int(f, kbflux);
884
885 if (i_flux_id > -1 && b_flux_id > -1) {
886 assert(idtvar != -1);
887 /* flux is non-conservative with the steady algorithm but flux field is of
888 dimension 1. Forbidden at parameters check */
889
890 cs_field_t *i_flux = cs_field_by_id(i_flux_id);
891 cs_field_t *b_flux = cs_field_by_id(b_flux_id);
892
893 /* rebuild before-last value of variable */
894 cs_real_t *prev_s_pvar;
895 BFT_MALLOC(prev_s_pvar, n_cells_ext, cs_real_t);
896 # pragma omp parallel for
897 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
898 prev_s_pvar[iel] = pvar[iel]-dpvar[iel];
899 }
900
901 cs_real_t *i_flux2;
902 BFT_MALLOC(i_flux2, 2*n_i_faces, cs_real_t);
903 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
904 i_flux2[2*face_id ] = 0.;
905 i_flux2[2*face_id+1] = 0.;
906 }
907
908 inc = 1;
909 iccocg = 1;
910 imasac = 0; /* mass accumluation not taken into account */
911
912 cs_face_convection_scalar(idtvar,
913 f_id,
914 *var_cal_opt,
915 icvflb,
916 inc,
917 iccocg,
918 imasac,
919 prev_s_pvar,
920 pvara,
921 icvfli,
922 coefap,
923 coefbp,
924 i_massflux,
925 b_massflux,
926 (cs_real_2_t *)i_flux2,
927 b_flux->val);
928 BFT_FREE(prev_s_pvar);
929
930 /* last increment in upwind to fulfill exactly the considered
931 balance equation */
932
933 inc = 0;
934 iccocg = 0;
935
936 cs_var_cal_opt_t var_cal_opt_loc;
937 int k_id = cs_field_key_id("var_cal_opt");
938 cs_field_get_key_struct(f, k_id, &var_cal_opt_loc);
939
940 var_cal_opt_loc.nswrgr = 0;
941 var_cal_opt_loc.blencv = 0.;
942
943 cs_face_convection_scalar(idtvar,
944 f_id,
945 var_cal_opt_loc,
946 icvflb,
947 inc,
948 iccocg,
949 imasac,
950 dpvar,
951 pvara,
952 icvfli,
953 coefap,
954 coefbp,
955 i_massflux,
956 b_massflux,
957 (cs_real_2_t *)i_flux2,
958 b_flux->val);
959
960 /* FIXME diffusion part */
961
962 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++)
963 i_flux->val[face_id] += i_flux2[2*face_id];
964 BFT_FREE(i_flux2);
965 }
966 }
967
968 /*==========================================================================
969 * 3. After having computed the new value, an estimator is computed for the
970 * prediction step of the velocity.
971 *==========================================================================*/
972
973 if (iescap > 0) {
974 /* Computation of the estimator of the current component */
975
976 /* smbini already contains unsteady terms and mass source terms
977 of the RHS updated at each sweep */
978
979 # pragma omp parallel for
980 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
981 smbrp[iel] = smbini[iel] - rovsdt[iel]*dpvar[iel];
982
983 inc = 1;
984 iccocg = 1;
985
986 /* Without relaxation even for a stationnary computation */
987
988 cs_balance_scalar(idtvar,
989 f_id,
990 imucpp,
991 imasac,
992 inc,
993 iccocg,
994 var_cal_opt,
995 pvar,
996 pvara,
997 coefap,
998 coefbp,
999 cofafp,
1000 cofbfp,
1001 i_massflux,
1002 b_massflux,
1003 i_visc,
1004 b_visc,
1005 viscel,
1006 xcpp,
1007 weighf,
1008 weighb,
1009 icvflb,
1010 icvfli,
1011 smbrp);
1012
1013 /* Contribution of the current component to the L2 norm stored in eswork */
1014
1015 # pragma omp parallel for
1016 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
1017 eswork[iel] = pow(smbrp[iel] / cell_vol[iel],2);
1018 }
1019
1020 /*==========================================================================
1021 * 4. Free solver setup
1022 *==========================================================================*/
1023
1024 cs_sles_free_native(f_id, var_name);
1025
1026 /* Free memory */
1027 BFT_FREE(dam);
1028 BFT_FREE(xam);
1029 if (conv_diff_mg) {
1030 BFT_FREE(dam_conv);
1031 BFT_FREE(xam_conv);
1032 BFT_FREE(dam_diff);
1033 BFT_FREE(xam_diff);
1034 }
1035
1036 BFT_FREE(smbini);
1037 if (iswdyp >= 1) {
1038 BFT_FREE(adxk);
1039 BFT_FREE(adxkm1);
1040 BFT_FREE(dpvarm1);
1041 BFT_FREE(rhs0);
1042 }
1043 }
1044
1045 /*----------------------------------------------------------------------------*/
1046 /*!
1047 * \brief This function solves an advection diffusion equation with source terms
1048 * for one time step for the vector variable \f$ \vect{a} \f$.
1049 *
1050 * The equation reads:
1051 *
1052 * \f[
1053 * \tens{f_s}^{imp}(\vect{a}^{n+1}-\vect{a}^n)
1054 * + \divv \left( \vect{a}^{n+1} \otimes \rho \vect {u}
1055 * - \mu \gradt \vect{a}^{n+1}\right)
1056 * = \vect{Rhs}
1057 * \f]
1058 *
1059 * This equation is rewritten as:
1060 *
1061 * \f[
1062 * \tens{f_s}^{imp} \delta \vect{a}
1063 * + \divv \left( \delta \vect{a} \otimes \rho \vect{u}
1064 * - \mu \gradt \delta \vect{a} \right)
1065 * = \vect{Rhs}^1
1066 * \f]
1067 *
1068 * where \f$ \delta \vect{a} = \vect{a}^{n+1} - \vect{a}^n\f$ and
1069 * \f$ \vect{Rhs}^1 = \vect{Rhs}
1070 * - \divv \left( \vect{a}^n \otimes \rho \vect{u}
1071 * - \mu \gradt \vect{a}^n \right)\f$
1072 *
1073 *
1074 * It is in fact solved with the following iterative process:
1075 *
1076 * \f[
1077 * \tens{f_s}^{imp} \delta \vect{a}^k
1078 * + \divv \left( \delta \vect{a}^k \otimes \rho \vect{u}
1079 * - \mu \gradt \delta \vect{a}^k \right)
1080 * = \vect{Rhs}^k
1081 * \f]
1082 *
1083 * where \f$ \vect{Rhs}^k = \vect{Rhs}
1084 * - \tens{f_s}^{imp} \left(\vect{a}^k-\vect{a}^n \right)
1085 * - \divv \left( \vect{a}^k \otimes \rho \vect{u}
1086 * - \mu \gradt \vect{a}^k \right)\f$
1087 *
1088 * Be careful, it is forbidden to modify \f$ \tens{f_s}^{imp} \f$ here!
1089 *
1090 * \param[in] idtvar indicator of the temporal scheme
1091 * \param[in] iterns external sub-iteration number
1092 * \param[in] f_id field id (or -1)
1093 * \param[in] name associated name if f_id < 0, or NULL
1094 * \param[in] ivisep indicator to take \f$ \divv
1095 * \left(\mu \gradt \transpose{\vect{a}} \right)
1096 * -2/3 \grad\left( \mu \dive \vect{a} \right)\f$
1097 * - 1 take into account,
1098 * - 0 otherwise
1099 * \param[in] iescap compute the predictor indicator if 1
1100 * \param[in] var_cal_opt pointer to a cs_var_cal_opt_t structure which
1101 * contains variable calculation options
1102 * \param[in] pvara variable at the previous time step
1103 * \f$ \vect{a}^n \f$
1104 * \param[in] pvark variable at the previous sub-iteration
1105 * \f$ \vect{a}^k \f$.
1106 * If you sub-iter on Navier-Stokes, then
1107 * it allows to initialize by something else than
1108 * \c pvara (usually \c pvar= \c pvara)
1109 * \param[in] coefav boundary condition array for the variable
1110 * (explicit part)
1111 * \param[in] coefbv boundary condition array for the variable
1112 * (implicit part)
1113 * \param[in] cofafv boundary condition array for the diffusion
1114 * of the variable (Explicit part)
1115 * \param[in] cofbfv boundary condition array for the diffusion
1116 * of the variable (Implicit part)
1117 * \param[in] i_massflux mass flux at interior faces
1118 * \param[in] b_massflux mass flux at boundary faces
1119 * \param[in] i_viscm \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1120 * at interior faces for the matrix
1121 * \param[in] b_viscm \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1122 * at boundary faces for the matrix
1123 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1124 * at interior faces for the r.h.s.
1125 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1126 * at boundary faces for the r.h.s.
1127 * \param[in] i_secvis secondary viscosity at interior faces
1128 * \param[in] b_secvis secondary viscosity at boundary faces
1129 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
1130 * \param[in] weighf internal face weight between cells i j in case
1131 * of tensor diffusion
1132 * \param[in] weighb boundary face weight for cells i in case
1133 * of tensor diffusion
1134 * \param[in] icvflb global indicator of boundary convection flux
1135 * - 0 upwind scheme at all boundary faces
1136 * - 1 imposed flux at some boundary faces
1137 * \param[in] icvfli boundary face indicator array of convection flux
1138 * - 0 upwind scheme
1139 * - 1 imposed flux
1140 * \param[in] fimp \f$ \tens{f_s}^{imp} \f$
1141 * \param[in] smbrp Right hand side \f$ \vect{Rhs}^k \f$
1142 * \param[in,out] pvar current variable
1143 * \param[out] eswork prediction-stage error estimator
1144 * (if iescap > 0)
1145 */
1146 /*----------------------------------------------------------------------------*/
1147
1148 void
cs_equation_iterative_solve_vector(int idtvar,int iterns,int f_id,const char * name,int ivisep,int iescap,cs_var_cal_opt_t * var_cal_opt,const cs_real_3_t pvara[],const cs_real_3_t pvark[],const cs_real_3_t coefav[],const cs_real_33_t coefbv[],const cs_real_3_t cofafv[],const cs_real_33_t cofbfv[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],cs_real_t i_viscm[],const cs_real_t b_viscm[],const cs_real_t i_visc[],const cs_real_t b_visc[],const cs_real_t i_secvis[],const cs_real_t b_secvis[],cs_real_6_t viscel[],const cs_real_2_t weighf[],const cs_real_t weighb[],int icvflb,const int icvfli[],cs_real_33_t fimp[],cs_real_3_t smbrp[],cs_real_3_t pvar[],cs_real_3_t eswork[])1149 cs_equation_iterative_solve_vector(int idtvar,
1150 int iterns,
1151 int f_id,
1152 const char *name,
1153 int ivisep,
1154 int iescap,
1155 cs_var_cal_opt_t *var_cal_opt,
1156 const cs_real_3_t pvara[],
1157 const cs_real_3_t pvark[],
1158 const cs_real_3_t coefav[],
1159 const cs_real_33_t coefbv[],
1160 const cs_real_3_t cofafv[],
1161 const cs_real_33_t cofbfv[],
1162 const cs_real_t i_massflux[],
1163 const cs_real_t b_massflux[],
1164 cs_real_t i_viscm[],
1165 const cs_real_t b_viscm[],
1166 const cs_real_t i_visc[],
1167 const cs_real_t b_visc[],
1168 const cs_real_t i_secvis[],
1169 const cs_real_t b_secvis[],
1170 cs_real_6_t viscel[],
1171 const cs_real_2_t weighf[],
1172 const cs_real_t weighb[],
1173 int icvflb,
1174 const int icvfli[],
1175 cs_real_33_t fimp[],
1176 cs_real_3_t smbrp[],
1177 cs_real_3_t pvar[],
1178 cs_real_3_t eswork[])
1179 {
1180 /* Local variables */
1181
1182 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1183
1184 int iconvp = var_cal_opt->iconv;
1185 int idiffp = var_cal_opt->idiff;
1186 int iwarnp = var_cal_opt->verbosity;
1187 int iswdyp = var_cal_opt->iswdyn;
1188 int idftnp = var_cal_opt->idften;
1189 int ndircp = var_cal_opt->ndircl;
1190 double epsrsp = var_cal_opt->epsrsm;
1191 double epsilp = var_cal_opt->epsilo;
1192 double relaxp = var_cal_opt->relaxv;
1193 double thetap = var_cal_opt->thetav;
1194
1195 const cs_real_t *cell_vol = cs_glob_mesh_quantities->cell_vol;
1196 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
1197 const cs_lnum_t n_faces = cs_glob_mesh->n_i_faces;
1198 const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
1199
1200 int isym, inc, isweep, niterf, nswmod, ibsize;
1201 int key_sinfo_id;
1202 int iesize, lvar, imasac;
1203 double residu, rnorm, ressol, rnorm2, thetex, alph, beta;
1204 double paxkrk, nadxk, paxm1rk, nadxkm1, paxm1ax;
1205
1206 cs_lnum_t eb_size[4], db_size[4];
1207
1208 cs_solving_info_t sinfo;
1209
1210 cs_field_t *f;
1211
1212 cs_real_t *xam;
1213 cs_real_33_t *dam;
1214 cs_real_3_t *dpvar, *smbini, *w1, *w2, *adxk, *adxkm1, *dpvarm1, *rhs0;
1215
1216 /*============================================================================
1217 * 0. Initialization
1218 *==========================================================================*/
1219
1220 /* Matrix block size */
1221 ibsize = 3;
1222 iesize = 1; /* CS_ISOTROPIC_DIFFUSION or CS_ANISOTROPIC_RIGHT_DIFFUSION */
1223 if (idftnp & CS_ANISOTROPIC_LEFT_DIFFUSION)
1224 iesize = 3;
1225
1226 if (cs_glob_porous_model == 3) { //FIXME iphydr + other option?
1227 if (iconvp > 0)
1228 iesize = 3;
1229 }
1230
1231 db_size[0] = ibsize;
1232 db_size[1] = ibsize;
1233 db_size[2] = ibsize;
1234 db_size[3] = ibsize*ibsize;
1235
1236 eb_size[0] = iesize;
1237 eb_size[1] = iesize;
1238 eb_size[2] = iesize;
1239 eb_size[3] = iesize*iesize;
1240
1241 /* Allocate temporary arrays */
1242 BFT_MALLOC(dam, n_cells_ext, cs_real_33_t);
1243 BFT_MALLOC(dpvar, n_cells_ext, cs_real_3_t);
1244 BFT_MALLOC(smbini, n_cells_ext, cs_real_3_t);
1245
1246 if (iswdyp >= 1) {
1247 BFT_MALLOC(adxk, n_cells_ext, cs_real_3_t);
1248 BFT_MALLOC(adxkm1, n_cells_ext, cs_real_3_t);
1249 BFT_MALLOC(dpvarm1, n_cells_ext, cs_real_3_t);
1250 BFT_MALLOC(rhs0, n_cells_ext, cs_real_3_t);
1251 }
1252
1253 /* solving info */
1254 key_sinfo_id = cs_field_key_id("solving_info");
1255 if (f_id > -1) {
1256 f = cs_field_by_id(f_id);
1257 cs_field_get_key_struct(f, key_sinfo_id, &sinfo);
1258 }
1259
1260 /* Name */
1261 const char *var_name = cs_sles_name(f_id, name);
1262
1263 /* Symmetric matrix, except if advection */
1264 isym = 1;
1265 if (iconvp > 0) isym = 2;
1266
1267 bool symmetric = (isym == 1) ? true : false;
1268
1269 /* be careful here, xam is interleaved*/
1270 if (iesize == 1)
1271 BFT_MALLOC(xam, isym*n_faces, cs_real_t);
1272 if (iesize == 3)
1273 BFT_MALLOC(xam, 3*3*isym*n_faces, cs_real_t);
1274
1275 /*============================================================================
1276 * 1. Building of the "simplified" matrix
1277 *==========================================================================*/
1278
1279 int tensorial_diffusion = 1;
1280
1281 if (idftnp & CS_ANISOTROPIC_LEFT_DIFFUSION)
1282 tensorial_diffusion = 2;
1283
1284 cs_matrix_wrapper_vector(iconvp,
1285 idiffp,
1286 tensorial_diffusion,
1287 ndircp,
1288 isym,
1289 eb_size,
1290 thetap,
1291 coefbv,
1292 cofbfv,
1293 fimp,
1294 i_massflux,
1295 b_massflux,
1296 i_viscm,
1297 b_viscm,
1298 dam,
1299 xam);
1300
1301 /* For steady computations, the diagonal is relaxed */
1302 if (idtvar < 0) {
1303 # pragma omp parallel for
1304 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1305 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1306 for (cs_lnum_t jsou = 0; jsou < 3; jsou++)
1307 dam[iel][isou][jsou] /= relaxp;
1308 }
1309 }
1310 }
1311
1312 /*============================================================================
1313 * 2. Iterative process to handle non orthogonlaities (starting from the
1314 * second iteration).
1315 *===========================================================================*/
1316
1317 /* Application du theta schema */
1318
1319 /* On calcule le bilan explicite total */
1320 thetex = 1. - thetap;
1321
1322 /* Si THETEX=0, ce n'est pas la peine d'en rajouter */
1323 if (fabs(thetex) > cs_math_epzero) {
1324 inc = 1;
1325
1326 /* The added convective scalar mass flux is:
1327 * (thetex*Y_\face-imasac*Y_\celli)*mf.
1328 * When building the explicit part of the rhs, one
1329 * has to impose 0 on mass accumulation. */
1330 imasac = 0;
1331
1332 var_cal_opt->thetav = thetex;
1333
1334 cs_balance_vector(idtvar,
1335 f_id,
1336 imasac,
1337 inc,
1338 ivisep,
1339 var_cal_opt,
1340 NULL, /* pvar == pvara */
1341 pvara,
1342 coefav,
1343 coefbv,
1344 cofafv,
1345 cofbfv,
1346 i_massflux,
1347 b_massflux,
1348 i_visc,
1349 b_visc,
1350 i_secvis,
1351 b_secvis,
1352 viscel,
1353 weighf,
1354 weighb,
1355 icvflb,
1356 icvfli,
1357 smbrp);
1358
1359 var_cal_opt->thetav = thetap;
1360 }
1361
1362 /* Before looping, the RHS without reconstruction is stored in smbini */
1363
1364 cs_lnum_t has_dc = mq->has_disable_flag;
1365 # pragma omp parallel if(n_cells > CS_THR_MIN)
1366 {
1367 # pragma omp for nowait
1368 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1369 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1370 smbini[cell_id][isou] = smbrp[cell_id][isou];
1371 smbrp[cell_id][isou] = 0.;
1372 }
1373 }
1374
1375 /* pvar is initialized on n_cells_ext to avoid a synchronization */
1376 # pragma omp for nowait
1377 for (cs_lnum_t iel = 0; iel < n_cells_ext; iel++) {
1378 for (cs_lnum_t isou = 0; isou < 3; isou++)
1379 pvar[iel][isou] = pvark[iel][isou];
1380 }
1381 }
1382
1383 /* In the following, cs_balance_vector is called with inc=1,
1384 * except for Weight Matrix (nswrsp=-1) */
1385 inc = 1;
1386
1387 if (var_cal_opt->nswrsm == -1) {
1388 var_cal_opt->nswrsm = 1;
1389 inc = 0;
1390 }
1391
1392 /* ---> INCREMENTATION ET RECONSTRUCTION DU SECOND MEMBRE */
1393
1394 /* On est entre avec un smb explicite base sur PVARA.
1395 * si on initialise avec PVAR avec autre chose que PVARA
1396 * on doit donc corriger SMBR (c'est le cas lorsqu'on itere sur navsto) */
1397
1398 /* The added convective scalar mass flux is:
1399 * (thetap*Y_\face-imasac*Y_\celli)*mf.
1400 * When building the implicit part of the rhs, one
1401 * has to impose 1 on mass accumulation. */
1402 imasac = 1;
1403
1404 cs_balance_vector(idtvar,
1405 f_id,
1406 imasac,
1407 inc,
1408 ivisep,
1409 var_cal_opt,
1410 pvar,
1411 pvara,
1412 coefav,
1413 coefbv,
1414 cofafv,
1415 cofbfv,
1416 i_massflux,
1417 b_massflux,
1418 i_visc,
1419 b_visc,
1420 i_secvis,
1421 b_secvis,
1422 viscel,
1423 weighf,
1424 weighb,
1425 icvflb,
1426 icvfli,
1427 smbrp);
1428
1429 if (CS_F_(vel)->id == f_id) {
1430 f = cs_field_by_name_try("velocity_explicit_balance");
1431
1432 if (f != NULL) {
1433 cs_real_3_t *cpro_cv_df_v = (cs_real_3_t *)f->val;
1434 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1435 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1436 for (cs_lnum_t isou = 0; isou < 3; isou++)
1437 cpro_cv_df_v[iel][isou] = smbrp[iel][isou];
1438 }
1439 }
1440 }
1441
1442 /* Dynamic relaxation*/
1443 if (iswdyp >= 1) {
1444 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1445 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1446 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1447 rhs0[iel][isou] = smbrp[iel][isou];
1448 smbini[iel][isou] = smbini[iel][isou]
1449 -fimp[iel][isou][0]*(pvar[iel][0] - pvara[iel][0])
1450 -fimp[iel][isou][1]*(pvar[iel][1] - pvara[iel][1])
1451 -fimp[iel][isou][2]*(pvar[iel][2] - pvara[iel][2]);
1452 smbrp[iel][isou] += smbini[iel][isou];
1453
1454 adxkm1[iel][isou] = 0.;
1455 adxk[iel][isou] = 0.;
1456 dpvar[iel][isou] = 0.;
1457 }
1458 }
1459
1460 /* ||A.dx^0||^2 = 0 */
1461 nadxk = 0.;
1462 }
1463 else {
1464 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1465 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1466 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1467 smbini[iel][isou] = smbini[iel][isou]
1468 -fimp[iel][isou][0]*(pvar[iel][0] - pvara[iel][0])
1469 -fimp[iel][isou][1]*(pvar[iel][1] - pvara[iel][1])
1470 -fimp[iel][isou][2]*(pvar[iel][2] - pvara[iel][2]);
1471 smbrp[iel][isou] += smbini[iel][isou];
1472 }
1473 }
1474 }
1475
1476 /* --- Convergence test */
1477 residu = sqrt(cs_gdot(3*n_cells, (cs_real_t *)smbrp, (cs_real_t *)smbrp));
1478
1479 /* ---> RESIDU DE NORMALISATION
1480 * (NORME C.L +TERMES SOURCES+ TERMES DE NON ORTHOGONALITE) */
1481
1482 /* Allocate a temporary array */
1483 BFT_MALLOC(w1, n_cells_ext, cs_real_3_t);
1484 BFT_MALLOC(w2, n_cells_ext, cs_real_3_t);
1485
1486 cs_real_t *pvar_i;
1487 BFT_MALLOC(pvar_i, n_cells_ext, cs_real_t);
1488
1489 /* Compute the L2 norm of the variable */
1490 for (cs_lnum_t i = 0; i < 3; i++) {
1491
1492 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
1493 pvar_i[cell_id] = pvar[cell_id][i];
1494
1495 cs_real_t p_mean = cs_gmean(n_cells, mq->cell_vol, pvar_i);
1496
1497 if (iwarnp >= 2)
1498 bft_printf("Spatial average of X_%d^n = %f\n", i, p_mean);
1499
1500 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
1501 w2[cell_id][i] = (pvar[cell_id][i] - p_mean);
1502
1503 }
1504 BFT_FREE(pvar_i);
1505
1506 cs_matrix_vector_native_multiply(symmetric,
1507 db_size,
1508 eb_size,
1509 f_id,
1510 (cs_real_t *)dam,
1511 xam,
1512 (cs_real_t *)w2,
1513 (cs_real_t *)w1);
1514
1515 BFT_FREE(w2);
1516
1517 if (iwarnp >= 2) {
1518 const cs_real_t *_w1 = (cs_real_t *)w1, *_smbrp = (cs_real_t *)smbrp;
1519 bft_printf("L2 norm ||AX^n|| = %f\n", sqrt(cs_gdot(3*n_cells, _w1, _w1)));
1520 bft_printf("L2 norm ||B^n|| = %f\n", sqrt(cs_gdot(3*n_cells, _smbrp, _smbrp)));
1521 }
1522
1523 # pragma omp parallel for
1524 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1525 for (int i = 0; i < 3; i++)
1526 w1[cell_id][i] += smbrp[cell_id][i];
1527 /* Remove contributions from penalized cells */
1528 if (has_dc * mq->c_disable_flag[has_dc * cell_id] != 0)
1529 for (int i = 0; i < 3; i++)
1530 w1[cell_id][i] = 0.;
1531 }
1532
1533 rnorm2 = cs_gdot(3*n_cells, (cs_real_t *)w1, (cs_real_t *)w1);
1534 rnorm = sqrt(rnorm2);
1535 sinfo.rhs_norm = rnorm;
1536
1537 /* Free memory */
1538 BFT_FREE(w1);
1539
1540 /* Warning: for Weight Matrix, one and only one sweep is done. */
1541 nswmod = CS_MAX(var_cal_opt->nswrsm, 1);
1542
1543 isweep = 1;
1544
1545 /* Reconstruction loop (beginning)
1546 *-------------------------------- */
1547 if (iterns <= 1)
1548 sinfo.n_it = 0;
1549
1550 while ((isweep <= nswmod && residu > epsrsp*rnorm) || isweep == 1) {
1551 /* --- Solving on the increment dpvar */
1552
1553 /* Dynamic relaxation of the system */
1554 if (iswdyp >= 1) {
1555 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1556 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1557 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1558 dpvarm1[iel][isou] = dpvar[iel][isou];
1559 dpvar[iel][isou] = 0.;
1560 }
1561 }
1562 }
1563 else {
1564 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1565 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1566 for (cs_lnum_t isou = 0; isou < 3; isou++)
1567 dpvar[iel][isou] = 0.;
1568 }
1569 }
1570
1571 /* Matrix block size */
1572 ibsize = 3;
1573
1574 db_size[0] = ibsize;
1575 db_size[1] = ibsize;
1576 db_size[2] = ibsize;
1577 db_size[3] = ibsize*ibsize;
1578
1579 /* Solver residual */
1580 ressol = residu;
1581
1582 cs_sles_solve_native(f_id,
1583 var_name,
1584 symmetric,
1585 db_size,
1586 eb_size,
1587 (cs_real_t *)dam,
1588 xam,
1589 epsilp,
1590 rnorm,
1591 &niterf,
1592 &ressol,
1593 (cs_real_t *)smbrp,
1594 (cs_real_t *)dpvar);
1595
1596 /* Dynamic relaxation of the system */
1597 if (iswdyp >= 1) {
1598
1599 /* Computation of the variable relaxation coefficient */
1600 lvar = -1;
1601
1602 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1603 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1604 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1605 adxkm1[iel][isou] = adxk[iel][isou];
1606 adxk[iel][isou] = - rhs0[iel][isou];
1607 }
1608 }
1609
1610 cs_balance_vector(idtvar,
1611 lvar,
1612 imasac,
1613 inc,
1614 ivisep,
1615 var_cal_opt,
1616 dpvar,
1617 NULL, /* dpvar */
1618 coefav,
1619 coefbv,
1620 cofafv,
1621 cofbfv,
1622 i_massflux,
1623 b_massflux,
1624 i_visc,
1625 b_visc,
1626 i_secvis,
1627 b_secvis,
1628 viscel,
1629 weighf,
1630 weighb,
1631 icvflb,
1632 icvfli,
1633 adxk);
1634
1635 /* ||E.dx^(k-1)-E.0||^2 */
1636 nadxkm1 = nadxk;
1637
1638 /* ||E.dx^k-E.0||^2 */
1639 nadxk = cs_gdot(3*n_cells, (cs_real_t *)adxk, (cs_real_t *)adxk);
1640
1641 /* < E.dx^k-E.0; r^k > */
1642 paxkrk = cs_gdot(3*n_cells, (cs_real_t *)smbrp, (cs_real_t *)adxk);
1643
1644 /* Relaxation with respect to dx^k and dx^(k-1) */
1645 if (iswdyp >= 2) {
1646 /* < E.dx^(k-1)-E.0; r^k > */
1647 paxm1rk = cs_gdot(3*n_cells, (cs_real_t *)smbrp, (cs_real_t *)adxkm1);
1648
1649 /* < E.dx^(k-1)-E.0; E.dx^k-E.0 > */
1650 paxm1ax = cs_gdot(3*n_cells, (cs_real_t *)adxk, (cs_real_t *)adxkm1);
1651
1652 if ((nadxkm1 > 1.e-30*rnorm2)
1653 && (nadxk*nadxkm1-pow(paxm1ax,2)) > 1.e-30*rnorm2)
1654 beta = (paxkrk*paxm1ax - nadxk*paxm1rk)/(nadxk*nadxkm1-pow(paxm1ax,2));
1655 else
1656 beta = 0.;
1657
1658 } else {
1659 beta = 0.;
1660 paxm1ax = 1.;
1661 paxm1rk = 0.;
1662 paxm1ax = 0.;
1663 }
1664
1665 /* The first sweep is not relaxed */
1666 if (isweep == 1) {
1667 alph = 1.;
1668 beta = 0.;
1669 } else if (isweep == 2) {
1670 beta = 0.;
1671 alph = -paxkrk/CS_MAX(nadxk, 1.e-30*rnorm2);
1672 } else {
1673 alph = -(paxkrk + beta*paxm1ax)/CS_MAX(nadxk, 1.e-30*rnorm2);
1674 }
1675
1676 /* Writing */
1677 if (iwarnp >= 3)
1678 bft_printf("%s Sweep: %d Dynamic relaxation: alpha = %12.5e, "
1679 "beta = %12.5e,\n< dI^k : R^k > = %12.5e, "
1680 "||dI^k ||^2 = %12.5e,\n< dI^k-1 : R^k > = %12.5e, "
1681 "||dI^k-1||^2 = %12.5e,\n< dI^k-1 : dI^k > = %12.5e\n",
1682 var_name, isweep, alph, beta, paxkrk, nadxk, paxm1rk,
1683 nadxkm1, paxm1ax);
1684 }
1685
1686 /* --- Update the solution with the increment */
1687
1688 if (iswdyp <= 0) {
1689 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1690 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1691 for (cs_lnum_t isou = 0; isou < 3; isou++)
1692 pvar[iel][isou] += dpvar[iel][isou];
1693 }
1694 }
1695 else if (iswdyp == 1) {
1696 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1697 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1698 for (cs_lnum_t isou = 0; isou < 3; isou++)
1699 pvar[iel][isou] += alph*dpvar[iel][isou];
1700 }
1701 }
1702 else if (iswdyp >= 2) {
1703 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1704 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1705 for (cs_lnum_t isou = 0; isou < 3; isou++)
1706 pvar[iel][isou] += alph*dpvar[iel][isou] + beta*dpvarm1[iel][isou];
1707 }
1708 }
1709
1710 /* --> Handle parallelism and periodicity */
1711
1712 if (cs_glob_rank_id >=0 || cs_glob_mesh->n_init_perio > 0)
1713 cs_mesh_sync_var_vect((cs_real_t *)pvar);
1714
1715 /* --- Update the right hand and compute the new residual */
1716
1717 if (iswdyp <= 0) {
1718 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1719 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1720 /* smbini already contains unsteady terms and mass source terms
1721 * of the RHS updated at each sweep */
1722 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1723 smbini[iel][isou] = smbini[iel][isou]
1724 - fimp[iel][isou][0]*dpvar[iel][0]
1725 - fimp[iel][isou][1]*dpvar[iel][1]
1726 - fimp[iel][isou][2]*dpvar[iel][2];
1727 smbrp[iel][isou] = smbini[iel][isou];
1728 }
1729 }
1730 }
1731 else if (iswdyp == 1) {
1732 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1733 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1734 /* smbini already contains unsteady terms and mass source terms
1735 * of the RHS updated at each sweep */
1736 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1737 smbini[iel][isou] = smbini[iel][isou]
1738 - fimp[iel][isou][0]*alph*dpvar[iel][0]
1739 - fimp[iel][isou][1]*alph*dpvar[iel][1]
1740 - fimp[iel][isou][2]*alph*dpvar[iel][2];
1741 smbrp[iel][isou] = smbini[iel][isou];
1742 }
1743 }
1744 }
1745 else if (iswdyp == 2) {
1746 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1747 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1748 /* smbini already contains unsteady terms and mass source terms
1749 * of the RHS updated at each sweep */
1750 for (cs_lnum_t isou = 0; isou < 3; isou++) {
1751 smbini[iel][isou] = smbini[iel][isou]
1752 - fimp[iel][isou][0]*(alph*dpvar[iel][0]
1753 + beta*dpvarm1[iel][0])
1754 - fimp[iel][isou][1]*(alph*dpvar[iel][1]
1755 + beta*dpvarm1[iel][1])
1756 - fimp[iel][isou][2]*(alph*dpvar[iel][2]
1757 + beta*dpvarm1[iel][2]);
1758 smbrp[iel][isou] = smbini[iel][isou];
1759 }
1760 }
1761 }
1762
1763 /* The added convective scalar mass flux is:
1764 * (thetex*Y_\face-imasac*Y_\celli)*mf.
1765 * When building the implicit part of the rhs, one
1766 * has to impose 1 on mass accumulation. */
1767 imasac = 1;
1768
1769 cs_balance_vector(idtvar,
1770 f_id,
1771 imasac,
1772 inc,
1773 ivisep,
1774 var_cal_opt,
1775 pvar,
1776 pvara,
1777 coefav,
1778 coefbv,
1779 cofafv,
1780 cofbfv,
1781 i_massflux,
1782 b_massflux,
1783 i_visc,
1784 b_visc,
1785 i_secvis,
1786 b_secvis,
1787 viscel,
1788 weighf,
1789 weighb,
1790 icvflb,
1791 icvfli,
1792 smbrp);
1793
1794 /* --- Convergence test */
1795 residu = sqrt(cs_gdot(3*n_cells, (cs_real_t *)smbrp, (cs_real_t *)smbrp));
1796
1797 /* Writing */
1798 sinfo.n_it = sinfo.n_it + niterf;
1799
1800 /* Writing */
1801 if (iwarnp >= 2) {
1802 bft_printf("%s: CV_DIF_TS, IT: %d, Res: %12.5e, Norm: %12.5e\n",
1803 var_name, isweep, residu, rnorm);
1804 bft_printf("%s: Current reconstruction sweep: %d, "
1805 "Iterations for solver: %d\n", var_name, isweep, niterf);
1806 }
1807
1808 isweep++;
1809 }
1810
1811 /* --- Reconstruction loop (end) */
1812
1813 /* Writing: convergence */
1814 if (fabs(rnorm)/sqrt(3.) > cs_math_epzero)
1815 sinfo.res_norm = residu/rnorm;
1816 else
1817 sinfo.res_norm = 0.;
1818
1819 if (iwarnp >= 1) {
1820 if (residu <= epsrsp*rnorm)
1821 bft_printf("%s : CV_DIF_TS, IT : %d, Res : %12.5e, Norm : %12.5e\n",
1822 var_name, isweep-1, residu, rnorm);
1823 /* Writing: non-convergence */
1824 else if (isweep > nswmod)
1825 bft_printf("@\n@ @@ WARNING: %s CONVECTION-DIFFUSION-SOURCE-TERMS\n@"
1826 "=========\n@ Maximum number of iterations %d reached\n@",
1827 var_name,nswmod);
1828 }
1829
1830 /* Save convergence info for fields */
1831 if (f_id > -1) {
1832 f = cs_field_by_id(f_id);
1833 cs_field_set_key_struct(f, key_sinfo_id, &sinfo);
1834 }
1835
1836 /*============================================================================
1837 * 3. After having computed the new value, an estimator is computed for the
1838 * prediction step of the velocity.
1839 *==========================================================================*/
1840
1841 if (iescap > 0) {
1842 /* ---> Computation of the estimator of the current component */
1843
1844 /* smbini already contains unsteady terms and mass source terms
1845 of the RHS updated at each sweep */
1846
1847 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1848 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1849 for (cs_lnum_t isou = 0; isou < 3; isou++)
1850 smbrp[iel][isou] = smbini[iel][isou] - fimp[iel][isou][0]*dpvar[iel][0]
1851 - fimp[iel][isou][1]*dpvar[iel][1]
1852 - fimp[iel][isou][2]*dpvar[iel][2];
1853 }
1854
1855 inc = 1;
1856
1857 /* Without relaxation even for a stationnary computation */
1858
1859 cs_balance_vector(idtvar,
1860 f_id,
1861 imasac,
1862 inc,
1863 ivisep,
1864 var_cal_opt,
1865 pvar,
1866 pvara,
1867 coefav,
1868 coefbv,
1869 cofafv,
1870 cofbfv,
1871 i_massflux,
1872 b_massflux,
1873 i_visc,
1874 b_visc,
1875 i_secvis,
1876 b_secvis,
1877 viscel,
1878 weighf,
1879 weighb,
1880 icvflb,
1881 icvfli,
1882 smbrp);
1883
1884 /* Contribution of the current component to the L2 norm stored in eswork */
1885
1886 # pragma omp parallel for if(n_cells > CS_THR_MIN)
1887 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
1888 for (cs_lnum_t isou = 0; isou < 3; isou++)
1889 eswork[iel][isou] = pow(smbrp[iel][isou] / cell_vol[iel],2);
1890 }
1891 }
1892
1893 /*==========================================================================
1894 * 4. Free solver setup
1895 *==========================================================================*/
1896
1897 cs_sles_free_native(f_id, var_name);
1898
1899 /* Save diagonal in case we want to use it */
1900 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1901 for (cs_lnum_t i = 0; i < 3; i++)
1902 for (cs_lnum_t j = 0; j < 3; j++)
1903 fimp[cell_id][i][j] = dam[cell_id][i][j];
1904 }
1905
1906 /* Free memory */
1907 BFT_FREE(dam);
1908 BFT_FREE(xam);
1909 BFT_FREE(smbini);
1910 BFT_FREE(dpvar);
1911 if (iswdyp >= 1) {
1912 BFT_FREE(adxk);
1913 BFT_FREE(adxkm1);
1914 BFT_FREE(dpvarm1);
1915 BFT_FREE(rhs0);
1916 }
1917 }
1918
1919 /*----------------------------------------------------------------------------*/
1920 /*! \brief This function solves an advection diffusion equation with source
1921 * terms for one time step for the symmetric tensor variable
1922 * \f$ \tens{\variat} \f$.
1923 *
1924 * The equation reads:
1925 *
1926 * \f[
1927 * \tens{f_s}^{imp}(\tens{\variat}^{n+1}-\tens{\variat}^n)
1928 * + \divt \left( \tens{\variat}^{n+1} \otimes \rho \vect {u}
1929 * - \mu \gradtt \tens{\variat}^{n+1}\right)
1930 * = \tens{Rhs}
1931 * \f]
1932 *
1933 * This equation is rewritten as:
1934 *
1935 * \f[
1936 * \tens{f_s}^{imp} \delta \tens{\variat}
1937 * + \divt \left( \delta \tens{\variat} \otimes \rho \vect{u}
1938 * - \mu \gradtt \delta \tens{\variat} \right)
1939 * = \tens{Rhs}^1
1940 * \f]
1941 *
1942 * where \f$ \delta \tens{\variat} = \tens{\variat}^{n+1} - \tens{\variat}^n\f$
1943 * and \f$ \tens{Rhs}^1 = \tens{Rhs}
1944 * - \divt \left( \tens{\variat}^n \otimes \rho \vect{u}
1945 * - \mu \gradtt \tens{\variat}^n \right)\f$
1946 *
1947 *
1948 * It is in fact solved with the following iterative process:
1949 *
1950 * \f[
1951 * \tens{f_s}^{imp} \delta \tens{\variat}^k
1952 * + \divt \left( \delta \tens{\variat}^k \otimes \rho \vect{u}
1953 * - \mu \gradtt \delta \tens{\variat}^k \right)
1954 * = \tens{Rhs}^k
1955 * \f]
1956 *
1957 * where \f$ \tens{Rhs}^k = \tens{Rhs}
1958 * - \tens{f_s}^{imp} \left(\tens{\variat}^k-\tens{\variat}^n \right)
1959 * - \divt \left( \tens{\variat}^k \otimes \rho \vect{u}
1960 * - \mu \gradtt \tens{\variat}^k \right)\f$
1961 *
1962 * Be careful, it is forbidden to modify \f$ \tens{f_s}^{imp} \f$ here!
1963 *
1964 * \param[in] idtvar indicator of the temporal scheme
1965 * \param[in] f_id field id (or -1)
1966 * \param[in] name associated name if f_id < 0, or NULL
1967 * \param[in] var_cal_opt pointer to a cs_var_cal_opt_t structure which
1968 * contains variable calculation options
1969 * \param[in] pvara variable at the previous time step
1970 * \f$ \vect{a}^n \f$
1971 * \param[in] pvark variable at the previous sub-iteration
1972 * \f$ \vect{a}^k \f$.
1973 * If you sub-iter on Navier-Stokes, then
1974 * it allows to initialize by something else than
1975 * pvara (usually pvar=pvara)
1976 * \param[in] coefats boundary condition array for the variable
1977 * (Explicit part)
1978 * \param[in] coefbts boundary condition array for the variable
1979 * (Impplicit part)
1980 * \param[in] cofafts boundary condition array for the diffusion
1981 * of the variable (Explicit part)
1982 * \param[in] cofbfts boundary condition array for the diffusion
1983 * of the variable (Implicit part)
1984 * \param[in] i_massflux mass flux at interior faces
1985 * \param[in] b_massflux mass flux at boundary faces
1986 * \param[in] i_viscm \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1987 * at interior faces for the matrix
1988 * \param[in] b_viscm \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1989 * at boundary faces for the matrix
1990 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1991 * at interior faces for the r.h.s.
1992 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1993 * at boundary faces for the r.h.s.
1994 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
1995 * \param[in] weighf internal face weight between cells i j in case
1996 * of tensor diffusion
1997 * \param[in] weighb boundary face weight for cells i in case
1998 * of tensor diffusion
1999 * \param[in] icvflb global indicator of boundary convection flux
2000 * - 0 upwind scheme at all boundary faces
2001 * - 1 imposed flux at some boundary faces
2002 * \param[in] icvfli boundary face indicator array of convection flux
2003 * - 0 upwind scheme
2004 * - 1 imposed flux
2005 * \param[in] fimp \f$ \tens{f_s}^{imp} \f$
2006 * \param[in] smbrp Right hand side \f$ \vect{Rhs}^k \f$
2007 * \param[in,out] pvar current variable
2008 */
2009 /*----------------------------------------------------------------------------*/
2010
2011 void
cs_equation_iterative_solve_tensor(int idtvar,int f_id,const char * name,cs_var_cal_opt_t * var_cal_opt,const cs_real_6_t pvara[],const cs_real_6_t pvark[],const cs_real_6_t coefats[],const cs_real_66_t coefbts[],const cs_real_6_t cofafts[],const cs_real_66_t cofbfts[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_viscm[],const cs_real_t b_viscm[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_6_t viscel[],const cs_real_2_t weighf[],const cs_real_t weighb[],int icvflb,const int icvfli[],const cs_real_66_t fimp[],cs_real_6_t smbrp[],cs_real_6_t pvar[])2012 cs_equation_iterative_solve_tensor(int idtvar,
2013 int f_id,
2014 const char *name,
2015 cs_var_cal_opt_t *var_cal_opt,
2016 const cs_real_6_t pvara[],
2017 const cs_real_6_t pvark[],
2018 const cs_real_6_t coefats[],
2019 const cs_real_66_t coefbts[],
2020 const cs_real_6_t cofafts[],
2021 const cs_real_66_t cofbfts[],
2022 const cs_real_t i_massflux[],
2023 const cs_real_t b_massflux[],
2024 const cs_real_t i_viscm[],
2025 const cs_real_t b_viscm[],
2026 const cs_real_t i_visc[],
2027 const cs_real_t b_visc[],
2028 cs_real_6_t viscel[],
2029 const cs_real_2_t weighf[],
2030 const cs_real_t weighb[],
2031 int icvflb,
2032 const int icvfli[],
2033 const cs_real_66_t fimp[],
2034 cs_real_6_t smbrp[],
2035 cs_real_6_t pvar[])
2036 {
2037 /* Local variables */
2038
2039 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
2040
2041 int iconvp = var_cal_opt->iconv;
2042 int idiffp = var_cal_opt->idiff;
2043 int iwarnp = var_cal_opt->verbosity;
2044 int iswdyp = var_cal_opt->iswdyn;
2045 int idftnp = var_cal_opt->idften;
2046 int ndircp = var_cal_opt->ndircl;
2047 double epsrsp = var_cal_opt->epsrsm;
2048 double epsilp = var_cal_opt->epsilo;
2049 double relaxp = var_cal_opt->relaxv;
2050 double thetap = var_cal_opt->thetav;
2051
2052 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
2053 const cs_lnum_t n_faces = cs_glob_mesh->n_i_faces;
2054 const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
2055
2056 int isym, inc, isweep, niterf, nswmod, ibsize;
2057 int key_sinfo_id;
2058 int iesize, lvar, imasac;
2059 double residu, rnorm, ressol, rnorm2, thetex, alph, beta;
2060 double paxkrk, nadxk, paxm1rk, nadxkm1, paxm1ax;
2061
2062 cs_lnum_t eb_size[4], db_size[4];
2063
2064 cs_solving_info_t sinfo;
2065
2066 cs_field_t *f;
2067
2068 cs_real_t *xam;
2069 cs_real_66_t *dam;
2070 cs_real_6_t *dpvar, *smbini, *w1, *w2, *adxk, *adxkm1, *dpvarm1, *rhs0;
2071
2072 /*============================================================================
2073 * 0. Initialization
2074 *==========================================================================*/
2075
2076 /* Matrix block size */
2077 ibsize = 6;
2078 iesize = 1; /* CS_ISOTROPIC_DIFFUSION or CS_ANISOTROPIC_RIGHT_DIFFUSION */
2079 if (idftnp & CS_ANISOTROPIC_LEFT_DIFFUSION) iesize = 6;
2080
2081 db_size[0] = ibsize;
2082 db_size[1] = ibsize;
2083 db_size[2] = ibsize;
2084 db_size[3] = ibsize*ibsize;
2085
2086 eb_size[0] = iesize;
2087 eb_size[1] = iesize;
2088 eb_size[2] = iesize;
2089 eb_size[3] = iesize*iesize;
2090
2091 /* Allocate temporary arrays */
2092 BFT_MALLOC(dam, n_cells_ext, cs_real_66_t);
2093 BFT_MALLOC(dpvar, n_cells_ext, cs_real_6_t);
2094 BFT_MALLOC(smbini, n_cells_ext, cs_real_6_t);
2095
2096 if (iswdyp >= 1) {
2097 BFT_MALLOC(adxk, n_cells_ext, cs_real_6_t);
2098 BFT_MALLOC(adxkm1, n_cells_ext, cs_real_6_t);
2099 BFT_MALLOC(dpvarm1, n_cells_ext, cs_real_6_t);
2100 BFT_MALLOC(rhs0, n_cells_ext, cs_real_6_t);
2101 }
2102
2103 /* solving info */
2104 key_sinfo_id = cs_field_key_id("solving_info");
2105 if (f_id > -1) {
2106 f = cs_field_by_id(f_id);
2107 cs_field_get_key_struct(f, key_sinfo_id, &sinfo);
2108 }
2109
2110 /* Name */
2111 const char *var_name = cs_sles_name(f_id, name);
2112
2113 /* Symmetric matrix, except if advection */
2114 isym = 1;
2115 if (iconvp > 0) isym = 2;
2116
2117 bool symmetric = (isym == 1) ? true : false;
2118
2119 /* be carefull here, xam is interleaved*/
2120 if (iesize == 1)
2121 BFT_MALLOC(xam, isym*n_faces, cs_real_t);
2122 if (iesize == 6)
2123 BFT_MALLOC(xam, 6*6*isym*n_faces, cs_real_t);
2124
2125 /*============================================================================
2126 * 1. Building of the "simplified" matrix
2127 *==========================================================================*/
2128
2129 int tensorial_diffusion = 1;
2130 if (iesize == 6)
2131 tensorial_diffusion = 2;
2132
2133 cs_matrix_wrapper_tensor(iconvp,
2134 idiffp,
2135 tensorial_diffusion,
2136 ndircp,
2137 isym,
2138 thetap,
2139 coefbts,
2140 cofbfts,
2141 fimp,
2142 i_massflux,
2143 b_massflux,
2144 i_viscm,
2145 b_viscm,
2146 dam,
2147 xam);
2148
2149 /* For steady computations, the diagonal is relaxed */
2150 if (idtvar < 0) {
2151 # pragma omp parallel for if(n_cells > CS_THR_MIN)
2152 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2153 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2154 for (cs_lnum_t jsou = 0; jsou < 6; jsou++)
2155 dam[iel][isou][jsou] /= relaxp;
2156 }
2157 }
2158 }
2159
2160 /*============================================================================
2161 * 2. Iterative process to handle non orthogonlaities (starting from the
2162 * second iteration).
2163 *===========================================================================*/
2164
2165 /* Application du theta schema */
2166
2167 /* On calcule le bilan explicite total */
2168 thetex = 1. - thetap;
2169
2170 /* Si THETEX=0, ce n'est pas la peine d'en rajouter */
2171 if (fabs(thetex) > cs_math_epzero) {
2172 inc = 1;
2173
2174 /* The added convective scalar mass flux is:
2175 * (thetex*Y_\face-imasac*Y_\celli)*mf.
2176 * When building the explicit part of the rhs, one
2177 * has to impose 0 on mass accumulation. */
2178 imasac = 0;
2179
2180 var_cal_opt->thetav = thetex;
2181
2182 cs_balance_tensor(idtvar,
2183 f_id,
2184 imasac,
2185 inc,
2186 var_cal_opt,
2187 NULL, /* pvar == pvara */
2188 pvara,
2189 coefats,
2190 coefbts,
2191 cofafts,
2192 cofbfts,
2193 i_massflux,
2194 b_massflux,
2195 i_visc,
2196 b_visc,
2197 viscel,
2198 weighf,
2199 weighb,
2200 icvflb,
2201 icvfli,
2202 smbrp);
2203
2204 var_cal_opt->thetav = thetap;
2205 }
2206
2207 /* Before looping, the RHS without reconstruction is stored in smbini */
2208
2209 cs_lnum_t has_dc = mq->has_disable_flag;
2210 # pragma omp parallel if(n_cells > CS_THR_MIN)
2211 {
2212 # pragma omp for nowait
2213 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
2214 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2215 smbini[cell_id][isou] = smbrp[cell_id][isou];
2216 smbrp[cell_id][isou] = 0.;
2217 }
2218 }
2219
2220 /* pvar is initialized on n_cells_ext to avoid a synchronization */
2221 # pragma omp for nowait
2222 for (cs_lnum_t iel = 0; iel < n_cells_ext; iel++) {
2223 for (cs_lnum_t isou = 0; isou < 6; isou++)
2224 pvar[iel][isou] = pvark[iel][isou];
2225 }
2226 }
2227
2228 /* In the following, cs_balance_vector is called with inc=1,
2229 * except for Weight Matrix (nswrsp=-1) */
2230 inc = 1;
2231
2232 if (var_cal_opt->nswrsm == -1) {
2233 var_cal_opt->nswrsm = 1;
2234 inc = 0;
2235 }
2236
2237 /* ---> INCREMENTATION ET RECONSTRUCTION DU SECOND MEMBRE */
2238
2239 /* On est entre avec un smb explicite base sur PVARA.
2240 * si on initialise avec PVAR avec autre chose que PVARA
2241 * on doit donc corriger SMBR (c'est le cas lorsqu'on itere sur navsto) */
2242
2243 /* The added convective scalar mass flux is:
2244 * (thetap*Y_\face-imasac*Y_\celli)*mf.
2245 * When building the implicit part of the rhs, one
2246 * has to impose 1 on mass accumulation. */
2247 imasac = 1;
2248
2249 cs_balance_tensor(idtvar,
2250 f_id,
2251 imasac,
2252 inc,
2253 var_cal_opt,
2254 pvar,
2255 pvara,
2256 coefats,
2257 coefbts,
2258 cofafts,
2259 cofbfts,
2260 i_massflux,
2261 b_massflux,
2262 i_visc,
2263 b_visc,
2264 viscel,
2265 weighf,
2266 weighb,
2267 icvflb,
2268 icvfli,
2269 smbrp);
2270
2271 /* Dynamic relaxation*/
2272 if (iswdyp >= 1) {
2273 # pragma omp parallel for
2274 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2275 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2276 rhs0[iel][isou] = smbrp[iel][isou];
2277 smbini[iel][isou] = smbini[iel][isou]
2278 -fimp[iel][isou][0]*(pvar[iel][0] - pvara[iel][0])
2279 -fimp[iel][isou][1]*(pvar[iel][1] - pvara[iel][1])
2280 -fimp[iel][isou][2]*(pvar[iel][2] - pvara[iel][2])
2281 -fimp[iel][isou][3]*(pvar[iel][3] - pvara[iel][3])
2282 -fimp[iel][isou][4]*(pvar[iel][4] - pvara[iel][4])
2283 -fimp[iel][isou][5]*(pvar[iel][5] - pvara[iel][5]);
2284 smbrp[iel][isou] += smbini[iel][isou];
2285
2286 adxkm1[iel][isou] = 0.;
2287 adxk[iel][isou] = 0.;
2288 dpvar[iel][isou] = 0.;
2289 }
2290
2291 /* ||A.dx^0||^2 = 0 */
2292 nadxk = 0.;
2293 }
2294 else {
2295 # pragma omp parallel for
2296 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2297 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2298 smbini[iel][isou] = smbini[iel][isou]
2299 - fimp[iel][isou][0]*(pvar[iel][0] - pvara[iel][0])
2300 - fimp[iel][isou][1]*(pvar[iel][1] - pvara[iel][1])
2301 - fimp[iel][isou][2]*(pvar[iel][2] - pvara[iel][2])
2302 - fimp[iel][isou][3]*(pvar[iel][3] - pvara[iel][3])
2303 - fimp[iel][isou][4]*(pvar[iel][4] - pvara[iel][4])
2304 - fimp[iel][isou][5]*(pvar[iel][5] - pvara[iel][5]);
2305 smbrp[iel][isou] += smbini[iel][isou];
2306 }
2307 }
2308 }
2309
2310 /* --- Convergence test */
2311 residu = sqrt(cs_gdot(6*n_cells, (cs_real_t *)smbrp, (cs_real_t *)smbrp));
2312
2313 /* ---> RESIDU DE NORMALISATION
2314 * (NORME C.L +TERMES SOURCES+ TERMES DE NON ORTHOGONALITE) */
2315
2316 /* Allocate a temporary array */
2317 BFT_MALLOC(w1, n_cells_ext, cs_real_6_t);
2318 BFT_MALLOC(w2, n_cells_ext, cs_real_6_t);
2319
2320 cs_real_t *pvar_i;
2321 BFT_MALLOC(pvar_i, n_cells_ext, cs_real_t);
2322
2323 /* Compute the L2 norm of the variable */
2324 for (cs_lnum_t i = 0; i < 6; i++) {
2325
2326 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
2327 pvar_i[cell_id] = pvar[cell_id][i];
2328
2329 cs_real_t p_mean = cs_gmean(n_cells, mq->cell_vol, pvar_i);
2330
2331 if (iwarnp >= 2)
2332 bft_printf("Spatial average of X_%d^n = %f\n", i, p_mean);
2333
2334 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++)
2335 w2[cell_id][i] = (pvar[cell_id][i] - p_mean);
2336
2337 }
2338 BFT_FREE(pvar_i);
2339
2340 cs_matrix_vector_native_multiply(symmetric,
2341 db_size,
2342 eb_size,
2343 f_id,
2344 (cs_real_t *)dam,
2345 xam,
2346 (cs_real_t *)w2,
2347 (cs_real_t *)w1);
2348
2349 BFT_FREE(w2);
2350
2351 if (iwarnp >= 2) {
2352 const cs_real_t *_w1 = (cs_real_t *)w1, *_smbrp = (cs_real_t *)smbrp;
2353 bft_printf("L2 norm ||AX^n|| = %f\n", sqrt(cs_gdot(3*n_cells, _w1, _w1)));
2354 bft_printf("L2 norm ||B^n|| = %f\n",
2355 sqrt(cs_gdot(3*n_cells, _smbrp, _smbrp)));
2356 }
2357
2358
2359 # pragma omp parallel for
2360 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
2361 for (cs_lnum_t i = 0; i < 6; i++)
2362 w1[cell_id][i] += smbrp[cell_id][i];
2363 /* Remove contributions from penalized cells */
2364 if (has_dc * mq->c_disable_flag[has_dc * cell_id] != 0) {
2365 for (cs_lnum_t i = 0; i < 6; i++)
2366 w1[cell_id][i] = 0.;
2367 }
2368 }
2369
2370 rnorm2 = cs_gdot(6*n_cells, (cs_real_t *)w1, (cs_real_t *)w1);
2371 rnorm = sqrt(rnorm2);
2372 sinfo.rhs_norm = rnorm;
2373
2374 /* Free memory */
2375 BFT_FREE(w1);
2376
2377 /* Warning: for Weight Matrix, one and only one sweep is done. */
2378 nswmod = CS_MAX(var_cal_opt->nswrsm, 1);
2379
2380 isweep = 1;
2381
2382 /* Reconstruction loop (beginning)
2383 *-------------------------------- */
2384 sinfo.n_it = 0;
2385
2386 while ((isweep <= nswmod && residu > epsrsp*rnorm) || isweep == 1) {
2387 /* --- Solving on the increment dpvar */
2388
2389 /* Dynamic relaxation of the system */
2390 if (iswdyp >= 1) {
2391 # pragma omp parallel for
2392 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2393 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2394 dpvarm1[iel][isou] = dpvar[iel][isou];
2395 dpvar[iel][isou] = 0.;
2396 }
2397 }
2398 }
2399 else {
2400 # pragma omp parallel for
2401 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2402 for (cs_lnum_t isou = 0; isou < 6; isou++)
2403 dpvar[iel][isou] = 0.;
2404 }
2405 }
2406
2407 /* Matrix block size */
2408 ibsize = 6;
2409
2410 db_size[0] = ibsize;
2411 db_size[1] = ibsize;
2412 db_size[2] = ibsize;
2413 db_size[3] = ibsize*ibsize;
2414
2415 /* Solver residual */
2416 ressol = residu;
2417
2418 cs_sles_solve_native(f_id,
2419 var_name,
2420 symmetric,
2421 db_size,
2422 eb_size,
2423 (cs_real_t *)dam,
2424 xam,
2425 epsilp,
2426 rnorm,
2427 &niterf,
2428 &ressol,
2429 (cs_real_t *)smbrp,
2430 (cs_real_t *)dpvar);
2431
2432 /* Dynamic relaxation of the system */
2433 if (iswdyp >= 1) {
2434
2435 /* Computation of the variable relaxation coefficient */
2436 lvar = -1;
2437
2438 # pragma omp parallel for
2439 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2440 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2441 adxkm1[iel][isou] = adxk[iel][isou];
2442 adxk[iel][isou] = - rhs0[iel][isou];
2443 }
2444 }
2445
2446 cs_balance_tensor(idtvar,
2447 lvar,
2448 imasac,
2449 inc,
2450 var_cal_opt,
2451 dpvar,
2452 NULL, /* dpvar */
2453 coefats,
2454 coefbts,
2455 cofafts,
2456 cofbfts,
2457 i_massflux,
2458 b_massflux,
2459 i_visc,
2460 b_visc,
2461 viscel,
2462 weighf,
2463 weighb,
2464 icvflb,
2465 icvfli,
2466 adxk);
2467
2468 /* ||E.dx^(k-1)-E.0||^2 */
2469 nadxkm1 = nadxk;
2470
2471 /* ||E.dx^k-E.0||^2 */
2472 nadxk = cs_gdot(6*n_cells, (cs_real_t *)adxk, (cs_real_t *)adxk);
2473
2474 /* < E.dx^k-E.0; r^k > */
2475 paxkrk = cs_gdot(6*n_cells, (cs_real_t *)smbrp, (cs_real_t *)adxk);
2476
2477 /* Relaxation with respect to dx^k and dx^(k-1) */
2478 if (iswdyp >= 2) {
2479 /* < E.dx^(k-1)-E.0; r^k > */
2480 paxm1rk = cs_gdot(6*n_cells, (cs_real_t *)smbrp, (cs_real_t *)adxkm1);
2481
2482 /* < E.dx^(k-1)-E.0; E.dx^k-E.0 > */
2483 paxm1ax = cs_gdot(6*n_cells, (cs_real_t *)adxk, (cs_real_t *)adxkm1);
2484
2485 if ((nadxkm1 > 1.e-30*rnorm2)
2486 && (nadxk*nadxkm1-pow(paxm1ax,2)) > 1.e-30*rnorm2)
2487 beta = (paxkrk*paxm1ax - nadxk*paxm1rk)/(nadxk*nadxkm1-pow(paxm1ax,2));
2488 else
2489 beta = 0.;
2490
2491 } else {
2492 beta = 0.;
2493 paxm1ax = 1.;
2494 paxm1rk = 0.;
2495 paxm1ax = 0.;
2496 }
2497
2498 /* The first sweep is not relaxed */
2499 if (isweep == 1) {
2500 alph = 1.;
2501 beta = 0.;
2502 } else if (isweep == 2) {
2503 beta = 0.;
2504 alph = -paxkrk/CS_MAX(nadxk, 1.e-30*rnorm2);
2505 } else {
2506 alph = -(paxkrk + beta*paxm1ax)/CS_MAX(nadxk, 1.e-30*rnorm2);
2507 }
2508
2509 /* Writing */
2510 if (iwarnp >= 3)
2511 bft_printf("%s Sweep: %d Dynamic relaxation: alpha = %12.5e, "
2512 "beta = %12.5e,\n< dI^k : R^k > = %12.5e, "
2513 "||dI^k ||^2 = %12.5e,\n< dI^k-1 : R^k > = %12.5e, "
2514 "||dI^k-1||^2 = %12.5e,\n< dI^k-1 : dI^k > = %12.5e\n",
2515 var_name, isweep, alph, beta, paxkrk, nadxk, paxm1rk,
2516 nadxkm1, paxm1ax);
2517 }
2518
2519 /* --- Update the solution with the increment */
2520
2521 if (iswdyp <= 0) {
2522 # pragma omp parallel for
2523 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2524 for (cs_lnum_t isou = 0; isou < 6; isou++)
2525 pvar[iel][isou] += dpvar[iel][isou];
2526 }
2527 else if (iswdyp == 1) {
2528 # pragma omp parallel for
2529 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2530 for (cs_lnum_t isou = 0; isou < 6; isou++)
2531 pvar[iel][isou] += alph*dpvar[iel][isou];
2532 }
2533 else if (iswdyp >= 2) {
2534 # pragma omp parallel for
2535 for (cs_lnum_t iel = 0; iel < n_cells; iel++)
2536 for (cs_lnum_t isou = 0; isou < 6; isou++)
2537 pvar[iel][isou] += alph*dpvar[iel][isou] + beta*dpvarm1[iel][isou];
2538 }
2539
2540 /* --> Handle parallelism and periodicity */
2541
2542 if (cs_glob_rank_id >=0 || cs_glob_mesh->n_init_perio > 0)
2543 cs_mesh_sync_var_sym_tens(pvar);
2544
2545 /* --- Update the right hand and compute the new residual */
2546
2547 if (iswdyp <= 0) {
2548 # pragma omp parallel for
2549 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2550 /* smbini already contains unsteady terms and mass source terms
2551 * of the RHS updated at each sweep */
2552 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2553 smbini[iel][isou] = smbini[iel][isou]
2554 - fimp[iel][isou][0]*dpvar[iel][0]
2555 - fimp[iel][isou][1]*dpvar[iel][1]
2556 - fimp[iel][isou][2]*dpvar[iel][2]
2557 - fimp[iel][isou][3]*dpvar[iel][3]
2558 - fimp[iel][isou][4]*dpvar[iel][4]
2559 - fimp[iel][isou][5]*dpvar[iel][5];
2560 smbrp[iel][isou] = smbini[iel][isou];
2561 }
2562 }
2563 }
2564 else if (iswdyp == 1) {
2565 # pragma omp parallel for
2566 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2567 /* smbini already contains unsteady terms and mass source terms
2568 * of the RHS updated at each sweep */
2569 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2570 smbini[iel][isou] = smbini[iel][isou]
2571 - fimp[iel][isou][0]*alph*dpvar[iel][0]
2572 - fimp[iel][isou][1]*alph*dpvar[iel][1]
2573 - fimp[iel][isou][2]*alph*dpvar[iel][2]
2574 - fimp[iel][isou][3]*alph*dpvar[iel][3]
2575 - fimp[iel][isou][4]*alph*dpvar[iel][4]
2576 - fimp[iel][isou][5]*alph*dpvar[iel][5];
2577 smbrp[iel][isou] = smbini[iel][isou];
2578 }
2579 }
2580 }
2581 else if (iswdyp == 2) {
2582 # pragma omp parallel for
2583 for (cs_lnum_t iel = 0; iel < n_cells; iel++) {
2584 /* smbini already contains unsteady terms and mass source terms
2585 * of the RHS updated at each sweep */
2586 for (cs_lnum_t isou = 0; isou < 6; isou++) {
2587 smbini[iel][isou] = smbini[iel][isou]
2588 - fimp[iel][isou][0]*(alph*dpvar[iel][0]
2589 + beta*dpvarm1[iel][0])
2590 - fimp[iel][isou][1]*(alph*dpvar[iel][1]
2591 + beta*dpvarm1[iel][1])
2592 - fimp[iel][isou][2]*(alph*dpvar[iel][2]
2593 + beta*dpvarm1[iel][2])
2594 - fimp[iel][isou][3]*(alph*dpvar[iel][3]
2595 + beta*dpvarm1[iel][3])
2596 - fimp[iel][isou][4]*(alph*dpvar[iel][4]
2597 + beta*dpvarm1[iel][4])
2598 - fimp[iel][isou][5]*(alph*dpvar[iel][5]
2599 + beta*dpvarm1[iel][5]);
2600 smbrp[iel][isou] = smbini[iel][isou];
2601 }
2602 }
2603 }
2604
2605 /* The added convective scalar mass flux is:
2606 * (thetex*Y_\face-imasac*Y_\celli)*mf.
2607 * When building the implicit part of the rhs, one
2608 * has to impose 1 on mass accumulation. */
2609 imasac = 1;
2610
2611 cs_balance_tensor(idtvar,
2612 f_id,
2613 imasac,
2614 inc,
2615 var_cal_opt,
2616 pvar,
2617 pvara,
2618 coefats,
2619 coefbts,
2620 cofafts,
2621 cofbfts,
2622 i_massflux,
2623 b_massflux,
2624 i_visc,
2625 b_visc,
2626 viscel,
2627 weighf,
2628 weighb,
2629 icvflb,
2630 icvfli,
2631 smbrp);
2632
2633 /* --- Convergence test */
2634 residu = sqrt(cs_gdot(6*n_cells, (cs_real_t *)smbrp, (cs_real_t *)smbrp));
2635
2636 /* Writing */
2637 sinfo.n_it = sinfo.n_it + niterf;
2638
2639 /* Writing */
2640 if (iwarnp >= 2) {
2641 bft_printf("%s: CV_DIF_TS, IT: %d, Res: %12.5e, Norm: %12.5e\n",
2642 var_name, isweep, residu, rnorm);
2643 bft_printf("%s: Current reconstruction sweep: %d, "
2644 "Iterations for solver: %d\n", var_name, isweep, niterf);
2645 }
2646
2647 isweep++;
2648 }
2649
2650 /* --- Reconstruction loop (end) */
2651
2652 /* Writing: convergence */
2653 if (fabs(rnorm)/sqrt(6.) > cs_math_epzero)
2654 sinfo.res_norm = residu/rnorm;
2655 else
2656 sinfo.res_norm = 0.;
2657
2658 if (iwarnp >= 1) {
2659 if (residu <= epsrsp*rnorm)
2660 bft_printf("%s : CV_DIF_TS, IT : %d, Res : %12.5e, Norm : %12.5e\n",
2661 var_name, isweep-1, residu, rnorm);
2662 /* Writing: non-convergence */
2663 else if (isweep > nswmod)
2664 bft_printf("@\n@ @@ WARNING: %s CONVECTION-DIFFUSION-SOURCE-TERMS\n@"
2665 "=========\n@ Maximum number of iterations %d reached\n@",
2666 var_name,nswmod);
2667 }
2668
2669 /* Save convergence info for fields */
2670 if (f_id > -1) {
2671 f = cs_field_by_id(f_id);
2672 cs_field_set_key_struct(f, key_sinfo_id, &sinfo);
2673 }
2674
2675 /*==========================================================================
2676 * 3. Free solver setup
2677 *==========================================================================*/
2678
2679 cs_sles_free_native(f_id, var_name);
2680
2681 /* Free memory */
2682 BFT_FREE(dam);
2683 BFT_FREE(xam);
2684 BFT_FREE(smbini);
2685 BFT_FREE(dpvar);
2686 if (iswdyp >= 1) {
2687 BFT_FREE(adxk);
2688 BFT_FREE(adxkm1);
2689 BFT_FREE(dpvarm1);
2690 BFT_FREE(rhs0);
2691 }
2692 }
2693
2694 /*----------------------------------------------------------------------------*/
2695
2696 END_C_DECLS
2697