Op_Conv_EF_Stab_PolyMAC_CDO_Face.cpp

/****************************************************************************
* Copyright (c) 2026, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#include <Op_Conv_EF_Stab_PolyMAC_CDO_Face.h>
#include <Discretisation_base.h>
#include <Dirichlet_homogene.h>
#include <Champ_Face_PolyMAC_CDO.h>
#include <Domaine_Cl_PolyMAC_family.h>
#include <Schema_Temps_base.h>
#include <Domaine_PolyMAC_CDO.h>
#include <Probleme_base.h>
#include <Matrix_tools.h>
#include <Array_tools.h>
#include <TRUSTLists.h>
#include <Dirichlet.h>
#include <Param.h>
#include <cmath>
 
Implemente_instanciable( Op_Conv_EF_Stab_PolyMAC_CDO_Face, "Op_Conv_EF_Stab_PolyMAC_CDO_Face", Op_Conv_PolyMAC_CDO_base );
Implemente_instanciable( Op_Conv_Amont_PolyMAC_CDO_Face, "Op_Conv_Amont_PolyMAC_CDO_Face", Op_Conv_EF_Stab_PolyMAC_CDO_Face );
Implemente_instanciable( Op_Conv_Centre_PolyMAC_CDO_Face, "Op_Conv_Centre_PolyMAC_CDO_Face", Op_Conv_EF_Stab_PolyMAC_CDO_Face );
 
// XD Op_Conv_EF_Stab_PolyMAC_CDO_Face interprete Op_Conv_EF_Stab_PolyMAC_CDO_Face BRACE Class
// XD_CONT Op_Conv_EF_Stab_PolyMAC_CDO_Face_PolyMAC_CDO
 
Sortie& Op_Conv_EF_Stab_PolyMAC_CDO_Face::printOn(Sortie& os) const { return Op_Conv_PolyMAC_CDO_base::printOn(os); }
Sortie& Op_Conv_Amont_PolyMAC_CDO_Face::printOn(Sortie& os) const { return Op_Conv_PolyMAC_CDO_base::printOn(os); }
Sortie& Op_Conv_Centre_PolyMAC_CDO_Face::printOn(Sortie& os) const { return Op_Conv_PolyMAC_CDO_base::printOn(os); }
 
Entree& Op_Conv_EF_Stab_PolyMAC_CDO_Face::readOn(Entree& is)
{
  Op_Conv_PolyMAC_CDO_base::readOn(is);
  Param param(que_suis_je());
  param.ajouter("alpha", &alpha_);            // XD_ADD_P double
  // XD_CONT parameter adjusting stabilization from 0 (centered scheme) to 1 (upwind scheme)
  param.lire_avec_accolades_depuis(is);
  return is;
}
 
Entree& Op_Conv_Amont_PolyMAC_CDO_Face::readOn(Entree& is)
{
  alpha_ = 1.0;
  return Op_Conv_PolyMAC_CDO_base::readOn(is);
}
 
Entree& Op_Conv_Centre_PolyMAC_CDO_Face::readOn(Entree& is)
{
  alpha_ = 0.0;
  return Op_Conv_PolyMAC_CDO_base::readOn(is);
}
 
double Op_Conv_EF_Stab_PolyMAC_CDO_Face::calculer_dt_stab() const
{
 
  double dt = 1e10;
  const Domaine_Poly_base& domaine = le_dom_poly_.valeur();
  const DoubleVect& fs = domaine.face_surfaces(), &pf = equation().milieu().porosite_face(), &ve = domaine.volumes(), &pe = equation().milieu().porosite_elem();
  const DoubleTab& vit = vitesse_->valeurs();
  const IntTab& e_f = domaine.elem_faces(), &f_e = domaine.face_voisins();
  const int N = vit.line_size();
  DoubleTrav flux(N); //sum of fluxes pf * |f| * vf, minimum volume of elements/faces affected by this flux
 
  for (int e = 0; e < domaine.nb_elem(); e++)
    {
      // Compute the effective volume of the element
      const double vol = pe(e) * ve(e);
      flux = 0.;
 
      // Loop over faces associated with the element
      for (int i = 0; i < e_f.dimension(1); i++)
        {
          int f = e_f(e, i);
          if (f < 0) continue; // non-existent face
 
          for (int n = 0; n < N; n++)
            {
              // Add incoming flux for component n: only incoming fluxes count
              double flux_f = pf(f) * fs(f) * std::max((e == f_e(f, 1) ? 1 : -1) * vit(f, n), 0.);
              flux(n) += flux_f;
            }
        }
 
      // Compute the time step for each component n
      for (int n = 0; n < N; n++)
        if (std::abs(flux(n)) > 1e-12)
          dt = std::min(dt, vol / flux(n));
    }
 
  return Process::mp_min(dt);
}
 
void Op_Conv_EF_Stab_PolyMAC_CDO_Face::completer()
{
  Op_Conv_PolyMAC_CDO_base::completer();
 
  /* just in case... */
  const Domaine_PolyMAC_CDO& domaine = le_dom_poly_.valeur();
  if (domaine.domaine().nb_joints() && domaine.domaine().joint(0).epaisseur() < 2)
    {
      Cerr << "Op_Conv_EF_Stab_PolyMAC_CDO_Face : largeur de joint insuffisante (minimum 2)!" << finl;
      Process::exit();
    }
  porosite_f.ref(mon_equation->milieu().porosite_face());
  porosite_e.ref(mon_equation->milieu().porosite_elem());
}
 
