/*============================================================================ * * Code_Saturne version 1.3 * ------------------------ * * * This file is part of the Code_Saturne Kernel, element of the * Code_Saturne CFD tool. * * Copyright (C) 1998-2007 EDF S.A., France * * contact: saturne-support@edf.fr * * The Code_Saturne Kernel is free software; you can redistribute it * and/or modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 of * the License, or (at your option) any later version. * * The Code_Saturne Kernel is distributed in the hope that it will be * useful, but WITHOUT ANY WARRANTY; without even the implied warranty * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with the Code_Saturne Kernel; if not, write to the * Free Software Foundation, Inc., * 51 Franklin St, Fifth Floor, * Boston, MA 02110-1301 USA * *============================================================================*/ /*============================================================================ * Définitions, variables globales, et fonctions de base *============================================================================*/ /* includes système */ #include #include #include #include #include #include /* Includes librairie BFT et FVM */ #include #include #include /* Includes librairie */ #include "cs_base.h" #include "cs_maillage.h" #include "cs_parallel.h" #include "cs_voiset.h" #ifdef __cplusplus extern "C" { #endif /* __cplusplus */ /*============================================================================ * Variables statiques locales *============================================================================*/ /*============================================================================ * Prototypes de fonctions privées *============================================================================*/ /*============================================================================ * Prototypes de fonctions FORTRAN appelées *============================================================================*/ /*---------------------------------------------------------------------------- * Calcul des gradients par moindres carrés (standard ou avec voisinage étendu) *----------------------------------------------------------------------------*/ void CS_PROCF (gradmc, GRADMC) ( const cs_int_t *const idbia0, /* --> numéro 1ère case libre dans IA */ const cs_int_t *const idbra0, /* --> numéro 1ère case libre dans RA */ const cs_int_t *const ncelet, /* --> nombre de cellules étendu */ const cs_int_t *const ncel, /* --> nombre de cellules */ const cs_int_t *const nfac, /* --> nombre de faces internes */ const cs_int_t *const nfabor, /* --> nombre de faces de bord */ const cs_int_t *const ncelbr, /* --> nombre de cellules de bord */ const cs_int_t *const nituse, /* --> longueur du tableau ituser[] */ const cs_int_t *const nrtuse, /* --> longueur du tableau rtuser[] */ const cs_int_t *const inc, /* --> 0 ou 1 : incrément ou non */ const cs_int_t *const iccocg, /* --> 1 ou 0 : recalcul de COCG ou non */ const cs_int_t *const nswrgp, /* --> >1 : avec reconstruction */ const cs_int_t *const idimte, /* --> 0, 1, 2 : variable scalaire, vecteur, tenseur */ const cs_int_t *const itenso, /* --> pour la périodicité de rotation */ const cs_int_t *const iphydp, /* --> prise en compte de P hydrostatiqu*/ const cs_int_t *const imrgra, /* --> mode de calcul du gradient */ const cs_int_t *const iwarnp, /* --> niveau d'impression */ const cs_int_t *const nfecra, /* --> unité sortie standard */ const cs_real_t *const epsrgp, /* --> précision pour calcul du gradient itératif */ const cs_real_t *const extrap, /* --> extrapolation des grad. au bord */ const cs_int_t ifacel[], /* --> liste des faces internes */ const cs_int_t ifabor[], /* --> liste des faces de bord */ const cs_int_t icelbr[], /* --> numéro des cellules ayant au moins une face de bord */ const cs_int_t ipcvse[], /* --> pour chaque cellule, position dans ielvse de la première cellule du voisinage étendu */ const cs_int_t ielvse[], /* --> numéro des cellules du voisinage étendu */ cs_int_t ituser[], /* <-> tab. complémentaire utilisateur */ cs_int_t ia[], /* <-> macro-tableau entier */ const cs_real_t volume[], /* --> volumes des cellules */ const cs_real_t surfac[], /* --> surfaces des faces internes */ const cs_real_t surfbo[], /* --> surfaces des faces de bord */ const cs_real_t surfbn[], /* --> norme de surfbo */ const cs_real_t pond[], /* --> pondération géométrique faces int*/ const cs_real_t dist[], /* --> distance entre I' et J' faces int*/ const cs_real_t distbr[], /* --> distance du cdgfbo à I' fac. brd */ const cs_real_t dijpf[], /* --> vecteur I'J' aux faces internes */ const cs_real_t diipb[], /* --> vecteur II' aux faces de bord */ const cs_real_t fextx[], /* --> composantes de la force */ const cs_real_t fexty[], /* extérieure générant la pression */ const cs_real_t fextz[], /* hydrostatique */ const cs_real_t xyzcen[], /* --> c.d.g. des cellules */ const cs_real_t cdgfac[], /* --> c.d.g. des faces internes */ const cs_real_t cdgfbo[], /* --> c.d.g. des faces de bord */ const cs_real_t coefap[], /* --> cond. limites aux faces de bord */ const cs_real_t coefbp[], /* --> cond. limites aux faces de bord */ const cs_real_t pvar[], /* --> variable dont on calcule le grad */ cs_real_t cocgb[], /* <-> contribution à COCG des faces internes des cellules de bord */ cs_real_t cocg[], /* <-> contribution à COCG des faces internes des cellules de bord */ cs_real_t rtuser[], /* <-> tab. complémentaire utilisateur */ cs_real_t dpdx[], /* <-- composante x du gradient */ cs_real_t dpdy[], /* <-- composante y du gradient */ cs_real_t dpdz[], /* <-- composante z du gradient */ cs_real_t bx[], /* - tab. de travail local */ cs_real_t by[], /* - tab. de travail local */ cs_real_t bz[], /* - tab. de travail local */ cs_real_t ra[] /* <-> macro-tableau réel */ ); /*---------------------------------------------------------------------------- * Limitation des gradients (standard ou avec voisinage étendu) *----------------------------------------------------------------------------*/ void CS_PROCF (limgrd, LIMGRD) ( const cs_int_t *const ndim, /* --> dimension de l'espace */ const cs_int_t *const ncelet, /* --> nombre de cellules étendu */ const cs_int_t *const ncel, /* --> nombre de cellules */ const cs_int_t *const nfac, /* --> nombre de faces internes */ const cs_int_t *const imrgra, /* --> mode de calcul du gradient */ const cs_int_t *const imligp, /* --> mode de limitation du gradient */ const cs_int_t *const idimte, /* --> 0, 1, 2 : variable scalaire, vecteur, tenseur */ const cs_int_t *const itenso, /* --> pour la périodicité de rotation */ const cs_int_t *const iwarnp, /* --> niveau d'impression */ const cs_int_t *const nfecra, /* --> unité sortie standard */ const cs_real_t *const climgp, /* --> coefficient de limitation du grad*/ const cs_int_t ifacel[], /* --> liste des faces internes */ const cs_int_t ipcvse[], /* --> pour chaque cellule, position dans ielvse de la première cellule du voisinage étendu */ const cs_int_t ielvse[], /* --> numéro des cellules du voisinage étendu */ const cs_real_t xyzcen[], /* --> c.d.g. des cellules */ const cs_real_t pvar[], /* --> variable dont on calcule le grad */ cs_real_t dpdx[], /* <-> composante x du gradient */ cs_real_t dpdy[], /* <-> composante y du gradient */ cs_real_t dpdz[], /* <-> composante z du gradient */ cs_real_t faclim[], /* - tab. de travail local */ cs_real_t denom[], /* - tab. de travail local */ cs_real_t denum[] /* - tab. de travail local */ ); void CS_PROCF (filtre, FILTRE) ( const cs_int_t *const ncelet, /* --> nombre de cellules étendu */ const cs_int_t *const ncel, /* --> nombre de cellules */ const cs_int_t *const nfac, /* --> nombre de faces internes */ const cs_int_t *const nfabor, /* --> nombre de faces de bord */ const cs_int_t ifacel[], /* --> liste des faces internes */ const cs_int_t ifabor[], /* --> liste des faces de bord */ const cs_int_t ipcvse[], /* --> pour chaque cellule, position dans ielvse de la première cellule du voisinage étendu */ const cs_int_t ielvse[], /* --> numéro des cellules du voisinage étendu */ const cs_real_t volume[], /* --> volumes des cellules */ const cs_real_t pvar[], /* --> variable a filtrer */ cs_real_t varfil[], /* --> resultat de la variable filtree */ cs_real_t w1[], /* - tab. de travail local */ cs_real_t w2[] /* - tab. de travail local */ ); /*============================================================================ * Fonctions publiques *============================================================================*/ /*---------------------------------------------------------------------------- * Definition du voisinage étendu pour le calcul des gradients par moindres * carrés : au premier passage, on choisit les cellules et on modifie la * structure maillage qui stocke le voisinage étendu. * ensuite, on se contente de retourner les pointeurs * * Attention : le voisinage étendu permet de toucher toutes les cellules * qui partagent un noeud avec une cellule donnée, mais il ne * comporte pas les cellules visibles d'une cellule au travers * des faces (elles sont accessibles par IFACEL) * * Interface Fortran : * * SUBROUTINE DEFVSE * ***************** * & ( NDIM , NCEL , NFAC , NNOD , * & IMRGRA , ANOMAX , * & IFACEL , IPNFAC , NODFAC , XYZCEN , SURFAC ) * *----------------------------------------------------------------------------*/ void CS_PROCF (defvse, DEFVSE) ( const cs_int_t *const ndim, /* --> Dimension = 3 */ const cs_int_t *const ncel, /* --> Nombre de cellules */ const cs_int_t *const nfac, /* --> Nombre de faces internes */ const cs_int_t *const nnod, /* --> Nombre de sommets */ const cs_int_t *const imrgra, /* --> Mode de reconstruction */ const cs_real_t *const anomax, /* --> Angle de non orthogonalité en radian au-delà duquel on retient dans le voisinage étendu de chaque cellule voisine les cellules dont un sommet est sur la face */ const cs_int_t *const ifacel, /* --> Connectivité face interne -> cellules*/ const cs_int_t *const ipnfac, /* --> Position du premier noeud d'une face dans nodfac */ const cs_int_t *const nodfac, /* --> Numéro des noeuds des faces internes */ const cs_real_t *const xyzcen, /* --> Coord. du centre des cellules */ const cs_real_t *const surfac /* --> Vecteur surface des faces internes */ ) { cs_maillage_t *maillage = cs_glob_maillage; cs_int_t ind_som; cs_int_t ind_fac; cs_int_t ind_cel; cs_int_t ind_cel_i; cs_int_t ind_cel_j; cs_int_t ind_bcl_cel; cs_int_t ind_bcl_cel_vse; cs_int_t ind_bcl_fac; cs_int_t num_cel; cs_int_t num_cel_vse; cs_int_t nbr_val_connect_som_cel; cs_int_t nbr_val_connect_som_cel_estime; cs_int_t nbr_cel_elimine; cs_int_t ind_val_last; cs_int_t ind_pos_last; cs_int_t icoo; cs_int_t *ielvse; cs_int_t *ipcvse; cs_int_t *pos_connect_som_fac; cs_int_t *val_connect_som_fac; cs_int_t *pos_connect_som_cel; cs_int_t *val_connect_som_cel; cs_real_t vij[CS_DIM_3]; cs_real_t vnf[CS_DIM_3]; cs_real_t normeij; cs_real_t normenf; cs_real_t cos_ij_nf; cs_real_t cos_ij_nf_min; cs_bool_t bool_cel_vu; const cs_real_t *xyzcen_avec_voiset; static cs_int_t cs_loc_premier_appel = 0; enum {X, Y, Z}; #define CS_LOC_MODULE(vect) \ sqrt(vect[X] * vect[X] + vect[Y] * vect[Y] + vect[Z] * vect[Z]) /* Au premier passage, on sélectionne les cellules */ if (cs_loc_premier_appel == 0) { /* Pour la suite, on note que l'on est passé une fois */ cs_loc_premier_appel = 1; /* Si on n'a pas de voisinage étendu, on s'arrête */ ipcvse = maillage->pos_connect_voiset_cel_cel_dom; ielvse = maillage->val_connect_voiset_cel_cel_dom; if ( ipcvse == NULL || ielvse == NULL || maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SANS) { bft_error (__FILE__, __LINE__, 0, _("Lorsque la méthode de calcul des gradients par moindres\n" "carrés sur support étendu est activée (IMRGRA = <%d>), \n" "il est indispensable de générer un voisinage étendu \n" "en utilisant l'Enveloppe en pré-traitement avec \n" "l'option -ve (COMMANDE_VOISET dans le lanceur)."), (int)(*imrgra)); } /* Selon imrgra, on sélectionne le voisinage de diverses façons */ /* Par défaut : on conserve tout ============================= */ /* imrgra = 3 : sélection sur critère de non orthogonalité ======================================================= pour chaque face interne, au-dessus d'un angle de non orthogonalité donné, on retient dans le voisinage étendu des deux cellules voisines de la face toutes les cellules qui s'appuient sur un des sommets de la face */ if ((*imrgra) == 3) { /* On construit d'abord la connectivité sommet -> cellules ------------------------------------------------------ On a déjà face -> sommets, il faut l'inverser on utilisera ensuite face -> cellules */ /* On construit la connectivité sommet -> faces en inversant face -> sommets */ /* Allocation de pos_connect_som_fac */ BFT_MALLOC(pos_connect_som_fac, (*nnod)+1, cs_int_t); /* On compte les faces de chaque sommet, on utilise pos_connect_som_fac*/ for (ind_som = 0 ; ind_som < (*nnod)+1 ; ind_som++) pos_connect_som_fac[ind_som] = 0; for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++) for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ; ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) { ind_som = nodfac[ind_bcl_fac] - 1; pos_connect_som_fac[ind_som+1]++; } /* Construction de pos_connect_som_fac */ pos_connect_som_fac[0] = 1; for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) pos_connect_som_fac[ind_som+1] += pos_connect_som_fac[ind_som]; /* Allocations de val_connect_som_fac */ BFT_MALLOC(val_connect_som_fac, pos_connect_som_fac[(*nnod)]-1, cs_int_t); /* Construction de val_connect_som_fac : pour cela, on utilise pos_connect_som_fac pour repérer la case à remplir dans val_connect_som_fac : on l'incrémente donc à chaque nouvelle valeur ajoutée. A la fin, on a donc décalé les valeurs de pos_connect_som_fac d'une case et il faut les remettre en place. */ for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++) for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ; ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) { ind_som = nodfac[ind_bcl_fac] - 1; val_connect_som_fac[pos_connect_som_fac[ind_som]-1] = ind_fac + 1; pos_connect_som_fac[ind_som]++; } for (ind_som = (*nnod) ; ind_som > 0 ; ind_som--) pos_connect_som_fac[ind_som] = pos_connect_som_fac[ind_som-1]; pos_connect_som_fac[0] = 1; #if 0 for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) { printf(" sommet %d :: ", ind_som+1) ; for (ind_bcl_cel = pos_connect_som_fac[ind_som ] - 1 ; ind_bcl_cel < pos_connect_som_fac[ind_som+1] - 1 ; ind_bcl_cel++) printf(" %d ", val_connect_som_fac[ind_bcl_cel]) ; printf("\n") ; } #endif /* On construit la connectivité sommet -> cellules avec sommet -> faces et face -> cellules */ /* Allocations */ nbr_val_connect_som_cel_estime = 3*(*nnod); BFT_MALLOC(val_connect_som_cel, nbr_val_connect_som_cel_estime, cs_int_t); BFT_MALLOC(pos_connect_som_cel, (*nnod)+1, cs_int_t); /* Initialisations */ nbr_val_connect_som_cel = 0; pos_connect_som_cel[0] = 1; /* Pour tous les sommets */ for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) { /* On balaye les faces */ for (ind_bcl_fac = pos_connect_som_fac[ind_som ] - 1 ; ind_bcl_fac < pos_connect_som_fac[ind_som+1] - 1 ; ind_bcl_fac++) { ind_fac = val_connect_som_fac[ind_bcl_fac] - 1; /* Et on note les cellules (il faut absolument qu'elles ne soient pas en plusieurs exemplaires */ /* Voisin 1 */ num_cel = ifacel[2*ind_fac ]; /* Test si deja vu */ bool_cel_vu = CS_FALSE; ind_bcl_cel = pos_connect_som_cel[ind_som] - 1; while ( (bool_cel_vu == CS_FALSE) && (ind_bcl_cel < nbr_val_connect_som_cel)) { if (num_cel == val_connect_som_cel[ind_bcl_cel]) bool_cel_vu = CS_TRUE; ind_bcl_cel++; } /* Si pas encore vu, on le note */ if (bool_cel_vu == CS_FALSE) { if (nbr_val_connect_som_cel + 1 > nbr_val_connect_som_cel_estime) { nbr_val_connect_som_cel_estime *= 2; BFT_REALLOC(val_connect_som_cel, nbr_val_connect_som_cel_estime, cs_int_t); } val_connect_som_cel[nbr_val_connect_som_cel] = num_cel; nbr_val_connect_som_cel++; } /* Voisin 2 */ num_cel = ifacel[2*ind_fac+1]; /* Test si deja vu */ bool_cel_vu = CS_FALSE; ind_bcl_cel = pos_connect_som_cel[ind_som] - 1; while ( (bool_cel_vu == CS_FALSE) && (ind_bcl_cel < nbr_val_connect_som_cel)) { if (num_cel == val_connect_som_cel[ind_bcl_cel]) bool_cel_vu = CS_TRUE; ind_bcl_cel++; } /* Si pas encore vu, on le note */ if (bool_cel_vu == CS_FALSE) { if (nbr_val_connect_som_cel + 1 > nbr_val_connect_som_cel_estime) { nbr_val_connect_som_cel_estime *=2; BFT_REALLOC(val_connect_som_cel, nbr_val_connect_som_cel_estime, cs_int_t); } val_connect_som_cel[nbr_val_connect_som_cel] = num_cel; nbr_val_connect_som_cel++; } } pos_connect_som_cel[ind_som+1] = nbr_val_connect_som_cel + 1; } /* Réallocation */ BFT_REALLOC(val_connect_som_cel, nbr_val_connect_som_cel, cs_int_t); /* Libération de sommet -> faces */ BFT_FREE(val_connect_som_fac); BFT_FREE(pos_connect_som_fac); #if 0 for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) { printf(" sommet %d :: ", ind_som+1) ; for (ind_bcl_cel = pos_connect_som_cel[ind_som ] - 1 ; ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ; ind_bcl_cel++) printf(" %d ", val_connect_som_cel[ind_bcl_cel]) ; printf("\n") ; } #endif /* On sélectionne enfin les cellules à conserver dans le voisinage --------------------------------------------------------------- On les marque d'abord, on élimine celles qui sont inutiles ensuite */ /* Marquage des cellules à éliminer : en négatif */ assert((*ndim) == CS_DIM_3); /* Coordonnées du centre */ if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE && maillage->nbr_cel_fac_par_voiset > 0) { xyzcen_avec_voiset = maillage->coord_cel_avec_voiset; } else xyzcen_avec_voiset = xyzcen; /* Sélection par les faces internes */ for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++) { /* On calcule le cosinus de l'angle de non orthogonalité des faces internes (angle entre la normale et la droite passant par le centre des cellules voisines */ /* Vecteur IJ */ ind_cel_i = ifacel[2*ind_fac ] - 1; ind_cel_j = ifacel[2*ind_fac+1] - 1; icoo = 0; vij[icoo] = xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo] - xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo]; icoo = 1; vij[icoo] = xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo] - xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo]; icoo = 2; vij[icoo] = xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo] - xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo]; normeij = CS_LOC_MODULE(vij); assert (normeij > 0.); /* Vecteur normal à la face */ icoo = 0; vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo]; icoo = 1; vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo]; icoo = 2; vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo]; normenf = CS_LOC_MODULE(vnf); assert (normenf > 0.); /* Produit scalaire vij.vnf */ cos_ij_nf = (vij[0]*vnf[0]+vij[1]*vnf[1]+vij[2]*vnf[2]) / (normeij*normenf); /* Comparaison à une limite prédéfinie, En dessous de cette limite, c'est non orthogonal et on garde la cellule dans le voisinage étendu des deux voisines de la face (on la marque en négatif, puis on changera tous les signes, pour que les cellules à éliminer soient en négatif */ cos_ij_nf_min = cos((*anomax)); /* Test sur l'angle de non orthogonalité */ if (cos_ij_nf < cos_ij_nf_min) { /* Pour chaque cellule voisine de la face : intersection du groupe des cellules du voisinage étendu que l'on n'a pas encore décidé de garder avec le groupe de cellules possédant un sommet sur la face */ /* Cellule voisine 1 (uniquement si elle est réelle) */ if (ind_cel_i < (*ncel)) { ind_cel = ind_cel_i; /* Voisinage étendu non encore retenu */ for (ind_bcl_cel_vse = ipcvse[ind_cel ] - 1 ; ind_bcl_cel_vse < ipcvse[ind_cel+1] - 1 ; ind_bcl_cel_vse++) { num_cel_vse = ielvse[ind_bcl_cel_vse]; if (num_cel_vse > 0) /* Cellules partageant un sommet avec la face */ for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ; ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) { ind_som = nodfac[ind_bcl_fac] - 1; for (ind_bcl_cel = pos_connect_som_cel[ind_som ] - 1 ; ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ; ind_bcl_cel++) { num_cel = val_connect_som_cel[ind_bcl_cel]; /* Comparaison et sélection */ if (num_cel == num_cel_vse && ielvse[ind_bcl_cel_vse] > 0) ielvse[ind_bcl_cel_vse] = - ielvse[ind_bcl_cel_vse]; } } } } /* Cellule voisine 2 (uniquement si elle est réelle) */ if (ind_cel_j < (*ncel)) { ind_cel = ind_cel_j; /* Voisinage étendu non encore retenu */ for (ind_bcl_cel_vse = ipcvse[ind_cel ] - 1 ; ind_bcl_cel_vse < ipcvse[ind_cel+1] - 1 ; ind_bcl_cel_vse++) { num_cel_vse = ielvse[ind_bcl_cel_vse]; if (num_cel_vse > 0) /* Cellules partageant un sommet avec la face */ for (ind_bcl_fac = ipnfac[ind_fac ] - 1 ; ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) { ind_som = nodfac[ind_bcl_fac] - 1; for (ind_bcl_cel = pos_connect_som_cel[ind_som ] - 1 ; ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ; ind_bcl_cel++) { num_cel = val_connect_som_cel[ind_bcl_cel]; /* Comparaison et sélection */ if (num_cel == num_cel_vse && ielvse[ind_bcl_cel_vse] > 0) ielvse[ind_bcl_cel_vse] = - ielvse[ind_bcl_cel_vse] ; } } } } } } /* Libération de sommet -> cellules */ BFT_FREE(val_connect_som_cel); BFT_FREE(pos_connect_som_cel); /* On change tous les signes de ielvse, pour que les cellules à éliminer soient en négatif */ for (ind_bcl_cel = 0 ; ind_bcl_cel < ipcvse[(*ncel)]-1 ; ind_bcl_cel++) ielvse[ind_bcl_cel] = - ielvse[ind_bcl_cel]; /* On élimine les cellules marquées négativement */ nbr_cel_elimine = 0; ind_val_last = -1; ind_pos_last = 0; for (ind_cel = 0 ; ind_cel < (*ncel) ; ind_cel++) { for (ind_bcl_cel = ind_pos_last ; ind_bcl_cel < ipcvse[ind_cel+1] - 1 ; ind_bcl_cel++) { if (ielvse[ind_bcl_cel] > 0) { ind_val_last++; ielvse[ind_val_last] = ielvse[ind_bcl_cel]; } else { nbr_cel_elimine++; } } ind_pos_last = ipcvse[ind_cel+1] - 1; ipcvse[ind_cel+1] = ipcvse[ind_cel+1] - nbr_cel_elimine; } maillage->nbr_val_connect_voiset_cel_cel_dom = maillage->nbr_val_connect_voiset_cel_cel_dom - nbr_cel_elimine; /* Réallocation avec la dimension exacte */ BFT_REALLOC(maillage->val_connect_voiset_cel_cel_dom, ipcvse[(*ncel)]-1, cs_int_t); bft_printf(_("\n %d cellules conservées dans le voisinage étendu\n"), maillage->nbr_val_connect_voiset_cel_cel_dom); #if 0 for (ind_cel = 0 ; ind_cel < (*ncel) ; ind_cel++) { printf(" cellule %d :: ", ind_cel+1); for (ind_bcl_cel = ipcvse[ind_cel ] - 1 ; ind_bcl_cel < ipcvse[ind_cel+1] - 1 ; ind_bcl_cel++) printf(" %d ", ielvse[ind_bcl_cel]); printf("\n"); } #endif } } } /*---------------------------------------------------------------------------- * Appel de la routine FORTRAN GRADMC pour le calcul des gradients par moindres * carrés. On ajoute aux arguments le voisinage étendu de la structure * maillage. * * Interface Fortran : * * SUBROUTINE CGRDMC * ***************** * * & ( IDBIA0 , IDBRA0 , * & NCELET , NCEL , NFAC , NFABOR , NCELBR , NITUSE , NRTUSE , * & INC , ICCOCG , NSWRGP , IDIMTE , ITENSO , IPHYDP , IMRGRA , * & IWARNP , NFECRA , EPSRGP , EXTRAP , * & IFACEL , IFABOR , ICELBR , ITUSER , IA , * & VOLUME , SURFAC , SURFBO , SURFBN , POND , * & DIST , DISTBR , DIJPF , DIIPB , * & FEXTX , FEXTY , FEXTZ , * & XYZCEN , CDGFAC , CDGFBO , COEFAP , COEFBP , PVAR , * & COCGB , COCG , RTUSER , * & DPDX , DPDY , DPDZ , * & BX , BY , BZ , RA ) * *----------------------------------------------------------------------------*/ void CS_PROCF (cgrdmc, CGRDMC) ( const cs_int_t *const idbia0, /* --> numéro 1ère case libre dans IA */ const cs_int_t *const idbra0, /* --> numéro 1ère case libre dans RA */ const cs_int_t *const ncelet, /* --> nombre de cellules étendu */ const cs_int_t *const ncel, /* --> nombre de cellules */ const cs_int_t *const nfac, /* --> nombre de faces internes */ const cs_int_t *const nfabor, /* --> nombre de faces de bord */ const cs_int_t *const ncelbr, /* --> nombre de cellules de bord */ const cs_int_t *const nituse, /* --> longueur du tableau ituser[] */ const cs_int_t *const nrtuse, /* --> longueur du tableau rtuser[] */ const cs_int_t *const inc, /* --> 0 ou 1 : incrément ou non */ const cs_int_t *const iccocg, /* --> 1 ou 0 : recalcul de COCG ou non */ const cs_int_t *const nswrgp, /* --> >1 : avec reconstruction */ const cs_int_t *const idimte, /* --> 0, 1, 2 : variable scalaire, vecteur, tenseur */ const cs_int_t *const itenso, /* --> pour la périodicité de rotation */ const cs_int_t *const iphydp, /* --> prise en compte de P hydrostatiqu*/ const cs_int_t *const imrgra, /* --> mode de calcul du gradient */ const cs_int_t *const iwarnp, /* --> niveau d'impression */ const cs_int_t *const nfecra, /* --> unité sortie standard */ const cs_real_t *const epsrgp, /* --> précision pour calcul du gradient itératif */ const cs_real_t *const extrap, /* --> extrapolation des grad. au bord */ const cs_int_t ifacel[], /* --> liste des faces internes */ const cs_int_t ifabor[], /* --> liste des faces de bord */ const cs_int_t icelbr[], /* --> numéro des cellules ayant au moins une face de bord */ cs_int_t ituser[], /* <-> tab. complémentaire utilisateur */ cs_int_t ia[], /* <-> macro-tableau entier */ const cs_real_t volume[], /* --> volumes des cellules */ const cs_real_t surfac[], /* --> surfaces des faces internes */ const cs_real_t surfbo[], /* --> surfaces des faces de bord */ const cs_real_t surfbn[], /* --> norme de surfbo */ const cs_real_t pond[], /* --> pondération géométrique faces int*/ const cs_real_t dist[], /* --> distance entre I' et J' faces int*/ const cs_real_t distbr[], /* --> distance du cdgfbo à I' fac. brd */ const cs_real_t dijpf[], /* --> vecteur I'J' aux faces internes */ const cs_real_t diipb[], /* --> vecteur II' aux faces de bord */ const cs_real_t fextx[], /* --> composantes de la force */ const cs_real_t fexty[], /* extérieure générant la pression */ const cs_real_t fextz[], /* hydrostatique */ const cs_real_t xyzcen[], /* --> c.d.g. des cellules */ const cs_real_t cdgfac[], /* --> c.d.g. des faces internes */ const cs_real_t cdgfbo[], /* --> c.d.g. des faces de bord */ const cs_real_t coefap[], /* --> cond. limites aux faces de bord */ const cs_real_t coefbp[], /* --> cond. limites aux faces de bord */ const cs_real_t pvar[], /* --> variable dont on calcule le grad */ cs_real_t cocgb[], /* <-> contribution à COCG des faces internes des cellules de bord */ cs_real_t cocg[], /* <-> contribution à COCG des faces internes des cellules de bord */ cs_real_t rtuser[], /* <-> tab. complémentaire utilisateur */ cs_real_t dpdx[], /* <-- composante x du gradient */ cs_real_t dpdy[], /* <-- composante y du gradient */ cs_real_t dpdz[], /* <-- composante z du gradient */ cs_real_t bx[], /* - tab. de travail local */ cs_real_t by[], /* - tab. de travail local */ cs_real_t bz[], /* - tab. de travail local */ cs_real_t ra[] /* <-> macro-tableau réel */ ) { #if defined(_CS_HAVE_MPI) cs_int_t iel; cs_int_t nceltot; #endif /*_CS_HAVE_MPI */ cs_int_t *ipcvse; cs_int_t *ielvse; const cs_real_t *xyzcen_avec_voiset; const cs_real_t *pvar_avec_voiset; const cs_real_t *fextx_avec_voiset; const cs_real_t *fexty_avec_voiset; const cs_real_t *fextz_avec_voiset; #if defined(_CS_HAVE_MPI) cs_real_t *pvar_voiset; cs_real_t *fextx_voiset; cs_real_t *fexty_voiset; cs_real_t *fextz_voiset; cs_real_t *fextx_avec_voiset_tmp; cs_real_t *fexty_avec_voiset_tmp; cs_real_t *fextz_avec_voiset_tmp; #endif /*_CS_HAVE_MPI */ cs_maillage_t *maillage = cs_glob_maillage; #if defined(_CS_HAVE_MPI) cs_maillage_tmp_t *maillage_tmp = cs_glob_maillage_tmp; #endif /*_CS_HAVE_MPI */ /* Remarque sur les tests faisant intervenir _CS_HAVE_MPI Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie definie que si on _CS_HAVE_MPI (si le parallélisme est possible) Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est correct car il ne peut y avoir un voisinage separé que en parallele etant donne le choix fait dans l'enveloppe (un voisinage separe en non parallele n'apporterait rien car il ne permettrait pas de reduire le volume des echanges). On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t par des tests sur _CS_HAVE_MPI pour permettre la compilation lorsque _CS_HAVE_MPI n'est pas defini Par coherence, on fait de meme avec les variables pvar_voiset fext?_voiset et fext?_avec_voiset_tmp, nceltot, iel */ /* La complexite du traitement dans cette fonction est liée à la prise en compte d'un voisinage étendu séparé du halo classique */ /* On va transmettre au gradient par moindre carrés une variable pvar_avec_voiset qui sera : maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe et qu'il est traité séparément pvar sinon (cas classiques) La variable maillage_tmp->var_cel_avec_voiset devra contenir pvar dans les ncelet premieres cases et les valeurs actualisees du voisinage étendu séparé dans le reste, noté pvar_voiset Il faut de même transmettre la variable allongée correcte pour les centres de gravité des cellules et, lorsque la pression hydrostatique est prise en compte, transmettre également des variables allongées pour la fext. L'utilisation des pointeurs fext?_avec_voiset_tmp au lieu de fext?_avec_voiset directement permet d'éviter un warning et de conserver le caractere const des fext? en argument. */ /* Echanges pour mise à jour du voisinage étendu séparé ---------------------------------------------------- */ #if defined(_CS_HAVE_MPI) if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) { /* Identification des pointeurs sur le voisinage étendu séparé */ if (maillage->nbr_cel_fac_par_voiset > 0) { pvar_voiset = maillage_tmp->var_cel_avec_voiset + (*ncelet); } else { pvar_voiset = NULL; } if (*iphydp != 0) { if (maillage->nbr_cel_fac_par_voiset > 0) { nceltot = (*ncelet)+maillage->nbr_cel_fac_par_voiset; BFT_MALLOC(fextx_avec_voiset_tmp,nceltot,cs_real_t); BFT_MALLOC(fexty_avec_voiset_tmp,nceltot,cs_real_t); BFT_MALLOC(fextz_avec_voiset_tmp,nceltot,cs_real_t); fextx_voiset = fextx_avec_voiset_tmp + (*ncelet); fexty_voiset = fexty_avec_voiset_tmp + (*ncelet); fextz_voiset = fextz_avec_voiset_tmp + (*ncelet); } else { fextx_avec_voiset_tmp = NULL; fexty_avec_voiset_tmp = NULL; fextz_avec_voiset_tmp = NULL; fextx_voiset = NULL; fexty_voiset = NULL; fextz_voiset = NULL; } } /* Echanges parallèles et périodiques Attention, un processeur peut n'avoir rien à recevoir (maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant des éléments à envoyer : on doit donc entrer dans parcve pour tous les processeurs lorsque l'on est en parallèle */ if (maillage->nbr_dom > 1) cs_parallel_parcve(pvar, pvar_voiset); if (maillage->nbr_per > 0) cs_perio_percve (pvar, pvar_voiset); if (*iphydp != 0) { if (maillage->nbr_dom > 1) { cs_parallel_parcve(fextx, fextx_voiset); cs_parallel_parcve(fexty, fexty_voiset); cs_parallel_parcve(fextz, fextz_voiset); } if (maillage->nbr_per > 0) { cs_perio_percve (fextx, fextx_voiset); cs_perio_percve (fexty, fexty_voiset); cs_perio_percve (fextz, fextz_voiset); } } } #endif /*_CS_HAVE_MPI */ /* Appel du gradient par moindres carres ------------------------------------- */ /* Identification du pointeur sur la variable allongée (de longueur ncelet + voisinage étendu séparé éventuel) et copie des valeurs 1-ncelet si nécessaire */ #if defined(_CS_HAVE_MPI) if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE && maillage->nbr_cel_fac_par_voiset > 0) { for (iel = 0 ; iel < (*ncelet) ; iel++) maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel]; pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset; if (*iphydp != 0) { for (iel = 0 ; iel < (*ncelet) ; iel++) fextx_avec_voiset_tmp[iel] = fextx[iel]; fextx_avec_voiset = fextx_avec_voiset_tmp; for (iel = 0 ; iel < (*ncelet) ; iel++) fexty_avec_voiset_tmp[iel] = fexty[iel]; fexty_avec_voiset = fexty_avec_voiset_tmp; for (iel = 0 ; iel < (*ncelet) ; iel++) fextz_avec_voiset_tmp[iel] = fextz[iel]; fextz_avec_voiset = fextz_avec_voiset_tmp; } else { fextx_avec_voiset = fextx; fexty_avec_voiset = fexty; fextz_avec_voiset = fextz; } } else { #endif /*_CS_HAVE_MPI */ pvar_avec_voiset = pvar; fextx_avec_voiset = fextx; fexty_avec_voiset = fexty; fextz_avec_voiset = fextz; #if defined(_CS_HAVE_MPI) } #endif /*_CS_HAVE_MPI */ /* Identification du pointeur sur les coordonnées allongées (de longueur ncelet + voisinage étendu séparé éventuel)*/ if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE && maillage->nbr_cel_fac_par_voiset > 0) { xyzcen_avec_voiset = maillage->coord_cel_avec_voiset; } else xyzcen_avec_voiset = xyzcen; /* Pointeurs sur le voisinage étendu */ ipcvse = maillage->pos_connect_voiset_cel_cel_dom; ielvse = maillage->val_connect_voiset_cel_cel_dom; /* Appel du calcul de gradient */ CS_PROCF(gradmc, GRADMC) (idbia0, idbra0, ncelet, ncel , nfac , nfabor, ncelbr, nituse, nrtuse, inc , iccocg, nswrgp, idimte, itenso, iphydp, imrgra, iwarnp, nfecra, epsrgp, extrap, ifacel, ifabor, icelbr, ipcvse, ielvse, ituser, ia , volume, surfac, surfbo, surfbn, pond , dist , distbr, dijpf , diipb , fextx_avec_voiset, fexty_avec_voiset, fextz_avec_voiset, xyzcen_avec_voiset, cdgfac, cdgfbo, coefap, coefbp, pvar_avec_voiset, cocgb , cocg , rtuser, dpdx , dpdy , dpdz , bx , by , bz , ra ); /* Libération si on a réservé quelque chose */ #if defined(_CS_HAVE_MPI) if ( (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) && (*iphydp != 0) && (maillage->nbr_cel_fac_par_voiset > 0)) { BFT_FREE(fextx_avec_voiset_tmp); BFT_FREE(fexty_avec_voiset_tmp); BFT_FREE(fextz_avec_voiset_tmp); } #endif /*_CS_HAVE_MPI */ } /*---------------------------------------------------------------------------- * Appel de la routine FORTRAN LIMGRD pour la limitation éventuelle des * gradients. On ajoute aux arguments le voisinage étendu de la structure * maillage. * * Interface Fortran : * * SUBROUTINE CLMGRD * ***************** * * & ( NDIM , NCELET , NCEL , NFAC , * & IMRGRA , IMLIGP , IDIMTE , ITENSO , IWARNP , NFECRA , CLIMGP , * & IFACEL , * & XYZCEN , PVAR , * & DPDX , DPDY , DPDZ , * & FACLIM , DENOM , DENUM ) * *----------------------------------------------------------------------------*/ void CS_PROCF (clmgrd, CLMGRD) ( const cs_int_t *const ndim, /* --> dimension de l'espace */ const cs_int_t *const ncelet, /* --> nombre de cellules étendu */ const cs_int_t *const ncel, /* --> nombre de cellules */ const cs_int_t *const nfac, /* --> nombre de faces internes */ const cs_int_t *const imrgra, /* --> mode de calcul du gradient */ const cs_int_t *const imligp, /* --> mode de limitation du gradient */ const cs_int_t *const idimte, /* --> 0, 1, 2 : variable scalaire, vecteur, tenseur */ const cs_int_t *const itenso, /* --> pour la périodicité de rotation */ const cs_int_t *const iwarnp, /* --> niveau d'impression */ const cs_int_t *const nfecra, /* --> unité sortie standard */ const cs_real_t *const climgp, /* --> coefficient de limitation du grad*/ const cs_int_t ifacel[], /* --> liste des faces internes */ const cs_real_t xyzcen[], /* --> c.d.g. des cellules */ const cs_real_t pvar[], /* --> variable dont on calcule le grad */ cs_real_t dpdx[], /* <-> composante x du gradient */ cs_real_t dpdy[], /* <-> composante y du gradient */ cs_real_t dpdz[], /* <-> composante z du gradient */ cs_real_t faclim[], /* - tab. de travail local */ cs_real_t denom[], /* - tab. de travail local */ cs_real_t denum[] /* - tab. de travail local */ ) { #if defined(_CS_HAVE_MPI) cs_int_t iel; cs_int_t nceltot; #endif /*_CS_HAVE_MPI */ cs_int_t *ipcvse; cs_int_t *ielvse; const cs_real_t *xyzcen_avec_voiset; const cs_real_t *pvar_avec_voiset; cs_real_t *dpdx_avec_voiset; cs_real_t *dpdy_avec_voiset; cs_real_t *dpdz_avec_voiset; cs_real_t *denum_avec_voiset; cs_real_t *denom_avec_voiset; #if defined(_CS_HAVE_MPI) cs_real_t *pvar_voiset; cs_real_t *dpdx_voiset; cs_real_t *dpdy_voiset; cs_real_t *dpdz_voiset; cs_real_t *dpdx_avec_voiset_tmp; cs_real_t *dpdy_avec_voiset_tmp; cs_real_t *dpdz_avec_voiset_tmp; cs_real_t *denum_avec_voiset_tmp; cs_real_t *denom_avec_voiset_tmp; #endif /*_CS_HAVE_MPI */ cs_maillage_t *maillage = cs_glob_maillage; #if defined(_CS_HAVE_MPI) cs_maillage_tmp_t *maillage_tmp = cs_glob_maillage_tmp; #endif /*_CS_HAVE_MPI */ /* Remarque sur les tests faisant intervenir _CS_HAVE_MPI Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie definie que si on _CS_HAVE_MPI (si le parallélisme est possible) Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est correct car il ne peut y avoir un voisinage separé que en parallele etant donne le choix fait dans l'enveloppe (un voisinage separe en non parallele n'apporterait rien car il ne permettrait pas de reduire le volume des echanges). On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t par des tests sur _CS_HAVE_MPI pour permettre la compilation lorsque _CS_HAVE_MPI n'est pas defini Par coherence, on fait de meme avec les variables pvar_voiset fext?_voiset et fext?_avec_voiset_tmp, nceltot, iel */ /* La complexite du traitement dans cette fonction est liée à la prise en compte d'un voisinage étendu séparé du halo classique */ /* On va transmettre au gradient par moindre carrés une variable pvar_avec_voiset qui sera : maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe et qu'il est traité séparément pvar sinon (cas classiques) La variable maillage_tmp->var_cel_avec_voiset devra contenir pvar dans les ncelet premieres cases et les valeurs actualisees du voisinage étendu séparé dans le reste, noté pvar_voiset Il faut de même transmettre la variable allongée correcte pour les centres de gravité des cellules et, pour le mode de limitation imligp = 1, transmettre également une variable allongée pour le gradient et les tableaux de travail denom et denum. */ /* Echanges pour mise à jour du voisinage étendu séparé ---------------------------------------------------- */ #if defined(_CS_HAVE_MPI) if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) { /* Identification des pointeurs sur le voisinage étendu séparé */ if (maillage->nbr_cel_fac_par_voiset > 0) { pvar_voiset = maillage_tmp->var_cel_avec_voiset + maillage->nbr_cel_etendu; } else { pvar_voiset = NULL; } /* Pour le gradient et deux tableaux de travail denom et denum. L'utilisation d'un pointeur _tmp pour le gradient (au lieu d'un pointeur dpdx_avec_voiset directement) provient du fait qu'on a recopié la fonction CGRDMC. Ici ce n'était pas indispensable, dpdx en entrée n'étant pas const. Pour denom et denum, pas nécessaire d'initialiser les tableaux on n'a donc pas besoin de denom_voiset et denum_voiset */ if (*imligp == 1) { if (maillage->nbr_cel_fac_par_voiset > 0) { nceltot = maillage->nbr_cel_etendu + maillage->nbr_cel_fac_par_voiset; BFT_MALLOC(dpdx_avec_voiset_tmp, nceltot, cs_real_t); BFT_MALLOC(dpdy_avec_voiset_tmp, nceltot, cs_real_t); BFT_MALLOC(dpdz_avec_voiset_tmp, nceltot, cs_real_t); BFT_MALLOC(denum_avec_voiset_tmp, nceltot, cs_real_t); BFT_MALLOC(denom_avec_voiset_tmp, nceltot, cs_real_t); dpdx_voiset = dpdx_avec_voiset_tmp + maillage->nbr_cel_etendu; dpdy_voiset = dpdy_avec_voiset_tmp + maillage->nbr_cel_etendu; dpdz_voiset = dpdz_avec_voiset_tmp + maillage->nbr_cel_etendu; } else { dpdx_avec_voiset_tmp = NULL; dpdy_avec_voiset_tmp = NULL; dpdz_avec_voiset_tmp = NULL; denum_avec_voiset_tmp = NULL; denom_avec_voiset_tmp = NULL; dpdx_voiset = NULL; dpdy_voiset = NULL; dpdz_voiset = NULL; } } /* Echanges parallèles et périodiques Attention, un processeur peut n'avoir rien à recevoir (maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant des éléments à envoyer : on doit donc entrer dans parcve pour tous les processeurs lorsque l'on est en parallèle */ if (maillage->nbr_dom > 1) cs_parallel_parcve(pvar, pvar_voiset); if (maillage->nbr_per > 0) cs_perio_percve(pvar, pvar_voiset); /* Echange pour les gradients (Pour les tableaux de travail, échange inutile) */ if (*imligp == 1) { if (maillage->nbr_dom > 1) { cs_parallel_parcve(dpdx, dpdx_voiset); cs_parallel_parcve(dpdy, dpdy_voiset); cs_parallel_parcve(dpdz, dpdz_voiset); } if (maillage->nbr_per > 0) { cs_perio_percve(dpdx, dpdx_voiset); cs_perio_percve(dpdy, dpdy_voiset); cs_perio_percve(dpdz, dpdz_voiset); } } } #endif /*_CS_HAVE_MPI */ /* Appel de la limitation du gradient ---------------------------------- */ /* Identification du pointeur sur la variable allongée (de longueur ncelet + voisinage étendu séparé éventuel) et copie des valeurs 1-ncelet si nécessaire */ #if defined(_CS_HAVE_MPI) if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE && maillage->nbr_cel_fac_par_voiset > 0) { for (iel = 0 ; iel < (*ncelet) ; iel++) maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel]; pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset; if (*imligp == 1) { for (iel = 0 ; iel < (*ncelet) ; iel++) dpdx_avec_voiset_tmp[iel] = dpdx[iel]; dpdx_avec_voiset = dpdx_avec_voiset_tmp; for (iel = 0 ; iel < (*ncelet) ; iel++) dpdy_avec_voiset_tmp[iel] = dpdy[iel]; dpdy_avec_voiset = dpdy_avec_voiset_tmp; for (iel = 0 ; iel < (*ncelet) ; iel++) dpdz_avec_voiset_tmp[iel] = dpdz[iel]; dpdz_avec_voiset = dpdz_avec_voiset_tmp; denum_avec_voiset = denum_avec_voiset_tmp; denom_avec_voiset = denom_avec_voiset_tmp; } else { dpdx_avec_voiset = dpdx; dpdy_avec_voiset = dpdy; dpdz_avec_voiset = dpdz; denum_avec_voiset = denum; denom_avec_voiset = denom; } } else { #endif /*_CS_HAVE_MPI */ pvar_avec_voiset = pvar; dpdx_avec_voiset = dpdx; dpdy_avec_voiset = dpdy; dpdz_avec_voiset = dpdz; denum_avec_voiset = denum; denom_avec_voiset = denom; #if defined(_CS_HAVE_MPI) } #endif /*_CS_HAVE_MPI */ /* Identification du pointeur sur les coordonnées allongées (de longueur ncelet + voisinage étendu séparé éventuel)*/ if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE && maillage->nbr_cel_fac_par_voiset > 0) { xyzcen_avec_voiset = maillage->coord_cel_avec_voiset; } else xyzcen_avec_voiset = xyzcen; /* Pointeurs sur le voisinage étendu */ ipcvse = maillage->pos_connect_voiset_cel_cel_dom; ielvse = maillage->val_connect_voiset_cel_cel_dom; /* Appel de la limitation du gradient */ CS_PROCF(limgrd, LIMGRD) (ndim, ncelet, ncel, nfac, imrgra, imligp, idimte, itenso, iwarnp, nfecra, climgp, ifacel, ipcvse, ielvse, xyzcen_avec_voiset, pvar_avec_voiset, dpdx_avec_voiset, dpdy_avec_voiset, dpdz_avec_voiset, faclim, denom_avec_voiset, denum_avec_voiset); /* Transfert du gradient limité si on a utilisé un tableau différent puis libération des tableaux locaux (voisinage étendu séparé, non vide et avec imligp == 1) */ #if defined(_CS_HAVE_MPI) if ( (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) && (maillage->nbr_cel_fac_par_voiset > 0) && (*imligp == 1)) { for (iel = 0 ; iel < (*ncelet) ; iel++) dpdx[iel] = dpdx_avec_voiset[iel]; for (iel = 0 ; iel < (*ncelet) ; iel++) dpdy[iel] = dpdy_avec_voiset[iel]; for (iel = 0 ; iel < (*ncelet) ; iel++) dpdz[iel] = dpdz_avec_voiset[iel]; BFT_FREE(dpdx_avec_voiset_tmp); BFT_FREE(dpdy_avec_voiset_tmp); BFT_FREE(dpdz_avec_voiset_tmp); BFT_FREE(denum_avec_voiset_tmp); BFT_FREE(denom_avec_voiset_tmp); } #endif /*_CS_HAVE_MPI */ } /*---------------------------------------------------------------------------- * Appel de la routine FORTRAN CFILTR pour le calcul des filtres avec voisinages * etendus pour le modele dynamique. * * Interface Fortran : * * SUBROUTINE CFILTR * ***************** * * & ( NCELET , NCEL , NFAC , NFABOR , * & IFACEL , IFABOR , * & VOLUME , * & PVAR , VARFIL , * & W1 , W2 ) * *----------------------------------------------------------------------------*/ void CS_PROCF (cfiltr, CFILTR) ( const cs_int_t *const ncelet, /* --> nombre de cellules étendu */ const cs_int_t *const ncel, /* --> nombre de cellules */ const cs_int_t *const nfac, /* --> nombre de faces internes */ const cs_int_t *const nfabor, /* --> nombre de faces de bord */ const cs_int_t ifacel[], /* --> liste des faces internes */ const cs_int_t ifabor[], /* --> liste des faces de bord */ const cs_real_t volume[], /* --> volumes des cellules */ const cs_real_t pvar[], /* --> variable a filter */ cs_real_t varfil[], /* --> resultat de la variable filtree */ cs_real_t w1[], /* - tab. de travail local */ cs_real_t w2[] /* - tab. de travail local */ ) { #if defined(_CS_HAVE_MPI) cs_int_t iel; cs_int_t nceltot; #endif /*_CS_HAVE_MPI */ cs_int_t *ipcvse; cs_int_t *ielvse; const cs_real_t *pvar_avec_voiset; const cs_real_t *volume_avec_voiset; #if defined(_CS_HAVE_MPI) cs_real_t *pvar_voiset; cs_real_t *volume_voiset; cs_real_t *volume_avec_voiset_tmp; #endif /*_CS_HAVE_MPI */ cs_maillage_t *maillage = cs_glob_maillage; #if defined(_CS_HAVE_MPI) cs_maillage_tmp_t *maillage_tmp = cs_glob_maillage_tmp; #endif /*_CS_HAVE_MPI */ /* Remarque sur les tests faisant intervenir _CS_HAVE_MPI Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie definie que si on _CS_HAVE_MPI (si le parallélisme est possible) Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est correct car il ne peut y avoir un voisinage separé que en parallele etant donne le choix fait dans l'enveloppe (un voisinage separe en non parallele n'apporterait rien car il ne permettrait pas de reduire le volume des echanges). On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t par des tests sur _CS_HAVE_MPI pour permettre la compilation lorsque _CS_HAVE_MPI n'est pas defini Par coherence, on fait de meme avec les variables pvar_voiset nceltot, iel */ /* La complexite du traitement dans cette fonction est liée à la prise en compte d'un voisinage étendu séparé du halo classique */ /* On va transmettre au gradient par moindre carrés une variable pvar_avec_voiset qui sera : maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe et qu'il est traité séparément pvar sinon (cas classiques) La variable maillage_tmp->var_cel_avec_voiset devra contenir pvar dans les ncelet premieres cases et les valeurs actualisees du voisinage étendu séparé dans le reste, noté pvar_voiset Il faut de même transmettre la variable allongée correcte pour les centres de gravité des cellules et, lorsque la pression hydrostatique est prise en compte, transmettre également des variables allongées pour la fext. L'utilisation des pointeurs fext?_avec_voiset_tmp au lieu de fext?_avec_voiset directement permet d'éviter un warning et de conserver le caractere const des fext? en argument. */ /* Echanges pour mise à jour du voisinage étendu séparé ---------------------------------------------------- */ #if defined(_CS_HAVE_MPI) if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) { /* Identification des pointeurs sur le voisinage étendu séparé */ if (maillage->nbr_cel_fac_par_voiset > 0) { pvar_voiset = maillage_tmp->var_cel_avec_voiset + (*ncelet); } else { pvar_voiset = NULL; } if (maillage->nbr_cel_fac_par_voiset > 0) { nceltot = (*ncelet)+maillage->nbr_cel_fac_par_voiset; BFT_MALLOC(volume_avec_voiset_tmp,nceltot,cs_real_t); volume_voiset = volume_avec_voiset_tmp + (*ncelet); } else { volume_avec_voiset_tmp = NULL; volume_voiset = NULL; } /* Echanges parallèles et périodiques Attention, un processeur peut n'avoir rien à recevoir (maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant des éléments à envoyer : on doit donc entrer dans parcve pour tous les processeurs lorsque l'on est en parallèle */ if (maillage->nbr_dom > 1) cs_parallel_parcve(pvar, pvar_voiset); cs_parallel_parcve(volume, volume_voiset); if (maillage->nbr_per > 0) cs_perio_percve (pvar, pvar_voiset); cs_perio_percve (volume, volume_voiset); } #endif /*_CS_HAVE_MPI */ #if defined(_CS_HAVE_MPI) if ( maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE && maillage->nbr_cel_fac_par_voiset > 0) { for (iel = 0 ; iel < (*ncelet) ; iel++) maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel]; pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset; for (iel = 0 ; iel < (*ncelet) ; iel++) volume_avec_voiset_tmp[iel] = volume[iel]; volume_avec_voiset = volume_avec_voiset_tmp; } else{ #endif /*_CS_HAVE_MPI */ pvar_avec_voiset = pvar; volume_avec_voiset = volume; #if defined(_CS_HAVE_MPI) } #endif /*_CS_HAVE_MPI */ /* Pointeurs sur le voisinage étendu */ ipcvse = maillage->pos_connect_voiset_cel_cel_dom; ielvse = maillage->val_connect_voiset_cel_cel_dom; /* Identification du pointeur sur la variable allongée (de longueur ncelet + voisinage étendu séparé éventuel) et copie des valeurs 1-ncelet si nécessaire */ CS_PROCF(filtre, FILTRE) (ncelet, ncel , nfac , nfabor, ifacel, ifabor, ipcvse, ielvse, volume_avec_voiset, pvar_avec_voiset, varfil, w1 , w2 ); /* Libération si on a réservé quelque chose */ #if defined(_CS_HAVE_MPI) if ( (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) && (maillage->nbr_cel_fac_par_voiset > 0)) { BFT_FREE(volume_avec_voiset_tmp); } #endif /*_CS_HAVE_MPI */ } /*----------------------------------------------------------------------------*/ #ifdef __cplusplus } #endif /* __cplusplus */