Equation_base.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 <Parametre_diffusion_implicite.h>
#include <Solveur_Implicite_base.h>
#include <Schema_Euler_explicite.h>
#include <Schema_Euler_Implicite.h>
#include <Source_dep_inco_base.h>
#include <Operateur_Conv_base.h>
#include <Op_Conv_negligeable.h>
#include <Operateur_Diff_base.h>
#include <Discretisation_base.h>
#include <Matrice_Morse_Diag.h>
#include <Frontiere_dis_base.h>
#include <Matrice_Morse_Sym.h>
#include <Operateur_base.h>
#include <TRUSTTab_parts.h>
#include <EcrFicPartage.h>
#include <Postraitement.h>
#include <Equation_base.h>
#include <Milieu_base.h>
#include <TRUST_2_PDI.h>
#include <Domaine_VF.h>
#include <SolveurSys.h>
#include <Domaine_VF.h>
#include <Operateur.h>
#include <EChaine.h>
#include <Avanc.h>
#include <Debog.h>
#include <Param.h>
#include <Perf_counters.h>
 
// XD condinit objet_lecture nul NO_BRACE Initial condition.
// XD attr nom chaine nom REQ Name of initial condition field.
// XD attr ch field_base ch REQ Type field and the initial values.
 
// XD condinits listobj condinits INHERITS_BRACE condinit NO_COMMA Initial conditions.
 
Implemente_base_sans_constructeur(Equation_base,"Equation_base",Objet_U);
// XD eqn_base mor_eqn eqn_base BRACE Basic class for equations.
/* Attributes further down in the cpp: */
 
Equation_base::Equation_base()
{
  sys_invariant_=1;
  implicite_=-1;
  has_time_factor_= false;
  Nom expr_equation_non_resolue="0";
  equation_non_resolue_.setNbVar(1);
  equation_non_resolue_.setString(expr_equation_non_resolue); // Default value, equation solved
  equation_non_resolue_.addVar("t");
  equation_non_resolue_.parseString();
  set_calculate_time_derivative(0);
}
 
int Equation_base::equation_non_resolue() const
{
  std::string str= equation_non_resolue_.getString();
  std::string str2("INPUT_INT_VALUE");
  if (str.compare(str2)==0)
    {
      Nom name = eq_non_resolue_input_.getName();
      int eq_non_resolue_int = mon_probleme->getOutputIntValue(name);
      return eq_non_resolue_int;
    }
  else
    {
      double t = le_schema_en_temps->temps_courant();
      equation_non_resolue_.setVar("t",t);
      return (int)equation_non_resolue_.eval();
    }
}
 
/*! @brief Returns the discretized domain associated with the equation.
 *
 * @return (Domaine_dis_base&) the discretized domain associated with the equation
 * @throws the discretized domain object (Domaine_dis) is invalid,
 * associated problem not discretized.
 */
Domaine_dis_base& Equation_base::domaine_dis()
{
  if (!le_dom_dis)
    {
      Cerr << "There is no object of type Domaine_dis yet associated to the equation " << que_suis_je() << finl;
      Cerr << "This means that the problem has not been discretized or" << finl;
      Cerr << "that instruction Discretiser is misplaced" << finl;
      exit();
    }
  return le_dom_dis.valeur();
}
 
/*! @brief Returns the discretized domain associated with the equation.
 *
 * (const version)
 *
 * @return (Domaine_dis_base&) the discretized domain associated with the equation
 * @throws the discretized domain object (Domaine_dis) is invalid,
 * associated problem not discretized.
 */
const Domaine_dis_base& Equation_base::domaine_dis() const
{
  if (!le_dom_dis)
    {
      Cerr << "There is no object of type Domaine_dis associated to the equation " << que_suis_je() << finl;
      Cerr << "This means that the problem has not been discretized" << finl;
      exit();
    }
  assert(le_dom_dis);
  return le_dom_dis.valeur();
}
 
/*! @brief Completes the construction (initialization) of objects associated with the equation.
 *
 *     Completes the sources, associates the equation with the unknown, completes the
 *     operators, completes the discretized boundary conditions.
 *     See the methods Source_base::completer(),
 *                       Operateur_base::completer()
 *                       Domaine_Cl_dis_base::completer()
 *                       Domaine_Cl_dis_base::completer(const Domaine_dis_base& )
 *
 */
void Equation_base::completer()
{
  inconnue().associer_eqn(*this);
  if (le_dom_Cl_dis)
    le_dom_Cl_dis->completer();
 
  inconnue().associer_domaine_cl_dis(le_dom_Cl_dis);
 
  if (mon_probleme->is_coupled())
    {
      bool err = false;
      int nb_cl = le_dom_Cl_dis->nb_cond_lim();
      for (int i=0; i<nb_cl; i++)
        {
          const Cond_lim_base& la_cl = le_dom_Cl_dis->les_conditions_limites(i);
          const Frontiere& la_frontiere = la_cl.frontiere_dis().frontiere();
 
          bool raccord_found = false;
 
          // VDF and PolyMAC_HFV
          if (la_cl.que_suis_je().debute_par("Paroi_Echange_contact") || la_cl.que_suis_je().debute_par("Echange_contact_Rayo"))
            raccord_found = true;
 
          // VEF
          if (la_cl.que_suis_je().debute_par("Scalaire_impose_paroi") || la_cl.que_suis_je().debute_par("paroi_temperature_imposee_rayo") || la_cl.que_suis_je().debute_par("temperature_imposee_paroi"))
            raccord_found = true;
 
          if ( raccord_found && !sub_type(Raccord_base,la_frontiere))
            {
              Cerr << "====================================================================" << finl;
              Cerr << "Boundary " << la_frontiere.le_nom() << " should be a raccord" << finl;
              Cerr << "Add in your data file between the definition and the partition of the domain " << mon_probleme->domaine().le_nom() << " : " << finl;
              Cerr << "Modif_bord_to_raccord " << mon_probleme->domaine().le_nom() << " " << la_frontiere.le_nom() << finl;
              err = true;
            }
        }
      if(err)
        exit();
    }
 
  //Addition of an empty element that will be returned if no model is found
  RefObjU mod;
  if (liste_modeles_.size()==0)
    liste_modeles_.add(mod);
 
  // for equations not calling preparer_calcul
  initialise_residu();
  int nb_op=nombre_d_operateurs();
  Nom msg="Equation_base::completer(), nb_op = " ;
  Debog::verifier(msg, nb_op);
  for(int i=0; i<nb_op; i++)
    operateur(i).completer();
 
  if (solveur_masse)  // [ABN]: In case of Front-Tracking, mass solver mass might be uninitialized ...
    solv_masse().completer();
 
  les_sources.completer();
  schema_temps().completer();
 
  // Determine implicite_ and sys_invariant_ variables:
  implicite_=0;
  for(int i=0; i<nombre_d_operateurs(); i++)
    {
      if(operateur(i).l_op_base().get_decal_temps()==1)
        {
          ++implicite_;
          sys_invariant_*=operateur(i).l_op_base().systeme_invariant();
        }
    }
}
 
 
/*! @brief Overrides Objet_U::printOn: prints the equation and its components to an output stream.
 *
 *     Prints the equation name, the mass solver, the source terms,
 *     the discretized boundary conditions, the unknowns, and the operators.
 *
 * @param (Sortie& os) an output stream
 * @return (Sortie&) the modified output stream
 */
Sortie& Equation_base::printOn(Sortie& os) const
{
  os << nom_ << finl;
  os << solveur_masse;
  os << les_sources;
  os << le_dom_Cl_dis;
  os << inconnue() ;
  for(int i=0; i<nombre_d_operateurs(); i++)
    {
      os << operateur(i).l_op_base();
    }
  return os;
}
 
/*! @brief Reads an equation from an input stream.
 *
 * The expected format is:
 *       {
 *         [Source { [sou_1], [sour_2], ...} ]
 *          Conditions_limites { [cl_1] [cl_2] ... }
 *          Conditions_initiales { [cl_init] }
 *       }
 *
 * @param (Entree& is) an input stream to read the equation from
 * @return (Entree&) the modified input stream
 * @throws bad read format, opening brace expected
 * @throws no initial conditions, you must initialize the unknown
 * @throws no boundary conditions, you must provide boundary conditions
 */
Entree& Equation_base::readOn(Entree& is)
{
  Cerr<<"Reading of data for a "<<que_suis_je()<<" equation"<<finl;
  Param param(que_suis_je());
  set_param(param);
  Nom expr_equation_non_resolue="0";
  param.ajouter("disable_equation_residual",&disable_equation_residual_); // XD attr disable_equation_residual entier disable_equation_residual OPT The equation residual will not be used for the problem residual used when checking time convergence or computing dynamic time-step
  equation_non_resolue_.setString(expr_equation_non_resolue);
  param.lire_avec_accolades_depuis(is);
  matrice_init = 0;
  return is;
}
 
// XD attr convection bloc_convection convection OPT Keyword to alter the convection scheme.
// XD attr diffusion bloc_diffusion diffusion OPT Keyword to specify the diffusion operator.
void Equation_base::set_param(Param& param) const
{
  param.ajouter_non_std("conditions_limites|boundary_conditions",(this),Param::REQUIRED);  // XD attr conditions_limites|boundary_conditions condlims conditions_limites OPT Boundary conditions.
  param.ajouter_non_std("conditions_initiales|initial_conditions",(this),Param::REQUIRED); // XD attr conditions_initiales|initial_conditions condinits conditions_initiales OPT Initial conditions.
  param.ajouter_non_std("sources",(this)); // XD attr sources sources sources OPT To introduce a source term into an
  // XD_CONT equation (in case of several source terms into the same equation, the blocks corresponding to the various
  // XD_CONT terms need to be separated by a comma)
  param.ajouter_non_std("ecrire_fichier_xyz_valeur",(this)); // XD attr ecrire_fichier_xyz_valeur ecrire_fichier_xyz_valeur ecrire_fichier_xyz_valeur OPT This keyword is used to write the values of a field only for some boundaries in a text file
  param.ajouter("parametre_equation",&parametre_equation_); // XD attr parametre_equation parametre_equation_base parametre_equation OPT Keyword used to specify additional parameters for the equation
  param.ajouter_non_std("equation_non_resolue",(this)); // XD attr equation_non_resolue chaine equation_non_resolue OPT
  // XD_CONT The equation will not be solved while condition(t) is verified if equation_non_resolue keyword is used.
  // XD_CONT Exemple: The Navier-Stokes equations are not solved between time t0 and t1. NL2 Navier_Sokes_Standard NL2 {
  // XD_CONT equation_non_resolue (t>t0)*(t<t1) }
  param.ajouter_non_std("rename_equation|renommer_equation",(this)); // XD attr renommer_equation chaine rename_equation OPT Rename the equation with a specific name.
}
 
