QDM_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 <Operateur_Diff_base.h>
#include <Discretisation_base.h>
#include <Schema_Temps_base.h>
#include <Pb_Multiphase_HEM.h>
#include <Champ_Composite.h>
#include <TRUSTTab_parts.h>
#include <QDM_Multiphase.h>
#include <Discret_Thyd.h>
#include <Fluide_base.h>
#include <EChaine.h>
#include <Param.h>
#include <vector>
#include <Nom.h>
 
Implemente_instanciable(QDM_Multiphase,"QDM_Multiphase",Navier_Stokes_std);
// XD QDM_Multiphase eqn_base QDM_Multiphase INHERITS_BRACE Momentum conservation equation for a multi-phase problem
// XD_CONT where the unknown is the velocity
// XD attr solveur_pression solveur_sys_base solveur_pression OPT Linear pressure system resolution method.
// XD attr evanescence bloc_lecture evanescence OPT Management of the vanishing phase (when alpha tends to 0 or 1)
 
// evanescence interprete evanescence 1 Management of the vanishing phase (when alpha tends to 0 or 1)
// attr homogene chaine homogene 1 Vanishing phases management : tends to homogeneous model when a phase vanishes (vl = vg)
// attr bloc_val bloc_lecture bloc_val 1 not set
 
// bloc_lecture interprete nul 1 not set
// attr alpha_res_min flottant alpha_res_min 0 Activation threshold for full replacement of vanishing phase equation (default value : 0)
// attr alpha_res flottant alpha_res 0 Activation threshold for gradual replacement of vanishing phase equation (tends to full replacement when alpha tends to alpha_res_min)
 
Sortie& QDM_Multiphase::printOn(Sortie& is) const
{
  return Equation_base::printOn(is);
}
 
/*! @brief Calls Equation_base::readOn(Entree& is). Verifies that the following have been read:
 *
 *         - the diffusive term,
 *         - the convective term,
 *         - the pressure solver
 *
 * @param is Input stream.
 * @return Modified input stream.
 * @throws diffusive term not specified in dataset; specify a negligible type for the operator if it is to be neglected
 * @throws convective term not specified in dataset; specify a negligible type for the operator if it is to be neglected
 * @throws pressure solver not defined in dataset
 */
 
Entree& QDM_Multiphase::readOn(Entree& is)
{
  evanescence_.associer_eqn(*this);
  Navier_Stokes_std::readOn(is);
  assert(le_fluide);
  if (!sub_type(Fluide_base,le_fluide.valeur()))
    {
      Cerr<<"ERROR : the QDM_Multiphase equation can be associated only to a fluid."<<finl;
      exit();
    }
  divergence.set_description((Nom)"Mass flow rate=Integral(rho*u*ndS) [kg.s-1]");
  terme_convectif->set_incompressible(0);
 
  const Pb_Multiphase& pb = ref_cast(Pb_Multiphase, probleme());
  if (!evanescence_ && pb.nb_phases() > 1)
    {
      // Special treatment for Pb_Multiphase_HEM : We enforce the evanescence to a specific value
      if (sub_type(Pb_Multiphase_HEM, probleme()))
        {
          EChaine eva("{ homogene { alpha_res 1 alpha_res_min 0.5 } }");
          eva >> evanescence_;
        }
      else
        {
          EChaine eva("{ homogene { alpha_res 1e-6 } }");
          eva >> evanescence_;
        }
    }
 
  /* velocity fields per phase for post-processing */
  noms_vit_phases_.dimensionner(pb.nb_phases()), vit_phases_.resize(pb.nb_phases());
  for (int i = 0; i < pb.nb_phases(); i++)
    noms_vit_phases_[i] = Nom("vitesse_") + pb.nom_phase(i);
 
  noms_grad_vit_phases_.dimensionner(pb.nb_phases()), grad_vit_phases_.resize(pb.nb_phases());
  for (int i = 0; i < pb.nb_phases(); i++)
    noms_grad_vit_phases_[i] = Nom("gradient_vitesse_") + pb.nom_phase(i);
 
  return is;
}
 
void QDM_Multiphase::set_param(Param& param) const
{
  Navier_Stokes_std::set_param(param);
  param.ajouter_non_std("evanescence|vanishing",(this));
}
 
int QDM_Multiphase::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  if (mot=="evanescence") is >> evanescence_;
  else return Navier_Stokes_std::lire_motcle_non_standard(mot, is);
  return 1;
}
 
