Milieu_Phase_field.cpp

/****************************************************************************
* 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.
*
*****************************************************************************/
 
#include <Milieu_Phase_field.h>
#include <Pb_Phase_field.h>
#include <Discret_Thyd.h>
#include <EChaine.h>
#include <Param.h>
 
Implemente_instanciable(Milieu_Phase_field, "Milieu_Phase_field", Fluide_Incompressible);
 
// XD approx_boussinesq objet_lecture nul NO_BRACE different mass density formulation are available depending if the
// XD_CONT Boussinesq approximation is made or not
// XD attr yes_or_no chaine(into=["oui","non"]) yes_or_no REQ To use or not the Boussinesq approximation.
// XD attr bloc_bouss bloc_boussinesq bloc_bouss REQ to choose the rho formulation
 
// XD bloc_boussinesq objet_lecture nul BRACE choice of rho formulation
// XD attr probleme ref_Pb_base probleme OPT Name of problem.
// XD attr rho_1 floattant rho_1 OPT value of rho
// XD attr rho_2 floattant rho_2 OPT value of rho
// XD attr rho_fonc_c bloc_rho_fonc_c rho_fonc_c_ OPT to use for define a general form for rho
 
// XD bloc_rho_fonc_c objet_lecture nul NO_BRACE if rho has a general form
// XD attr Champ_Fonc_Fonction chaine(into=["Champ_Fonc_Fonction"]) Champ_Fonc_Fonction OPT Champ_Fonc_Fonction
// XD attr problem_name ref_Pb_base problem_name OPT Name of problem.
// XD attr concentration chaine(into=["concentration"]) concentration OPT concentration
// XD attr dim entier dim OPT dimension of the problem
// XD attr val chaine val OPT function of rho
// XD attr Champ_Uniforme chaine(into=["Champ_Uniforme"]) Champ_Uniforme OPT Champ_Uniforme
// XD attr fielddim entier fielddim OPT dimension of the problem
// XD attr val2 chaine val2 OPT function of rho
 
// XD visco_dyn_cons objet_lecture nul NO_BRACE different treatment of the kinematic viscosity could be done depending
// XD_CONT of the use of the Boussinesq approximation or the constant dynamic viscosity approximation
// XD attr yes_or_no chaine(into=["oui","non"]) yes_or_no REQ To use or not the constant dynamic viscosity
// XD attr bloc_visco bloc_visco2 bloc_visco REQ to choose the mu formulation
 
// XD bloc_visco2 objet_lecture nul BRACE choice of mu formulation
// XD attr probleme ref_Pb_base probleme OPT Name of problem.
// XD attr mu_1 floattant mu_1 OPT value of mu
// XD attr mu_2 floattant mu_2 OPT value of mu
// XD attr mu_fonc_c bloc_mu_fonc_c mu_fonc_c_ OPT to use for define a general form for mu
 
// XD bloc_mu_fonc_c objet_lecture nul NO_BRACE if mu has a general form
// XD attr Champ_Fonc_Fonction chaine(into=["Champ_Fonc_Fonction"]) Champ_Fonc_Fonction OPT Champ_Fonc_Fonction
// XD attr problem_name ref_Pb_base problem_name OPT Name of problem.
// XD attr concentration chaine(into=["concentration"]) concentration OPT concentration
// XD attr dim entier dim OPT dimension of the problem
// XD attr val chaine val OPT function of mu
 
Sortie& Milieu_Phase_field::printOn(Sortie& os) const
{
  return os << que_suis_je() << " " << le_nom() << finl;
}
 
Entree& Milieu_Phase_field::readOn(Entree& is)
{
  return Fluide_Incompressible::readOn(is);
}
 
void Milieu_Phase_field::set_param(Param& param) const
{
  // XD milieu_phase_field fluide_incompressible milieu_phase_field INHERITS_BRACE Class for Phase_field fluid.
  Fluide_Incompressible::set_param(param);
  param.ajouter("Fermeture|Closure", &fermeture_, Param::REQUIRED);
  param.ajouter_non_std("approximation_de_boussinesq", (this), Param::REQUIRED); // XD_ADD_P approx_boussinesq
  // XD_CONT To use or not the Boussinesq approximation.
  param.ajouter_non_std("viscosite_dynamique_constante", (this), Param::OPTIONAL); // XD_ADD_P visco_dyn_cons
  // XD_CONT To use or not a viscosity which will depends on concentration C (in fact, C is the unknown of Cahn-Hilliard
  // XD_CONT equation).
}
 