int Equation_base::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  if (mot=="sources")
    {
      lire_sources(is);
      return 1;
    }
  else if (mot=="conditions_limites|boundary_conditions")
    {
      lire_cl(is);
      verif_Cl();
      return 1;
    }
  else if (mot=="conditions_initiales|initial_conditions")
    {
      lire_cond_init(is);
      return 1;
    }
  else if (mot=="ecrire_fichier_xyz_valeur")
    {
      xyz_field_values_file_.associer_eqn(*this);
      is >> xyz_field_values_file_;
      return 1;
    }
  else if (mot=="rename_equation|renommer_equation")
    {
      Cerr << "Equation " << nom_ << " renamed to : ";
      is >> nom_;
      Cerr << nom_ << finl;
      return 1;
    }
  else if (mot=="equation_non_resolue")
    {
      // differentiate an entry with an expression from an entry with the keyword eq_non_resoue_input_
      // On complete:
      Motcle expr_equation_non_resolue;
      is >> expr_equation_non_resolue;
      equation_non_resolue_.setString(expr_equation_non_resolue);
      if (expr_equation_non_resolue=="INPUT_INT_VALUE")
        is >> eq_non_resolue_input_;
      else
        equation_non_resolue_.parseString();
      return 1;
    }
  return -1;
}
 
 
/*! @brief Reading of source terms in an input stream.
 *
 * @param (Entree& is) input stream to read source terms
 * @return (Entree&) the modified input stream
 */
Entree& Equation_base::lire_sources(Entree& is)
{
  static bool already_read=false;
  if (already_read)
    {
      Cerr << "Error: the 'sources' keyword should only appear once per equation." << finl;
      Cerr << "To declare several source terms, list them inside a single block separated by commas:" << finl;
      Cerr << "    sources { TermA { ... } , TermB { ... } }" << finl;
      // Process::exit(); // comment that out since it appears pb_multiphase is doing shadowy stuff...
    }
 
  already_read=true;
  Cerr << "Reading of source terms" << finl ;
  is >> les_sources;
  return is;
}
 
/*! @brief Returns the source terms associated with the equation
 *
 * @return (Sources&) the list of source terms associated with the equation
 */
Sources& Equation_base::sources()
{
  return les_sources;
}
 
/*! @brief Returns the source terms associated with the equation (const version)
 *
 * @return (Sources&) the list of source terms associated with the equation
 */
const Sources& Equation_base::sources() const
{
  return les_sources;
}
 
/*! @brief Reading of initial conditions in an input stream.
 *
 * The reading format is as follows:
 *      {
 *       Name [MUST BE THE NAME OF THE UNKNOWN]
 *       [READ A GIVEN FIELD]
 *      }
 *
 * @param (Entree& is) input stream
 * @return (Entree&) the modified input stream
 * @throws format error, opening brace expected
 * @throws bad name for the unknown
 * @throws format error, closing brace expected
 */
Entree& Equation_base::lire_cond_init(Entree& is)
{
  Cerr << "Reading of initial conditions\n";
  Nom nom;
  Motcle motlu;
  is >> nom;
  motlu = nom;
  if(motlu!=Motcle("{"))
    {
      Cerr << "We expected a { while reading " << que_suis_je() << finl;
      Cerr << "and not : " << nom << finl;
      exit();
    }
  is >> nom;
  motlu = nom;
  if (motlu != Motcle(inconnue().le_nom()))
    {
      Cerr << nom << " is not the name of the unknown "
           << inconnue().le_nom() << finl;
      exit();
    }
  OWN_PTR(Champ_Don_base) ch_init;
  is >> ch_init;
  const int nb_comp = ch_init->nb_comp();
  verifie_ch_init_nb_comp(inconnue(),nb_comp);
 
  //Cerr<<"inconnue().que_suis_je() = "<<inconnue().que_suis_je()<<finl;
 
  inconnue().affecter(ch_init.valeur());
  is >> nom;
  motlu = nom;
  if(motlu!=Motcle("}"))
    {
      Cerr << "We expected a } while reading " << que_suis_je() << finl;
      Cerr << "and not : " << nom << finl;
      exit();
    }
  return is;
}
 
/*! @brief Reading of boundary conditions in an input stream.
 *
 * see Domaine_Cl_dis_base::readOn
 *
 * @param (Entree& is) input stream
 * @return (Entree&) the modified input stream
 * @throws the discretized boundary condition domain is empty
 */
Entree& Equation_base::lire_cl(Entree& is)
{
  Cerr << "Reading of boundaries conditions\n";
  if(!le_dom_Cl_dis)
    {
      Cerr << "Error while reading boundaries conditions : " <<
           que_suis_je() << finl;
      Cerr << "The Domaine_Cl_dis_base is empty ..." << finl;
      exit();
    }
  is >> le_dom_Cl_dis.valeur();
  return is;
}
 
/*! @brief We save the unknown, then the source terms to an output stream.
 *
 * @param (Sortie& os)
 * @return (int) the return code of Champ_Inc::sauvegarder()
 */
int Equation_base::sauvegarder(Sortie& os) const
{
  int ret = inconnue().sauvegarder(os);
  for (int i = 0; i < les_sources.size(); i++)
    ret += les_sources(i)->sauvegarder(os);
  return ret;
}
 
 
/*! @brief for PDI IO: retrieve name, type and dimensions of the data to save/restore.
 * This has to be overrode for all the equations that either:
 * - have extra fields (ie in addition to the unknown) or extra scalars to save/restore.
 * - want to save the unknown but with a different name
 * These data will then be written in a YAML file, to initialize PDI.
 * They have to be shared with PDI afterwards, when they need to be read/written (via TRUST_2_PDI)
 *
 */
std::vector<YAML_data> Equation_base::data_a_sauvegarder() const
{
  std::vector<YAML_data> data = inconnue().data_a_sauvegarder();
  return data;
}
 
/*! @brief We resume the unknown from an input stream.
 *
 * [WE SEARCH FOR THE UNKNOWN BY ITS NAME]
 *      [WE READ THE UNKNOWN]
 *      See Champ_Inc::reprendre(Entree&)
 *
 * @param (Entree& fich) input stream (file) to read
 * @return (int) always returns 1
 * @throws restart error, end of file reached without finding the unknown
 */
int Equation_base::reprendre(Entree& fich)
{
  if(!TRUST_2_PDI::is_PDI_restart())
    {
      double temps = schema_temps().temps_courant();
      Nom field_tag(inconnue().le_nom());
      field_tag += inconnue().que_suis_je();
      field_tag += probleme().domaine().le_nom();
      field_tag += Nom(temps,probleme().reprise_format_temps());
      // This part ensure the backwards compatibility of TRUST regarding the renaming of the PolyMAC family of discretisation. It will be deprecated in the TRUST 2.0
 
      if (probleme().discretisation().is_poly_family())
        {
          Nom field_tag_syno = create_polymacfamily_syno(field_tag);
          avancer_fichier_with_syno(fich,field_tag,field_tag_syno);
        }
      // end of the backward compatibility
      else
        avancer_fichier(fich,field_tag);
    }
  int ret = inconnue().reprendre(fich);
  for (int i = 0; i < les_sources.size(); i++)
    ret = ret && les_sources(i)->reprendre(fich);
  return ret;
}
 
/*! @brief Create a synonym of a field name in order to ensure backward compatibility with old names of the PolyMAC discretisation family
 *
 * @param field_tag
 * @return synonym of the field tag
 */
Nom Equation_base::create_polymacfamily_syno(const Nom& field_tag) const
{
  Nom field_tag_syno=field_tag;
  auto create_syno= [field_tag](std::string pattern, std::string replace)
  {
    std::string s = field_tag.getString();
    std::size_t pos = s.find(pattern);
    while (pos != std::string::npos)
      {
        s.replace(pos, pattern.length(), replace);
        pos = s.find(pattern, pos + replace.length());
      }
    return Nom(s);
  };
 
  if (probleme().discretisation().is_PolyMAC_CDO())
    field_tag_syno=create_syno("PolyMAC_CDO","PolyMAC");
  else if (probleme().discretisation().is_PolyMAC_MPFA())
    field_tag_syno=create_syno("PolyMAC_MPFA","PolyMAC_P0");
  else if (probleme().discretisation().is_PolyMAC_HFV())
    field_tag_syno=create_syno("PolyMAC_HFV","PolyMAC_P0P1NC");
  return field_tag_syno;
}
/*! @brief Ask the time scheme if an output is needed.
 *
 * Returns 1 if an output is needed.
 *      Simple call to Schema_Temps_base::limpr()
 *
 * @return (int) 1 if an output is needed, 0 otherwise.
 */
int Equation_base::limpr() const
{
  return le_schema_en_temps->limpr();
}
 
/*! @brief Prints the equation operators to an output stream, unconditionally.
 *
 *     calls Operateur_base::impr(os)
 *
 * @param (Sortie& os) output stream
 * @return (int) always returns 1
 */
int Equation_base::impr(Sortie& os) const
{
  for(int i=0; i<nombre_d_operateurs(); i++)
    operateur(i).impr(os);
  //if (je_suis_maitre() && les_sources.size()>0) os << finl << "Printing source terms for equation " << que_suis_je() << " : " << finl;
  les_sources.impr(os);
  return 1;
}
 
/*! @brief Prints the equation operators if the time scheme indicates it is necessary.
 *
 *      [IF limpr() THEN impr(os)]
 *
 * @param (Sortie& os) output stream
 */
void Equation_base::imprimer(Sortie& os) const
{
  xyz_field_values_file_.write_fields();
  if (limpr() )
    impr(os);
}
 
/*! @brief Returns the time derivative of the unknown I of the equation: dI/dt = M-1*(sum(operators(I) + sources))
 *
 *  In : derivee contains I (but immediatly set to 0)
 *  Out: derivee contains dI/dt
 *
 */
