Debog_Pb.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 <Octree_Double.h>
#include <Probleme_base.h>
#include <PE_Groups.h>
#include <Debog_Pb.h>
#include <Equation_base.h>
#include <EChaine.h>
#include <Domaine_VF.h>
#include <Domaine.h>
#include <Param.h>
 
OBS_PTR(Debog_Pb) Debog_Pb::instance_debog_;
 
Implemente_instanciable(Debog_Pb,"Debog_pb",Objet_U);
Implemente_instanciable(Debog_Pb_Wrapper,"Debog",Interprete);
// XD debog interprete debog NO_BRACE Class to debug some differences between two TRUST versions on a same data file.
// XD_CONT NL2 If you want to compare the results of the same code in sequential and parallel calculation, first run
// XD_CONT (mode=0) in sequential mode (the files fichier1 and fichier2 will be written first) then the second run in
// XD_CONT parallel calculation (mode=1). NL2 During the first run (mode=0), it prints into the file DEBOG, values at
// XD_CONT different points of the code thanks to the C++ instruction call. see for example in
// XD_CONT Kernel/Framework/Resoudre.cpp file the instruction: Debog::verifier(msg,value); Where msg is a string and
// XD_CONT value may be a double, an integer or an array. NL2 During the second run (mode=1), it prints into a file
// XD_CONT Err_Debog.dbg the same messages than in the DEBOG file and checks if the differences between results from
// XD_CONT both codes are less than a given value (error). If not, it prints Ok else show the differences and the lines
// XD_CONT where it occured.
// XD attr pb ref_pb_gen_base pb REQ Name of the problem to debug.
// XD attr fichier1 chaine file1 REQ Name of the file where domain will be written in sequential calculation.
// XD attr fichier2 chaine file2 REQ Name of the file where faces will be written in sequential calculation.
// XD attr seuil floattant seuil REQ Minimal value (by default 1.e-20) for the differences between the two codes.
// XD attr mode entier mode REQ By default -1 (nothing is written in the different files), you will set 0 for the
// XD_CONT sequential run, and 1 for the parallel run.
 
 
Sortie& Debog_Pb::printOn(Sortie& os) const
{
  exit();
  return os;
}
 
/*! @brief Method called as soon as an error is found in real space.
 *
 * (useful for setting a breakpoint in gdb)
 *
 */
void Debog_Pb::error_function()
{
  if (exit_on_error_)
    {
      Cerr << "Error found in Debog_Pb (see debog.log files)" << finl;
      Process::exit();
    }
}
 
Entree& Debog_Pb::readOn(Entree& is)
{
  fichier_domaine_ = "dom.debog";
  fichier_faces_ = "faces.debog";
  fichier_debog_ = "DEBOG";
 
  Param param(que_suis_je());
  param.ajouter("nom_pb", &nom_pb_, Param::REQUIRED);
  param.ajouter("fichier_domaine", &fichier_domaine_);
  param.ajouter("fichier_faces", &fichier_faces_);
  param.ajouter("fichier_debog", &fichier_debog_);
  param.ajouter("seuil_absolu", &seuil_absolu_);
  param.ajouter("seuil_relatif", &seuil_relatif_);
  param.ajouter("ignorer_messages", &msg_a_ignorer_);
  param.ajouter("noeuds_doubles_ignores",&noeuds_doubles_ignores_);
  param.ajouter("mode_db", &mode_db_, Param::REQUIRED);
  param.dictionnaire("write", 0);
  param.dictionnaire("read", 1);
  param.ajouter("exit_on_error", &exit_on_error_);
  param.lire_avec_accolades(is);
 
  Objet_U& obj = Interprete::objet(nom_pb_);
  if (!sub_type(Probleme_base, obj))
    {
      Cerr << "Error in Debog::interpreter: " << obj << " is not a Problem_base object" << finl;
      Process::exit();
    }
 
  Probleme_base& pb = ref_cast(Probleme_base, obj);
  ref_pb_ = pb;
 
  if (mode_db_ == 0)
    {
      // Writing:
      write_geometry_data();
      write_debog_data_file_.ouvrir(fichier_debog_);
      write_debog_data_file_.setf(ios::scientific);
      write_debog_data_file_.precision(20);
    }
  else
    {
      read_geometry_data();
      debog_data_file_.ouvrir(fichier_debog_);
      debog_data_file_.set_error_action(Entree::ERROR_EXIT);
    }
 
  if (Process::je_suis_maitre()) log_file_.ouvrir("debog.log");
 
  detailed_log_file_.ouvrir("debog_detail.log");
 
  instance_debog_ = *this;
  return is;
}
 
 
int Debog_Pb::test_ignore_msg(const char* const msg)
{
  return msg_a_ignorer_.rang(msg) >= 0;
}
 