int Milieu_Phase_field::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  Motcle motlu;
  if (mot == "approximation_de_boussinesq")
    {
      Motcle temp_boussi_;
      is >> temp_boussi_;
      bool fonctionsALire = false;
      if (temp_boussi_ == "oui")
        {
          boussi_ = 1;
          diff_boussi_ = 0;
          if (!has_beta_c())
            {
              fonctionsALire = true;
              Cerr << "Approximation de Boussinesq utilisee mais beta_co non defini comme propriete du fluide" << finl;
              Cerr << "On s'attend donc a lire les fonctions rho(c) et drho/dc" << finl;
            }
        }
      else
        {
          boussi_ = 0;
          fonctionsALire = true;
        }
      is >> motlu;
      if (motlu != "{")
        {
          Cerr << "Error while reading approximation_de_boussinesq" << finl;
          Cerr << "We expected a { instead of " << motlu << finl;
          Process::exit();
        }
      is >> motlu;
      if (fonctionsALire)
        {
          if (motlu == "rho_fonc_c")
            {
              Cerr << motlu << finl;
              is >> rho_;
              is >> drhodc_;
              is >> motlu;
            }
          else
            {
              double rho1, rho2;
              Nom prob;
              while (motlu != "}")
                {
                  if (motlu == "probleme")
                    is >> prob;
                  if (motlu == "rho_1")
                    is >> rho1;
                  if (motlu == "rho_2")
                    is >> rho2;
                  is >> motlu;
                }
              Nom chaine("Champ_Fonc_Fonction ");
              chaine += prob;
              chaine += " concentration 1 ";
              Nom rhoM(0.5 * (rho1 + rho2));
              Nom drho(rho2 - rho1);
              chaine += rhoM;
              chaine += "+val*(";
              chaine += drho;
              chaine += ")";
              EChaine echaine(chaine);
              echaine >> rho_;
              Nom chaine2("Champ_Uniforme 1 ");
              chaine2 += drho;
              EChaine echaine2(chaine2);
              Cerr << "@@@ " << chaine2 << finl;
              echaine2 >> drhodc_;
            }
        }
      if (motlu != "}")
        {
          Cerr << "Error while reading approximation_de_boussinesq" << finl;
          Cerr << "We expected a } instead of " << motlu << finl;
          Process::exit();
        }
      return 1;
    }
  else if (mot == "viscosite_dynamique_constante")
    {
      Motcle temp_diffbou_;
      is >> temp_diffbou_;
      if (temp_diffbou_ == "oui")
        {
          diff_boussi_ = 1;
        }
      else
        {
          diff_boussi_ = 0;
          is >> motlu;
          assert(motlu == "{");
          is >> motlu;
          if (motlu == "mu_fonc_c")
            {
              Cerr << motlu << finl;
              is >> mu_;
              is >> motlu;
            }
          else
            {
              double mu1, mu2;
              Nom prob;
              while (motlu != "}")
                {
                  if (motlu == "probleme")
                    is >> prob;
                  if (motlu == "mu_1")
                    is >> mu1;
                  if (motlu == "mu_2")
                    is >> mu2;
                  is >> motlu;
                }
              Nom chaine("Champ_Fonc_Fonction ");
              chaine += prob;
              chaine += " concentration 1 ";
              Nom muM(0.5 * (mu1 + mu2));
              Nom dmu(mu2 - mu1);
              chaine += muM;
              chaine += "+val*(";
              chaine += dmu;
              chaine += ")";
              EChaine echaine(chaine);
              echaine >> mu_;
            }
          assert(motlu == "}");
        }
      if (diff_boussi_ == 1)
        {
          // RLT: pas sur de tout a fait comprendre si c'est encore d'actualite mais, en tout cas, en lien avec modification
          // commentee dans Navier_Stokes_phase_field::calculer_pas_de_temps
          Cerr << "LE PAS DE TEMPS DE STABILITE (DE DIFFUSION) DOIT ETRE FORCE DANS VDF/Operateurs/OpDifVDFFaCs.cpp : ligne 115" << finl;
          Cerr << "ATTENTION : LE PAS DE TEMPS DE DIFFUSION SERA ALORS FORCE A LA VALEUR DE 1 SECONDE !!!" << finl;
          // Ceci est utile car normalement dans ce cas le pas de temps de diffusion ne depend pas de rho
          // car la viscosite dynamique est constante (on prend alors mu/rho0)
          // Or dans ce cas TRUST semble calculer ce pas temps avec une influence venant de rho (Debug ?)
        }
      return 1;
    }
  else
    return Fluide_Incompressible::lire_motcle_non_standard(mot, is);
}
 