DoubleTab& Equation_base::derivee_en_temps_inco(DoubleTab& derivee)
{
  derivee = 0.;
 
  DoubleTrav secmem(derivee);
  // secmem = sum(operators) + sources + equation specific terms
 
 
  derivee_en_temps_inco_sources(secmem);
  const double time_factor = get_time_factor();
 
  bool calcul_explicite = false;
  if (parametre_equation_ && sub_type(Parametre_implicite, parametre_equation_.valeur()))
    {
      Parametre_implicite& param2 = ref_cast(Parametre_implicite, parametre_equation_.valeur());
      calcul_explicite = param2.calcul_explicite();
    }
 
  if (schema_temps().diffusion_implicite() && !calcul_explicite)
    {
      // Add convection operator only if equation has one
      derivee=inconnue().valeurs();
      if (nombre_d_operateurs()>1)
        {
          derivee_en_temps_conv(secmem, derivee);
          if (has_time_factor_)
            {
              secmem *= time_factor;
            }
        }
    }
  else
    {
      verify_scheme();
 
      // Add all explicit operators
      for(int i=0; i<nombre_d_operateurs(); i++)
        if(operateur(i).l_op_base().get_decal_temps()!=1)
          {
            if (has_time_factor_)
              {
                DoubleTrav secmem_tmp(secmem);
                operateur(i).ajouter(secmem_tmp);
                if (i == 1) secmem_tmp *= time_factor;
                secmem += secmem_tmp;
              }
            else
              {
                operateur(i).ajouter(secmem);
              }
          }
    }
  les_sources.ajouter(secmem);
  if (calculate_time_derivative())
    {
      // Store dI/dt(n) = M-1 secmem :
      derivee_en_temps().valeurs()=secmem;
      solveur_masse->appliquer(derivee_en_temps().valeurs());
      schema_temps().modifier_second_membre((*this),secmem); // Change secmem for some schemes (eg: Adams_Bashforth)
    }
 
  corriger_derivee_expl(secmem); // Add specific term for an equation (eg: -gradP for Navier Stokes)
 
  if (implicite_==0)
    {
      solveur_masse->appliquer(secmem); // M-1 * secmem
      if (schema_temps().diffusion_implicite() && !calcul_explicite)
        {
          // Solve: (1/dt + M-1*A)*dI = M-1 * secmem
          // where A is the diffusion
          Equation_base::Gradient_conjugue_diff_impl(secmem, derivee);
        }
      else
        {
          derivee = secmem;
          derivee.echange_espace_virtuel();
        }
      corriger_derivee_impl(derivee);  // Solve specific implicit term for an equation (eg: pressure for Navier Stokes)
 
    }
  else if (implicite_>0)
    {
      // TRUST support notices that this part has never been covered...
      //implicit
      // M dU/dt + AU* = f -BUn;
      // U* = Un+dt dU/dt
      // (M/dt + A) U* = f -BUn + M/dt Un
      //
      double dt=schema_temps().pas_de_temps();
      for(int i=0; i<nombre_d_operateurs(); i++)
        {
          //loop over operators
          Operateur_base& op=operateur(i).l_op_base();
          if(!op.get_matrice())
            op.set_matrice().typer("Matrice_Morse");
          if(op.get_decal_temps()==1)
            {
              //if (op.set_matrice()->nb_lignes()<2)
              {
                Matrice_Morse& matrice=ref_cast(Matrice_Morse,op.set_matrice().valeur());
                op.dimensionner(matrice);
                sys_invariant_=0;
              }
              if(!sys_invariant_)
                {
                  Matrice_Morse& matrice=ref_cast(Matrice_Morse, op.set_matrice().valeur());
                  op.contribuer_a_avec(inconnue().valeurs(), matrice);
                  solv_masse().ajouter_masse(dt, op.set_matrice().valeur());
 
                  if(
                    (op.get_solveur()->que_suis_je()=="Solv_Cholesky")
                    ||
                    (op.get_solveur()->que_suis_je()=="Solv_GCP")
                  )
                    {
                      Matrice_Morse_Sym new_mat(matrice);
                      new_mat.set_est_definie(1);
                      op.set_matrice()=new_mat;
                      ref_cast_non_const(SolveurSys,op.get_solveur())->reinit();
                    }
                }
            }
        }
      if(implicite_==1)
        {
          // Single implicit operator.
          // Assumed to be the first one (diffusion!!)
          Operateur_base& op=operateur(0).l_op_base();
          Matrice_Base& matrice=op.set_matrice().valeur();
          // DoubleTrav secmem(derivee);
          secmem=derivee;
          solv_masse().ajouter_masse(dt, secmem, inconnue().valeurs());
          op.contribuer_au_second_membre(secmem );
          op.set_solveur().resoudre_systeme(matrice,
                                            secmem,
                                            derivee
                                           );
          solv_masse().corriger_solution(derivee,inconnue().valeurs());
 
          derivee-=inconnue().valeurs();
          derivee/=dt;
 
          //Used only to compute boundary fluxes when diffusion is made implicit!
          DoubleTab resu;
          resu=derivee;
          operateur(0).calculer(inconnue().valeurs(), resu);
        }
      else
        {
          // multiple implicit operators ...
          Cerr << "Must be coded ... " << finl;
          exit();
        }
    }
  else
    {
      Cerr << "Error in Equation_base::derivee_en_temps_inco" << finl;
      Cerr << "implicite_ = " << implicite_ << " has not been initialized!" << finl;
      Cerr << "May be " << que_suis_je() << "::completer() method doesn't call Equation_base::completer()" << finl;
      Cerr << "Contact TRUST support." << finl;
      Process::exit();
    }
 
  return derivee;
}
 
/*! @brief Returns the problem associated with the equation.
 *
 * @return (Probleme_base&) the problem associated with the equation
 */
Probleme_base& Equation_base::probleme()
{
  return mon_probleme.valeur();
}
 
/*! @brief Returns the problem associated with the equation.
 *
 * (const version)
 *
 * @return (Probleme_base&) the problem associated with the equation
 */
const Probleme_base& Equation_base::probleme() const
{
  return mon_probleme.valeur();
}
 
/*! @brief Method called when the object instance is created in the data file (Interprete::ajouter)
 *
 */
void Equation_base::nommer(const Nom& un_nom)
{
  nom_ = un_nom;
}
 
/*! @brief Associates with the Problem passed as parameter.
 *
 * Also associates the sources, operators and mass
 *     solver with the equation.
 *
 * @param (Probleme_base& pb) the problem with which the equation must be associated
 */
void Equation_base::associer_pb_base(const Probleme_base& pb)
{
  mon_probleme=pb;
  // Modif B. Mathieu: for front_tracking, the equation name
  // is given in the data file during object instantiation
  // (the "nommer" method is called at that moment).
  // If the equation is not instantiated, the name is empty, we give
  // the default name which is defined below:
  Nom nom_vide;
  if (nom_ == nom_vide)
    {
      nom_ = pb.le_nom();
      nom_ += que_suis_je();
    }
  // end modif B.M.
  les_sources.associer_eqn(*this);
  int nb_op = nombre_d_operateurs();
  for(int i=0; i<nb_op; i++)
    operateur(i).associer_eqn(*this);
 
}
 
/*! @brief Discretizes the equation.
 *
 * Type the domaine_Cl_dis, the format, associate it with the equation.
 *     Type the mass solver, associate the discretized domain and
 *     the discretized boundary condition domain with it.
 *
 */
void Equation_base::discretiser()
{
  Cerr << "Discretizing of the boundary conditions... ";
  le_dom_Cl_dis.typer(discretisation().domaine_cl_dis_type());
  le_dom_Cl_dis->associer(domaine_dis());
  Cerr << " OK " << finl;
 
  le_dom_Cl_dis->associer_eqn(*this);
  le_dom_Cl_dis->associer_inconnue(inconnue());
 
  /*
   * XXX : Elie Saikali : to type correctly the mass solver ...
   */
  Nom typ = discretisation().get_name_of_type_for("Solveur_Masse", "??" /* rien */, *this);
  solveur_masse.typer(typ);
  solveur_masse->associer_eqn(*this);
  solveur_masse->associer_domaine_dis_base(domaine_dis());
  solveur_masse->associer_domaine_cl_dis_base(le_dom_Cl_dis.valeur());
 
  if (calculate_time_derivative())
    {
      Motcle directive("temperature");
      if (inconnue().is_vectorial())
        directive="vitesse";
 
      Nom nom("derivee_en_temps_");
      nom += inconnue().le_nom();
 
      Nom unite(inconnue().unites()[0]);
      unite += "/s";
 
      discretisation().discretiser_champ(directive,domaine_dis(),nom,unite,
                                         inconnue().nb_comp(),
                                         schema_temps().nb_valeurs_temporelles(),
                                         schema_temps().temps_courant(),derivee_en_temps_);
 
      champs_compris_.ajoute_champ(derivee_en_temps_);
    }
}
void Equation_base::associer_milieu_equation()
{
  milieu().associer_equation(this);
}
 
/*! @brief Associates the time scheme with the equation.
 *
 * @param (Schema_Temps_base& un_schema_en_temps) the time scheme to associate with the equation
 */
void Equation_base::associer_sch_tps_base(const Schema_Temps_base& un_schema_en_temps)
{
  le_schema_en_temps=un_schema_en_temps;
}
 
/*! @brief Returns the time scheme associated with the equation.
 *
 * @return (Schema_Temps_base&) the time scheme associated with the equation
 * @throws no time scheme associated with the equation
 */
Schema_Temps_base& Equation_base::schema_temps()
{
  if(!le_schema_en_temps)
    {
      Cerr << "Error : " << que_suis_je()
           << "has not been associated to a time scheme " << finl;
      exit();
    }
  return le_schema_en_temps.valeur();
}
 
/*! @brief Returns the time scheme associated with the equation.
 *
 * @return (Schema_Temps_base&) the time scheme associated with the equation
 * @throws no time scheme associated with the equation
 */
const Schema_Temps_base& Equation_base::schema_temps() const
{
  if(!le_schema_en_temps)
    {
      Cerr << "Error : " << que_suis_je()
           << "has not been associated to a time scheme " << finl;
      exit();
    }
  return le_schema_en_temps.valeur();
}
 
/*! @brief Associates the discretized domain with the equation.
 *
 * @param (Domaine_dis_base& z) the discretized domain to associate
 */
void Equation_base::associer_domaine_dis(const Domaine_dis_base& z)
{
  le_dom_dis=z;
}
 
/*! @brief The value of the unknown at the time step has been calculated.
 *
 * This method advances the present until the time passed as parameter.
 *    It also updates the properties of the medium.
 *
 * @param (double temps) the time step for update
 */
void Equation_base::mettre_a_jour(double temps)
{
  // We turn the unknown wheel
  // Update the unknown:
  inconnue().mettre_a_jour(temps);
  if (calculate_time_derivative()) derivee_en_temps().mettre_a_jour(temps);
 
  les_sources.mettre_a_jour(temps);
 
  // We turn the boundary condition wheel
  // Update the boundary condition:
  if (temps > schema_temps().temps_courant()) domaine_Cl_dis().avancer(temps);
}
 
//update of champ_conserve / champ_convecte : called by Probleme_base::mettre_a_jour() after updating the medium
void Equation_base::mettre_a_jour_champs_conserves(double temps, int reset)
{
  if (reset && champ_conserve_) champ_conserve_->reset_champ_calcule(); //force calculation of all cells
  if (reset && champ_convecte_) champ_convecte_->reset_champ_calcule();
  if (champ_conserve_) champ_conserve_->mettre_a_jour(temps);
  if (champ_convecte_) champ_convecte_->mettre_a_jour(temps);
}
 
/*! @brief Reinitialize what must be.
 *
 * This method is called when a time step is abandoned,
 *     for example because the calculation has diverged.
 *     So you must clean up what could pollute the new resolution.
 *
 * @throws WrongContext
 */
void Equation_base::abortTimeStep()
{
  for (int i=0; i<nombre_d_operateurs(); i++)
    operateur(i).l_op_base().abortTimeStep();
  inconnue().abortTimeStep();
  if (champ_conserve_) champ_conserve_->abortTimeStep();
  if (champ_convecte_) champ_convecte_->abortTimeStep();
}
 
/*! @brief Reset current time of the equation. Used from ICoCo.
 * See documentation of Problem_base::resetTime()
 */
void Equation_base::resetTime(double time)
{
  if(solveur_masse) solveur_masse->resetTime(time);
  les_sources.resetTime(time);
  le_dom_Cl_dis->resetTime(time);
  for (int i=0; i<nombre_d_operateurs(); i++)
    operateur(i).l_op_base().resetTime(time);
  inconnue().resetTime(time);
  if (champ_conserve_) champ_conserve_->resetTime(time);
  if (champ_convecte_) champ_convecte_->resetTime(time);
}
 
 
/*! @brief virtual method to correct the unknown during implicit iterations for example K-eps must remain positive
 *
 *  mass fractions between 0 and 1
 *
 */
void Equation_base::valider_iteration()
{
}
 
