Masse_Multiphase.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 <EcritureLectureSpecial.h>
#include <Pb_Multiphase_HEM.h>
#include <Masse_Multiphase.h>
#include <Champ_Uniforme.h>
#include <Matrice_Morse.h>
#include <Discret_Thyd.h>
#include <Fluide_base.h>
#include <Domaine_VF.h>
#include <TRUSTTrav.h>
#include <Domaine.h>
#include <EChaine.h>
#include <Param.h>
 
Implemente_instanciable(Masse_Multiphase, "Masse_Multiphase", Convection_Diffusion_std);
// XD Masse_Multiphase eqn_base Masse_Multiphase INHERITS_BRACE Mass consevation equation for a multi-phase problem
// XD_CONT where the unknown is the alpha (void fraction)
 
Sortie& Masse_Multiphase::printOn(Sortie& is) const { return Equation_base::printOn(is); }
 
Entree& Masse_Multiphase::readOn(Entree& is)
{
  assert(l_inco_ch_);
  assert(le_fluide_);
  terme_convectif.associer_eqn(*this), evanescence_.associer_eqn(*this);
 
  Equation_base::readOn(is);
 
  terme_convectif.set_fichier("Debit");
  terme_convectif.set_description((Nom)"Mass flow rate=Integral(-alpha*rho*u*ndS) [kg/s] if SI units used");
 
  champs_compris_.ajoute_champ(l_inco_ch_);
 
  // Special treatment for Pb_Multiphase_HEM
  // We enforce the presence of a source term related to the interfacial flux automatically.
  if (sub_type(Pb_Multiphase_HEM, probleme()))
    {
      bool check_source_FICC = false;
 
      for (int ii = 0; ii < sources().size(); ii++)
        if (sources()(ii)->que_suis_je().debute_par("Flux_interfacial"))
          check_source_FICC = true;
 
      if (!check_source_FICC)
        {
          EChaine source_FI("{ flux_interfacial }");
          lire_sources(source_FI);
        }
    }
 
  return is;
}
 
void Masse_Multiphase::set_param(Param& param) const
{
  Equation_base::set_param(param);
  param.ajouter_non_std("convection",(this));
  param.ajouter_non_std("evanescence|vanishing",(this));
  param.ajouter_condition("is_read_convection","The convection operator must be read, select negligeable type if you want to neglect it.");
}
 
int Masse_Multiphase::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  if (mot=="evanescence") is >> evanescence_;
  else return Convection_Diffusion_std::lire_motcle_non_standard(mot, is);
  return 1;
}
 
const Operateur& Masse_Multiphase::operateur(int i) const
{
  if (i)
    {
      Cerr << "Masse_Multiphase : wrong operator number " << i << finl;
      Process::exit();
    }
  return terme_convectif;
}
 
Operateur& Masse_Multiphase::operateur(int i)
{
  if (i)
    {
      Cerr << "Masse_Multiphase : wrong operator number " << i << finl;
      Process::exit();
    }
  return terme_convectif;
}
 
void Masse_Multiphase::dimensionner_matrice_sans_mem(Matrice_Morse& matrice)
{
  Equation_base::dimensionner_matrice_sans_mem(matrice);
  evanescence_->dimensionner(matrice);
}
 
int Masse_Multiphase::has_interface_blocs() const
{
  int ok = Convection_Diffusion_std::has_interface_blocs();
  if (evanescence_) ok &= evanescence_->has_interface_blocs();
  return ok;
}
 
/* evanescence is processed last */
void Masse_Multiphase::dimensionner_blocs(matrices_t matrices, const tabs_t& semi_impl) const
{
  Equation_base::dimensionner_blocs(matrices, semi_impl);
  if (evanescence_) evanescence_->dimensionner_blocs(matrices, semi_impl);
}
 
void Masse_Multiphase::assembler_blocs_avec_inertie(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl)
{
  Equation_base::assembler_blocs_avec_inertie(matrices, secmem, semi_impl);
  if (evanescence_) evanescence_->ajouter_blocs(matrices, secmem, semi_impl);
}
 
/*! @brief Associates a physical medium to the equation; the medium is cast to Fluide_base.
 *
 * @param un_milieu Physical medium to associate.
 * @throws not all physical properties of the fluid are specified
 */
void Masse_Multiphase::associer_milieu_base(const Milieu_base& un_milieu)
{
  const Fluide_base& un_fluide = ref_cast(Fluide_base,un_milieu);
  associer_fluide(un_fluide);
}
 
/*! @brief Associates the unknown of the equation to the state law.
 *
 */
void Masse_Multiphase::completer()
{
  Equation_base::completer();
  terme_convectif->check_multiphase_compatibility();
 
  const Domaine_dis_base& zdis = domaine_dis();
  if (discretisation().is_vdf())
    {
      // initialize the grad operator if VDF
      Op_Grad_.associer_eqn(*this);
      Op_Grad_.typer();
      Op_Grad_.l_op_base().associer_eqn(*this);
      const Domaine_Cl_dis_base& zcl = domaine_Cl_dis();
      const Champ_Inc_base& inco = inconnue();
      Op_Grad_->associer(zdis, zcl, inco);
    }
 
  verifier_somme_alpha();
}
 
