Equation_IBM_proto.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 <Discretisation_base.h>
#include <Equation_IBM_proto.h>
#include <Schema_Temps_base.h>
#include <Source_PDF_base.h>
#include <Operateur_base.h>
#include <Probleme_base.h>
#include <Equation_base.h>
#include <Operateur.h>
#include <Param.h>
#include <Nom.h>
#include <Perf_counters.h>
 
void Equation_IBM_proto::set_param_ibm_proto(Param& param) const
{
  param.ajouter("correction_variable_initiale",&correction_variable_initiale_,Param::OPTIONAL);
}
 
Entree& Equation_IBM_proto::readOn_ibm_proto(Entree& is, Equation_base& eq)
{
  eq_IBM_ = eq;
 
  int i_source = -1;
  int nb_source_pdf = 0;
  int nb_source  = eq.sources().size();
  Cerr << "nb_source : " << nb_source << finl;
  for(int i = 0; i<nb_source; i++)
    {
      //Cerr << (sources())[i]->que_suis_je() << ", " << (sources())[i]->que_suis_je().find("gne") << finl;
      if ((eq.sources())[i]->que_suis_je().find("Source_PDF") > -1)
        {
          i_source = i;
          nb_source_pdf++;
        }
    }
  if (nb_source_pdf>1)
    {
      Cerr<<"(IBM) More than one Source_PDF_base detected."<<finl;
      Cerr<<"(IBM) Case not supported."<<finl;
      Cerr<<"Aborting..."<<finl;
      abort();
    }
  i_source_pdf_ = i_source;
  if (i_source_pdf_ != -1)
    {
      Cerr<<"(IBM) Immersed Interface with Source_PDF_base for equation : "<< eq_IBM_->le_nom()<<finl;
      eq.sources()[i_source_pdf_]->set_description((Nom)"Sum of the PDF source term on the obstacle (~= force/power induced by the obstacle under some assumptions)");
      Nom the_suffix = eq_IBM_->le_nom().getSuffix(eq_IBM_->probleme().le_nom());
      eq.sources()[i_source_pdf_]->set_fichier((Nom)the_suffix+"_Bilan_term_induced_by_obstacle");
    }
 
  // XXX Elie Saikali
  // ok for temperature for now (test case Comte_Bellot_EF_IBC in trio ! to be corrected later Michel !)
  if (i_source_pdf_ == -1 && eq_IBM_->que_suis_je().debute_par("Navier_Stokes"))
    {
      Cerr << "No PDF source term read in your equation " << eq_IBM_->que_suis_je() << " !!!"<< finl;
      Cerr << "Either add this source term or use a non_IBM problem !!!"<< finl;
      Process::exit();
    }
 
  return is;
}
 
bool Equation_IBM_proto::initTimeStep_ibm_proto(double ddt)
{
  if (i_source_pdf_ == -1) return false;
  Cerr<<"(IBM) initTimeStep_IBM: update of immersed values for equation : "<< eq_IBM_->le_nom() <<finl;
  Source_PDF_base& src = dynamic_cast<Source_PDF_base&>((eq_IBM_->sources())[i_source_pdf_].valeur());
  src.calculer_variable_imposee();
 
  if ( src.get_modele().get_PDF_mobile() && src.get_modele().is_vitesse_PDF_donnee() )
    {
      Cerr<<"(IBM) PDF_mobile : updating the shape following shape velocity"<<finl;
      const Probleme_base& pb = eq_IBM_->probleme();
      const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
      const DoubleTab& coords = le_dom_dis.domaine().coord_sommets();
      Domaine_VF& le_dom_VF = ref_cast_non_const(Domaine_VF, pb.domaine_dis());
      double delta_t = eq_IBM_->schema_temps().pas_de_temps();
      double temps= eq_IBM_->schema_temps().temps_courant();
 
      // Boundary displacement
      PDF_model& my_model = ref_cast_non_const(PDF_model, src.get_modele());
      my_model.affecter_vitesse_shape_IBM(le_dom_VF, coords, temps);
      DoubleTab& vitesse = ref_cast_non_const(DoubleTab, src.get_modele().get_vitesse_shape_IBM()) ;
      double raid = my_model.raid();
      src.update_elem_IBM(vitesse, delta_t,raid);
    }
  return true;
}
 
