Milieu_composite.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 <Fluide_Incompressible.h>
#include <Discretisation_base.h>
#include <Schema_Temps_base.h>
#include <Milieu_composite.h>
#include <Champ_Uniforme.h>
#include <Equation_base.h>
#include <Probleme_base.h>
#include <Pb_Multiphase.h>
#include <Interprete.h>
#include <Domaine_VF.h>
 
Implemente_instanciable(Milieu_composite, "Milieu_composite", Fluide_base);
// XD Milieu_composite listobj Milieu_composite INHERITS_BRACE milieu_base NO_COMMA Composite medium made of several sub
// XD_CONT mediums.
 
Sortie& Milieu_composite::printOn(Sortie& os) const { return os; }
 
Entree& Milieu_composite::readOn(Entree& is)
{
  int i = 0;
  std::vector<std::pair<std::string, int>> especes; // string for phase_name. int: 0 for liquid and 1 for gas
  Cerr << "Reading Milieu_composite..." << finl;
  Nom mot;
  is >> mot;
  if (mot != "{")
    {
      Cerr << "Milieu_composite : { expected instead of " << mot << finl;
      Process::exit();
    }
 
  for (is >> mot; mot != "}"; is >> mot)
    {
      if (Motcle(mot) == "POROSITES_CHAMP") is >> ch_porosites_;
      else if (Motcle(mot) == "DIAMETRE_HYD_CHAMP") is >> ch_diametre_hyd_;
      else if (Motcle(mot) == "POROSITES")
        {
          Cerr << "You should use porosites_champ and not porosites ! Call the 911 !" << finl;
          Process::exit();
        }
      else if (Motcle(mot) == "GRAVITE")
        {
          Cerr << que_suis_je() << " : gravity should not be defined in Pb_Multiphase ! Use source_qdm if you want gravity in QDM equation !" << finl;
          Process::exit();
        }
      else if (!mot.debute_par("saturation") && !mot.debute_par("interface")) // add phases
        {
          noms_phases_.add(mot);
          Cerr << "Milieu_composite: add phase " << mot << " ... " << finl;
          OWN_PTR(Fluide_base) fluide;
          fluide.typer_lire_simple(is, "Typing the fluid medium ...");
 
          if (fluide->has_porosites())
            {
              Cerr << que_suis_je() + " : porosity should be defined only once in the milieu_composite block, not in " + fluide->que_suis_je() << finl;
              Process::exit();
            }
          if (fluide->a_gravite())
            {
              Cerr << que_suis_je() + " : gravity should be defined only once in the milieu_composite block, not in " + fluide->que_suis_je() << finl;
              Process::exit();
            }
          if (fluide->has_hydr_diam())
            {
              Cerr << que_suis_je() + " : diametre_hyd_champ should be defined only once in the milieu_composite block, not in " + fluide->que_suis_je() << finl;
              Process::exit();
            }
          if (!sub_type(Fluide_base,fluide.valeur()))
            {
              Cerr << que_suis_je() + " : only real fluids should be read and not a fluid of type " + fluide->que_suis_je() << finl;
              Process::exit();
            }
 
          fluide->set_id_composite(i++);
          fluide->nommer(mot); // XXX
          fluides_.push_back(fluide);
          especes.push_back(check_fluid_name(fluide->le_nom()));
        }
      else if (mot.debute_par("saturation")) // add saturation
        {
          has_saturation_ = true;
          Cerr << "Milieu_composite: add saturation " << mot << " ... " << finl;
          sat_lu_.typer_lire_simple(is, "Typing the saturation ...");
        }
      else // add interface
        {
          has_interface_ = true;
          Cerr << "Milieu_composite: add interface " << mot << " ... " << finl;
          inter_lu_.typer_lire_simple(is, "Typing the interface ...");
        }
    }
 
  if (has_saturation() && has_interface())
    {
      Cerr << "You define both interface and saturation in Milieu_composite ???" << finl;
      Process::exit();
    }
 
  // Interface processing
  const int N = (int)fluides_.size();
  for (int n = 0; n < N; n++)
    {
      std::vector<Interface_base *> inter;
      for (int m = 0; m < N; m++)
        {
          const std::string espn = especes[n].first, espm = especes[m].first;
          const int pn = especes[n].second, pm = especes[m].second;
 
          if (pn != pm && (has_interface() || (espn == espm && has_saturation())))
            {
              Cerr << "Interface between fluid " << n << " : " << fluides_[n]->le_nom() << " and " << m << " : " << fluides_[m]->le_nom() << finl;
              inter.push_back(&ref_cast(Interface_base, has_saturation_ ? sat_lu_.valeur() : inter_lu_.valeur()));
              const Saturation_base *sat = sub_type(Saturation_base, *inter.back()) ? &ref_cast(Saturation_base, *inter.back()) : nullptr;
              if (sat && sat->get_Pref() > 0) // for the law e = e0 + cp * (T - T0)
                {
                  const double hn = pn ? sat->Hvs(sat->get_Pref()) : sat->Hls(sat->get_Pref()),
                               hm = pm ? sat->Hvs(sat->get_Pref()) : sat->Hls(sat->get_Pref()),
                               T0 = sat->Tsat(sat->get_Pref());
                  fluides_[n]->set_h0_T0(hn, T0), fluides_[m]->set_h0_T0(hm, T0);
                }
            }
          else inter.push_back(nullptr);
        }
      tab_interface_.push_back(inter);
    }
  return is;
}
 