void Masse_Multiphase::verifier_somme_alpha()
{
  const DoubleTab& vals = l_inco_ch_->valeurs();
  const int ne = vals.dimension(0), nl = vals.line_size();
  DoubleVect vals_somme(ne);
  vals_somme = 0.;
 
  for (int i = 0; i < ne; i++)
    for (int j = 0; j < nl; j++)
      vals_somme(i) += vals(i, j);
 
  const double min_a = mp_min_vect(vals_somme), max_a = mp_max_vect(vals_somme);
 
  if (min_a < 1. - 1.e-12 || max_a > 1. + 1.e-12)
    {
      Cerr << "WHAT ?? The sum of the void fraction (per cell) is not 1 !!!! You should do something !" << finl;
      Process::exit();
    }
}
 
/*! @brief Discretizes the equation.
 *
 */
void Masse_Multiphase::discretiser()
{
  int nb_valeurs_temp = schema_temps().nb_valeurs_temporelles();
  double temps = schema_temps().temps_courant();
  const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
  Cerr << "Volume fraction discretization" << finl;
  // Use temperature directive since the discretization is identical
  const Pb_Multiphase& pb = ref_cast(Pb_Multiphase, probleme());
  dis.discretiser_champ("temperature",domaine_dis(),"alpha","sans_dimension", pb.nb_phases(),nb_valeurs_temp,temps,l_inco_ch_);
  l_inco_ch_->fixer_nature_du_champ(pb.nb_phases() == 1 ? scalaire : pb.nb_phases() == dimension ? vectoriel : multi_scalaire); //pfft
  for (int i = 0; i < pb.nb_phases(); i++)
    l_inco_ch_->fixer_nom_compo(i, Nom("alpha_") + pb.nom_phase(i));
  champs_compris_.ajoute_champ(l_inco_ch_);
  Equation_base::discretiser();
 
  Cerr << "Masse_Multiphase::discretiser() ok" << finl;
}
 
/*! @brief Returns the physical medium associated with the equation.
 *
 * (a Fluide_base upcasted to Milieu_base)
 *     (const version)
 *
 * @return The Fluide_base upcasted to Milieu_base.
 */
const Milieu_base& Masse_Multiphase::milieu() const
{
  return fluide();
}
 
 
/*! @brief Returns the physical medium associated with the equation.
 *
 * (a Fluide_base upcasted to Milieu_base)
 *
 * @return The Fluide_base upcasted to Milieu_base.
 */
Milieu_base& Masse_Multiphase::milieu()
{
  return fluide();
}
 
 
/*! @brief Returns the incompressible fluid associated with the equation (const version).
 *
 * @return The incompressible fluid associated with the equation.
 * @throws If no fluid is associated with the equation.
 */
const Fluide_base& Masse_Multiphase::fluide() const
{
  if (!le_fluide_)
    {
      Cerr << "You forgot to associate the fluid to the problem named " << probleme().le_nom() << finl;
      Process::exit();
    }
  return le_fluide_.valeur();
}
 
 
/*! @brief Returns the incompressible fluid associated with the equation.
 *
 * @return The incompressible fluid associated with the equation.
 * @throws If no fluid is associated with the equation.
 */
Fluide_base& Masse_Multiphase::fluide()
{
  assert(le_fluide_);
  return le_fluide_.valeur();
}
 
/*! @brief Returns the name of the application domain of the equation.
 *
 * Here "Fraction_volumique".
 *
 * @return The name of the application domain.
 */
const Motcle& Masse_Multiphase::domaine_application() const
{
  static Motcle domaine ="Fraction_volumique";
  return domaine;
}
 
/*! @brief Associates a fluid of type Fluide_base with the equation.
 *
 * @param un_fluide The fluid medium to associate with the equation.
 */
void Masse_Multiphase::associer_fluide(const Fluide_base& un_fluide)
{
  assert(sub_type(Fluide_base,un_fluide));
  le_fluide_ = ref_cast(Fluide_base,un_fluide);
}
 