void QDM_Multiphase::dimensionner_matrice_sans_mem(Matrice_Morse& matrice)
{
  Navier_Stokes_std::dimensionner_matrice_sans_mem(matrice);
  if (evanescence_) evanescence_->dimensionner(matrice);
}
 
int QDM_Multiphase::has_interface_blocs() const
{
  int ok = Navier_Stokes_std::has_interface_blocs();
  if (evanescence_) ok &= evanescence_->has_interface_blocs();
  return ok;
}
 
/* evanescence is processed last */
void QDM_Multiphase::dimensionner_blocs(matrices_t matrices, const tabs_t& semi_impl) const
{
  Navier_Stokes_std::dimensionner_blocs(matrices, semi_impl);
  if (evanescence_) evanescence_->dimensionner_blocs(matrices, semi_impl);
}
 
void QDM_Multiphase::assembler_blocs_avec_inertie(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl)
{
  Navier_Stokes_std::assembler_blocs_avec_inertie(matrices, secmem, semi_impl);
  if (evanescence_) evanescence_->ajouter_blocs(matrices, secmem, semi_impl);
}
 
void QDM_Multiphase::mettre_a_jour(double temps)
{
  Equation_base::mettre_a_jour(temps);  // skip Navier_Stokes_std::mettre_a_jour
  pression().mettre_a_jour(temps);
  pression_pa().mettre_a_jour(temps);
 
  int i, j, n, N = ref_cast(Pb_Multiphase, probleme()).nb_phases(), d, D = dimension;
  for (n = 0; n < N; n++)
    if (vit_phases_[n])
      {
        vit_phases_[n]->mettre_a_jour(temps);
        DoubleTab_parts psrc(inconnue().valeurs()), pdst(vit_phases_[n]->valeurs());
        for (i = 0; i < std::min(psrc.size(), pdst.size()); i++)
          {
            DoubleTab& src = psrc[i], &dst = pdst[i];
            if (src.line_size() == N) /* one column per component */
              for (j = 0; j < src.dimension_tot(0); j++) dst(j) = src(j, n);
            else if (src.line_size() == N * D) /* storage N * d + n */
              for (j = 0; j < src.dimension_tot(0); j++)
                for (d = 0; d < D; d++)
                  dst(j, d) = src(j, N * d + n);
            else abort(); // unknown storage layout
          }
      }
 
  if (grad_u) grad_u->mettre_a_jour(temps);
  if (la_vorticite) la_vorticite->mettre_a_jour(temps);
  if (Taux_cisaillement) Taux_cisaillement->mettre_a_jour(temps);
 
  const bool is_vdf = discretisation().is_vdf();
  for (n = 0; n < N; n++)
    if (grad_vit_phases_[n])
      {
        if (is_vdf)
          {
            for (int e = 0; e < domaine_dis().nb_elem(); e++)
              for (i = 0; i < D; i++)
                for (j = 0; j < D; j++)
                  grad_vit_phases_[n]->valeurs()(e, D * i + j) = grad_u->valeurs()(e, N * (D * i + j) + n);
          }
        else
          {
            DoubleTab_parts psrc(grad_u->valeurs()), pdst(grad_vit_phases_[n]->valeurs());
            for (i = 0; i < psrc.size(); i++)
              for (j = 0; j < pdst.size(); j++)
                {
                  DoubleTab& src = psrc[i], &dst = pdst[j];
                  if (src.line_size() == N * D * D && dst.line_size() == D * D) /* une colonne par composante */
                    {
                      for (int k = 0; k < src.dimension_tot(0); k++)
                        for (int dU = 0; dU < D; dU++)
                          for (int dX = 0; dX < D; dX++)
                            dst(k, dX + D * dU) = src(k, dX, dU + n * D); // rows and columns are swapped when going through DoubleTab_parts
                    }
                }
          }
        grad_vit_phases_[n]->mettre_a_jour(temps);
      }
  if (gradient_P)
    {
      gradient_P->valeurs() = 0;
      gradient->ajouter(la_pression->valeurs(), gradient_P->valeurs());
      solv_masse().appliquer_impl(gradient_P->valeurs());
    }
}
 
