Fermeture_Thermo_base.h

/****************************************************************************
* Copyright (c) 2025, 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.
*
*****************************************************************************/
 
#ifndef Fermeture_Thermo_base_included
#define Fermeture_Thermo_base_included
 
#include <Objet_U.h>
#include <Param.h>
#include <Champ_Fonc_base.h>
#include <Potentiel_Chimique_base.h>
#include <Table.h>
 
/*! @brief classe Fermeture_Thermo_base
 *
 *      Classe de base des fermetures thermodynamiques, c'est-à-dire des modèles de
 *      potentiels chimiques choisis pour la forme de l'énergie libre de Gibbs.
 *
 * @sa Equation_base
 */
class Fermeture_Thermo_base : public Objet_U
{
 
  Declare_base_sans_constructeur(Fermeture_Thermo_base);
 
public:
 
  Fermeture_Thermo_base();
 
  // Mobilité (dégénérée ou non)
  DoubleTab mobilite_func_c_cst(const DoubleTab&) const;
  DoubleTab mobilite_func_c_defaut(const DoubleTab&) const;
  DoubleTab mobilite_func_c_general(const DoubleTab&) const;
  DoubleTab mobilite_func_auto_diffusion(const DoubleTab&) const;
 
  DoubleTab mobilite(const DoubleTab& d) const { return (this->*mobilite_func_c)(d); } // Renvoie le champ entier sur le domaine sous forme de DoubleTab
  DoubleTab (Fermeture_Thermo_base::*mobilite_func_c)(const DoubleTab&) const;
 
  // coefficient Kappa
  inline DoubleVect kappa() const { return kappa_i_; }
 
  // Potentiel chimique
  inline Potentiel_Chimique_base& potentiel_chimique() { return mu_.valeur(); }
  inline const Potentiel_Chimique_base& potentiel_chimique() const { return mu_.valeur(); }
 
  // coefficient Beta
  inline DoubleVect beta() const { return beta_; }
 
  inline int nb_parametres_d_ordre() const { return nb_parametres_d_ordre_; } // On appelle ce paramètre nb_equations_ et non nb_phases, car il correspond à
  inline int type_mobilite_auto_diffusion() const { return type_mobilite_auto_diffusion_; }
  inline double mult_mobilite() const { return mult_mobilite_; }
 
  inline double molarVolume() const { return molarVolume_; }
  inline double temperature() const { return temperature_; }
 
protected :
 
  int nb_parametres_d_ordre_ = -123;
 
  DoubleVect mob_i_, kappa_i_, beta_;
 
  OWN_PTR(Potentiel_Chimique_base) mu_;
 
  Table mobilite_forme_expr_;
 
  int type_mobilite_auto_diffusion_ = -123;
  int type_kappa_rotation_ = -123;
  double mult_mobilite_ = -123.;
 
  double kappa_ref_ = -123., angle_kappa_ = -123.;
  double diagonal_coeff_ = -123.;
  double temperature_ = -123., molarVolume_ = -123.;
 
  DoubleVect coeff_auto_diffusion_;
 
};
 
inline DoubleTab Fermeture_Thermo_base::mobilite_func_c_cst(const DoubleTab& c) const
{
  DoubleTab mobiliteMatrix_tab(c.dimension(0), nb_parametres_d_ordre_*nb_parametres_d_ordre_);
  // Cerr << "[Closure] Mobility: Constant function is used." << finl;
  // Cerr << "[Closure] Concentration field is sized " << c.dimension(0) << " x " << c.line_size() << finl;
  for (int i = 0; i < c.dimension(0); i++)
    for (int j = 0; j < nb_parametres_d_ordre_*nb_parametres_d_ordre_; j++)
      mobiliteMatrix_tab(i, j) = mult_mobilite_ * mob_i_(j);
  return mobiliteMatrix_tab;
}
 