void Debog_Pb::goto_msg(const char *const message)
{
  assert(mode_db_ == 1);
  Nom dummy;
  Nom n;
  int num;
  Nom msg(message); // copy without leading/trailing spaces
  while (msg.debute_par(" "))
    msg.suffix(" ");
  while (msg.finit_par(" "))
    msg.prefix(" ");
  while (1)
    {
      // read until we find a "msg" word:
      do
        {
          debog_data_file_ >> dummy;
        }
      while (dummy != "msg");
 
      debog_data_file_ >> dummy;
      if (dummy != ":")
        {
          Cerr << "Error in Debog_Pb::goto_msg(" << msg << ")\n invalid data in DEBOG file" << finl;
          exit();
        }
      num = -1;
      debog_data_file_ >> num;
      debog_data_file_ >> dummy;
      if (dummy != ":")
        {
          Cerr << "Error in Debog_Pb::goto_msg(" << msg << ")\n invalid data in DEBOG file" << finl;
          exit();
        }
      n = "";
      while (1)
        {
          debog_data_file_ >> dummy;
          if (dummy == "FinMsg")
            break;
          n += dummy;
          n += " ";
        }
      if (n.debute_par(msg))
        break;
      if (Process::je_suis_maitre())
        {
          log_file_ << "Skipping message " << n << " (looking for " << msg << ")" << finl;
        }
    }
  if (Process::je_suis_maitre())
    log_file_ << "Reading message " << num << "(file) " << debog_msg_count_ << "(current) " << msg << finl;
 
  detailed_log_file_ << "Reading message " << num << "(file) " << debog_msg_count_ << "(current) " << msg << finl;
}
 
void Debog_Pb::write_geometry_data()
{
  if (Process::is_parallel())
    {
      Cerr << "Error in Debog.cpp: cannot write geometry data in parallel." << finl;
      Process::exit();
    }
  const Domaine& dom = ref_pb_->domaine();
  const Domaine_dis_base& zd = ref_pb_->domaine_dis();
  const Domaine_VF& zvf = ref_cast(Domaine_VF, zd);
  {
    SFichier f(fichier_domaine_);
    f.precision(20);
    f << dom;
  }
  {
    SFichier f(fichier_faces_);
    f.precision(20);
    f << zvf.elem_faces();
    f << zvf.xv();
    f << zvf.face_voisins();
    if (zvf.xa().dimension(0) > 0)
      f << zvf.xa();
  }
  register_item(dom.md_vector_sommets(), "SOM");
  register_item(dom.md_vector_elements(), "ELEM");
  register_item(zvf.md_vector_faces(), "FACE");
  if (zvf.md_vector_aretes())
    register_item(zvf.md_vector_aretes(), "ARETES");
}
 
void Debog_Pb::register_item(const MD_Vector& md, const Nom& id)
{
  known_md_.add(md);
  renum_id_.add(id);
}
 
void Debog_Pb::add_renum_item(const DoubleTab& coord_ref, const DoubleTab& coord_par, const MD_Vector& md, const Nom& id)
{
  register_item(md, id);
 
  const double epsilon = Objet_U::precision_geom;
 
  // We create an octree containing the points of the reference domain
  // and we search for each point of the local domain in this octree.
  Octree_Double octree;
  octree.build_nodes(coord_ref, 1 /* include_virtual_items */);
 
  // Add a new renum array:
  renum_array_.add(IntVect());
  IntVect& renum = renum_array_[renum_array_.size()-1];
  MD_Vector_tools::creer_tableau_distribue(md, renum, RESIZE_OPTIONS::NOCOPY_NOINIT);
 
  if (renum.size_totale() != coord_par.dimension_tot(0))
    {
      Cerr << "Internal error in Debog::add_renum_item: coordinate array does not match md descriptor" << finl;
      Process::exit();
    }
  const int dim = coord_par.dimension(1);
  ArrOfDouble center(dim);
  ArrOfInt elements;
 
  const int n = renum.size_totale();
  for (int i = 0; i < n; i++)
    {
      for (int j = 0; j < dim; j++)
        center[j] = coord_par(i,j);
      octree.search_elements_box(center, epsilon, elements);
      // the first search returns all potentially close vertices.
      // a test must then be performed on each vertex:
      octree.search_nodes_close_to(center, coord_ref, elements, epsilon);
      const int k = elements.size_array();
      if (k != 1)
        {
          Cerr << "Debog::add_renum_item: Error. Id=" << id
               << "\n Item " << i << " with following coordinates was found " << k << " times within epsilon=" << epsilon << " in the reference geometry: " << center << finl;
          if (noeuds_doubles_ignores_==0)
            {
              Cerr << "If you want to discard nodes comparison in the mesh, add 'noeuds_doubles_ignores 1' option in debog_pb keyword." << finl;
              exit();
            }
          else
            {
              Cerr << "Warning, this item will be ignored during the comparison..." << finl;
              renum[i] = -1;
            }
        }
      else
        {
          // renum[i] is the index of the i-th local item in the reference geometry
          renum[i] = elements[0];
        }
    }
}
 