void Op_Conv_EF_Stab_PolyMAC_CDO_Face::dimensionner(Matrice_Morse& mat) const
{
  if (has_interface_blocs())
    {
      Operateur_base::dimensionner(mat);
      return;
    }
 
  const Domaine_PolyMAC_CDO& domaine = le_dom_poly_.valeur();
  const Champ_Face_PolyMAC_CDO& ch = ref_cast(Champ_Face_PolyMAC_CDO, equation().inconnue());
  const IntTab& e_f = domaine.elem_faces(), &f_e = domaine.face_voisins(), &equiv = domaine.equiv();
  int i, j, k, l, e, eb, f, fb, fd, fc;
 
  ch.fcl();
 
  Stencil stencil(0, 2);
 
  for (f = 0; f < domaine.nb_faces_tot(); f++)
    if (f_e(f, 0) >= 0 && (f_e(f, 1) >= 0 ||  ch.fcl()(f, 0) == 1 ||  ch.fcl()(f, 0) == 3))
      {
        for (i = 0; i < 2; i++)
          if ((e = f_e(f, i)) >= 0)
            {
              for (k = 0; k < e_f.dimension(1) && (fb = e_f(e, k)) >= 0; k++)
                if (fb < domaine.nb_faces() &&  ch.fcl()(fb, 0) < 2) //"faces" part
                  {
                    if ((fc = equiv(f, i, k)) >= 0 || f_e(f, 1) < 0)
                      for (j = 0; j < 2; j++) //equivalence: face fd -> face fb
                        {
                          fd = (j == i ? fb : fc); //element/face sources
                          if (fd >= 0) stencil.append_line(fb,fd);
                        }
                    else for (j = 0; j < 2; j++)  //no equivalence: n_f * element operator
                        {
                          for (eb = f_e(f, j), l = 0; l < e_f.dimension(1) && (fc = e_f(eb, l)) >= 0; l++)
                            stencil.append_line(fb,fc);
                        }
                  }
            }
      }
 
  // Sort and remove duplicates
  tableau_trier_retirer_doublons(stencil);
 
  // Allocate the matrix
  int taille = domaine.nb_faces_tot() + (dimension < 3 ? domaine.domaine().nb_som_tot() : domaine.domaine().nb_aretes_tot());
  Matrix_tools::allocate_morse_matrix(taille, taille, stencil, mat);
}
 
// adds the convection contribution to the right-hand side resu
// returns resu
 