std::pair<std::string, int> Milieu_composite::check_fluid_name(const Nom& name)
{
  int phase = name.debute_par("gaz");
  Nom espece = phase ? name.getSuffix("gaz_") : name.getSuffix("liquide_");
 
  // Remove suffixes "_group1" and "_group2" for IATE 2-group model
  std::string nomjdd = espece.getString();
  std::string suffix1 = "_group1";
  std::string suffix2 = "_group2";
  if (nomjdd.size() >= suffix1.size() && nomjdd.compare(nomjdd.size() - suffix1.size(), suffix1.size(), suffix1) == 0)
    {
      nomjdd.erase(nomjdd.size() - suffix1.size(), suffix1.size());
    }
  else if (nomjdd.size() >= suffix2.size() && nomjdd.compare(nomjdd.size() - suffix2.size(), suffix2.size(), suffix2) == 0)
    {
      nomjdd.erase(nomjdd.size() - suffix2.size(), suffix2.size());
    }
 
  return std::make_pair(nomjdd, phase);
}
 
bool Milieu_composite::has_interface(int k, int l) const
{
  if (tab_interface_[k][l]) return true;
  return false;
}
 
Interface_base& Milieu_composite::get_interface(int k, int l) const
{
  return *(tab_interface_[k][l]);
}
 
bool Milieu_composite::has_saturation(int k, int l) const
{
  if (tab_interface_[k][l] && sub_type(Saturation_base, *tab_interface_[k][l])) return true;
  return false;
}
 
Saturation_base& Milieu_composite::get_saturation(int k, int l) const
{
  return ref_cast(Saturation_base, *(tab_interface_[k][l]));
}
 
const Fluide_base& Milieu_composite::get_fluid(const int i) const
{
  assert(i >= 0 && i < (int )fluides_.size());
  return fluides_[i].valeur();
}
 
Fluide_base& Milieu_composite::get_fluid(const int i)
{
  assert(i >= 0 && i < (int )fluides_.size());
  return fluides_[i].valeur();
}
 
int Milieu_composite::initialiser(const double temps)
{
  for (auto &itr : fluides_) itr->initialiser(temps);
 
  const bool res_en_T = equation_.count("temperature") ? true : false;
 
  const Equation_base& eqn = res_en_T ? equation("temperature") : equation("enthalpie");
  Champ_Inc_base& ch_rho = ref_cast(Champ_Inc_base, ch_rho_.valeur()),
                  &ch_e = ref_cast(Champ_Inc_base, ch_e_int_.valeur()),
                   &ch_h_ou_T = ref_cast(Champ_Inc_base, ch_h_ou_T_.valeur());
 
  ch_rho.associer_eqn(eqn);
  ch_rho.init_champ_calcule(*this, calculer_masse_volumique);
 
  ch_e.associer_eqn(eqn);
  ch_e.init_champ_calcule(*this, calculer_energie_interne);
 
  ch_h_ou_T.associer_eqn(eqn);
  res_en_T ? ch_h_ou_T.init_champ_calcule(*this, calculer_enthalpie) : ch_h_ou_T.init_champ_calcule(*this, calculer_temperature_multiphase);
 
  t_init_ = temps;
 
  // XXX Elie Saikali: useful for restart cases!
  ch_rho_->changer_temps(temps);
  ch_e_int_->changer_temps(temps);
  ch_h_ou_T_->changer_temps(temps);
 
  return Milieu_base::initialiser_porosite(temps);
}
 