bool QDM_Multiphase::initTimeStep(double dt)
{
  Schema_Temps_base& sch=schema_temps();
  ConstDoubleTab_parts ppart(pression().valeurs());
  /* if pression_pa() is smaller than pression() (e.g., auxiliary variables PolyMAC_HFV), copy only the first part */
  const DoubleTab& p_red = pression_pa().valeurs().dimension_tot(0) < pression().valeurs().dimension_tot(0) ? ppart[0] : pression().valeurs();
  for (int i=1; i<=sch.nb_valeurs_futures(); i++)
    {
      // Update time in the pressure field
      pression().changer_temps_futur(sch.temps_futur(i),i);
      pression().futur(i)=pression().valeurs();
      pression_pa().changer_temps_futur(sch.temps_futur(i),i);
      pression_pa().futur(i) = p_red;
    }
  return Equation_base::initTimeStep(dt);
}
 
void QDM_Multiphase::abortTimeStep()
{
  Equation_base::abortTimeStep();
}
 
void QDM_Multiphase::discretiser_vitesse()
{
  const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
  dis.vitesse(schema_temps(), domaine_dis(), la_vitesse, ref_cast(Pb_Multiphase, probleme()).nb_phases());
}
 
void QDM_Multiphase::discretiser_grad_p()
{
// Does nothing! Called by default in Navier_Stokes_std.discretiser() but not required in Pb_Multiphase.
// The discretization happens in QDM_Multiphase.creer_champ()
}
 
const Champ_Don_base& QDM_Multiphase::diffusivite_pour_transport() const
{
  return le_fluide->viscosite_dynamique();
}
 
const Champ_base& QDM_Multiphase::diffusivite_pour_pas_de_temps() const
{
  return le_fluide->viscosite_cinematique();
}
 
const Champ_base& QDM_Multiphase::vitesse_pour_transport() const
{
  return la_vitesse;
}
 
/*! @brief Completes the base equation, associates the pressure to the equation,
 *
 *     completes the divergence, gradient and pressure solver.
 *     Adds 2 source terms: one representing the centrifugal force
 *     in the axisymmetric case, the other involved in the
 *     2D axisymmetric resolution.
 *
 */
void QDM_Multiphase::completer()
{
  Cerr<<" Navier_Stokes_std::completer_deb"<<finl;
  Navier_Stokes_std::completer();
  Cerr<<" Navier_Stokes_std::completer_fin"<<finl;
  Cerr<<"unknow field type  "<<inconnue().que_suis_je()<<finl;
  Cerr<<"unknow field name  "<<inconnue().le_nom()<<finl;
  Cerr<<"equation type "<<inconnue().equation().que_suis_je()<<finl;
 
  /* list of things that must be compatible with multiphase */
  std::vector<const MorEqn*> morceaux = { &solveur_masse.valeur(), &les_sources, &terme_convectif.valeur(), &terme_diffusif.valeur(), &gradient.valeur() };
  for (auto mor : morceaux) mor->check_multiphase_compatibility();
}
 
void QDM_Multiphase::get_noms_champs_postraitables(Noms& noms,Option opt) const
{
  Navier_Stokes_std::get_noms_champs_postraitables(noms,opt);
 
  Noms noms_compris;
  const Pb_Multiphase& pb = ref_cast(Pb_Multiphase, probleme());
  for (int i = 0; i < pb.nb_phases(); i++)
    {
      noms_compris.add(noms_grad_vit_phases_[i]);
      noms_compris.add(noms_vit_phases_[i]);
    }
  if (opt==DESCRIPTION)
    Cerr<<" QDM_Multiphase : "<< noms_compris <<finl;
  else
    noms.add(noms_compris);
}
 