DoubleTab& Op_Conv_EF_Stab_PolyMAC_CDO_Face::ajouter(const DoubleTab& inco, DoubleTab& secmem) const
{
  if (has_interface_blocs())
    return Operateur_base::ajouter(inco, secmem);
 
  const Domaine_PolyMAC_CDO& domaine = le_dom_poly_.valeur();
  const Champ_Face_PolyMAC_CDO& ch = ref_cast(Champ_Face_PolyMAC_CDO, equation().inconnue());
  const Conds_lim& cls = la_zcl_poly_->les_conditions_limites();
 
  const IntTab& f_e = domaine.face_voisins(), &e_f = domaine.elem_faces(), &equiv = domaine.equiv();
  const DoubleTab& xp = domaine.xp(), &xv = domaine.xv(), &vfd = domaine.volumes_entrelaces_dir(), &vit = vitesse_->valeurs();
  const DoubleVect& fs = domaine.face_surfaces(), &ve = domaine.volumes(),  &pf = porosite_f, &pe = porosite_e;
  const DoubleTab& nf = domaine.face_normales();
 
  int i, j, k, l, e, eb, f, fb, fc, fd, m, n, N = inco.line_size(), d, D = dimension, comp = !incompressible_;
  double mult;
 
  assert(N == 1);
  DoubleTrav dfac(2, N, N);
  double sum_dfac[2];
  const int nb_faces_tot = domaine.nb_faces_tot();
  const int nb_faces = domaine.nb_faces();
  int nb_face_elem = e_f.dimension(1);
  for (f = 0; f < nb_faces_tot; f++)
    {
      const int elem0 = f_e(f, 0);
      const int elem1 = f_e(f, 1);
      if (elem0 >= 0 && (elem1 >= 0 || ch.fcl()(f, 0) == 1 || ch.fcl()(f, 0) == 3))
        {
          //mass: diagonal + added mass if correlation
          const double inv_masse = 1.0 / (std::fabs(vit[f]) > 1e-10 ? inco(f) / vit[f] : 1.0);
          for (i = 0, dfac = 0; i < 2; i++)
            {
              //contribution a dfac
              for (eb = f_e(f, i), n = 0; n < N; n++)
                {
                  for (m = 0; m < N; m++)
                    dfac(ch.fcl()(f, 0) == 1 ? 0 : i, n, m) += fs(f) * inco[f] * pe(eb >= 0 ? eb : elem0)
                                                               * (1. + (vit[f] * (i ? -1 : 1) >= 0 ? 1. : vit[f] ? -1.
                                                                        : 0.) *
                                                                  alpha_) * 0.5;
                }
              sum_dfac[i] = 0;
              for (n = 0; n < N; n++)
                for (m = 0; m < N; m++)
                  sum_dfac[i] += dfac(i, n, m);
            }
          for (i = 0; i < 2; i++)
            if ((e = f_e(f, i)) >= 0)
              {
                double inv_ve = 1.0 / ve(e);
                for (k = 0; k < nb_face_elem && (fb = e_f(e, k)) >= 0; k++)
                  if (fb < nb_faces && ch.fcl()(fb, 0) < 2) //"faces" part
                    {
                      if ((fc = equiv(f, i, k)) >= 0 || elem1 < 0)
                        for (j = 0; j < 2; j++) //equivalence: face fd -> face fb
                          {
                            eb = f_e(f, j), fd = (j == i ? fb : fc); //element/face sources
                            mult = (fd < 0 || domaine.dot(&nf(fb, 0), &nf(fd, 0)) > 0 ? 1 : -1) *
                                   (fd >= 0 ? pf(fd) / pe(eb) : 1); //multiplier to convert from vf to ve
                            for (n = 0; n < N; n++)
                              for (m = 0; m < N; m++)
                                if (dfac(j, n, m))
                                  {
                                    double fac = (e == elem0 ? 1 : -1) * vfd(fb, e != f_e(fb, 0)) *
                                                 dfac(j, n, m) * inv_ve;
                                    if (fd >= 0)
                                      secmem[fb] -= fac * mult * vit[fd]; //other computed face
                                    else
                                      {
                                        const Cond_lim_base& my_cl = cls[ch.fcl()(f, 1)].valeur();
                                        if (sub_type(Dirichlet, my_cl)) // otherwise: wall -> no contribution
                                          for (d = 0; d < D; d++)  //Dirichlet BC
                                            secmem[fb] -= fac * nf(fb, d) / fs(fb) *
                                                          ref_cast(Dirichlet, my_cl).val_imp(
                                                            ch.fcl()(f, 2), N * d + m) * inv_masse;
                                      }
                                    if (comp) secmem[fb] += fac * vit[fb]; //v div(alpha rho v) part
                                  }
                          }
                      else
                        {
                          for (j = 0; j < 2; j++)  //no equivalence: n_f * element operator
                            {
                              for (eb = f_e(f, j), l = 0; l < nb_face_elem && (fc = e_f(eb, l)) >= 0; l++)
                                {
                                  double num =
                                    (e == f_e(fb, 0) ? 1 : -1) * (e == elem0 ? 1 : -1) * fs(fc) * fs(fb) *
                                    domaine.dot(&xv(fc, 0), &xv(fb, 0), &xp(eb, 0), &xp(e, 0)) *
                                    (eb == f_e(fc, 0) ? 1 : -1);
                                  double den = ve(eb) * ve(e);
                                  if (std::fabs(num) > 1e-9 * den)
                                    {
                                      double num_den = num / den;
                                      secmem[fb] -= sum_dfac[j] * num_den * vit[fc];
                                    }
                                }
                              if (comp)
                                for (l = 0; l < nb_face_elem && (fc = e_f(e, l)) >= 0; l++)
                                  {
                                    double num = (e == f_e(fb, 0) ? 1 : -1) * (e == elem0 ? 1 : -1) * fs(fc) *
                                                 fs(fb) *
                                                 domaine.dot(&xv(fc, 0), &xv(fb, 0), &xp(e, 0), &xp(e, 0)) *
                                                 (e == f_e(fc, 0) ? 1 : -1);
                                    double den = ve(e) * ve(e);
                                    if (std::fabs(num) > 1e-9 * den)
                                      {
                                        double num_den = num / den;
                                        secmem[fb] += sum_dfac[j] * num_den * vit[fc];
                                      }
                                  }
                            }
                        }
                    }
              }
        }
    }
 
  return secmem;
}
 