DoubleTab& Equation_IBM_proto::derivee_en_temps_inco_ibm_proto(DoubleTab& derivee)
{
  if (i_source_pdf_ == -1) return derivee;
  Cerr<<"(IBM) Immersed Interface: Dirichlet value in Equation for PDF (if any)."<<finl;
  Source_PDF_base& src = dynamic_cast<Source_PDF_base&>((eq_IBM_->sources())[i_source_pdf_].valeur());
 
  // PDF IB value term : -rho/delta_t  ksi_gamma/epsilon U_gamma
  DoubleTab secmem_pdf(derivee);
  src.calculer_pdf(secmem_pdf);
  derivee -= secmem_pdf;
  derivee.echange_espace_virtuel();
 
  // Time term : -rho/delta_t ksi_gamma Un
  int pdf_bilan = src.get_modele().pdf_bilan();
  if (pdf_bilan == 1|| pdf_bilan == 2)
    {
      int i_traitement_special = 101;
      int nb_comp = derivee.dimension(1);
      if (nb_comp == Objet_U::dimension)
        {
          //vector; NS
          // Navier_Stokes_std& eq_NS = ref_cast_non_const(Navier_Stokes_std, eq_IBM);
          if (eq_IBM_->nombre_d_operateurs() > 1)
            {
              // if (eq_NS.vitesse_pour_transport().le_nom()=="rho_u") i_traitement_special = 1;
              i_traitement_special = 101;
            }
        }
      else if (nb_comp == 1)
        {
          //scalar; time term in rho
          i_traitement_special = 1;
        }
      else
        {
          Cerr << "Equation_IBM_proto: for scalar or vector only; dim = " <<nb_comp<< finl;
          Process::exit();
        }
 
      DoubleTrav secmem_pdf_time(derivee);
      src.calculer(secmem_pdf_time, i_traitement_special);
      secmem_pdf += secmem_pdf_time;
    }
  else if(pdf_bilan != 0 )
    {
      Cerr<<"Source_PDF_EF: Model pdf_bilan must be 0; 1 or 2 only"<<finl;
      Process::exit();
    }
 
  // Save secmem_pdf
  secmem_pdf.echange_espace_virtuel();
  src.set_sec_mem_pdf(secmem_pdf);
 
  return derivee;
}
 
void Equation_IBM_proto::preparer_calcul_ibm_proto()
{
  if (i_source_pdf_ == -1) return;
 
  Source_PDF_base& src = dynamic_cast<Source_PDF_base&>((eq_IBM_->sources())[i_source_pdf_].valeur());
  src.updateChampRho();
  DoubleTab coeff;
  coeff = src.compute_coeff_matrice();
  coeff.echange_espace_virtuel();
  Cerr<<"(IBM) preparer_calcul_IBM: Min/max coeff : "<< mp_min_vect(coeff) << " " << mp_max_vect(coeff) <<finl;
  champ_coeff_pdf_som_ = coeff;
  modify_initial_variable_ibm_proto(eq_IBM_->inconnue().valeurs());
}
 
void Equation_IBM_proto::modify_initial_variable_ibm_proto(DoubleTab& variable)
{
  if (correction_variable_initiale_==1)
    {
      Cerr<<"(IBM) Immersed Interface: modified initial variable."<<finl;
      const Source_PDF_base& src = dynamic_cast<Source_PDF_base&>((eq_IBM_->sources())[i_source_pdf_].valeur());
      src.correct_variable(champ_coeff_pdf_som_, variable);
      variable.echange_espace_virtuel();
    }
}
 