void Milieu_composite::preparer_calcul()
{
  mettre_a_jour(t_init_);
  fluid_properties_initialised_ = true;
}
 
void Milieu_composite::discretiser(const Probleme_base& pb, const  Discretisation_base& dis)
{
  Cerr << "Composite medium discretization." << finl;
 
  const int N = (int)fluides_.size(), nc = pb.equation(0).inconnue().nb_valeurs_temporelles();
  const Domaine_dis_base& domaine_dis = pb.equation(0).domaine_dis();
  const double temps = pb.schema_temps().temps_courant();
 
  res_en_T_ = sub_type(Pb_Multiphase,pb) ? ref_cast(Pb_Multiphase,pb).resolution_en_T() : true;
 
  for (auto& itr : fluides_) itr->discretiser(pb, dis);
 
  /* masse volumique, energie interne, enthalpie : champ_inc */
  OWN_PTR(Champ_Inc_base) rho_inc, ei_inc, h_ou_T_inc;
  dis.discretiser_champ("champ_elem", domaine_dis, "masse_volumique", "kg/m^3", N, nc, temps, rho_inc);
  dis.discretiser_champ("champ_elem", domaine_dis, "energie_interne", "J/m^3", N, nc, temps, ei_inc);
  if (res_en_T_)
    dis.discretiser_champ("champ_elem", domaine_dis, "enthalpie", "J/m^3", N, nc, temps, h_ou_T_inc);
  else
    dis.discretiser_champ("champ_elem", domaine_dis, "temperature", "C", N, nc, temps, h_ou_T_inc);
 
  ch_rho_ = rho_inc, ch_e_int_ = ei_inc, ch_h_ou_T_ = h_ou_T_inc;
 
  /* autres champs : champ_fonc */
  dis.discretiser_champ("champ_elem", domaine_dis, "viscosite_dynamique", "kg/m/s", N, temps, ch_mu_);
  dis.discretiser_champ("champ_elem", domaine_dis, "viscosite_cinematique", "m2/s", N, temps, ch_nu_);
  dis.discretiser_champ("champ_elem", domaine_dis, "diffusivite", "m2/s", N, temps, ch_alpha_);
  dis.discretiser_champ("champ_elem", domaine_dis, "alpha_fois_rho", "kg/m/s", N, temps, ch_alpha_fois_rho_);
  dis.discretiser_champ("champ_elem", domaine_dis, "conductivite", "W/m/K", N, temps, ch_lambda_);
  dis.discretiser_champ("champ_elem", domaine_dis, "capacite_calorifique", "J/kg/K", N, temps, ch_Cp_);
  dis.discretiser_champ("champ_elem", domaine_dis, "masse_volumique_melange", "kg/m^3", 1, temps, rho_m_);
  dis.discretiser_champ("champ_elem", domaine_dis, "enthalpie_melange", "J/m^3", 1, temps, h_m_);
 
  champs_compris_.ajoute_champ(ch_rho_);
  champs_compris_.ajoute_champ(ch_e_int_);
  champs_compris_.ajoute_champ(ch_h_ou_T_);
 
  std::vector<OWN_PTR(Champ_Don_base)* > fields = {&ch_mu_, &ch_nu_, &ch_lambda_, &ch_alpha_, &ch_alpha_fois_rho_, &ch_Cp_, &rho_m_, &h_m_};
  for (auto && f: fields) champs_compris_.ajoute_champ((*f).valeur());
 
  // discretize sigma / Tsat fields if needed ...
  if ((int) tab_interface_.size() > 0)
    {
      Cerr << "Milieu_composite::" << __func__ << " ==> Surface tension discretization ..." << finl;
      if (has_saturation_) Cerr << "Milieu_composite::" << __func__ << " ==> Saturation temperature discretization ..." << finl;
 
      for (int k = 0; k < N; k++)
        for (int l = k + 1; l < N; l++)
          {
            int phase = fluides_[k]->le_nom().debute_par("gaz");
            Nom espece = phase ? fluides_[k]->le_nom().getSuffix("gaz_") : fluides_[k]->le_nom().getSuffix("liquide_");
            if (has_interface(k, l)) // OK if interface/saturation
              {
                Interface_base& inter = get_interface(k, l);
                inter.assoscier_pb(pb);
                Nom sig_nom = Nom("surface_tension_") + espece;
                inter.discretiser_sigma(sig_nom, temps);
                champs_compris_.ajoute_champ(inter.get_sigma_champ());
              }
 
            if (has_saturation(k, l)) // OK if saturation only
              {
                Saturation_base& sat = get_saturation(k, l);
                Nom Tsat_nom = Nom("Tsat_") + espece;
                sat.discretiser_Tsat(Tsat_nom, temps);
                champs_compris_.ajoute_champ(sat.get_Tsat_champ());
              }
          }
    }
 
  // Finally, discretize porosity + hydraulic diameter
  Milieu_base::discretiser_porosite(pb,dis);
  Milieu_base::discretiser_diametre_hydro(pb,dis);
}
 