/*! @brief Everything that does not depend on other possible problems.
 *
 * @return (int) always returns 1
 */
int Equation_base::preparer_calcul()
{
  solveur_masse->preparer_calcul();
  int nb_op=nombre_d_operateurs();
  for(int i=0; i<nb_op; i++)
    {
      operateur(i).l_op_base().preparer_calcul();
      dt_op_bak.add(0.);
    }
  initialise_residu();
  double temps=schema_temps().temps_courant();
  inconnue().verifie_valeurs_cl();
  inconnue().changer_temps(temps);
  if (calculate_time_derivative()) derivee_en_temps().changer_temps(temps);
  le_dom_Cl_dis->initialiser(temps);
 
  Nom msg=que_suis_je();
  msg+=" dans Equation_base::preparer_calcul 1";
  Debog::verifier(msg ,inconnue());
 
  les_sources.initialiser(temps);
  domaine_Cl_dis().imposer_cond_lim(inconnue(),temps);
  inconnue().valeurs().echange_espace_virtuel();
 
  /* initialization of parametre_equation() by the time scheme if it allows it */
  if (sub_type(Schema_Implicite_base, schema_temps()))
    ref_cast(Schema_Implicite_base, schema_temps()).solveur()->get_and_set_parametre_equation(*this);
 
  msg=que_suis_je();
  msg+=" dans Equation_base::preparer_calcul ";
  Debog::verifier(msg ,inconnue());
 
  return 1;
}
 
 
/*! @brief Allocation and initialization of the unknown and boundary conditions until present+dt.
 *
 *   + other initializations for calculations on the next
 *     time step : operators, mass solver->
 *
 * @return (0 in case of error, 1 otherwise)
 */
bool Equation_base::initTimeStep(double dt)
{
  Schema_Temps_base& sch=schema_temps();
  double temps_present=sch.temps_courant();
 
  // The time step must be that of the scheme !!
  assert(sch.pas_de_temps()==dt);
 
  // For each future time
  for (int i=1; i<=sch.nb_valeurs_futures(); i++)
    {
      double tps=sch.temps_futur(i);
      // Update of the time in the unknown
      inconnue().changer_temps_futur(tps,i);
      if (champ_conserve_) champ_conserve_->changer_temps_futur(tps,i);
      if (champ_convecte_) champ_convecte_->changer_temps_futur(tps,i);
      if (calculate_time_derivative()) derivee_en_temps().changer_temps_futur(tps,i);
 
      inconnue().futur(i)=inconnue().valeurs();
      if (champ_conserve_) champ_conserve_->futur(i) = champ_conserve().valeurs();
      if (champ_convecte_) champ_convecte_->futur(i) = champ_convecte().valeurs();
      if (calculate_time_derivative()) derivee_en_temps().futur(i)=derivee_en_temps().valeurs();
 
      // Update of the time in the boundary conditions
      domaine_Cl_dis().changer_temps_futur(tps,i);
    }
 
  // Update of the default time of the boundary conditions
  domaine_Cl_dis().set_temps_defaut(sch.temps_defaut());
 
  // Update of the operators
  for(int i=0; i<nombre_d_operateurs(); i++)
    operateur(i).mettre_a_jour(temps_present);
 
  // Update of the mass solver at the present time
  if (solveur_masse)
    solveur_masse->mettre_a_jour(temps_present);
 
  return true;
}
 
bool Equation_base::updateGivenFields()
{
  Schema_Temps_base& sch=schema_temps();
  double temps_present=sch.temps_courant();
  double temps_futur=temps_present+sch.pas_de_temps();
 
  // For each future time
  for (int i=1; i<=sch.nb_valeurs_futures(); i++)
    {
      double tps=sch.temps_futur(i);
      // Calculation of boundary conditions at this time.
      domaine_Cl_dis().mettre_a_jour(tps);
    }
  // Calculation of the rate of change of boundary conditions between present and future times.
  domaine_Cl_dis().calculer_derivee_en_temps(temps_present,temps_futur);
 
  //Update of operators
  for (int i = 0; i < nombre_d_operateurs(); i++)
    if (operateur(i).op_non_nul())
      operateur(i).l_op_base().mettre_a_jour(temps_present);
 
  // Update of source terms at the present time
  les_sources.mettre_a_jour(temps_present);
 
  return true;
}
 
/*! @brief Returns the discretization associated with the equation.
 *
 * @return (Discretisation_base&) a discretization associated with the equation
 * @throws no problem associated
 */
const Discretisation_base& Equation_base::discretisation() const
{
  if(!mon_probleme)
    {
      Cerr << "Error : " << que_suis_je() << " has not been associated to a problem ! " << finl;
      exit();
    }
  return mon_probleme->discretisation();
}
 
void Equation_base::creer_champ(const Motcle& motlu)
{
  const Equation_base& me_const = (*this); // to get a const equation !!!!!
  const Nom& nom_inco = me_const.inconnue().le_nom();
  Nom inco(nom_inco);
  inco += "_residu";
  if (motlu == Motcle(inco))
    if (!field_residu_)
      {
        discretisation().residu(domaine_dis(), inconnue(), field_residu_);
        champs_compris_.ajoute_champ(field_residu_);
      }
 
  for (int i = 0; i < nombre_d_operateurs(); i++)
    if (operateur(i).op_non_nul())
      operateur(i).l_op_base().creer_champ(motlu);
 
  for (auto &itr : les_sources)
    if (itr)
      itr->creer_champ(motlu);
}
 
bool Equation_base::has_champ(const Motcle& nom, OBS_PTR(Champ_base) &ref_champ) const
{
  Nom inco_residu(inconnue().le_nom());
  inco_residu += "_residu";
 
  if (nom == Motcle(inco_residu))
    {
      ref_champ = Equation_base::get_champ(nom);
      return true;
    }
 
  for (const auto &itr : list_champ_combi)
    if (itr->le_nom() == nom)
      {
        ref_champ = Equation_base::get_champ(nom);
        return true;
      }
 
  if (champs_compris_.has_champ(nom, ref_champ))
    return true;
 
  if (milieu().has_champ(nom, ref_champ))
    return true;
 
  for (int i = 0; i < nombre_d_operateurs(); i++)
    if (operateur(i).op_non_nul())
      if (operateur(i).l_op_base().has_champ(nom, ref_champ))
        return true;
 
  for (const auto &itr : les_sources)
    if (itr)
      if (itr->has_champ(nom, ref_champ))
        return true;
 
  return false; /* nothing found */
}
 
bool Equation_base::has_champ(const Motcle& nom) const
{
  Nom inco_residu(inconnue().le_nom());
  inco_residu += "_residu";
 
  if (nom == Motcle(inco_residu))
    return true;
 
  for (const auto &itr : list_champ_combi)
    if (itr->le_nom() == nom)
      return true;
 
  if (champs_compris_.has_champ(nom))
    return true;
 
  if (milieu().has_champ(nom))
    return true;
 
  for (int i = 0; i < nombre_d_operateurs(); i++)
    if (operateur(i).op_non_nul())
      if (operateur(i).l_op_base().has_champ(nom))
        return true;
 
  for (const auto &itr : les_sources)
    if (itr)
      if (itr->has_champ(nom))
        return true;
 
  return false; /* nothing found */
}
 
const Champ_base& Equation_base::get_champ(const Motcle& nom) const
{
  Nom inco_residu(inconnue().le_nom());
  inco_residu += "_residu";
  if (nom == Motcle(inco_residu))
    {
      Champ_Fonc_base& ch = ref_cast_non_const(Champ_Fonc_base, field_residu_.valeur());
      double temps_init = schema_temps().temps_init();
      if (((ch.temps() != inconnue().temps()) || (ch.temps() == temps_init)) && (inconnue().mon_equation_non_nul()))
        ch.mettre_a_jour(inconnue().temps());
    }
 
  for (const auto &itr : list_champ_combi)
    if (itr->le_nom() == nom && itr->temps() != inconnue().temps())
      ref_cast_non_const(Champ_Fonc_base,itr.valeur()).mettre_a_jour(inconnue().temps());
 
  OBS_PTR(Champ_base) ref_champ;
 
  if (champs_compris_.has_champ(nom, ref_champ))
    return ref_champ;
 
  if (milieu().has_champ(nom, ref_champ))
    return ref_champ;
 
  for (int i = 0; i < nombre_d_operateurs(); i++)
    if (operateur(i).op_non_nul())
      if (operateur(i).l_op_base().has_champ(nom, ref_champ))
        return ref_champ;
 
  for (const auto &itr : les_sources)
    if (itr)
      if (itr->has_champ(nom, ref_champ))
        return ref_champ;
 
  throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
}
 
void Equation_base::get_noms_champs_postraitables(Noms& noms, Option opt) const
{
  if (opt == DESCRIPTION)
    Cerr << que_suis_je() << " : " << champs_compris_.liste_noms_compris() << finl;
  else
    noms.add(champs_compris_.liste_noms_compris());
 
  milieu().get_noms_champs_postraitables(noms, opt);
  for (int i = 0; i < nombre_d_operateurs(); i++)
    if (operateur(i).op_non_nul())
      operateur(i).l_op_base().get_noms_champs_postraitables(noms, opt);
 
  for (const auto &itr : les_sources)
    if (itr)
      itr->get_noms_champs_postraitables(noms, opt);
}
 
/*! @brief Calculation of the next time step.
 *
 * Returns the inverse of the sum of the inverses of the
 *     time steps calculated by the operators.
 *     These time steps are those for the explicit Euler scheme.
 *     The time step is not limited by dt_max, this is done in corriger_dt_calcule
 *
 * @return (double) inverse of the sum of the inverses of the time steps calculated by the operators
 */