void Debog_Pb::read_geometry_data()
{
  const Domaine& dom = ref_pb_->domaine();
  const Domaine_dis_base& zd = ref_pb_->domaine_dis();
  const Domaine_VF& zvf = ref_cast(Domaine_VF, zd);
  {
    DoubleTab coord_som_seq; // vertices
    DoubleTab xp_seq; // element centers
    // We must enter a single-processor group for Domaine::readOn:
    {
      OWN_PTR(Comm_Group) group;
      ArrOfInt liste_procs(1); // List of 1 processor containing proc 0
      PE_Groups::create_group(liste_procs, group, 1);
      if (PE_Groups::enter_group(group.valeur()))
        {
          EFichier f(fichier_domaine_);
          Domaine dom_seq;
          f >> dom_seq;
          coord_som_seq = dom_seq.coord_sommets();
          const Elem_geom_base& elem = dom_seq.type_elem().valeur();
          elem.calculer_centres_gravite(xp_seq);
          PE_Groups::exit_group();
        }
    }
    // All processors receive the coordinates of the reference vertices
    // (destroy the descriptor otherwise printOn refuses to send the vector)
    coord_som_seq.set_md_vector(MD_Vector());
    envoyer_broadcast(coord_som_seq, 0);
    // Same for element centers
    xp_seq.set_md_vector(MD_Vector());
    envoyer_broadcast(xp_seq, 0);
 
    // Compute the renumberings:
    const DoubleTab& coord_som = dom.coord_sommets();
    const MD_Vector& md_som = dom.coord_sommets().get_md_vector();
    add_renum_item(coord_som_seq, coord_som, md_som, "SOM");
    const DoubleTab& xp = zvf.xp();
    // The xp array has no parallel structure, we take it from the elements...
    const MD_Vector& md_elem = dom.les_elems().get_md_vector();
    add_renum_item(xp_seq, xp, md_elem, "ELEM");
  }
  {
    DoubleTab xv_seq;
    DoubleTab xa_seq;
    int avec_aretes = 1;
    if (Process::je_suis_maitre())
      {
        EFichier f(fichier_faces_);
        f.precision(20);
        {
          IntTab elem_faces;
          f >> elem_faces; // not used
        }
        f >> xv_seq;
        {
          IntTab face_voisins;
          f >> face_voisins; // not used
        }
        f.set_error_action(Entree::ERROR_EXCEPTION);
        try
          {
            f >> xa_seq;
          }
        catch (Entree_Sortie_Error&)
          {
            avec_aretes = 0;
          }
      }
    envoyer_broadcast(xv_seq, 0);
    const DoubleTab& xv = zvf.xv();
    const MD_Vector& md_face = zvf.md_vector_faces();
    add_renum_item(xv_seq, xv, md_face, "FACE");
 
    envoyer_broadcast(avec_aretes, 0);
    if (avec_aretes)
      {
        envoyer_broadcast(xa_seq, 0);
        const DoubleTab& xa = zvf.xa();
        const MD_Vector& md_arete = zvf.md_vector_aretes();
        add_renum_item(xa_seq, xa, md_arete, "ARETE");
      }
  }
}
 
enum Debog_Exceptions { RENUM_ARRAY_NOT_FOUND };
 
const IntVect& Debog_Pb::find_renum_vector(const MD_Vector& mdv, Nom& id) const
{
  const int n = renum_array_.size();
  for (int i = 0; i < n; i++)
    {
      const IntVect& renum = renum_array_[i];
      if (renum.get_md_vector() == mdv)
        {
          id = renum_id_[i];
          return renum;
        }
    }
  Cerr << "Error, we don't find renum_vector for the array on item: "<< id << finl;
  throw RENUM_ARRAY_NOT_FOUND;
}
 
void Debog_Pb::set_nom_pb_actuel(const Nom& nom)
{
  nom_pb_actuel_ = nom;
}
 
