Eq_rayo_semi_transp.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 <Eq_rayo_semi_transp.h>
#include <Pb_rayo_semi_transp.h>
#include <Matrice_Morse_Sym.h>
#include <Champ_Uniforme.h>
#include <Matrice_Bloc.h>
#include <Discret_Thyd.h>
#include <Fluide_base.h>
#include <EChaine.h>
#include <Param.h>
 
Implemente_instanciable(Eq_rayo_semi_transp, "Eq_rayo_semi_transp", Equation_base);
 
bool Eq_rayo_semi_transp::initTimeStep(double dt)
{
  schema_temps().set_dt() = dt;
  return Equation_base::initTimeStep(dt);
}
 
void Eq_rayo_semi_transp::resoudre(double temps)
{
  rayo_solv_->resoudre(temps);
}
 
bool Eq_rayo_semi_transp::resoudre()
{
  rayo_solv_->resoudre(schema_temps().temps_courant() + schema_temps().pas_de_temps());
  return true;
}
 
Sortie& Eq_rayo_semi_transp::printOn(Sortie& s) const { return s << que_suis_je() << finl; }
 
Entree& Eq_rayo_semi_transp::readOn(Entree& is) { return Equation_base::readOn(is); }
 
void Eq_rayo_semi_transp::set_param(Param& param) const
{
  param.ajouter_non_std("conditions_limites|boundary_conditions", (this), Param::REQUIRED);
  param.ajouter_non_std("solveur", (this), Param::REQUIRED);
}
 
int Eq_rayo_semi_transp::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  int retval = 1;
  if (mot == "conditions_limites|boundary_conditions")
    {
      lire_cl(is);
      verif_Cl();
    }
  else if (mot == "solveur")
    {
      Cerr << "Reading and typing of the radiation equation solver :" << finl;
      Nom nom_solveur("Solv_");
      Nom type_solv_sys;
      is >> type_solv_sys;
      nom_solveur += type_solv_sys;
      Cerr << "Name of the radiation equation solver : " << nom_solveur << finl;
      solveur_.typer(nom_solveur);
      is >> solveur_.valeur();
      solveur_.nommer("solveur_irradiance");
    }
  else
    retval = -1;
 
  return retval;
}
 
void Eq_rayo_semi_transp::completer()
{
  const Domaine_dis_base& dom_dis = domaine_dis();
  for (int i = 0; i < dom_dis.nb_front_Cl(); i++)
    {
      const Frontiere_dis_base& la_fr_dis = dom_dis.frontiere_dis(i);
      le_dom_Cl_dis->les_conditions_limites(i)->associer_fr_dis_base(la_fr_dis);
    }
 
  // typing of the diffusion operator
  Cerr << "Reading and typing of the diffusion operator of equation " << que_suis_je() << finl;
 
  if (sub_type(Fluide_base, fluide()))
    if (fluide().is_rayo_semi_transp())
      {
        if (fluide().is_longueur_rayo_discretised())
          terme_diffusif_.associer_diffusivite(fluide().longueur_rayo());
        else
          {
            Cerr << "Error in Eq_rayo_semi_transp::completer." << finl;
            Cerr << "You may not have discretized the problem of type Pb_Couple_rayo_semi_transp." << finl;
            Process::exit();
          }
      }
    else
      {
        Cerr << "Error : the radiative properties of the incompressible fluid have not" << finl;
        Cerr << "been defined while a semi transparent radiation problem is used." << finl;
        Cerr << "The fields kappa and indice which respectively define the absoption coefficient" << finl;
        Cerr << "and the refraction index must be added to the fluid properties." << finl;
        Process::exit();
      }
  else
    {
      Cerr << "Error while reading the Radiation equation. Your fluid is of type " << fluide().que_suis_je() << finl;
      Cerr << "Currently only fluid of type Fluide_base can be considered with the semi transparent radiation model." << finl;
      Process::exit();
    }
 
  EChaine diff("{ }");
  diff >> terme_diffusif_;
  terme_diffusif_->associer_diffusivite(fluide().longueur_rayo());
  terme_diffusif_.completer();
  terme_diffusif_->dimensionner(la_matrice_);
 
  Equation_base::completer();
 
  // The matrix is assembled once and for all at the start of the computation
  // XXX Warning: this is only valid if kappa is constant over time
  rayo_solv_->assembler_matrice();
}
 
/*! @brief Associates a physical medium with the equation.
 *
 * @param un_milieu The physical medium to associate with the equation.
 */
void Eq_rayo_semi_transp::associer_milieu_base(const Milieu_base& un_milieu)
{
  if (sub_type(Fluide_base, un_milieu))
    if (fluide().is_rayo_semi_transp())
      {
        const Fluide_base& un_fluide = ref_cast(Fluide_base, un_milieu);
        associer_fluide(un_fluide);
        le_fluide_ = un_fluide;
      }
    else
      {
        Cerr << "Error : the radiative properties of the incompressible fluid have not" << finl;
        Cerr << "been defined while a semi transparent radiation problem is used." << finl;
        Cerr << "The fields kappa and indice which respectively define the absoption coefficient" << finl;
        Cerr << "and the refraction index must be added to the fluid properties." << finl;
        Process::exit();
      }
  else
    {
      Cerr << "Error while reading the Radiation equation. Your fluid is of type " << fluide().que_suis_je() << finl;
      Cerr << "Currently only fluid of type Fluide_base can be considered with the semi transparent radiation model." << finl;
      Process::exit();
    }
}
 