double Equation_base::calculer_pas_de_temps() const
{
  double dt = 0;
  int nb_op = nombre_d_operateurs();
  const Schema_Temps_base& sch = le_schema_en_temps.valeur();
  bool cfl_based = sub_type(Schema_Euler_Implicite, sch) && ref_cast(Schema_Euler_Implicite, sch).facsec_cfl();
  for(int i=0; i<nb_op; i++)
    {
      const Operateur_base& op=operateur(i).l_op_base();
      bool diff_impl = sub_type(Operateur_Diff_base,op) && (sch.diffusion_implicite() || cfl_based);
      double dt_op = op.get_decal_temps()==1 ? DMAXFLOAT : operateur(i).calculer_pas_de_temps();
      Debog::verifier("Equation_base::calculer_pas_de_temps dt ",dt);
      if (dt_op>0 && !diff_impl)
        {
          // A half harmonic mean is justified by the fact that diffusion and convection are two phenomena that accumulate
          // The information at each time step cannot cross more than one computational cell so dt*U + dt*vitesse_diffusion(~alpha/dx) < dx
          // so dt < dx/(U+alpha/dx) = 1/(1/(dx/U)+1/(dx^2/alpha)) : it is indeed a half harmonic mean...
          dt = dt + 1./dt_op;
        }
      if (sch.limpr())
        {
          if (i == 0) Cout << " " << finl << "Printing of the next provisional time steps for the equation: " << que_suis_je() << finl;
          if (sub_type(Operateur_Conv_base,op))
            Cout << "   convective";
          else if (sub_type(Operateur_Diff_base,op))
            Cout << "   diffusive";
          else
            Cout << "   operator ";
          Cout<<" time step : "<< dt_op_bak[i] << finl;
        }
      dt_op_bak[i]=dt_op;
    }
  if (dt==0.)
    dt = DMAXFLOAT;
  else
    dt = 1./dt;
 
  // Implicit case
  if (sub_type(Schema_Euler_Implicite,sch) && nb_op>1)
    {
      if (cfl_based)
        // Case of zero initial velocity for example (dt=DMAFLOAT), we use the diffusion time step to start the calculation:
        if (dt>1e10) dt = dt_op_bak[0];
    }
  // Case of implicit diffusion and several operators:
  if (sub_type(Schema_Euler_explicite,sch) && nb_op>1)
    {
      if (sch.diffusion_implicite())
        {
          // Case of zero initial velocity for example (dt=DMAFLOAT), we use the diffusion time step to start the calculation:
          if (dt>1e10) dt = dt_op_bak[0];
        }
      else
        {
          // Warning on the first time steps if impliciting the diffusion
          // in an explicit scheme is worth it (dt_diff<<min(dt_max,dt_conv))
          // We do the test after the first 10 time steps in case of startup with zero velocity
          int nw = 100;
          if (sch.nb_pas_dt() > 10 && sch.nb_pas_dt() < nw && sch.limpr())
            if (dt_op_bak[0] < 0.01 * std::min(sch.pas_temps_max(), dt_op_bak[1]))
              {
                Cerr << finl << "**** Advice (printed only on the first " << nw << " time steps) ****" << finl;
                Cerr << "You could use diffusion_implicite option into the Euler scheme to increase" << finl;
                Cerr << "the stability time step of the " << this->que_suis_je() << " equation by impliciting the diffusive operator." << finl;
              }
        }
    }
  Debog::verifier("Equation_base::calculer_pas_de_temps dt ",dt);
  return dt;
}
 
// Calculation of local time: Vect of size number of faces of the domain
// Calculate whether the "steady" option is used in the "Euler implicit"
// The local time step is a convection temp step
void Equation_base::calculer_pas_de_temps_locaux(DoubleTab& dt_op) const
{
  int nb_op=nombre_d_operateurs();
  for(int i=0; i<nb_op; i++)
    {
      const Operateur_base& op=operateur(i).l_op_base();
      if (sub_type(Operateur_Conv_base,op))
        {
          if(sub_type(Op_Conv_negligeable, op))
            {
              Cerr<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<finl;
              Cerr<<"STEADY option is not compatible with the 'CONVECTION { NEGLIGEABLE }' model!"<<finl;
              Cerr << "Please, contact TRUST support." << finl;
              Cerr<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<finl;
              exit();
            }
          else
            operateur(i).calculer_pas_de_temps_locaux(dt_op);
 
        }
    }
 
}
 
/*! @brief Verifies the compatibility of boundary conditions with the equation.
 *
 * see Conds_lim::compatible_avec_eqn().
 *
 * @return (int)
 */
int Equation_base::verif_Cl() const
{
  const Conds_lim& les_cl=le_dom_Cl_dis->les_conditions_limites();
  les_cl.compatible_avec_eqn(*this);
  return les_cl.compatible_avec_discr(discretisation());
}
 
/*! @brief Returns "indeterminate" Navier_Stokes_standard for example overrides this method
 *
 *     to return "Hydraulic"
 *
 * @return (Motcle&) the name of the application domain
 */
const Motcle& Equation_base::domaine_application() const
{
  static Motcle defaut="indetermine";
  return defaut;
}
 
/*! @brief Verification of the number of components read for the specification of a field.
 *
 * Currently used for reading an initial or boundary condition.
 *
 * @param (ch_ref  : an unknown field of the considered equation)
 */
void Equation_base::verifie_ch_init_nb_comp(const Champ_Inc_base& ch_ref, const int nb_comp) const
{
  const Nature_du_champ nature = ch_ref.nature_du_champ();
  const int nb_composantes = ch_ref.nb_comp();
  const Nom& nom = ch_ref.le_nom();
 
  if (nature==scalaire)
    {
      if (nb_comp!=1)
        {
          Cerr<<"The nature of the field "<<nom<<" unknown to the equation name "<<le_nom()<<" is scalar."<<finl;
          Cerr<<"The number of components readed for this field "<<nb_comp<<" is not compatible with its nature."<<finl;
          Cerr<<"We must read only one component for this field."<<finl;
          exit();
        }
    }
  else if (nature==quadrature_scalaire)
    {
      //DG size depends on the geometry
    }
  else if (nature==vectoriel)
    {
      if (nb_comp!=Objet_U::dimension)
        {
          Cerr<<"The nature of the field "<<nom<<" unknown to the equation name "<<le_nom()<<" is vector."<<finl;
          Cerr<<"The number of components readed for this field "<<nb_comp<<" is not compatible with its nature."<<finl;
          Cerr<<"It should read "<<Objet_U::dimension<<" components for this field."<<finl;
          exit();
        }
    }
  else if (nature==quadrature_vectoriel)
    {
      //DG size depends on the geometry
    }
  else if (nature==multi_scalaire)
    {
      if (nb_comp!=nb_composantes)
        {
          Cerr<<"The nature of the field "<<nom<<" unknown to the equation name "<<le_nom()<<" is multiscalar."<<finl;
          Cerr<<"The number of components readed for this field "<<nb_comp<<" does not match the number expected."<<finl;
          Cerr<<"It should read "<<nb_composantes<<" components for this field."<<finl;
          exit();
        }
    }
  else if (nature==basis_function_order_1_scalar)
    {
      /* if (nb_comp!=3)
        {
          Cerr<<"The nature of the field "<<nom<<" unknown to the equation name "<<le_nom()<<" is a basis_function_order_1_scalar."<<finl;
          Cerr<<"The number of components readed for this field "<<nb_comp<<" does not match the number expected."<<finl;
          Cerr<<"It should read "<<3<<" components for this field."<<finl;
          exit();
        } */
    }
  else if (nature==basis_function_order_2_scalar)
    {
      /* if (nb_comp!=6)
        {
          Cerr<<"The nature of the field "<<nom<<" unknown to the equation name "<<le_nom()<<" is a basis_function_order_2_scalar."<<finl;
          Cerr<<"The number of components readed for this field "<<nb_comp<<" does not match the number expected."<<finl;
          Cerr<<"It should read "<<6<<" components for this field."<<finl;
          exit();
        } */
    }
  else if (nature==basis_function_order_1_vectorial)
    {
      /* if (nb_comp!=3*Objet_U::dimension)
        {
          Cerr<<"The nature of the field "<<nom<<" unknown to the equation name "<<le_nom()<<" is a basis_function_order_1_vectorial."<<finl;
          Cerr<<"The number of components readed for this field "<<nb_comp<<" does not match the number expected."<<finl;
          Cerr<<"It should read "<<3*Objet_U::dimension<<" components for this field."<<finl;
          exit();
        } */
    }
  else if (nature==basis_function_order_2_vectorial)
    {
      /* if (nb_comp!=6*Objet_U::dimension)
        {
          Cerr<<"The nature of the field "<<nom<<" unknown to the equation name "<<le_nom()<<" is a basis_function_order_2_vectorial."<<finl;
          Cerr<<"The number of components readed for this field "<<nb_comp<<" does not match the number expected."<<finl;
          Cerr<<"It should read "<<6*Objet_U::dimension<<" components for this field."<<finl;
          exit();
        } */
    }
  else
    {
      Cerr<<"The nature of the field is not recognized"<<finl;
      exit();
    }
}
 
/*! @brief Add convection term In:  solution is the unknown I
 *
 *  Out: secmem is increased by convection(I)
 *
 */
DoubleTab& Equation_base::derivee_en_temps_conv(DoubleTab& secmem, const DoubleTab& solution)
{
  double dt = le_schema_en_temps->pas_de_temps();
  double dt_convection;
  if (le_schema_en_temps->no_conv_subiteration_diffusion_implicite())
    {
      // false but allows not to sub-iterate on the convection
      dt_convection=1e38;
    }
  else
    dt_convection= operateur(1).calculer_pas_de_temps();
  if(inf_strict(dt_convection,dt))
    {
      // Several steps (nstep) are done if the convective time step
      // is smaller than the time step (dt) :
      double epsilon = 1.e-10;
      int nstep = (int)((dt / dt_convection)*(1+epsilon)) + 1 ;
      if (( nstep % 2 ) != 0 ) nstep += 1 ;
      double dt_loc = dt/double(nstep) ;
      Cout << "Convection Semi Implicite: Number of Sub-Cycles : " << nstep << " with dt_loc = " << dt_loc << " dt =" << dt << finl ;
 
      secmem = 0;
      DoubleTrav solution_loc(solution);
      solution_loc = solution;
      DoubleTrav derivee;
      derivee.copy(secmem, RESIZE_OPTIONS::NOCOPY_NOINIT);
      for (int i=0; i<nstep; i++)
        {
          derivee = 0. ;
          operateur(1).ajouter(solution_loc, derivee); // derivee=conv(solution_loc)
          derivee.echange_espace_virtuel();
          secmem += derivee ;
 
          // If it is not the last iteration:
          if (i+1 < nstep)
            {
              solveur_masse->appliquer(derivee);
              solution_loc.ajoute(dt_loc, derivee, VECT_REAL_ITEMS);
              solution_loc.echange_espace_virtuel();
            }
        }
      secmem /= double(nstep) ;
      return secmem;
    }
  else
    {
      operateur(1).ajouter(solution, secmem); // secmem+=conv(solution)
      return secmem;
    }
}
 
/*! @brief Solve: (1/dt + M-1*L)*dI = M-1 * secmem with a Conjugate Gradient matrix-free algorithm by default
 *
 *  L is the diffusion
 *  M is the mass
 *  In : solution=I(n)
 *  Out: solution=dI/dt
 *
 */