void Milieu_composite::mettre_a_jour(double temps)
{
  for (auto& itr : fluides_) itr->mettre_a_jour(temps);
 
  ch_rho_->mettre_a_jour(temps);
  ch_e_int_->mettre_a_jour(temps);
  ch_h_ou_T_->mettre_a_jour(temps);
 
  std::vector<OWN_PTR(Champ_Don_base)* > fields = {&ch_mu_, &ch_nu_, &ch_lambda_, &ch_alpha_, &ch_alpha_fois_rho_, &ch_Cp_, &rho_m_, &h_m_};
  for (auto && f: fields) (*f)->mettre_a_jour(temps);
 
  mettre_a_jour_tabs();
 
  // MAJ sigma & Tsat
  if ((int) tab_interface_.size() > 0)
    {
      const int N = (int) fluides_.size();
      for (int k = 0; k < N; k++)
        for (int l = k + 1; l < N; l++)
          if (has_interface(k, l)) // OK if interface/saturation
            get_interface(k, l).mettre_a_jour(temps);
    }
 
  Milieu_base::mettre_a_jour_porosite(temps);
}
 
void Milieu_composite::mettre_a_jour_tabs()
{
  const int N = (int)fluides_.size();
 
  /* viscosite dynamique */
  {
    DoubleTab& tab = ch_mu_->valeurs();
    const int Nl = ch_mu_->valeurs().dimension_tot(0);
    for (int n = 0; n < N; n++)
      {
        const Champ_base& ch_n = fluides_[n]->viscosite_dynamique();
        const int cch = sub_type(Champ_Uniforme, ch_n);
        for (int i = 0; i < Nl; i++) tab(i, n) = ch_n.valeurs()(!cch * i, 0);
      }
  }
 
  /* viscosite cinematique */
  {
    DoubleTab& tab = ch_nu_->valeurs();
    const int Nl = ch_nu_->valeurs().dimension_tot(0);
    for (int n = 0; n < N; n++)
      {
        const Champ_base& ch_n = fluides_[n]->viscosite_cinematique();
        const int cch = sub_type(Champ_Uniforme, ch_n);
        for (int i = 0; i < Nl; i++) tab(i, n) = ch_n.valeurs()(!cch * i, 0);
      }
  }
 
  /* diffusivite */
  {
    DoubleTab& tab = ch_alpha_->valeurs();
    const int Nl = ch_alpha_->valeurs().dimension_tot(0);
    for (int n = 0; n < N; n++)
      {
        const Champ_base& ch_n = fluides_[n]->diffusivite();
        const int cch = sub_type(Champ_Uniforme, ch_n);
        for (int i = 0; i < Nl; i++) tab(i, n) = ch_n.valeurs()(!cch * i, 0);
      }
  }
 
  /* alpha fois rho */
  {
    DoubleTab& tab = ch_alpha_fois_rho_->valeurs();
    const int Nl = ch_alpha_fois_rho_->valeurs().dimension_tot(0);
    for (int n = 0; n < N; n++)
      {
        const Champ_base& ch_n = fluides_[n]->diffusivite_fois_rho();
        const int cch = sub_type(Champ_Uniforme, ch_n);
        for (int i = 0; i < Nl; i++) tab(i, n) = ch_n.valeurs()(!cch * i, 0);
      }
  }
 
  /* conductivite */
  {
    DoubleTab& tab = ch_lambda_->valeurs();
    const int Nl = ch_lambda_->valeurs().dimension_tot(0);
    for (int n = 0; n < N; n++)
      {
        const Champ_base& ch_n = fluides_[n]->conductivite();
        const int cch = sub_type(Champ_Uniforme, ch_n);
        for (int i = 0; i < Nl; i++) tab(i, n) = ch_n.valeurs()(!cch * i, 0);
      }
  }
 
  /* capacite calorifique */
  {
    DoubleTab& tab = ch_Cp_->valeurs();
    const int Nl = ch_Cp_->valeurs().dimension_tot(0);
    for (int n = 0; n < N; n++)
      {
        const Champ_base& ch_n = fluides_[n]->capacite_calorifique();
        const int cch = sub_type(Champ_Uniforme, ch_n);
        for (int i = 0; i < Nl; i++) tab(i, n) = ch_n.valeurs()(!cch * i, 0);
      }
  }
 
  /* calcul des quantites melange */
  DoubleTab& trm = rho_m_->valeurs(), &thm = h_m_->valeurs();
  trm = 0, thm = 0;
 
  const DoubleTab& a = equation("alpha").inconnue().valeurs(), &r = ch_rho_->valeurs(),
                   &ent = res_en_T_ ? ch_h_ou_T_->valeurs() : equation("enthalpie").inconnue().valeurs();
 
  const int Nl = rho_m_->valeurs().dimension_tot(0);
 
  // density
  for (int n = 0; n < N; n++)
    for (int i = 0; i < Nl; i++) trm(i) += a(i, n) * r(i, n);
 
  // enthalpy
  for (int n = 0; n < N; n++)
    for (int i = 0; i < Nl; i++) thm(i) += a(i, n) * r(i, n) * ent(i, n);
  for (int i = 0; i < Nl; i++) thm(i) /= trm(i);
}
 