void Milieu_Phase_field::discretiser(const Probleme_base& pb, const Discretisation_base& dis)
{
  pb_pf_ = ref_cast(Pb_Phase_field, pb);
  Fluide_Incompressible::discretiser(pb, dis);
 
  if (mu_)
    {
      dis.nommer_completer_champ_physique(pb.domaine_dis(), "viscosite_dynamique", "Pa.s", mu_, pb);
      champs_compris_.ajoute_champ(mu_);
    }
 
  const DoubleTab& rho0Tab = masse_volumique().valeurs();
  int dim = rho0Tab.nb_dim();
  switch(dim)
    {
    case 2:
      rho0_ = rho0Tab(0, 0);
      break;
    case 3:
      rho0_ = rho0Tab(0, 0, 0);
      break;
    default:
      Cerr << "Milieu_Phase_field : Problem with rho0 dimension:" << dim << finl;
      Process::exit();
      break;
    }
 
//  if (has_beta_c())
//    betac_ = beta_c();
 
  if (!rho_)
    {
      if (boussi_ != 1)
        Process::exit("Milieu_Phase_field::discretiser : should not be here (1)");
 
      dis.discretiser_champ("CHAMP_ELEM", pb.domaine_dis(), "masse_volumique", "kg/m3", 1, 0., rho_);
      champs_compris_.ajoute_champ(rho_);
      dis.discretiser_champ("CHAMP_ELEM", pb.domaine_dis(), "derivee_masse_volumique", "kg/m3", 1, 0., drhodc_);
      champs_compris_.ajoute_champ(drhodc_);
    }
  if (rho_->le_nom() == "anonyme")
    {
      if (boussi_ == 1 && has_beta_c())
        Process::exit("Milieu_Phase_field::discretiser : should not be here (2)");
 
      dis.nommer_completer_champ_physique(pb.domaine_dis(), "masse_volumique", "kg/m3", rho_, pb);
      champs_compris_.ajoute_champ(rho_);
      dis.nommer_completer_champ_physique(pb.domaine_dis(), "derivee_masse_volumique", "kg/m3", drhodc_, pb);
      champs_compris_.ajoute_champ(drhodc_);
    }
 
}
 
/*! @brief Calcule rho_n+1 pour utilisation dans la pression
 *
 */