void Equation_base::Gradient_conjugue_diff_impl(DoubleTrav& secmem, DoubleTab& solution, int size_terme_mul, const DoubleTab& terme_mul)
{
  if (le_schema_en_temps->impr_diffusion_implicite())
    Cout << "Implicited diffusion algorithm applied on " << que_suis_je() << " equation:" << finl;
  int marq_tot = 0;
  int size_s = sources().size();
  ArrOfInt marq(size_s);
  for (int i = 0; i < size_s; i++)
    {
      if (sub_type(Source_dep_inco_base, sources()(i).valeur()))
        {
          marq[i] = 1;
          marq_tot = 1;
        }
    }
  bool has_diffusion_implicit_solver = parametre_equation() &&
                                       sub_type(Parametre_diffusion_implicite, parametre_equation().valeur()) &&
                                       ref_cast(Parametre_diffusion_implicite, parametre_equation().valeur()).solveur();
  if (has_diffusion_implicit_solver)
    {
      // PL: Solve (M/dt + L)*dI = Secmem(sans diffusion) with a matrix build to use more solvers
      if (marq_tot)
        {
          Cerr << "solveur_dffusion_implicite can't be used yet with source terms depending from the unknowns. " << finl;
          Process::exit();
        }
      if (size_terme_mul)
        {
          Cerr << "solveur_dffusion_implicite can't be used yet with penalization. " << finl;
          Process::exit();
        }
      statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
      // We multiply secmem by M (which was divided by M before the call...)
      {
        // Scope to release Trav quickly
        DoubleTrav copie(secmem);
        copie = secmem;
        secmem = 0;
        solveur_masse->ajouter_masse(1, secmem, copie);
      }
      // Build matrix A:
      Matrice_Morse& matrice = ref_cast(Parametre_diffusion_implicite, parametre_equation().valeur()).matrice();
      if (matrice.ordre()==0) dimensionner_matrice(matrice);
      matrice.clean(); // A=0
      // Add diffusion matrix L into matrix
      {
        DoubleTrav present(solution); // I(n)
        present = solution;
        operateur(0).l_op_base().contribuer_a_avec(present, matrice); // A=L
        statistics().end_count(STD_COUNTERS::matrix_assembly);
        operateur(0).ajouter(present, secmem);
        statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
        matrice.ajouter_multvect(present, secmem);
        // Add M/dt into matrix
        schema_temps().ajouter_inertie(matrice,secmem,(*this)); // A=M/dt+L
        modifier_pour_Cl(matrice,secmem);
        statistics().end_count(STD_COUNTERS::matrix_assembly);
        // Solve to get I(n+1):
        SolveurSys& solveur = ref_cast(Parametre_diffusion_implicite, parametre_equation().valeur()).solveur();
        solveur->reinit(); // New matrix but same sparsity
        int niter = solveur.resoudre_systeme(matrice, secmem, solution);
        Cout << "Diffusion operator implicited for the equation " << que_suis_je()
             << " : Conjugate gradient converged in " << niter << " iterations." << finl;
        // CHD 230501 : Call to diffusive operator to update flux_bords (boundary fluxes):
        operateur(0).ajouter(inconnue(), secmem);
 
        solution-=present; // dI=I(n+1)-I(n)
        solution/=schema_temps().pas_de_temps(); // dI/dt
      }
    }
  else
    {
      statistics().begin_count(STD_COUNTERS::implicit_diffusion,statistics().get_last_opened_counter_level()+1);
      int n = secmem.size_totale();
      if (solution.size_totale() != n)
        {
          Cerr << "the size of the solution and the second member does not match";
          exit();
        }
 
      // The maximum number of iterations can now be limited by niter_max_diff_impl
      int nmax = le_schema_en_temps->niter_max_diffusion_implicite();
 
      double aCKN = 1; // Crank - Nicholson -> 0.5
 
      // Calculation of the Diagonal matrix
 
      int precond_diag = 0;
      // We precondition by the diagonal if we penalize
      // because otherwise the initial residual is too large
      if (size_terme_mul) precond_diag = 1;
 
      double seuil_diffusion_implicite = le_schema_en_temps->seuil_diffusion_implicite();
      // Possible recovery of Parametre_diffusion_implicite options
      if (parametre_equation() && (sub_type(Parametre_diffusion_implicite, parametre_equation().valeur())))
        {
          const Parametre_diffusion_implicite& param = ref_cast(Parametre_diffusion_implicite,
                                                                parametre_equation().valeur());
          if (param.crank())
            aCKN = 0.5;
          precond_diag = param.precoditionnement_diag();
          if (param.seuil_diffusion_implicite() > 0)
            seuil_diffusion_implicite = param.seuil_diffusion_implicite();
          if (param.nb_it_max() > 0)
            nmax = param.nb_it_max();
        }
      statistics().end_count(STD_COUNTERS::implicit_diffusion,0,0);
      /////////////////////
      // Preconditioning
      /////////////////////
      double dt = le_schema_en_temps->pas_de_temps();
      if (precond_diag)
        {
          statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
          const int nb_case = inconnue().valeurs().dimension_tot(0);
          const int nb_comp = inconnue().valeurs().line_size();
          if (nb_comp * nb_case != n)
            {
              Cerr << "the size of the unknown and the second member does not match" << finl;
              Cerr << "dimension_tot nbdim1 nb_comp = " << nb_case << " " << nb_comp << finl;
              Cerr << "size_totale = " << n << finl;
              exit();
            }
          if (diag_.ordre()==1) diag_.dimensionne_diag(n);
          diag_.clean(); // Set to 0
          operateur(0).l_op_base().contribuer_a_avec(inconnue().valeurs(), diag_);
          for (int i = 0; i < size_s; i++)
            if (marq[i])
              sources()(i)->contribuer_a_avec(inconnue().valeurs(), diag_);
          // The diagonal is proportional to the control volume....
          // We must apply the mass solver
          DoubleTrav tab_tempo(inconnue().valeurs());
          {
            Matrice_Morse_View matrice; // ToDo Kokkos CMatrice_Morse_View diag = diag_.view_ro();
            matrice.set(diag_);
            DoubleTabView tempo = tab_tempo.view_wo();
            Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                                 Kokkos::RangePolicy<>(0, nb_case), KOKKOS_LAMBDA(
                                   const int ca)
            {
              for (int ncp = 0; ncp < nb_comp; ncp++)
                tempo(ca, ncp) = matrice.diag(ca * nb_comp + ncp);
            });
            end_gpu_timer(__KERNEL_NAME__);
          }
          solveur_masse->appliquer(tab_tempo);
          tab_tempo.echange_espace_virtuel();
          {
            // We invert... // Crank - Nicholson
            // The matrix corresponds to - the jacobian (to have a plus precisely, GF)
            CDoubleTabView tempo = tab_tempo.view_ro();
            CDoubleTabView terme_mul_v = terme_mul.view_ro();
            Matrice_Morse_View matrice;  // ToDo Kokkos Matrice_Morse_View diag = diag_.view_rw();
            matrice.set(diag_);
            Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                                 Kokkos::RangePolicy<>(0, nb_case), KOKKOS_LAMBDA(
                                   const int ca)
            {
              double tmp = (size_terme_mul ? terme_mul_v(ca, 0) : 1) / dt;
              for (int ncpa = 0; ncpa < nb_comp; ncpa++)
                matrice.diag(ca * nb_comp + ncpa) = 1. / (tmp + tempo(ca, ncpa) * aCKN);
            });
            end_gpu_timer(__KERNEL_NAME__);
          }
          statistics().end_count(STD_COUNTERS::matrix_assembly);
        }
      // Lambda function to apply diffusive operator and source terms to avoid code repetition
      auto matvec = [&](const DoubleTab& input_field, DoubleTab& output_field)
      {
        // Stop the counter because operator diffusion is also counted
        statistics().end_count(STD_COUNTERS::implicit_diffusion,0,0);
        operateur(0).ajouter(input_field, output_field);
        if (marq_tot)
          {
            for (int i = 0; i < size_s; i++)
              if (marq[i])
                ref_cast(Source_dep_inco_base, sources()(i).valeur()).ajouter_(input_field, output_field);
          }
        statistics().begin_count(STD_COUNTERS::implicit_diffusion,statistics().get_last_opened_counter_level()+1);
        output_field*=-1;
      };
      statistics().begin_count(STD_COUNTERS::implicit_diffusion,statistics().get_last_opened_counter_level()+1);
      // We use p to calculate phiB :
      DoubleTrav p(solution);
      DoubleTrav moins_phiB(solution); // the Boundary part of the operator with p=0
      matvec(p, moins_phiB);
 
      DoubleTrav resu(solution);
      matvec(solution, resu);
      solveur_masse->appliquer(resu);
      resu.echange_espace_virtuel();
      {
        // Create scope to release DoubleTrav quickly
        DoubleTrav sol;
        if (size_terme_mul)
          {
            sol = solution;
            ToDo_Kokkos("critical for IBC");
            //operator_multiply(sol, term_mul);
            // tab_multiply_any_shape(sol, terme_mul);
            for (int i = 0; i < size_terme_mul; i++)
              sol(i) *= terme_mul(i);
          }
        else
          sol.ref(solution);
 
        // we remove the source terms depending on the unknown ; we will add them back later
        if (marq_tot)
          {
            DoubleTrav sum_sources(secmem);
            statistics().end_count(STD_COUNTERS::implicit_diffusion,0,0);
            for (int i = 0; i < size_s; i++)
              if (marq[i])
                sources()(i).ajouter(sum_sources);
            statistics().begin_count(STD_COUNTERS::implicit_diffusion,statistics().get_last_opened_counter_level()+1);
            solv_masse().appliquer(sum_sources);
            secmem.ajoute(-1., sum_sources); // ,VECT_REAL_ITEMS);
          }
        secmem.ajoute(1. / dt, sol, VECT_REAL_ITEMS);
        resu.ajoute(1. / dt, sol, VECT_REAL_ITEMS);
      }
 
      DoubleTrav residu; // residu = Ax
      residu = resu;
      residu -= secmem;      // residu = Ax-B
      DoubleTrav z;
      if (precond_diag)
        {
          z = residu;
          diag_.multvect(residu, z); // z = D-1(Ax-B)
        }
      else
        z.ref(residu); // z = Ax-B
      p -= z;
 
      double initial_residual = mp_carre_norme_vect(z); // =||Ax(0)-B|| ou ||D-1.(Ax(0)-B)||
      /* Alternative computation of initial_residual:
         DoubleTab tmp(secmem);
         if (precond_diag)
         diag_.multvect(secmem, tmp);
         else
         tmp.ref(secmem);
         double initial_residual = mp_carre_norme_vect(tmp) ; // =||B|| ou ||D-1.B|| */
 
      if (le_schema_en_temps->impr_diffusion_implicite()) Cout << "Residu(0)=" << initial_residual << finl;
      double seuil = seuil_diffusion_implicite * seuil_diffusion_implicite;
      // We calculate a relative threshold (PL: 18/07/12, I have the impression there's an error, it's seuil*=initial_residual;
      // Moreover, it would be good to do residual = sqrt(mp_carre_norme_vect(tmp)) at several places...
      if (initial_residual > seuil)
        seuil = seuil_diffusion_implicite * initial_residual;
 
      double residual = seuil;
      double prodrz_old = mp_prodscal(residu, z);
      DoubleTrav pp;
      if (size_terme_mul) pp=p;
      DoubleTrav merk;
      merk = solution;
      int niter = 0;
      if (initial_residual != 0)
        while (niter++ <= nmax)
          {
            resu = moins_phiB; // We remove the contribution of the boundaries.
            p.echange_espace_virtuel();
            matvec(p, resu);
            solveur_masse->appliquer(resu);
 
            if (aCKN!=1) resu *= aCKN; // Crank - Nicholson
            if (size_terme_mul)
              {
                for (int i = 0; i < size_terme_mul; i++)
                  pp(i) = p(i) * terme_mul(i);
                resu.ajoute(1. / dt, pp, VECT_REAL_ITEMS);
              }
            else
              resu.ajoute(1. / dt, p, VECT_REAL_ITEMS);
            double alfa = prodrz_old / mp_prodscal(resu, p);
            // Pointless to do a virtual space exchange:
            solution.ajoute(alfa, p, VECT_REAL_ITEMS);
            residu.ajoute(alfa, resu, VECT_REAL_ITEMS);
 
            if (precond_diag)
              diag_.multvect(residu, z); // preconditioning by diag^(-1)
            // otherwise z=residu
 
            // Optimization: combine 2 MPI reductions into 1
            residual = local_carre_norme_vect(z);
            double prodrz_new = local_prodscal(residu, z);
            Process::mp_sum_for_each(residual, prodrz_new);
            if (le_schema_en_temps->impr_diffusion_implicite())
              Cout << "Iteration n=" << niter << " Residu(n)/Residu(0)=" << residual / initial_residual << finl;
 
            if (residual<seuil)
              break;
            else
              {
                p *= (prodrz_new / prodrz_old);
                p -= z;
                prodrz_old = prodrz_new;
                // We warn if convergence is abnormally long on very large cases
                if (niter == 100)
                  {
                    Cout
                        << "Already 100 iterations of the conjugate gradient solver used for diffusion implicited algorithm."
                        << finl;
                    Cout
                        << "The calculation is diverging, may be because of incorrect initial and/or boundary conditions for the equation "
                        << que_suis_je() << finl;
                    Cout
                        << "You may also try to stop the calculation and rerun it with a lower facsec, say 0.9 or may be less: 0.5."
                        << finl;
                  }
              }
          }
      if ((le_schema_en_temps->no_error_if_not_converged_diffusion_implicite() == 0) && (niter > nmax))
        {
          Cerr << "No convergence of the implicited diffusion algorithm for the equation " << que_suis_je() << " in "
               << nmax << " iterations." << finl;
          Cerr << "Residue : " << residual << " Threshold : " << seuil << finl;
          Cerr << "======================================================================================" << finl;
          Cerr << "The problem is post processed to help you to see where the non convergence is located." << finl;
          schema_temps().set_stationnaire_atteint()=0;
          probleme().postraiter(1);
          probleme().mettre_a_jour(schema_temps().temps_courant());
          probleme().finir();
          Cerr << "The problem " << probleme().le_nom() << " has been saved too." << finl;
          exit();
        }
      else
        Cout << "Diffusion operator implicited for the equation " << que_suis_je()
             << " : Conjugate gradient converged in " << niter << " iterations." << finl;
      solution -= merk;
      solution.echange_espace_virtuel();
 
      // CHD 230501 : last call to update flux_bords (boundary fluxes):
      resu = 0.;
      matvec(inconnue().valeurs(), resu);
 
      // Since 1.6.8 returns dI/dt:
      solution/=dt;
 
      // End the counter
      statistics().end_count(STD_COUNTERS::implicit_diffusion);
    }
}
 