inline DoubleTab Fermeture_Thermo_base::mobilite_func_c_defaut(const DoubleTab& c) const
{
  assert(c.line_size() == nb_parametres_d_ordre_);
  DoubleTab mobiliteMatrix_tab(c.dimension(0), nb_parametres_d_ordre_*nb_parametres_d_ordre_);
  // Cerr << "Mobility: Constant function is used." << finl;
  for (int i = 0; i < c.dimension(0); i++)
    for (int j = 0; j < nb_parametres_d_ordre_; j++)
      for (int k = 0; k < nb_parametres_d_ordre_; k++)
        mobiliteMatrix_tab(i, k+nb_parametres_d_ordre_*j) = std::max(mult_mobilite_ * mob_i_(k+nb_parametres_d_ordre_*j) * c(i, k) * (1.0 - c(i, k)), 0.);
  return mobiliteMatrix_tab;
}
 
inline DoubleTab Fermeture_Thermo_base::mobilite_func_c_general(const DoubleTab& c) const
{
  assert(c.line_size() == 1);
  DoubleTab mobiliteMatrix_tab(c.dimension(0), 1);
  Cerr << "Mobility: General form is used." << finl;
  for (int i = 0; i < c.dimension(0); i++)
    mobiliteMatrix_tab(i, 0) = std::max(mult_mobilite_ * mob_i_(0) * mobilite_forme_expr_.val(c(i, 0), 0), 0.);
  return mobiliteMatrix_tab;
}
 
inline DoubleTab Fermeture_Thermo_base::mobilite_func_auto_diffusion(const DoubleTab& c) const
{
  Cerr << "Mobility: Auto-diffusion is used." << finl;
  DoubleTab temp_c(c.dimension(0), c.line_size() + 1);
  DoubleTab temp_mobilite(c.dimension(0), nb_parametres_d_ordre_, nb_parametres_d_ordre_); //define temp to avoid resize here
  DoubleTab mobiliteMat = temp_mobilite;
  mobiliteMat = 0;
  DoubleTab mobilite_ij(temp_mobilite);
  DoubleTab somme_c(c.dimension(0));
  somme_c = 0;
  DoubleTab mobiliteMatrix_tab(c.dimension(0), nb_parametres_d_ordre_ * nb_parametres_d_ordre_);
 
  // creation d'un tableau temp_c qui contiendra toutes les concentrations des n composants non pas les n-1 seulements (n-1 eqs).
  //La derniere colonne correspondra a 1-somme sur (n-1) des concentrations
 
  for (int l = 0; l < c.dimension(0); l++)
    for (int k = 0; k < c.line_size(); k++)
      {
        temp_c(l, k) = c(l, k);
        somme_c(l) += c(l, k);
        temp_c(l, c.line_size()) = 1.0 - somme_c(l);
      }
 
  // calcul de la mobilite mobilite dans un doubleTab (n1,n2,n3)
  for (int k = 0; k < nb_parametres_d_ordre_ + 1; k++)
    {
      for (int j = 0; j < nb_parametres_d_ordre_; j++)
        for (int i = 0; i < nb_parametres_d_ordre_; i++)
          for (int l = 0; l < c.dimension(0); l++)
            {
              if (j == k)
                mobilite_ij(l, i, j) = 1 - temp_c(l, j);
              else
                mobilite_ij(l, i, j) = -temp_c(l, j);
              if (i == k)
                mobilite_ij(l, i, j) *= 1 - temp_c(l, i);
              else
                mobilite_ij(l, i, j) *= -temp_c(l, i);
 
              mobilite_ij(l, i, j) *= temp_c(l, k) * coeff_auto_diffusion_(k);
            }
      mobiliteMat += mobilite_ij;
    }
 
  //retour vers un DoubleTab (n1,n2)
  for (int l = 0; l < c.dimension(0); l++)
    for (int j = 0; j < nb_parametres_d_ordre_; j++)
      for (int i = 0; i < nb_parametres_d_ordre_; i++)
        mobiliteMatrix_tab(l, i + nb_parametres_d_ordre_ * j) = mobiliteMat(l, i, j);
 
  mobiliteMatrix_tab *= molarVolume_ * 8.314 * temperature_;
 
  return mobiliteMatrix_tab;
}
 
#endif /* Fermeture_Thermo_base_included */