1 #ifndef __CS_GRADIENT_H__ 2 #define __CS_GRADIENT_H__ 3 4 /*============================================================================ 5 * Gradient reconstruction. 6 *============================================================================*/ 7 8 /* 9 This file is part of Code_Saturne, a general-purpose CFD tool. 10 11 Copyright (C) 1998-2021 EDF S.A. 12 13 This program is free software; you can redistribute it and/or modify it under 14 the terms of the GNU General Public License as published by the Free Software 15 Foundation; either version 2 of the License, or (at your option) any later 16 version. 17 18 This program is distributed in the hope that it will be useful, but WITHOUT 19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 21 details. 22 23 You should have received a copy of the GNU General Public License along with 24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 25 Street, Fifth Floor, Boston, MA 02110-1301, USA. 26 */ 27 28 /*----------------------------------------------------------------------------*/ 29 30 /*---------------------------------------------------------------------------- 31 * Local headers 32 *----------------------------------------------------------------------------*/ 33 34 #include "cs_base.h" 35 #include "cs_halo.h" 36 #include "cs_internal_coupling.h" 37 #include "cs_mesh.h" 38 #include "cs_mesh_quantities.h" 39 40 /*----------------------------------------------------------------------------*/ 41 42 BEGIN_C_DECLS 43 44 /*============================================================================= 45 * Local Macro definitions 46 *============================================================================*/ 47 48 /*============================================================================ 49 * Type definition 50 *============================================================================*/ 51 52 /*---------------------------------------------------------------------------- 53 * Gradient reconstruction method 54 *----------------------------------------------------------------------------*/ 55 56 typedef enum { 57 58 CS_GRADIENT_GREEN_ITER, /*!< Iterative */ 59 CS_GRADIENT_LSQ, /*!< Least-squares */ 60 CS_GRADIENT_GREEN_LSQ, /*!< Green-Gauss reconstruction with 61 least squares gradient face values */ 62 CS_GRADIENT_GREEN_VTX /*!< Green-Gauss with vertex interpolation */ 63 64 } cs_gradient_type_t; 65 66 /*---------------------------------------------------------------------------- 67 * Gradient limiter mode 68 *----------------------------------------------------------------------------*/ 69 70 typedef enum { 71 72 CS_GRADIENT_LIMIT_NONE = -1, /*!< no limiter */ 73 CS_GRADIENT_LIMIT_CELL = 0, /*!< cell values extrapolated by cell gradient 74 to neighbor cells can not exceed the 75 limiting factor times the actual value */ 76 CS_GRADIENT_LIMIT_FACE = 1, /*!< cell values extrapolated by face gradient 77 (mean of cell gradients) extrapolated to 78 neighbor cells can not exceed the limiting 79 factor times the actual value */ 80 81 } cs_gradient_limit_t; 82 83 /*============================================================================ 84 * Global variables 85 *============================================================================*/ 86 87 /* Short names for gradient types */ 88 89 extern const char *cs_gradient_type_name[]; 90 91 /*============================================================================ 92 * Public function prototypes for Fortran API 93 *============================================================================*/ 94 95 /*---------------------------------------------------------------------------- 96 * Compute the steady balance due to porous modelling for the pressure 97 * gradient 98 *----------------------------------------------------------------------------*/ 99 100 void CS_PROCF (grdpor, GRDPOR) 101 ( 102 const int *const inc /* <-- 0 or 1: increment or not */ 103 ); 104 105 /*---------------------------------------------------------------------------- 106 * Compute cell gradient of vector field. 107 *----------------------------------------------------------------------------*/ 108 109 void CS_PROCF (cgdvec, CGDVEC) 110 ( 111 const int *const f_id, /* <-- field id, or -1 */ 112 const int *const imrgra, /* <-- gradient computation mode */ 113 const int *const inc, /* <-- 0 or 1: increment or not */ 114 const int *const n_r_sweeps,/* <-- >1: with reconstruction */ 115 const int *const iwarnp, /* <-- verbosity level */ 116 const int *const imligp, /* <-- type of clipping */ 117 const cs_real_t *const epsrgp, /* <-- precision for iterative 118 gradient calculation */ 119 const cs_real_t *const climgp, /* <-- clipping coefficient */ 120 const cs_real_3_t coefav[], /* <-- boundary condition term */ 121 const cs_real_33_t coefbv[], /* <-- boundary condition term */ 122 cs_real_3_t pvar[], /* <-- gradient's base variable */ 123 cs_real_33_t grad[] /* <-> gradient of the variable 124 (du_i/dx_j : gradv[][i][j]) */ 125 ); 126 127 /*============================================================================= 128 * Public function prototypes 129 *============================================================================*/ 130 131 /*---------------------------------------------------------------------------- 132 * Initialize gradient computation API. 133 *----------------------------------------------------------------------------*/ 134 135 void 136 cs_gradient_initialize(void); 137 138 /*---------------------------------------------------------------------------- 139 * Finalize gradient computation API. 140 *----------------------------------------------------------------------------*/ 141 142 void 143 cs_gradient_finalize(void); 144 145 /*----------------------------------------------------------------------------*/ 146 /*! 147 * \brief Free saved gradient quantities. 148 * 149 * This is required when the mesh changes, so that the on-demand computation 150 * will be updated. 151 */ 152 /*----------------------------------------------------------------------------*/ 153 154 void 155 cs_gradient_free_quantities(void); 156 157 /*----------------------------------------------------------------------------*/ 158 /*! 159 * \brief Compute cell gradient of scalar field or component of vector or 160 * tensor field. 161 * 162 * \param[in] var_name variable name 163 * \param[in] gradient_type gradient type 164 * \param[in] halo_type halo type 165 * \param[in] inc if 0, solve on increment; 1 otherwise 166 * \param[in] recompute_cocg should COCG FV quantities be recomputed ? 167 * \param[in] n_r_sweeps if > 1, number of reconstruction sweeps 168 * (only used by CS_GRADIENT_GREEN_ITER) 169 * \param[in] tr_dim ignored 170 * \param[in] hyd_p_flag flag for hydrostatic pressure 171 * \param[in] w_stride stride for weighting coefficient 172 * \param[in] verbosity verbosity level 173 * \param[in] clip_mode clipping mode 174 * \param[in] epsilon precision for iterative gradient calculation 175 * \param[in] clip_coeff clipping coefficient 176 * \param[in] f_ext exterior force generating the 177 * hydrostatic pressure 178 * \param[in] bc_coeff_a boundary condition term a 179 * \param[in] bc_coeff_b boundary condition term b 180 * \param[in, out] var gradient's base variable 181 * \param[in, out] c_weight cell variable weight, or NULL 182 * \param[in] cpl associated internal coupling, or NULL 183 * \param[out] grad gradient 184 */ 185 /*----------------------------------------------------------------------------*/ 186 187 void 188 cs_gradient_scalar(const char *var_name, 189 cs_gradient_type_t gradient_type, 190 cs_halo_type_t halo_type, 191 int inc, 192 bool recompute_cocg, 193 int n_r_sweeps, 194 int tr_dim, 195 int hyd_p_flag, 196 int w_stride, 197 int verbosity, 198 cs_gradient_limit_t clip_mode, 199 double epsilon, 200 double clip_coeff, 201 cs_real_3_t f_ext[], 202 const cs_real_t bc_coeff_a[], 203 const cs_real_t bc_coeff_b[], 204 cs_real_t var[restrict], 205 cs_real_t *restrict c_weight, 206 const cs_internal_coupling_t *cpl, 207 cs_real_t grad[restrict][3]); 208 209 /*----------------------------------------------------------------------------*/ 210 /*! 211 * \brief Compute cell gradient of vector field. 212 * 213 * \param[in] var_name variable name 214 * \param[in] gradient_type gradient type 215 * \param[in] halo_type halo type 216 * \param[in] inc if 0, solve on increment; 1 otherwise 217 * \param[in] n_r_sweeps if > 1, number of reconstruction sweeps 218 * (only used by CS_GRADIENT_GREEN_ITER) 219 * \param[in] verbosity verbosity level 220 * \param[in] clip_mode clipping mode 221 * \param[in] epsilon precision for iterative gradient calculation 222 * \param[in] clip_coeff clipping coefficient 223 * \param[in] bc_coeff_a boundary condition term a 224 * \param[in] bc_coeff_b boundary condition term b 225 * \param[in, out] var gradient's base variable 226 * \param[in, out] c_weight cell variable weight, or NULL 227 * \param[in] cpl associated internal coupling, or NULL 228 * \param[out] gradv gradient 229 (\f$ \der{u_i}{x_j} \f$ is gradv[][i][j]) 230 */ 231 /*----------------------------------------------------------------------------*/ 232 233 void 234 cs_gradient_vector(const char *var_name, 235 cs_gradient_type_t gradient_type, 236 cs_halo_type_t halo_type, 237 int inc, 238 int n_r_sweeps, 239 int verbosity, 240 cs_gradient_limit_t clip_mode, 241 double epsilon, 242 double clip_coeff, 243 const cs_real_t bc_coeff_a[][3], 244 const cs_real_t bc_coeff_b[][3][3], 245 cs_real_t var[restrict][3], 246 cs_real_t *restrict c_weight, 247 const cs_internal_coupling_t *cpl, 248 cs_real_t gradv[restrict][3][3]); 249 250 /*----------------------------------------------------------------------------*/ 251 /*! 252 * \brief Compute cell gradient of tensor. 253 * 254 * \param[in] var_name variable name 255 * \param[in] gradient_type gradient type 256 * \param[in] halo_type halo type 257 * \param[in] inc if 0, solve on increment; 1 otherwise 258 * \param[in] n_r_sweeps if > 1, number of reconstruction sweeps 259 * (only used by CS_GRADIENT_GREEN_ITER) 260 * \param[in] verbosity verbosity level 261 * \param[in] clip_mode clipping mode 262 * \param[in] epsilon precision for iterative gradient calculation 263 * \param[in] clip_coeff clipping coefficient 264 * \param[in] bc_coeff_a boundary condition term a 265 * \param[in] bc_coeff_b boundary condition term b 266 * \param[in, out] var gradient's base variable 267 * \param[out] grad gradient 268 (\f$ \der{t_ij}{x_k} \f$ is grad[][ij][k]) 269 */ 270 /*----------------------------------------------------------------------------*/ 271 272 void 273 cs_gradient_tensor(const char *var_name, 274 cs_gradient_type_t gradient_type, 275 cs_halo_type_t halo_type, 276 int inc, 277 int n_r_sweeps, 278 int verbosity, 279 cs_gradient_limit_t clip_mode, 280 double epsilon, 281 double clip_coeff, 282 const cs_real_6_t bc_coeff_a[], 283 const cs_real_66_t bc_coeff_b[], 284 cs_real_6_t *restrict var, 285 cs_real_63_t *restrict grad); 286 287 /*----------------------------------------------------------------------------*/ 288 /*! 289 * \brief Compute cell gradient of scalar field or component of vector or 290 * tensor field. 291 * 292 * This variant of the \ref cs_gradient_scalar function assumes ghost cell 293 * values for input arrays (var and optionally c_weight) 294 * have already been synchronized. 295 * 296 * \param[in] var_name variable name 297 * \param[in] gradient_type gradient type 298 * \param[in] halo_type halo type 299 * \param[in] inc if 0, solve on increment; 1 otherwise 300 * \param[in] recompute_cocg should COCG FV quantities be recomputed ? 301 * \param[in] n_r_sweeps if > 1, number of reconstruction sweeps 302 * (only used by CS_GRADIENT_GREEN_ITER) 303 * \param[in] hyd_p_flag flag for hydrostatic pressure 304 * \param[in] w_stride stride for weighting coefficient 305 * \param[in] verbosity verbosity level 306 * \param[in] clip_mode clipping mode 307 * \param[in] epsilon precision for iterative gradient calculation 308 * \param[in] clip_coeff clipping coefficient 309 * \param[in] f_ext exterior force generating the 310 * hydrostatic pressure 311 * \param[in] bc_coeff_a boundary condition term a 312 * \param[in] bc_coeff_b boundary condition term b 313 * \param[in] var gradient's base variable 314 * \param[in] c_weight cell variable weight, or NULL 315 * \param[in] cpl associated internal coupling, or NULL 316 * \param[out] grad gradient 317 */ 318 /*----------------------------------------------------------------------------*/ 319 320 void 321 cs_gradient_scalar_synced_input(const char *var_name, 322 cs_gradient_type_t gradient_type, 323 cs_halo_type_t halo_type, 324 int inc, 325 bool recompute_cocg, 326 int n_r_sweeps, 327 int hyd_p_flag, 328 int w_stride, 329 int verbosity, 330 cs_gradient_limit_t clip_mode, 331 double epsilon, 332 double clip_coeff, 333 cs_real_t f_ext[][3], 334 const cs_real_t bc_coeff_a[], 335 const cs_real_t bc_coeff_b[], 336 const cs_real_t var[restrict], 337 const cs_real_t c_weight[restrict], 338 const cs_internal_coupling_t *cpl, 339 cs_real_t grad[restrict][3]); 340 341 /*----------------------------------------------------------------------------*/ 342 /*! 343 * \brief Compute cell gradient of vector field. 344 * 345 * This variant of the \ref cs_gradient_vector function assumes ghost cell 346 * values for input arrays (var and optionally c_weight) 347 * have already been synchronized. 348 * 349 * \param[in] var_name variable name 350 * \param[in] gradient_type gradient type 351 * \param[in] halo_type halo type 352 * \param[in] inc if 0, solve on increment; 1 otherwise 353 * \param[in] n_r_sweeps if > 1, number of reconstruction sweeps 354 * (only used by CS_GRADIENT_GREEN_ITER) 355 * \param[in] verbosity verbosity level 356 * \param[in] clip_mode clipping mode 357 * \param[in] epsilon precision for iterative gradient calculation 358 * \param[in] clip_coeff clipping coefficient 359 * \param[in] bc_coeff_a boundary condition term a 360 * \param[in] bc_coeff_b boundary condition term b 361 * \param[in] var gradient's base variable 362 * \param[in] c_weight cell variable weight, or NULL 363 * \param[in] cpl associated internal coupling, or NULL 364 * \param[out] grad gradient 365 (\f$ \der{u_i}{x_j} \f$ is gradv[][i][j]) 366 */ 367 /*----------------------------------------------------------------------------*/ 368 369 void 370 cs_gradient_vector_synced_input(const char *var_name, 371 cs_gradient_type_t gradient_type, 372 cs_halo_type_t halo_type, 373 int inc, 374 int n_r_sweeps, 375 int verbosity, 376 cs_gradient_limit_t clip_mode, 377 double epsilon, 378 double clip_coeff, 379 const cs_real_t bc_coeff_a[][3], 380 const cs_real_t bc_coeff_b[][3][3], 381 const cs_real_t var[restrict][3], 382 const cs_real_t c_weight[restrict], 383 const cs_internal_coupling_t *cpl, 384 cs_real_t grad[restrict][3][3]); 385 386 /*----------------------------------------------------------------------------*/ 387 /*! 388 * \brief Compute cell gradient of tensor. 389 * 390 * This variant of the \ref cs_gradient_tensor function assumes ghost cell 391 * values for input arrays (var and optionally c_weight) 392 * have already been synchronized. 393 * 394 * \param[in] var_name variable name 395 * \param[in] gradient_type gradient type 396 * \param[in] halo_type halo type 397 * \param[in] inc if 0, solve on increment; 1 otherwise 398 * \param[in] n_r_sweeps if > 1, number of reconstruction sweeps 399 * (only used by CS_GRADIENT_GREEN_ITER) 400 * \param[in] verbosity verbosity level 401 * \param[in] clip_mode clipping mode 402 * \param[in] epsilon precision for iterative gradient calculation 403 * \param[in] clip_coeff clipping coefficient 404 * \param[in] bc_coeff_a boundary condition term a 405 * \param[in] bc_coeff_b boundary condition term b 406 * \param[in, out] var gradient's base variable 407 * \param[out] grad gradient 408 (\f$ \der{t_ij}{x_k} \f$ is grad[][ij][k]) 409 */ 410 /*----------------------------------------------------------------------------*/ 411 412 void 413 cs_gradient_tensor_synced_input(const char *var_name, 414 cs_gradient_type_t gradient_type, 415 cs_halo_type_t halo_type, 416 int inc, 417 int n_r_sweeps, 418 int verbosity, 419 cs_gradient_limit_t clip_mode, 420 double epsilon, 421 double clip_coeff, 422 const cs_real_t bc_coeff_a[][6], 423 const cs_real_t bc_coeff_b[][6][6], 424 const cs_real_t var[restrict][6], 425 cs_real_63_t *restrict grad); 426 427 /*----------------------------------------------------------------------------*/ 428 /*! 429 * \brief Compute the gradient of a scalar field at a given cell 430 * using least-squares reconstruction. 431 * 432 * This assumes ghost cell values which might be used are already 433 * synchronized. 434 * 435 * When boundary conditions are provided, both the bc_coeff_a and bc_coeff_b 436 * arrays must be given. If boundary values are known, bc_coeff_a 437 * can point to the boundary values array, and bc_coeff_b set to NULL. 438 * If bc_coeff_a is NULL, bc_coeff_b is ignored. 439 * 440 * \param[in] m pointer to associated mesh structure 441 * \param[in] fvq pointer to associated finite volume quantities 442 * \param[in] c_id cell id 443 * \param[in] halo_type halo type 444 * \param[in] bc_coeff_a boundary condition term a, or NULL 445 * \param[in] bc_coeff_b boundary condition term b, or NULL 446 * \param[in] var gradient's base variable 447 * \param[in] c_weight cell variable weight, or NULL 448 * \param[out] grad gradient 449 */ 450 /*----------------------------------------------------------------------------*/ 451 452 void 453 cs_gradient_scalar_cell(const cs_mesh_t *m, 454 const cs_mesh_quantities_t *fvq, 455 cs_lnum_t c_id, 456 cs_halo_type_t halo_type, 457 const cs_real_t bc_coeff_a[], 458 const cs_real_t bc_coeff_b[], 459 const cs_real_t var[], 460 const cs_real_t c_weight[], 461 cs_real_t grad[3]); 462 463 /*----------------------------------------------------------------------------*/ 464 /*! 465 * \brief Compute the gradient of a vector field at a given cell 466 * using least-squares reconstruction. 467 * 468 * This assumes ghost cell values which might be used are already 469 * synchronized. 470 * 471 * When boundary conditions are provided, both the bc_coeff_a and bc_coeff_b 472 * arrays must be given. If boundary values are known, bc_coeff_a 473 * can point to the boundary values array, and bc_coeff_b set to NULL. 474 * If bc_coeff_a is NULL, bc_coeff_b is ignored. 475 * 476 * \param[in] m pointer to associated mesh structure 477 * \param[in] fvq pointer to associated finite volume quantities 478 * \param[in] c_id cell id 479 * \param[in] halo_type halo type 480 * \param[in] bc_coeff_a boundary condition term a, or NULL 481 * \param[in] bc_coeff_b boundary condition term b, or NULL 482 * \param[in] var gradient's base variable 483 * \param[in] c_weight cell variable weight, or NULL 484 * \param[out] grad gradient 485 */ 486 /*----------------------------------------------------------------------------*/ 487 488 void 489 cs_gradient_vector_cell(const cs_mesh_t *m, 490 const cs_mesh_quantities_t *fvq, 491 cs_lnum_t c_id, 492 cs_halo_type_t halo_type, 493 const cs_real_t bc_coeff_a[][3], 494 const cs_real_t bc_coeff_b[][3][3], 495 const cs_real_t var[restrict][3], 496 const cs_real_t c_weight[], 497 cs_real_t grad[3][3]); 498 499 /*----------------------------------------------------------------------------*/ 500 /*! 501 * \brief Compute the gradient of a tensor field at a given cell 502 * using least-squares reconstruction. 503 * 504 * This assumes ghost cell values which might be used are already 505 * synchronized. 506 * 507 * When boundary conditions are provided, both the bc_coeff_a and bc_coeff_b 508 * arrays must be given. If boundary values are known, bc_coeff_a 509 * can point to the boundary values array, and bc_coeff_b set to NULL. 510 * If bc_coeff_a is NULL, bc_coeff_b is ignored. 511 * 512 * \param[in] m pointer to associated mesh structure 513 * \param[in] fvq pointer to associated finite volume quantities 514 * \param[in] c_id cell id 515 * \param[in] halo_type halo type 516 * \param[in] bc_coeff_a boundary condition term a, or NULL 517 * \param[in] bc_coeff_b boundary condition term b, or NULL 518 * \param[in] var gradient's base variable 519 * \param[in] c_weight cell variable weight, or NULL 520 * \param[out] grad gradient 521 */ 522 /*----------------------------------------------------------------------------*/ 523 524 void 525 cs_gradient_tensor_cell(const cs_mesh_t *m, 526 const cs_mesh_quantities_t *fvq, 527 cs_lnum_t c_id, 528 cs_halo_type_t halo_type, 529 const cs_real_t bc_coeff_a[][6], 530 const cs_real_t bc_coeff_b[][6][6], 531 const cs_real_t var[restrict][6], 532 const cs_real_t c_weight[], 533 cs_real_t grad[6][3]); 534 535 /*---------------------------------------------------------------------------- 536 * Determine gradient type by Fortran "imrgra" value 537 * 538 * parameters: 539 * imrgra <-- Fortran gradient option 540 * gradient_type --> gradient type 541 * halo_type --> halo type 542 *----------------------------------------------------------------------------*/ 543 544 void 545 cs_gradient_type_by_imrgra(int imrgra, 546 cs_gradient_type_t *gradient_type, 547 cs_halo_type_t *halo_type); 548 549 /*----------------------------------------------------------------------------*/ 550 /*! 551 * \brief compute the steady balance due to porous modelling for the pressure 552 * gradient. 553 * 554 * \param[in] inc if 0, solve on increment; 1 otherwise 555 */ 556 /*----------------------------------------------------------------------------*/ 557 558 void 559 cs_gradient_porosity_balance(int inc); 560 561 /*----------------------------------------------------------------------------*/ 562 563 END_C_DECLS 564 565 #endif /* __CS_GRADIENT__ */ 566