void Milieu_composite::associer_equation(const Equation_base *eqn) const
{
  Fluide_base::associer_equation(eqn);
  // propagate to underlying media (real fluids need this)
  for (const auto& itr : fluides_) itr->associer_equation(eqn);
}
 
int Milieu_composite::check_unknown_range() const
{
  int ok = 1;
  for (int i = 0; ok && i < (int)fluides_.size(); i++) ok &= fluides_[i]->check_unknown_range(); // each fluid is checked
  return ok;
}
 
void Milieu_composite::calculer_masse_volumique(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Milieu_composite& mil = ref_cast(Milieu_composite, obj);
  const Domaine_VF& zvf = ref_cast(Domaine_VF, mil.equation_.begin()->second->domaine_dis());
  int i, Ni = val.dimension_tot(0), Nb = bval.dimension_tot(0), n, N = (int)mil.fluides_.size();
  std::vector<const DoubleTab *> split(N);
  for (n = 0; n < N; n++) split[n] = &mil.fluides_[n]->masse_volumique().valeurs();
  for (i = 0; i < Ni; i++)
    for (n = 0; n < N; n++) val(i, n) = (*split[n])(i * (split[n]->dimension(0) > 1), 0);
 
  std::vector<DoubleTab> bsplit(N);
  for (n = 0; n < N; n++)
    if (mil.fluides_[n]->masse_volumique().a_un_domaine_dis_base())
      bsplit[n] = mil.fluides_[n]->masse_volumique().valeur_aux_bords();
    else bsplit[n].resize(bval.dimension_tot(0), 1), mil.fluides_[n]->masse_volumique().valeur_aux(zvf.xv_bord(), bsplit[n]);
  for (i = 0; i < Nb; i++)
    for (n = 0; n < N; n++) bval(i, n) = bsplit[n](i * (split[n]->dimension(0) > 1), 0);
 
  /* derivees */
  std::vector<const tabs_t *> split_der(N);
  for (n = 0; n < N; n++) split_der[n] = sub_type(Champ_Inc_base, mil.fluides_[n]->masse_volumique()) ? &ref_cast(Champ_Inc_base, mil.fluides_[n]->masse_volumique()).derivees() : nullptr;
  std::set<std::string> noms_der;
  for (n = 0; n < N; n++)
    if (split_der[n])
      for (auto &&n_d : *split_der[n]) noms_der.insert(n_d.first);
  for (auto &&nom : noms_der)
    {
      for (n = 0; n < N; n++) split[n] = split_der[n] && split_der[n]->count(nom) ? &split_der[n]->at(nom) : nullptr;
      DoubleTab& der = deriv[nom];
      for (der.resize(Ni, N), i = 0; i < Ni; i++)
        for (n = 0; n < N; n++) der(i, n) = split[n] ? (*split[n])(i * (split[n]->dimension(0) > 1)) : 0;
    }
}
 