/*! @brief Returns the physical medium of the equation (the Fluide_base upcast to Milieu_base).
 *
 * @return The Fluide_base of the equation upcast to Milieu_base.
 */
const Milieu_base& Eq_rayo_semi_transp::milieu() const
{
  if (!le_fluide_)
    {
      Cerr << "You forgot to associate the fluid to the problem named " << probleme().le_nom() << finl;
      Process::exit();
    }
  return le_fluide_.valeur();
}
 
/*! @brief Returns the physical medium of the equation (the Fluide_base upcast to Milieu_base).
 *
 *     (non-const version)
 *
 * @return The Fluide_base of the equation upcast to Milieu_base.
 */
Milieu_base& Eq_rayo_semi_transp::milieu()
{
  if (!le_fluide_)
    {
      Cerr << "You forgot to associate the fluid to the problem named " << probleme().le_nom() << finl;
      Process::exit();
    }
  return le_fluide_.valeur();
}
 
/*! @brief Returns the operator at the specified index: returns terme_diffusif if i = 0, exits if i > 0.
 *
 *     (const version)
 *
 * @param i Index of the operator to return.
 * @return The specified operator.
 * @throws The equation has at most 1 operator.
 */
const Operateur& Eq_rayo_semi_transp::operateur(int i) const
{
  switch(i)
    {
    case 0:
      return terme_diffusif_;
    default:
      Cerr << "Error for Eq_rayo_semi_transp::operateur(int i)" << finl;
      Cerr << "Eq_rayo_semi_transp has " << nombre_d_operateurs() << " operators " << finl;
      Cerr << "and you are trying to access the " << i << " th one." << finl;
      Process::exit();
    }
  return terme_diffusif_;
}
 
/*! @brief Returns the operator at the specified index: returns terme_diffusif if i = 0, exits if i > 0.
 *
 *     (non-const version)
 *
 * @param i Index of the operator to return.
 * @return The specified operator.
 * @throws The equation has at most 1 operator.
 */
Operateur& Eq_rayo_semi_transp::operateur(int i)
{
  switch(i)
    {
    case 0:
      return terme_diffusif_;
    default:
      Cerr << "Error for Eq_rayo_semi_transp::operateur(int i)" << finl;
      Cerr << "Eq_rayo_semi_transp has " << nombre_d_operateurs() << " operators " << finl;
      Cerr << "and you are trying to access the " << i << " th one." << finl;
      Process::exit();
    }
  return terme_diffusif_;
}
 
void Eq_rayo_semi_transp::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());
}
 
void Eq_rayo_semi_transp::discretiser()
{
  // Discretisation of the radiation equation
  const Discret_Thyd& dis = ref_cast(Discret_Thyd, discretisation());
  Cerr << "Radiation equation discretisation" << finl;
  dis.discretiser_champ("temperature", domaine_dis(), "irradiance", "w/m2", 1, 1 /* one slot */, schema_temps().temps_courant(), irradiance_);
  champs_compris_.ajoute_champ(irradiance_);
 
  Equation_base::discretiser();
 
  // typing Rayo_semi_transp_solver_base
  Nom discr = dis.que_suis_je(), type = "Rayo_semi_transp_solver_";
  if (discr == "VEFPreP1B") discr = "VEF";
  type += discr;
  rayo_solv_.typer(type);
  rayo_solv_->associer_equation_rayo(*this);
}
 
/*! @brief Returns the discretisation associated with the equation.
 *
 * @return The discretisation associated with the equation.
 * @throws If no problem is associated.
 */
const Discretisation_base& Eq_rayo_semi_transp::discretisation() const
{
  return pb_rayo_semi_transp_->discretisation();
}
 
void Eq_rayo_semi_transp::associer_pb_base(const Probleme_base& pb)
{
  Equation_base::associer_pb_base(pb);
  pb_rayo_semi_transp_ = ref_cast(Pb_rayo_semi_transp, pb);
  associer_sch_tps_base(pb.schema_temps());
}
 