// adds the contributions of operators and sources
void Equation_IBM_proto::assembler_ibm_proto(Matrice_Morse& matrice, const DoubleTab& inco, DoubleTab& resu)
{
  const double rhoCp = eq_IBM_->get_time_factor();
  int size_s = eq_IBM_->sources().size();
 
  // Verification test of the contribuer_a_avec method
  for (int op=0; op<eq_IBM_->nombre_d_operateurs(); op++)
    eq_IBM_->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: convective flux Op in VDF:
  // Op=T*u*S and A=-d(T*u*S)/dT=-u*S
  const Discretisation_base::type_calcul_du_residu& type_codage=eq_IBM_->probleme().discretisation().codage_du_calcul_du_residu();
  if (type_codage==Discretisation_base::VIA_CONTRIBUER_AU_SECOND_MEMBRE)
    {
      if ( eq_IBM_->probleme().discretisation().que_suis_je() == "EF")
        {
          // We compute sum(residual) via contribuer_au_second_membre (typically non-implicitized BCs)
          // This approach requires coding 3 methods (contribuer_a_avec, contribuer_au_second_membre and ajouter for the explicit case)
          if (!is_IBM())
            eq_IBM_->sources().contribuer_a_avec(inco, matrice);
          else
            {
              for (int i = 0; i < size_s; i++)
                {
                  if (eq_IBM_->sources()(i).valeur().que_suis_je().find("Source_PDF") <= -1)
                    {
                      const Source_base& src_base = eq_IBM_->sources()(i).valeur();
                      src_base.contribuer_a_avec(inco, matrice);
                    }
                }
            }
          statistics().end_count(STD_COUNTERS::matrix_assembly);
          if (!is_IBM())
            eq_IBM_->sources().ajouter(resu);
          else
            {
              for (int i = 0; i < size_s; i++)
                if (eq_IBM_->sources()(i).valeur().que_suis_je().find("Source_PDF") <= -1)
                  eq_IBM_->sources()(i).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 < eq_IBM_->nombre_d_operateurs(); op++)
            {
              eq_IBM_->operateur(op).l_op_base().contribuer_a_avec(inco, matrice);
              eq_IBM_->operateur(op).l_op_base().contribuer_au_second_membre(resu);
            }
        }
      else
        {
          Cerr << "VIA_CONTRIBUER_AU_SECOND_MEMBRE not coded for " << eq_IBM_->que_suis_je() << ":assembler" << finl;
          Cerr << "with discretisation " <<  eq_IBM_->probleme().discretisation().que_suis_je() << "" << finl;
          Process::exit();
        }
      // // We compute sum(residual) via contribuer_au_second_membre (typically non-implicitized BCs)
      // // This approach requires coding 3 methods (contribuer_a_avec, contribuer_au_second_membre and ajouter for the explicit case)
      // sources().contribuer_a_avec(inco,matrice);
      // sources().ajouter(resu);
      // 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 compute sum(residual) via sum(operator)+sources+A*Inco(n)
      // This approach requires coding only two methods (contribuer_a_avec and ajouter)
      // So slightly more expensive in computation time but less code to write/maintain
      for (int op=0; op<eq_IBM_->nombre_d_operateurs(); op++)
        {
          Matrice_Morse mat(matrice);
          mat.get_set_coeff() = 0.0;
          eq_IBM_->operateur(op).l_op_base().contribuer_a_avec(inco, mat);
          if (op == 1) mat *= rhoCp; // the derivative is multiplied by rhoCp for convection
          matrice += mat;
          statistics().end_count(STD_COUNTERS::matrix_assembly);
          {
            DoubleTab resu_tmp(resu);
            resu_tmp = 0.;
            eq_IBM_->operateur(op).ajouter(inco, resu_tmp);
            if (op == 1) resu_tmp *= rhoCp;
            resu += resu_tmp;
          }
          statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
        }
      if (! is_IBM() )
        eq_IBM_->sources().contribuer_a_avec(inco,matrice);
      else
        {
          for (int i = 0; i < size_s; i++)
            if (eq_IBM_->sources()(i).valeur().que_suis_je().find("Source_PDF") <= -1)
              {
                const Source_base& src_base = eq_IBM_->sources()(i).valeur();
                src_base.contribuer_a_avec(inco, matrice);
              }
        }
      statistics().end_count(STD_COUNTERS::matrix_assembly);
      if (!is_IBM())
        eq_IBM_->sources().ajouter(resu);
      else
        {
          for (int i = 0; i < size_s; i++)
            if (eq_IBM_->sources()(i).valeur().que_suis_je().find("Source_PDF") <= -1)
              eq_IBM_->sources()(i).ajouter(resu);
        }
      statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
      matrice.ajouter_multvect(inco, resu); // Adding A*Inco(n)
      // PL (11/04/2018): We would like to compute the source contribution 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 Convection_Diffusion_Temperature_IBM::assembler for " << (int)type_codage << finl;
      Process::exit();
    }
 
  // to avoid infinite PDF terms when adding A*Inco(n)
  if ( is_IBM() )
    for (int i = 0; i < size_s; i++)
      {
        if (eq_IBM_->sources()(i).valeur().que_suis_je().find("Source_PDF") > -1)
          {
            const Source_base& src_base = eq_IBM_->sources()(i).valeur();
            src_base.contribuer_a_avec(inco,matrice);
          }
      }
  // add PDF source with the correct sign
  derivee_en_temps_inco_ibm_proto(resu);
}