void Milieu_composite::calculer_energie_interne(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Milieu_composite& mil = ref_cast(Milieu_composite, obj);
  int i, Ni = val.dimension_tot(0), Nb = bval.dimension_tot(0), n, N = (int)mil.fluides_.size();
  std::vector<const DoubleTab *> split(N);
  for (n = 0; n < N; n++) split[n] = &mil.fluides_[n]->energie_interne().valeurs();
  for (i = 0; i < Ni; i++)
    for (n = 0; n < N; n++) val(i, n) = (*split[n])(i * (split[n]->dimension(0) > 1), 0);
 
  std::vector<DoubleTab> bsplit(N);
  for (n = 0; n < N; n++) bsplit[n] = mil.fluides_[n]->energie_interne().valeur_aux_bords();
  for (i = 0; i < Nb; i++)
    for (n = 0; n < N; n++) bval(i, n) = bsplit[n](i * (split[n]->dimension(0) > 1), 0);
 
  /* derivees */
  std::vector<const tabs_t *> split_der(N);
  for (n = 0; n < N; n++) split_der[n] = sub_type(Champ_Inc_base, mil.fluides_[n]->energie_interne()) ? &ref_cast(Champ_Inc_base, mil.fluides_[n]->energie_interne()).derivees() : nullptr;
  std::set<std::string> noms_der;
  for (n = 0; n < N; n++)
    if (split_der[n])
      for (auto &&n_d : *split_der[n]) noms_der.insert(n_d.first);
  for (auto &&nom : noms_der)
    {
      for (n = 0; n < N; n++) split[n] = split_der[n] && split_der[n]->count(nom) ? &split_der[n]->at(nom) : nullptr;
      DoubleTab& der = deriv[nom];
      for (der.resize(Ni, N), i = 0; i < Ni; i++)
        for (n = 0; n < N; n++) der(i, n) = split[n] ? (*split[n])(i * (split[n]->dimension(0) > 1)) : 0;
    }
}
 
void Milieu_composite::calculer_enthalpie(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Milieu_composite& mil = ref_cast(Milieu_composite, obj);
  int i, Ni = val.dimension_tot(0), Nb = bval.dimension_tot(0), n, N = (int)mil.fluides_.size();
  std::vector<const DoubleTab *> split(N);
  for (n = 0; n < N; n++) split[n] = &mil.fluides_[n]->enthalpie().valeurs();
  for (i = 0; i < Ni; i++)
    for (n = 0; n < N; n++) val(i, n) = (*split[n])(i * (split[n]->dimension(0) > 1), 0);
 
  std::vector<DoubleTab> bsplit(N);
  for (n = 0; n < N; n++) bsplit[n] = mil.fluides_[n]->enthalpie().valeur_aux_bords();
  for (i = 0; i < Nb; i++)
    for (n = 0; n < N; n++) bval(i, n) = bsplit[n](i * (split[n]->dimension(0) > 1), 0);
 
  /* derivees */
  std::vector<const tabs_t *> split_der(N);
  for (n = 0; n < N; n++) split_der[n] = sub_type(Champ_Inc_base, mil.fluides_[n]->enthalpie()) ? &ref_cast(Champ_Inc_base, mil.fluides_[n]->enthalpie()).derivees() : nullptr;
  std::set<std::string> noms_der;
  for (n = 0; n < N; n++)
    if (split_der[n])
      for (auto &&n_d : *split_der[n]) noms_der.insert(n_d.first);
  for (auto &&nom : noms_der)
    {
      for (n = 0; n < N; n++) split[n] = split_der[n] && split_der[n]->count(nom) ? &split_der[n]->at(nom) : nullptr;
      DoubleTab& der = deriv[nom];
      for (der.resize(Ni, N), i = 0; i < Ni; i++)
        for (n = 0; n < N; n++) der(i, n) = split[n] ? (*split[n])(i * (split[n]->dimension(0) > 1)) : 0;
    }
}
 