/*! @brief Assemble the matrix.
 *
 */
inline void Op_Conv_EF_Stab_PolyMAC_CDO_Face::contribuer_a_avec(const DoubleTab& inco, Matrice_Morse& matrice) const
{
  if (has_interface_blocs())
    {
      Operateur_base::contribuer_a_avec(inco, matrice);
      return;
    }
 
  const Domaine_PolyMAC_CDO& domaine = le_dom_poly_.valeur();
  const Champ_Face_PolyMAC_CDO& ch = ref_cast(Champ_Face_PolyMAC_CDO, equation().inconnue());
  const IntTab& f_e = domaine.face_voisins(), &e_f = domaine.elem_faces(), &equiv = domaine.equiv();
  const DoubleTab& xp = domaine.xp(), &xv = domaine.xv(), &vfd = domaine.volumes_entrelaces_dir(), &vit = vitesse_->valeurs();
  const DoubleVect& fs = domaine.face_surfaces(), &ve = domaine.volumes(), &pf = porosite_f, &pe = porosite_e;
  const DoubleTab& nf = domaine.face_normales();
 
  int i, j, k, l, e, eb, f, fb, fc, fd, m, n, N = inco.line_size(), comp = !incompressible_;
  double mult;
 
  assert(N == 1);
  DoubleTrav dfac(2, N, N), masse(N, N);
  for (f = 0; f < domaine.nb_faces_tot(); f++)
    if (f_e(f, 0) >= 0 && (f_e(f, 1) >= 0 ||  ch.fcl()(f, 0) == 1 ||  ch.fcl()(f, 0) == 3))
      {
        for (i = 0, dfac = 0; i < 2; i++)
          {
            //mass: diagonal + added mass if correlation
            masse(0, 0) = std::fabs(vit[f]) > 1e-10 ? inco(f) / vit[f] : 1.0;
            //contribution a dfac
            for (eb = f_e(f, i), n = 0; n < N; n++)
              for (m = 0; m < N; m++)
                dfac( ch.fcl()(f, 0) == 1 ? 0 : i, n, m) += fs(f) * inco[f] * pe(eb >= 0 ? eb : f_e(f, 0))
                                                            * (1. + (vit[f] * (i ? -1 : 1) >= 0 ? 1. : vit[f] ? -1. : 0.) * alpha_) / 2;
          }
        for (i = 0; i < 2; i++)
          if ((e = f_e(f, i)) >= 0)
            {
              for (k = 0; k < e_f.dimension(1) && (fb = e_f(e, k)) >= 0; k++)
                if (fb < domaine.nb_faces() &&  ch.fcl()(fb, 0) < 2) //"faces" part
                  {
                    if ((fc = equiv(f, i, k)) >= 0 || f_e(f, 1) < 0)
                      for (j = 0; j < 2; j++) //equivalence: face fd -> face fb
                        {
                          eb = f_e(f, j), fd = (j == i ? fb : fc); //element/face sources
                          mult = (fd < 0 || domaine.dot(&nf(fb, 0), &nf(fd, 0)) > 0 ? 1 : -1) * (fd >= 0 ? pf(fd) / pe(eb) : 1); //multiplier to convert from vf to ve
                          for (n = 0; n < N; n++)
                            for (m = 0; m < N; m++)
                              if (dfac(j, n, m))
                                {
                                  double fac = (e == f_e(f, 0) ? 1 : -1) * vfd(fb, e != f_e(fb, 0)) * dfac(j, n, m) / ve(e);
                                  if (fd >= 0) matrice(fb,fd) += fac * mult; //other computed face
                                  if (comp) matrice(fb,fb) -= fac; //v div(alpha rho v) part
                                }
                        }
                    else for (j = 0; j < 2; j++)  //no equivalence: n_f * element operator
                        {
                          for (eb = f_e(f, j), l = 0; l < e_f.dimension(1) && (fc = e_f(eb, l)) >= 0; l++)
                            {
                              double num = (e == f_e(fb, 0) ? 1 : -1) * (e == f_e(f, 0) ? 1 : -1) * fs(fc) * fs(fb) * domaine.dot(&xv(fc, 0), &xv(fb, 0), &xp(eb, 0), &xp(e, 0)) * (eb == f_e(fc, 0) ? 1 : -1);
                              double den = ve(eb) * ve(e);
                              if (std::fabs(num) > 1e-9 * den)
                                {
                                  double num_den = num/den;
                                  for (n = 0; n < N; n++)
                                    for (m = 0; m < N; m++)
                                      if (dfac(j, n, m))
                                        {
                                          double fac = dfac(j, n, m) * num_den;
                                          matrice(fb,fc) += fac;
                                        }
                                }
                            }
                          if (comp)
                            for (l = 0; l < e_f.dimension(1) && (fc = e_f(e, l)) >= 0; l++)
                              {
                                double num = (e == f_e(fb, 0) ? 1 : -1) * (e == f_e(f, 0) ? 1 : -1) * fs(fc) * fs(fb) * domaine.dot(&xv(fc, 0), &xv(fb, 0), &xp(e, 0), &xp(e, 0)) * (e == f_e(fc, 0) ? 1 : -1);
                                double den = ve(e) * ve(e);
                                if (std::fabs(num) > 1e-9 * den)
                                  {
                                    double num_den = num/den;
                                    for (n = 0; n < N; n++)
                                      for (m = 0; m < N; m++)
                                        if (dfac(j, n, m))
                                          {
                                            double fac = dfac(j, n, m) * num_den;
                                            matrice(fb,fc) -= fac;
                                          }
                                  }
                              }
                        }
                  }
            }
      }
}
 
void Op_Conv_EF_Stab_PolyMAC_CDO_Face::set_incompressible(const int flag)
{
  if (flag == 0)
    {
      Cerr << "Compressible form of operator \"" << que_suis_je() << "\" :" << finl;
      Cerr << "Discretization of \u2207(inco \u2297 v) - v \u2207.(inco)" << finl;
    }
  incompressible_ = flag;
}