void Milieu_Phase_field::calculer_rho(const bool init)
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_pf_->equation(1));
  Sources& list_sources = ref_cast_non_const(Sources, eq_c.sources());
  int type_systeme_naire = fermeture_->get_type_systeme_naire();
 
  int i = 1;
  double temps = pb_pf_->schema_temps().temps_futur(i);
  if (init)
    {
      i = 0;
      temps = pb_pf_->schema_temps().temps_courant();
    }
  if (boussi_ == 0 || !has_beta_c())
    {
      // normalement, quelque chose comme 'rho_->mettre_a_jour(schema_temps().temps_futur(i))' devrait faire l'affaire
      // probleme, cette evaluation se fait avec c(n) alors que l'on veut evaluer avec c(n+1)
      // on introduit une methode specifique (qui depend de la discretisation) que l'on met arbitrairement dans Source_Con_Phase_field_base
      Source_Con_Phase_field_base& source_pf_base = ref_cast(Source_Con_Phase_field_base, list_sources(0).valeur());
      source_pf_base.calculer_champ_fonc_c(temps, rho_, eq_c.inconnue().futur(i));
      source_pf_base.calculer_champ_fonc_c(temps, drhodc_, eq_c.inconnue().futur(i));
    }
  else
    {
      const DoubleTab& c = eq_c.inconnue().futur(i);
      //Mirantsoa 264902
      DoubleTab& rhoTab = rho_->valeurs(); /**/
      DoubleTab& drhodcTab = drhodc_->valeurs();/**/
      DoubleTab rho0Tab(rhoTab);
 
      if (type_systeme_naire == 0)
        {
          drhodcTab = rho0_;
          tab_multiply_any_shape(drhodcTab, beta_c().valeurs());
          rhoTab = c;
          tab_multiply_any_shape(rhoTab, drhodcTab);
          rhoTab += rho0_;
          //Cerr << "c = "<< c<<finl;
 
        }
      else if (type_systeme_naire == 1)
        {
          DoubleTab beta_(rhoTab);
          DoubleTab betacTab(c);
 
          //calcul de drhodcTab comme rho0*somme(betac)
          drhodcTab = rho0_;
          for (i = 0; i < beta_.dimension(0); i++)
            {
              for (int j = 0; j < betacTab.line_size(); j++)
                {
                  betacTab(i, j) = beta_c().valeurs()(0, j);
                  beta_(i) += betacTab(i, j);
                }
            }
          tab_multiply_any_shape(drhodcTab, beta_);
 
          //calcul de rhoTab comme rho0+rho0*somme(betac*c)
          rhoTab = 0;
          tab_multiply_any_shape(betacTab, c);
          for (i = 0; i < betacTab.dimension(0); i++)
            {
              for (int j = 0; j < betacTab.line_size(); j++)
                {
                  rhoTab(i) += betacTab(i, j);
                }
            }
          rho0Tab = rho0_;
          tab_multiply_any_shape(rhoTab, rho0Tab);
          rhoTab += rho0_;
        }
    }
  rho_->valeurs().echange_espace_virtuel();
  drhodc_->valeurs().echange_espace_virtuel();
}
 
void Milieu_Phase_field::calculer_mu(const bool init)
{
  if (boussi_ == 0 && diff_boussi_ == 0)
    {
      const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_pf_->equation(1));
      int i = 1;
      double temps = pb_pf_->schema_temps().temps_futur(i);
      if (init)
        {
          i = 0;
          temps = pb_pf_->schema_temps().temps_courant();
        }
      // normalement, quelque chose comme 'mu_->mettre_a_jour(schema_temps().temps_futur(i))' devrait faire l'affaire
      // probleme, cette evaluation se fait avec c(n) alors que l'on veut evaluer avec c(n+1)
      // on introduit une methode specifique (qui depend de la discretisation) que l'on met arbitrairement dans Source_Con_Phase_field_base
      Sources& list_sources = ref_cast_non_const(Sources, eq_c.sources());
      Source_Con_Phase_field_base& source_pf = ref_cast(Source_Con_Phase_field_base, list_sources(0).valeur());
      source_pf.calculer_champ_fonc_c(temps, mu_, eq_c.inconnue().futur(i));
      mu_->valeurs().echange_espace_virtuel();
    }
}
 
int Milieu_Phase_field::initialiser(const double temps)
{
  Fluide_Incompressible::initialiser(temps);
  Cerr << "Milieu_Phase_field::initialiser()" << finl;
  // Calcul de la masse volumique
  rho_->initialiser(temps);
 
  drhodc_->initialiser(temps);
 
  calculer_rho(true); // RLT: cette minitialisation etait incorrecte dans TrioCFD 1.8.0 car on utilisait c.futur()
  // Calcul de la viscosite dynamique dans le cas boussi_==0 et diff_boussi_==0
  if (mu_)
    {
      mu_->initialiser(temps);
      calculer_mu(true); // RLT: cette minitialisation etait incorrecte dans dans TrioCFD 1.8.0 car on utilisait c.futur()
    }
 
  rho_->changer_temps(temps);
  drhodc_->changer_temps(temps);
 
  if (mu_)
    mu_->changer_temps(temps);
 
  return 1;
}
 
void Milieu_Phase_field::mettre_a_jour(double temps)
{
  Fluide_Incompressible::mettre_a_jour(temps);
}
 
/*
 * A voir si l'ordre est pas important, on fait ca dans mettre_a_jour ...
 */
void Milieu_Phase_field::update_rho_mu(double temps)
{
  calculer_rho();
  rho_->changer_temps(temps);
  drhodc_->changer_temps(temps);
  if (mu_)
    {
      calculer_mu(); // RLT: cette mise a jour etait manquante dans TrioCFD 1.8.0
      mu_->changer_temps(temps);
    }
}