void Milieu_composite::calculer_temperature_multiphase(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Milieu_composite& mil = ref_cast(Milieu_composite, obj);
  int i, Ni = val.dimension_tot(0), Nb = bval.dimension_tot(0), n, N = (int)mil.fluides_.size();
  std::vector<const DoubleTab *> split(N);
  for (n = 0; n < N; n++) split[n] = &mil.fluides_[n]->temperature_multiphase().valeurs();
  for (i = 0; i < Ni; i++)
    for (n = 0; n < N; n++) val(i, n) = (*split[n])(i * (split[n]->dimension(0) > 1), 0);
 
  std::vector<DoubleTab> bsplit(N);
  for (n = 0; n < N; n++) bsplit[n] = mil.fluides_[n]->temperature_multiphase().valeur_aux_bords();
  for (i = 0; i < Nb; i++)
    for (n = 0; n < N; n++) bval(i, n) = bsplit[n](i * (split[n]->dimension(0) > 1), 0);
 
  /* derivees */
  std::vector<const tabs_t *> split_der(N);
  for (n = 0; n < N; n++) split_der[n] = sub_type(Champ_Inc_base, mil.fluides_[n]->temperature_multiphase()) ? &ref_cast(Champ_Inc_base, mil.fluides_[n]->temperature_multiphase()).derivees() : nullptr;
  std::set<std::string> noms_der;
  for (n = 0; n < N; n++)
    if (split_der[n])
      for (auto &&n_d : *split_der[n]) noms_der.insert(n_d.first);
  for (auto &&nom : noms_der)
    {
      for (n = 0; n < N; n++) split[n] = split_der[n] && split_der[n]->count(nom) ? &split_der[n]->at(nom) : nullptr;
      DoubleTab& der = deriv[nom];
      for (der.resize(Ni, N), i = 0; i < Ni; i++)
        for (n = 0; n < N; n++) der(i, n) = split[n] ? (*split[n])(i * (split[n]->dimension(0) > 1)) : 0;
    }
}
 
void Milieu_composite::abortTimeStep()
{
  Fluide_base::abortTimeStep();
  if (ch_e_int_) ch_e_int_->abortTimeStep();
  if (ch_h_ou_T_) ch_h_ou_T_->abortTimeStep();
}
 
bool Milieu_composite::initTimeStep(double dt)
{
  for (auto& itr : fluides_) itr->initTimeStep(dt);
 
  if (!equation_.size()) return true; // no associated equation -> ???
  const Schema_Temps_base& sch = equation_.begin()->second->schema_temps(); // retrieve the time scheme from the first equation
 
  /* fields for which time levels must be created */
  std::vector<Champ_Inc_base *> vch;
  if (ch_rho_ && sub_type(Champ_Inc_base, ch_rho_.valeur()))
    vch.push_back(&ref_cast(Champ_Inc_base, ch_rho_.valeur()));
  if (ch_e_int_ && sub_type(Champ_Inc_base, ch_e_int_.valeur()))
    vch.push_back(&ref_cast(Champ_Inc_base, ch_e_int_.valeur()));
  if (ch_h_ou_T_ && sub_type(Champ_Inc_base, ch_h_ou_T_.valeur()))
    vch.push_back(&ref_cast(Champ_Inc_base, ch_h_ou_T_.valeur()));
 
  for (auto &pch : vch)
    for (int i = 1; i <= sch.nb_valeurs_futures(); i++)
      {
        pch->changer_temps_futur(sch.temps_futur(i), i);
        pch->futur(i) = pch->valeurs();
      }
  return true;
}