void Eq_rayo_semi_transp::Mat_Morse_to_Mat_Bloc(Matrice& matrice_tmp)
{
  const int n1 = rayo_solv_->nb_colonnes_tot();
  const int n2 = rayo_solv_->nb_colonnes();
 
  Matrice_Bloc& matrice = ref_cast(Matrice_Bloc, matrice_tmp.valeur());
  Matrice_Morse& MBrr = ref_cast(Matrice_Morse, matrice.get_bloc(0, 0).valeur());
  Matrice_Morse& MBrv = ref_cast(Matrice_Morse, matrice.get_bloc(0, 1).valeur());
 
  auto& tab1RR = MBrr.get_set_tab1();
  auto& tab2RR = MBrr.get_set_tab2();
  auto& coeffRR = MBrr.get_set_coeff();
  auto& tab1RV = MBrv.get_set_tab1();
  auto& tab2RV = MBrv.get_set_tab2();
  auto& coeffRV = MBrv.get_set_coeff();
 
  DoubleTab ligne_tmp(n1);
  for (int i = 0; i < n2; i++)
    {
      int k;
      // Copy the first block of the matrix into a temporary array:
      //      ligne_tmp = 0;
      for (k = la_matrice_.get_tab1()(i) - 1; k < la_matrice_.get_tab1()(i + 1) - 1; k++)
        ligne_tmp(la_matrice_.get_tab2()(k) - 1) = la_matrice_.get_coeff()(k);
 
      // Fill in the real part of the matrix
      for (k = tab1RR(i) - 1; k < tab1RR(i + 1) - 1; k++)
        coeffRR[k] = ligne_tmp(tab2RR[k] - 1);
 
      // Fill in the virtual part
      for (k = tab1RV(i) - 1; k < tab1RV(i + 1) - 1; k++)
        coeffRV[k] = ligne_tmp(n2 + tab2RV[k] - 1);
    }
}
 
void Eq_rayo_semi_transp::dimensionner_Mat_Bloc_Morse_Sym(Matrice& matrice_tmp)
{
  const int n1 = rayo_solv_->nb_colonnes_tot();
  const int n2 = rayo_solv_->nb_colonnes();
 
  int iligne;
  const auto& tab1 = la_matrice_.get_set_tab1();
  const auto& tab2 = la_matrice_.get_set_tab2();
 
  matrice_tmp.typer("Matrice_Bloc");
  Matrice_Bloc& matrice = ref_cast(Matrice_Bloc, matrice_tmp.valeur());
  matrice.dimensionner(1, 2);
  matrice.get_bloc(0, 0).typer("Matrice_Morse_Sym");
  matrice.get_bloc(0, 1).typer("Matrice_Morse");
 
  Matrice_Morse_Sym& MBrr = ref_cast(Matrice_Morse_Sym, matrice.get_bloc(0, 0).valeur());
  Matrice_Morse& MBrv = ref_cast(Matrice_Morse, matrice.get_bloc(0, 1).valeur());
  MBrr.dimensionner(n2, 0);
  MBrv.dimensionner(n2, 0);
 
  auto& tab1RR = MBrr.get_set_tab1();
  auto& tab2RR = MBrr.get_set_tab2();
  auto& tab1RV = MBrv.get_set_tab1();
  auto& tab2RV = MBrv.get_set_tab2();
 
  IntVect compteur_MBrr(n2);
  IntVect compteur_MBrv(n2);
  compteur_MBrr = 0;
  compteur_MBrv = 0;
 
  // Iterate over the rows of la_matrice to count the non-zero
  // elements of each row
  int jcolonne;
  for (iligne = 0; iligne < n2; iligne++)
    {
      int k;
      for (k = tab1(iligne) - 1; k < tab1(iligne + 1) - 1; k++)
        {
          jcolonne = tab2(k) - 1;
          if (jcolonne < n2)
            {
              // the corresponding element is in the RR part of la_matrice
              if ((jcolonne >= iligne) && (jcolonne < n2))
                {
                  // the corresponding element is located above the diagonal of la_matrice
                  compteur_MBrr(iligne)++;
                }
            }
          else
            {
              // the corresponding element is in the RV part of la_matrice
              compteur_MBrv(iligne)++;
            }
        }
    }
 
  // Fill tab1RR and tab1RV
  tab1RR(0) = 1;
  tab1RV(0) = 1;
  for (int i = 0; i < n2; i++)
    {
      tab1RR(i + 1) = compteur_MBrr(i) + tab1RR(i);
      tab1RV(i + 1) = compteur_MBrv(i) + tab1RV(i);
    }
  // Size tab2RR and tab2RV
  MBrr.dimensionner(n2, tab1RR(n2) - 1);
  MBrv.dimensionner(n2, n1 - n2, tab1RV(n2) - 1);
 
  // Fill tab2RR and tab2RV
  int compteurRR, compteurRV;
  for (iligne = 0; iligne < n2; iligne++)
    {
      int k;
      compteurRR = tab1RR(iligne) - 1;
      compteurRV = tab1RV(iligne) - 1;
      for (k = tab1(iligne) - 1; k < tab1(iligne + 1) - 1; k++)
        {
          jcolonne = tab2(k) - 1;
          if (jcolonne < n2)
            {
              // the corresponding element is in the RR part of la_matrice
              if ((jcolonne >= iligne) && (jcolonne < n2))
                {
                  // the corresponding element is located above the diagonal of la_matrice
                  tab2RR(compteurRR) = tab2(k);
                  compteurRR++;
                }
            }
          else
            {
              // the corresponding element is in the RV part of la_matrice
              tab2RV(compteurRV) = tab2(k) - n2;
              compteurRV++;
            }
        }
    }
}