DoubleTab& Equation_base::corriger_derivee_expl(DoubleTab& d)
{
  return d;
}
DoubleTab& Equation_base::corriger_derivee_impl(DoubleTab& d)
{
  return d;
}
 
const RefObjU& Equation_base::get_modele(Type_modele type) const
{
  return liste_modeles_(0);
 
}
// MODIF ELI LAUCOIN (22/11/2007) : I add an advance and a backward
// by default, this only advances or goes back the unknown (and champ_conserve_ if initialized)
void Equation_base::avancer(int i)
{
  inconnue().avancer(i);
  if (champ_conserve_) champ_conserve_->avancer(i);
  if (champ_convecte_) champ_convecte_->avancer(i);
}
 
void Equation_base::reculer(int i)
{
  inconnue().reculer(i);
  if (champ_conserve_) champ_conserve_->reculer(i);
  if (champ_convecte_) champ_convecte_->reculer(i);
}
// END MODIF ELI LAUCOIN (22/11/2007)
 
 
#define BLOQUE Cerr<<__PRETTY_FUNCTION__<< " "<<__FILE__<<":"<<(int)__LINE__<<" not coded, retrieves coding simpler" <<finl;exit()
 
// methods for the implicit
 
/* can use a memoization (discretisations PolyMAC_HFV)*/
void Equation_base::dimensionner_matrice(Matrice_Morse& matrice)
{
  if (matrice_init)
    {
      matrice.get_set_tab1().ref_array(matrice_stockee.get_set_tab1());
      matrice.get_set_tab2().ref_array(matrice_stockee.get_set_tab2());
      matrice.get_set_coeff().ref_array(matrice_stockee.get_set_coeff());
      matrice.set_nb_columns(matrice_stockee.nb_colonnes());
      matrice.sorted_ = matrice_stockee.sorted_;
      matrice.constant_stencil() = matrice_stockee.constant_stencil();
      matrice.get_set_coeff() = 0.0;
      return;
    }
 
  dimensionner_matrice_sans_mem(matrice); //le vrai appel
 
  matrice.get_set_coeff() = 0.0;  // just to be sure ...
 
  if (probleme().discretisation().is_poly_family() || probleme().discretisation().is_vef())
    {
      // PL: ToDo why sorting leads to issue in VEF ???
      if (probleme().discretisation().is_poly_family())
        {
          matrice.sort_stencil();
          matrice_stockee.sorted_ = 1;
        }
      matrice.constant_stencil()=true; // stencil will not change anymore
      matrice_stockee.get_set_tab1().ref_array(matrice.get_set_tab1());
      matrice_stockee.get_set_tab2().ref_array(matrice.get_set_tab2());
      matrice_stockee.get_set_coeff().ref_array(matrice.get_set_coeff());
      matrice_stockee.set_nb_columns(matrice.nb_colonnes());
      matrice_stockee.constant_stencil()=matrice.constant_stencil();
      matrice_init = 1;
    }
}
 
void Equation_base::dimensionner_matrice_sans_mem(Matrice_Morse& matrice)
{
  // [ABN]: dimensioning of the implicit matrix: it can receive input from:
  //  - the operators
  //  - the mass solver
  //  - the sources
  // We proceed in this order, with the idea that operators should take most of the slots in the Morse
  // structure. Having the mass solver first would lead (in VDF, EF and VEF) to diagonal slots
  // being reserved before sub-diagonal terms. In the Morse structure coefficients would then
  // not be ordered by increasing columns number. This is not desirable, notably for EF which seem to
  // be sensitive to this ordering (although they should not ....). Some Baltiks too (FLICA5 - thermic part)
  // rely on this.
  if (matrice_init) // memoization
    {
      matrice.get_set_tab1().ref_array(matrice_stockee.get_set_tab1());
      matrice.get_set_tab2().ref_array(matrice_stockee.get_set_tab2());
      matrice.get_set_coeff().ref_array(matrice_stockee.get_set_coeff());
      matrice.set_nb_columns(matrice_stockee.nb_colonnes());
      matrice.constant_stencil() = matrice_stockee.constant_stencil();
      return;
    }
 
  bool isInit = false;
  // Operators first (they're likely to take most of the slots in the Morse structure)
  for(int opp=0; opp < nombre_d_operateurs(); opp++)
    {
      if (!isInit)
        operateur(opp).l_op_base().dimensionner(matrice);
      else
        {
          Matrice_Morse mat2;
          operateur(opp).l_op_base().dimensionner(mat2);
          if (mat2.nb_colonnes())
            matrice += mat2;  // this only works if the matrix has been given its overall size first
        }
      isInit = isInit || (matrice.nb_colonnes() != 0);
 
    }
 
  //  ... then the mass solver ...
  if (!isInit)
    solv_masse().dimensionner(matrice);
  else
    {
      Matrice_Morse mat2;
      solv_masse().dimensionner(mat2);
      matrice += mat2; // this only works if the matrix has been given its overall size first
    }
 
  // and finally sources (mass solver has surely done something already so no need for "isInit" test)
  les_sources.dimensionner(matrice);
}
 
void Equation_base::dimensionner_termes_croises(Matrice_Morse& matrice, const Probleme_base& autre_pb, int nl, int nc)
{
  matrice.dimensionner(nl, nc, 0);
  for(int i = 0; i < nombre_d_operateurs(); i++)
    {
      Matrice_Morse mat2;
      operateur(i).l_op_base().dimensionner_termes_croises(mat2, autre_pb, nl, nc);
      if (mat2.nb_colonnes()) matrice += mat2;
    }
}
 
void Equation_base::ajouter_termes_croises(const DoubleTab& inco, const Probleme_base& autre_pb, const DoubleTab& autre_inco, DoubleTab& resu) const
{
  for (int i = 0; i < nombre_d_operateurs(); i++)
    operateur(i).l_op_base().ajouter_termes_croises(inco, autre_pb, autre_inco, resu);
}
 
void Equation_base::contribuer_termes_croises(const DoubleTab& inco, const Probleme_base& autre_pb, const DoubleTab& autre_inco, Matrice_Morse& matrice) const
{
  for (int i = 0; i < nombre_d_operateurs(); i++)
    operateur(i).l_op_base().contribuer_termes_croises(inco, autre_pb, autre_inco, matrice);
}
 
// adds the contributions of operators and sources
void Equation_base::assembler(Matrice_Morse& matrice, const DoubleTab& inco, DoubleTab& resu)
{
  // Test to verify the contribuer_a_avec method
  for (int op=0; op<nombre_d_operateurs(); op++)
    operateur(op).l_op_base().tester_contribuer_a_avec(inco, matrice);
 
  // Contribution of operators and sources:
  // [Vol/dt+A]Inco(n+1)=sum(residual)+Vol/dt*Inco(n)
  // Typically: if Op=flux(Inco) then the implicit matrix A contains a contribution -dflux/dInco
  // Example: Op convective flux in VDF:
  // Op=T*u*S and A=-d(T*u*S)/dT=-u*S
  const Discretisation_base::type_calcul_du_residu& type_codage=probleme().discretisation().codage_du_calcul_du_residu();
  if (type_codage==Discretisation_base::VIA_CONTRIBUER_AU_SECOND_MEMBRE)
    {
      // We calculate sum(residual) by contribuer_au_second_membre (typically non-implicit boundary conditions)
      // This approach requires coding 3 methods (contribuer_a_avec, contribuer_au_second_membre and add for explicit)
      sources().contribuer_a_avec(inco,matrice);
      statistics().end_count(STD_COUNTERS::matrix_assembly,0,0);
      sources().ajouter(resu);
      statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
      matrice.ajouter_multvect(inco, resu); // Add source residual first
      for (int op = 0; op < nombre_d_operateurs(); op++)
        {
          operateur(op).l_op_base().contribuer_a_avec(inco, matrice);
          operateur(op).l_op_base().contribuer_au_second_membre(resu);
        }
    }
  else if (type_codage==Discretisation_base::VIA_AJOUTER)
    {
      // We calculate sum(residual) by sum(operator)+sources+A*Inco(n)
      // This approach requires coding only two methods (contribuer_a_avec and add)
      // So a bit more costly in calculation time but less code to write/maintain
 
      /*
       * XXX XXX XXX Elie Saikali, Gauthier Fauchet, aout 2025 :
       *
       * We need to calculate : [ I^(n+1) - I^n ] / dt = f( I^(n+1) )
       *
       * But we dont know now what is f( I^(n+1) )
       *
       * 2 possibilities
       *
       *  1- [ I^(k+1) - I^k ] / dt = f( I^(k) )
       *    then use fixed point algo to have I^k = I^(k+1) = I^(n+1)
       *
       *
       *  2- Use Newton algorithm (we use this here)
       *
       *    f( I^(n+1) ) = f ( I^n ) + df/dI ( I^(n+1) - I^n )
       *
       *    df/dI is the matrix -A obtained by the method contribuer_a_avec
       *
       *    so : f( I^(n+1) ) = f ( I^n ) - A ( I^(n+1) - I^n ) = f ( I^n ) - A * I^(n+1) + A * I^n
       *
       *    Thus, 1/dt I^(n+1) + A * I^(n+1) = 1/dt I^n + f ( I^n ) + A * I^n
       *
       *    I^(n+1) [ Identity / dt + A ] =  1/dt I^n + f ( I^n ) + A * I^n
       *
       *    and f ( I^n ) is obtained by the method ajouter
       */
      for (int op=0; op<nombre_d_operateurs(); op++)
        {
          operateur(op).l_op_base().contribuer_a_avec(inco, matrice);
          statistics().end_count(STD_COUNTERS::matrix_assembly,0,0);
          operateur(op).ajouter(resu);
          statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
        }
      sources().contribuer_a_avec(inco,matrice);
      statistics().end_count(STD_COUNTERS::matrix_assembly,0,0);
      sources().ajouter(resu);
      statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
      matrice.ajouter_multvect(inco, resu); // Addition of A*Inco(n)
      // PL (11/04/2018): We would like to calculate the contribution of sources first
      // as in the VIA_CONTRIBUER_AU_SECOND_MEMBRE case but the Canal_perio_3D case (keps
      // periodic crashes: there is a periodicity error in the source terms of the KEps model...
    }
  else
    {
      Cerr << "Unknown value in Equation_base::assembler for " << (int)type_codage << finl;
      Process::exit();
    }
}
 