void QDM_Multiphase::creer_champ(const Motcle& motlu)
{
  Navier_Stokes_std::creer_champ(motlu);
  if (la_vorticite)
    if (!grad_u) creer_champ("gradient_vitesse");
  int i = noms_vit_phases_.rang(motlu);
  if (i >= 0 && !vit_phases_[i])
    {
      discretisation().discretiser_champ("vitesse",domaine_dis(), noms_vit_phases_[i], "m/s",dimension, 1, 0, vit_phases_[i]);
      champs_compris_.ajoute_champ(vit_phases_[i]);
    }
  i = noms_grad_vit_phases_.rang(motlu);
  if (i >= 0 && !grad_vit_phases_[i])
    {
      int D = dimension ;
      Noms noms(D * D), unites(D * D);
      std::vector<Nom> composantsVitesse({Nom("dU"), Nom("dV"), Nom("dW")});
      std::vector<Nom> composantsDerivee({Nom("dx"), Nom("dy"), Nom("dz")});
      for (int dU = 0 ; dU< D ; dU++)
        for (int dX = 0 ; dX < D ; dX++)
          {
            noms[ D * dU + dX]=Nom(composantsVitesse[dU] + composantsDerivee[dX]);
            unites[ D * dU + dX] = Nom("m2/s");
          }
      noms[0] = noms_grad_vit_phases_[i]; // To give it the correct name in discretiser_champ; consequence: the first output coordinate does not have the right name
      Motcle typeChamp = "champ_elem" ;
      discretisation().discretiser_champ(typeChamp, domaine_dis(), multi_scalaire, noms , unites, D*D, 0, grad_vit_phases_[0]);
      champs_compris_.ajoute_champ(grad_vit_phases_[i]);
    }
 
  if (motlu == "gradient_pression")
    {
      if (!gradient_P)
        {
          const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
          dis.gradient_P(schema_temps(), domaine_dis(), gradient_P, ref_cast(Pb_Multiphase, probleme()).nb_phases());
          champs_compris_.ajoute_champ(gradient_P);
        }
    }
 
}
 
Entree& QDM_Multiphase::lire_cond_init(Entree& is)
{
  Cerr << "Reading of initial conditions\n";
  Nom nom;
  is >> nom;
  if (nom != "{")
    {
      Cerr << que_suis_je() << ": expected { instead of " << nom << finl;
      Process::exit();
    }
  int vit_lu = 0, press_lu = 0;
  for (is >> nom; nom != "}"; is >> nom)
    if (nom == "vitesse" || nom == "velocity")
      {
        OWN_PTR(Champ_Don_base) src;
        is >> src;
 
        if (src->que_suis_je() == "Champ_Composite")
          {
            const int nb_phases = ref_cast(Pb_Multiphase, probleme()).nb_phases(), nb_dim = ref_cast(Champ_Composite,src.valeur()).get_champ_composite_dim();
            if ( nb_dim != nb_phases)
              {
                Cerr << que_suis_je() << ": velocity initial condition Champ_Composite should have "<<  nb_phases << " fields and not " << nb_dim << " !" << finl;
                Process::exit();
              }
          }
 
        verifie_ch_init_nb_comp(la_vitesse, src->nb_comp());
        la_vitesse->affecter(src), vit_lu = 1;
        la_vitesse->passe() = la_vitesse->valeurs();
      }
    else if (nom == "pression" || nom == "pressure")
      {
        OWN_PTR(Champ_Don_base) src;
        is >> src, verifie_ch_init_nb_comp(la_pression, src->nb_comp());
        la_pression->affecter(src);
        la_pression_en_pa->passe() = la_pression_en_pa->valeurs() = la_pression->passe() = la_pression->valeurs();
        press_lu = 1;
      }
    else
      {
        Cerr << que_suis_je() << ": expected vitesse|velocity|pression|pressure instead of " << nom << finl;
        Process::exit();
      }
 
  if (!vit_lu)
    {
      Cerr << que_suis_je() << ": velocity initial condition not found." << finl;
      Process::exit();
    }
  if (!press_lu)
    {
      Cerr << que_suis_je() << ": pressure initial condition not found." << finl;
      Process::exit();
    }
 
  return is;
}
 
int QDM_Multiphase::preparer_calcul()
{
  Equation_base::preparer_calcul(); // to avoid calling Navier_Stokes_std::preparer_calcul() !
 
  // XXX Elie Saikali : useful for restart case!
  const double temps = schema_temps().temps_courant();
  pression().changer_temps(temps);
  pression_pa().changer_temps(temps);
 
  if (evanescence_)
    evanescence_->preparer_calcul();
 
  return 1;
}
 
double QDM_Multiphase::alpha_res() const
{
  const Pb_Multiphase& pb = ref_cast(Pb_Multiphase, probleme());
  if (pb.nb_phases() == 1) return 0.;
  if (!evanescence_) Process::exit( "QDM_Multiphase::alpha_res : the evanescence operator should have been created already !" );
  if sub_type(Operateur_Evanescence_base, evanescence_.valeur()) return ref_cast(Operateur_Evanescence_base, evanescence_.valeur()).alpha_res();
  return -1.;
}