/*! @brief md_lignes: descriptor of row indices (i.e. descriptor of vector b in A*x=b) md_colonnes: same, column indices (i.e. descriptor of vector x in A*x=b)
 *
 */
void Debog_Pb::verifier_matrice(const char *const msg, const Matrice_Base& matrice, const MD_Vector& md_lignes, const MD_Vector& md_colonnes)
{
  // Warning: this implementation only works for md_colonnes types registered in renum_array_
  //  (i.e. simple types: vertices, elements, faces, edges, and not MD_Vector_composite)
  // and only for linesize==1 (one row and one column of the matrix per item)
 
  // For each column i of the matrix, writing the product matrix * vector_i
  //  with vector_i[j] = (i==j);
  DoubleVect vecteur_i;
  DoubleVect tmp;
  MD_Vector_tools::creer_tableau_distribue(md_colonnes, vecteur_i);
  MD_Vector_tools::creer_tableau_distribue(md_lignes, tmp);
  Nom id;
 
  // Loop over the sequential items of the vector x in A*x=b
  const trustIdType nbc0 = md_colonnes->nb_items_seq_tot();
  if (nbc0 > std::numeric_limits<int>::max())
    Process::exit("Debog_Pb::verifier_matrice() - total number of items too big!");
  const int nb_colonnes = static_cast<int>(nbc0);
 
  for (int i = 0; i < nb_colonnes; i++)
    {
      int renum_i;
      if (mode_db_ == 0)
        {
          renum_i = i;
        }
      else
        {
          // Search for the local index corresponding to the sequential index i
          // If i is a virtual item we take it anyway: vecteur_i will have its virtual space up to date
          // If i does not exist on this processor, we set -1
          renum_i = -1;
          // Renum_colonnes[j] is the index in the sequential reference vector of item j on this processor
          const IntVect& renum_colonnes = find_renum_vector(md_colonnes, id);
          const int sz = renum_colonnes.size_totale();
          for (int j = 0; j < sz; j++)
            {
              if (renum_colonnes[j] == i)
                {
                  renum_i = j;
                  break;
                }
            }
        }
      if (renum_i >= 0)
        vecteur_i[renum_i] = 1.;
      matrice.multvect(vecteur_i, tmp);
      Nom msg2(msg);
      msg2 += " produit avec base ";
      msg2 += Nom(i);
      verifier(msg2, tmp);
      if (renum_i >= 0)
        vecteur_i[renum_i] = 0.;
    }
}
 
void Debog_Pb::verifier_Mat_elems(const char* const msg, const Matrice_Base& la_matrice)
{
  // Search for the md_vector of elements
 
  const int n = renum_id_.size();
  int i;
  for (i = 0; i < n; i++)
    {
      if (renum_id_[i] == "ELEM")
        break;
    }
  if (i == n)
    {
      Cerr << "Error in Debog_Pb::verifier_Mat_elems no registered MD_Vector for ELEM type" << finl;
      exit();
    }
  const MD_Vector& md = known_md_[i];
  verifier_matrice(msg, la_matrice, md, md);
}
 
////////////////////////////////////////////////////////////////
 
Sortie& Debog_Pb_Wrapper::printOn(Sortie& os) const
{
  Process::exit();
  return os;
}
 
Entree& Debog_Pb_Wrapper::readOn(Entree& is)
{
  Process::exit();
  return is;
}
 
Entree& Debog_Pb_Wrapper::interpreter(Entree& is)
{
  Nom nom_pb, nom_fic, nom_fic_faces, seuil;
  is >> nom_pb;
  int mode_db;
  is >> nom_fic >> nom_fic_faces >> seuil >> mode_db;
 
  Nom n("Debog_pb debog Lire debog { nom_pb ");
  n += nom_pb;
  n += " fichier_domaine ";
  n += nom_fic;
  n += " fichier_faces ";
  n += nom_fic_faces;
  n += " fichier_debog DEBOG";
  n += " seuil_absolu ";
  n += Nom(seuil);
  n += " seuil_relatif ";
  n += Nom(seuil);
  n += " mode_db ";
  if (mode_db == 0) n += "write";
  else n += "read";
  n += "  }";
 
  if (je_suis_maitre()) Cerr << "Debog_Pb_Wrapper: convert to new syntax: " << finl << n << finl;
 
  EChaine echaine(n);
  Interprete_bloc::interprete_courant().interpreter_bloc(echaine, Interprete_bloc::BLOC_EOF, 0 /* verifie_sans_interpreter */);
  return is;
}