// modifies the matrix and the second member according to the boundary conditions
void Equation_base::modifier_pour_Cl(Matrice_Morse& mat_morse, DoubleTab& secmem) const
{
  for (int op=0; op<nombre_d_operateurs(); op++)
    {
      operateur(op).l_op_base().modifier_pour_Cl(mat_morse,secmem);
    }
}
// assembles, adds inertia, and modifies for boundary conditions.
void Equation_base::assembler_avec_inertie( Matrice_Morse& mat_morse,const DoubleTab& present, DoubleTab& secmem)
{
  statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
  assembler(mat_morse,present,secmem);
  schema_temps().ajouter_inertie(mat_morse,secmem,(*this));
  modifier_pour_Cl(mat_morse,secmem);
  statistics().end_count(STD_COUNTERS::matrix_assembly);
}
 
/* verifies that all terms are compatible */
int Equation_base::has_interface_blocs() const
{
  int ok = 1;
  /* operators, mass, sources */
  for (int op = 0; op < nombre_d_operateurs(); op++) ok &= operateur(op).l_op_base().has_interface_blocs();
  ok &= solv_masse().has_interface_blocs();
  for (int i = 0; i < les_sources.size(); i++) ok &= les_sources(i)->has_interface_blocs();
  return ok;
}
 
/* like dimensionner_matrice(), but without memoization */
void Equation_base::dimensionner_blocs(matrices_t matrices, const tabs_t& semi_impl) const
{
  /* operators, mass, sources */
  int nb = nombre_d_operateurs();
  for (int op = 0; op < nb; op++) operateur(op).l_op_base().dimensionner_blocs(matrices, semi_impl);
  solv_masse().dimensionner_blocs(matrices);
  for (int i = 0; i < les_sources.size(); i++) les_sources(i)->dimensionner_blocs(matrices, semi_impl);
}
 
void Equation_base::assembler_blocs(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl) const
{
  /* setting to zero */
  secmem = 0;
  for (auto && i_m : matrices) i_m.second->get_set_coeff() = 0;
  /* operators, sources, mass */
  for (int i = 0; i < nombre_d_operateurs(); i++)
    operateur(i).l_op_base().ajouter_blocs(matrices, secmem, semi_impl);
 
  statistics().end_count(STD_COUNTERS::ajouter_blocs,0,0);
 
  statistics().begin_count(STD_COUNTERS::source_terms,statistics().get_last_opened_counter_level()+1);
  for (int i = 0; i < les_sources.size(); i++)
    les_sources(i)->ajouter_blocs(matrices, secmem, semi_impl);
 
  statistics().end_count(STD_COUNTERS::source_terms);
 
 
  statistics().begin_count(STD_COUNTERS::ajouter_blocs,statistics().get_last_opened_counter_level()+1);
  if (!(discretisation().is_poly_family() || probleme().que_suis_je().debute_par("Pb_Multiphase") || que_suis_je().debute_par("Equation_flux")))
    {
      const std::string& nom_inco = inconnue().le_nom().getString();
      Matrice_Morse *mat = matrices.count(nom_inco) ? matrices.at(nom_inco) : nullptr;
      if(mat) mat->ajouter_multvect(inconnue().valeurs(), secmem);
    }
}
 
void Equation_base::assembler_blocs_avec_inertie(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl)
{
  statistics().begin_count(STD_COUNTERS::ajouter_blocs,statistics().get_last_opened_counter_level()+1);
  assembler_blocs(matrices, secmem, semi_impl);
  solv_masse().set_penalisation_flag(0);
  schema_temps().ajouter_blocs(matrices, secmem, *this);
 
  if (!discretisation().is_poly_family())
    {
      const std::string& nom_inco = inconnue().le_nom().getString();
      Matrice_Morse *mat = matrices.count(nom_inco) ? matrices.at(nom_inco) : nullptr;
      modifier_pour_Cl(*mat,secmem);
    }
  statistics().end_count(STD_COUNTERS::ajouter_blocs);
}
 
/* creation of champ_conserve_, cl_champ_conserve_ */
void Equation_base::init_champ_conserve() const
{
  if (champ_conserve_) return; //already done
  int Nt = inconnue().nb_valeurs_temporelles(),
      Nl = inconnue().valeurs().size_reelle_ok() ? inconnue().valeurs().dimension(0) : -1,
      Nc = inconnue().valeurs().line_size();
  //champ_conserve_ : same type / support as the unknown
  discretisation().creer_champ(champ_conserve_, domaine_dis(), inconnue().que_suis_je(), "N/A", "N/A", Nc, Nl, Nt, schema_temps().temps_courant());
  champ_conserve_->associer_eqn(*this);
  auto nom_fonc = get_fonc_champ_conserve();
  champ_conserve_->nommer(nom_fonc.first);
  champ_conserve_->init_champ_calcule(*this, nom_fonc.second);
}
 
/* default calculation method of champ_conserve : product of coefficient_temporel * unknown */
void Equation_base::calculer_champ_conserve(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Equation_base& eqn = ref_cast(Equation_base, obj);
  const Champ_base *coeff = eqn.solv_masse().has_coefficient_temporel() ? &eqn.get_champ(eqn.solv_masse().get_name_of_coefficient_temporel()) : nullptr; //temporal coeff
  const Champ_Inc_base& inco = eqn.inconnue();
  ConstDoubleTab_parts part(inco.valeurs());
  //value of the field itself
  val = part[0];
  if (coeff) tab_multiply_any_shape(val, coeff->valeurs(), VECT_ALL_ITEMS);
 
  bval = inco.valeur_aux_bords();
  if (coeff) tab_multiply_any_shape(bval, coeff->valeur_aux_bords(), VECT_ALL_ITEMS);
 
  DoubleTab& der = deriv[inco.le_nom().getString()];
  if (coeff) der = coeff->valeurs();
  else der.resize(val.dimension_tot(0), val.line_size()), der = 1;
}
 
// Printing of the residual in file (usually dt_ev)
// This method can be overridden by equations
// printing particular residuals (K-Eps, Concentrations,...)
void Equation_base::imprime_residu(SFichier& fic)
{
  int size = residu_.size_array();
  for (int i=0; i<size; i++)
    {
      Nom tab = schema_temps().norm_residu() == "max" ? "\t " : "\t\t ";
      fic << residu_(i) << tab;
      int nb_unknowns = inconnue().valeurs().dimension_tot(0);
      double residual_limit = ( nb_unknowns == 0 ? DMAXFLOAT : DMAXFLOAT/nb_unknowns);
      if (residu_(i)>residual_limit)
        {
          fic << finl;
          Cerr << "Equation " << le_nom() << " is diverging:" << finl;
          Cerr << "A residual ||dI/dt|| for this equation is bigger than " << residual_limit << " !" << finl;
          Cerr << "Have a look at the " << nom_du_cas() << ".dt_ev file." << finl;
          Process::exit();
        }
    }
  // Display min/max
  if (schema_temps().impr_extremums() && limpr())
    {
      double vmax = mp_max_vect(inconnue().valeurs());
      double vmin = mp_min_vect(inconnue().valeurs());
      Cout << finl << inconnue().le_nom() << " field [min/max]: ";
      if (je_suis_maitre()) printf("[ %.2e / %.2e ]",vmin,vmax);
      Cout << finl;
    }
}
 
// Returns the expression of the residual (likewise can be overridden)
Nom Equation_base::expression_residu()
{
  int size = residu_.size_array();
  Nom tmp(""),ajout("");
  if (probleme().is_coupled())
    {
      ajout="_";
      ajout+=probleme().le_nom();
    }
  for (int i=0; i<size; i++)
    {
      Nom norm = schema_temps().norm_residu();
      tmp+="Ri=";
      tmp+=norm;
      if(norm != "max") tmp+="-norm";
      tmp+=ajout;
      tmp+="|d";
      // How to name clearly AND concisely
      // the field of the expression:
      if (size==1)
        {
          // Scalar field
          tmp+=(Motcle)inconnue().le_nom()[0];
        }
      else
        {
          // Multi-scalar field
          if (inconnue().nom_compo(0)[0]!=inconnue().nom_compo(size-1)[0])
            {
              // Example, K_Eps field
              tmp+=(Motcle)inconnue().nom_compo(i)[0];
            }
          else
            {
              // Example, concentration field
              tmp+=(Motcle)inconnue().le_nom()[0];
              tmp+=(Nom)i;
            }
        }
      tmp+="/dt|\t ";
    }
  return tmp;
}
 
// Dimension the residual arrays
void Equation_base::initialise_residu(int size)
{
  if (size==0)
    size = inconnue().valeurs().line_size();
  residu_.resize(size);
  residu_initial_.resize(size);
  residu_=0;
  residu_initial_=0;
}
 
// Fills the field_residu_ field if it exists
void Equation_base::set_residuals(const DoubleTab& residual)
{
  if(field_residu_)
    {
      DoubleTab& tab = field_residu_->valeurs();
      if (tab.dimension_tot(0) == residual.dimension_tot(0))
        tab = residual;
      else
        {
          ConstDoubleTab_parts parts(residual);
          if (parts[0].dimension_tot(0) == tab.dimension_tot(0))
            tab = parts[0];
        }
    }
}
 
const Operateur& Equation_base::operateur_fonctionnel(int i) const
{
  Cerr<<"The method operateur_fonctionnel is not coded for equation "<<que_suis_je()<<finl;
  exit();
  throw;
}
 
Operateur& Equation_base::operateur_fonctionnel(int i)
{
  Cerr<<"The method operateur_fonctionnel is not coded for equation "<<que_suis_je()<<finl;
  exit();
  throw;
}
 
int Equation_base::nombre_d_operateurs_tot() const
{
  return nombre_d_operateurs();
}