void Masse_Multiphase::calculer_alpha_rho(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Equation_base& eqn = ref_cast(Equation_base, obj);
  const Champ_base& ch_rho = eqn.milieu().masse_volumique();
  const Champ_Inc_base& ch_alpha = eqn.inconnue(),
                        *pch_rho = sub_type(Champ_Inc_base, ch_rho) ? &ref_cast(Champ_Inc_base, ch_rho) : nullptr;
 
  const DoubleTab& alpha = ch_alpha.valeurs(),
                   &rho = ch_rho.valeurs();
 
  const int nl = val.dimension_tot(0),
            N = val.line_size(),
            cR = sub_type(Champ_Uniforme, ch_rho);
 
  /* field values */
  for (int i = 0; i < nl; i++)
    for (int n = 0; n < N; n++)
      val(i, n) = alpha(i, n) * rho(!cR * i, n);
 
  /* boundary values */
  /* valeur_aux_bords can only be used if ch_rho has a domaine_dis_base */
  ch_rho.a_un_domaine_dis_base() ? bval = ch_rho.valeur_aux_bords() :
                                          ch_rho.valeur_aux(ref_cast(Domaine_VF, eqn.domaine_dis()).xv_bord(), bval);
 
  tab_multiply_any_shape(bval, ch_alpha.valeur_aux_bords(), VECT_ALL_ITEMS);
 
  /* derivatives */
  DoubleTab& d_a = deriv["alpha"]; // derivative w.r.t. alpha: rho
  d_a.resize(nl, N);
 
  for (int i = 0; i < nl; i++)
    for (int n = 0; n < N; n++)
      d_a(i, n) = rho(!cR * i, n);
 
  if (pch_rho)
    for (auto &&d_c : pch_rho->derivees()) // derivatives w.r.t. the dependencies of rho
      {
        DoubleTab& der = deriv[d_c.first];
        der.resize(nl, N);
 
        for (int i = 0; i < nl; i++)
          for (int n = 0; n < N; n++)
            der(i, n) = d_c.second(i, n) * alpha(i, n);
      }
}
 
void Masse_Multiphase::calculer_alpha_rho_conv(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Equation_base& eqn = ref_cast(Equation_base, obj);
  const Pb_Multiphase& pbm = ref_cast(Pb_Multiphase, eqn.probleme());
  const Champ_base& ch_rho = eqn.milieu().masse_volumique();
  const Champ_Inc_base& ch_alpha = eqn.inconnue(), *pch_rho = sub_type(Champ_Inc_base, ch_rho) ? &ref_cast(Champ_Inc_base, ch_rho) : nullptr;
  const DoubleTab& alpha = ch_alpha.valeurs(), &rho = ch_rho.valeurs();
  const int  nl = val.dimension_tot(0),
             N = val.line_size(),
             cR = sub_type(Champ_Uniforme, ch_rho);
 
  /* field values */
  for (int i = 0; i < nl; i++)
    for (int n = 0; n < N; n++)
      val(i, n) = (alpha(i, n) - pbm.alpha_inf_phase(n)) * rho(!cR * i, n);
 
  /* boundary values */
  /* valeur_aux_bords can only be used if ch_rho has a domaine_dis_base */
  ch_rho.a_un_domaine_dis_base() ? bval = ch_rho.valeur_aux_bords() :
                                          ch_rho.valeur_aux(ref_cast(Domaine_VF, eqn.domaine_dis()).xv_bord(), bval);
 
  // tab_multiply_any_shape(bval, ch_alpha.valeur_aux_bords(), VECT_ALL_ITEMS);
  DoubleTrav calpha;
  calpha = ch_alpha.valeur_aux_bords();
 
  const int n2 = calpha.dimension_tot(0);
  for (int i = 0; i < n2; i++)
    for (int n = 0; n < N; n++)
      calpha(i, n) = calpha(i, n) - pbm.alpha_inf_phase(n);
 
  tab_multiply_any_shape(bval, calpha, VECT_ALL_ITEMS);
 
  /* derivatives */
  DoubleTab& d_a = deriv["alpha"]; // derivative w.r.t. alpha: rho
  d_a.resize(nl, N);
 
  for (int i = 0; i < nl; i++)
    for (int n = 0; n < N; n++)
      d_a(i, n) = rho(!cR * i, n);
 
  if (pch_rho)
    for (auto &&d_c : pch_rho->derivees()) // derivatives w.r.t. the dependencies of rho
      {
        DoubleTab& der = deriv[d_c.first];
        der.resize(nl, N);
 
        for (int i = 0; i < nl; i++)
          for (int n = 0; n < N; n++)
            der(i, n) = d_c.second(i, n) * alpha(i, n);
      }
}
 
void Masse_Multiphase::init_champ_convecte() const
{
  if (champ_convecte_)
    return; // already done
 
  const int Nt = inconnue().nb_valeurs_temporelles(),
            Nl = inconnue().valeurs().size_reelle_ok() ? inconnue().valeurs().dimension(0) : -1,
            Nc = inconnue().valeurs().line_size();
 
  // champ_convecte_: same type / support as the unknown
  discretisation().creer_champ(champ_convecte_, domaine_dis(), inconnue().que_suis_je(), "N/A", "N/A", Nc, Nl, Nt, schema_temps().temps_courant());
 
  champ_convecte_->associer_eqn(*this);
  auto nom_fonc = get_fonc_champ_convecte();
  champ_convecte_->nommer(nom_fonc.first);
  champ_convecte_->init_champ_calcule(*this, nom_fonc.second);
}